double-math.c

     
   1  //! @file double-math.c
   2  //! @author J. Marcel van der Veer
   3  //!
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  //!
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! LONG REAL, LONG COMPLEX routines. 
  25  
  26  // References:
  27  //
  28  //   Milton Abramowitz and Irene Stegun, Handbook of Mathematical Functions,
  29  //   Dover Publications, New York [1970]
  30  //   https://en.wikipedia.org/wiki/Abramowitz_and_Stegun
  31  
  32  #include "a68g.h"
  33  #include "a68g-genie.h"
  34  #include "a68g-prelude.h"
  35  #include "a68g-double.h"
  36  #include "a68g-numbers.h"
  37  #include "a68g-math.h"
  38  
  39  #if (A68_LEVEL >= 3)
  40  
  41  DOUBLE_T a68_beta_inc_double_real (DOUBLE_T s, DOUBLE_T t, DOUBLE_T x)
  42  {
  43  // Incomplete beta function I{x}(s, t).
  44  // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
  45    if (x < 0.0q || x > 1.0q) {
  46      errno = ERANGE;
  47      return -1.0q;
  48    } else {
  49      const INT_T lim = 16 * sizeof (DOUBLE_T);
  50  // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
  51      if (x > (s + 1.0q) / (s + t + 2.0q)) {
  52  // B{x}(s, t) = 1 - B{1-x}(t, s)
  53        return 1.0q - a68_beta_inc_double_real (s, t, 1.0q - x);
  54      }
  55  // Lentz's algorithm for continued fraction.
  56      DOUBLE_T W = 1.0q, F = 1.0q, c = 1.0q, d = 0.0q;
  57      BOOL_T cont = A68_TRUE;
  58      for (int N = 0, m = 0; cont && N < lim; N++) {
  59        DOUBLE_T T;
  60        if (N == 0) {
  61          T = 1.0q;
  62        } else if (N % 2 == 0) {
  63  // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
  64          T = x * m * (t - m) / (s + 2.0q * m - 1.0q) / (s + 2.0q * m);
  65        } else {
  66  // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
  67          T = -x * (s + m) * (s + t + m) / (s + 2.0q * m + 1.0q) / (s + 2.0q * m);
  68          m++;
  69        }
  70        d = 1.0q / (T * d + 1.0q);
  71        c = T / c + 1.0q;
  72        F *= c * d;
  73        if (F == W) {
  74          cont = A68_FALSE;
  75        } else {
  76          W = F;
  77        }
  78      }
  79  // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
  80      DOUBLE_T beta = expq (lgammaq (s) + lgammaq (t) - lgammaq (s + t));
  81      return powq (x, s) * powq (1.0q - x, t) / s / beta * (F - 1.0q);
  82    }
  83  }
  84  
  85  //! @brief PROC (LONG REAL) LONG REAL csc
  86  
  87  DOUBLE_T csc_double_real (DOUBLE_T x)
  88  {
  89    DOUBLE_T z = sinq (x);
  90    A68_OVERFLOW (z == 0.0q);
  91    return 1.0q / z;
  92  }
  93  
  94  //! @brief PROC (LONG REAL) LONG REAL acsc
  95  
  96  DOUBLE_T acsc_double_real (DOUBLE_T x)
  97  {
  98    A68_OVERFLOW (x == 0.0q);
  99    return asinq (1.0q / x);
 100  }
 101  
 102  //! @brief PROC (LONG REAL) LONG REAL sec
 103  
 104  DOUBLE_T sec_double_real (DOUBLE_T x)
 105  {
 106    DOUBLE_T z = cosq (x);
 107    A68_OVERFLOW (z == 0.0q);
 108    return 1.0q / z;
 109  }
 110  
 111  //! @brief PROC (LONG REAL) LONG REAL asec
 112  
 113  DOUBLE_T asec_double_real (DOUBLE_T x)
 114  {
 115    A68_OVERFLOW (x == 0.0q);
 116    return acosq (1.0q / x);
 117  }
 118  
 119  //! @brief PROC (LONG REAL) LONG REAL cot
 120  
 121  DOUBLE_T cot_double_real (DOUBLE_T x)
 122  {
 123    DOUBLE_T z = sinq (x);
 124    A68_OVERFLOW (z == 0.0q);
 125    return cosq (x) / z;
 126  }
 127  
 128  //! @brief PROC (LONG REAL) LONG REAL acot
 129  
 130  DOUBLE_T acot_double_real (DOUBLE_T x)
 131  {
 132    A68_OVERFLOW (x == 0.0q);
 133    return atanq (1 / x);
 134  }
 135  
 136  //! brief PROC (LONG REAL) LONG REAL sindg
 137  
 138  DOUBLE_T sindg_double_real (DOUBLE_T x)
 139  {
 140    return sinq (x * CONST_PI_OVER_180_Q);
 141  }
 142  
 143  //! brief PROC (LONG REAL) LONG REAL cosdg
 144  
 145  DOUBLE_T cosdg_double_real (DOUBLE_T x)
 146  {
 147    return cosq (x * CONST_PI_OVER_180_Q);
 148  }
 149  
 150  //! brief PROC (LONG REAL) LONG REAL tandg
 151  
 152  DOUBLE_T tandg_double_real (DOUBLE_T x)
 153  {
 154    return tanq (x * CONST_PI_OVER_180_Q);
 155  }
 156  
 157  //! brief PROC (LONG REAL) LONG REAL asindg
 158  
 159  DOUBLE_T asindg_double_real (DOUBLE_T x)
 160  {
 161    return asinq (x) * CONST_180_OVER_PI_Q;
 162  }
 163  
 164  //! brief PROC (LONG REAL) LONG REAL acosdg
 165  
 166  DOUBLE_T acosdg_double_real (DOUBLE_T x)
 167  {
 168    return acosq (x) * CONST_180_OVER_PI_Q;
 169  }
 170  
 171  //! brief PROC (LONG REAL) LONG REAL atandg
 172  
 173  DOUBLE_T atandg_double_real (DOUBLE_T x)
 174  {
 175    return atanq (x) * CONST_180_OVER_PI_Q;
 176  }
 177  
 178  // PROC (LONG REAL) LONG REAL cotdg
 179  
 180  DOUBLE_T cotdg_double_real (DOUBLE_T x)
 181  {
 182    DOUBLE_T z = a68_sindg (x);
 183    A68_OVERFLOW (z == 0);
 184    return a68_cosdg (x) / z;
 185  }
 186  
 187  // PROC (LONG REAL) LONG REAL acotdg
 188  
 189  DOUBLE_T acotdg_double_real (DOUBLE_T z)
 190  {
 191    A68_OVERFLOW (z == 0);
 192    return a68_atandg (1 / z);
 193  }
 194  
 195  // @brief PROC (LONG REAL) LONG REAL sinpi
 196  
 197  DOUBLE_T sinpi_double_real (DOUBLE_T x)
 198  {
 199    x = fmodq (x, 2.0q);
 200    if (x <= -1.0q) {
 201      x += 2.0q;
 202    } else if (x > 1.0q) {
 203      x -= 2.0q;
 204    }
 205  // x in <-1, 1].
 206    if (x == 0.0q || x == 1.0q) {
 207      return 0.0q;
 208    } else if (x == 0.5q) {
 209      return 1.0q;
 210    }
 211    if (x == -0.5q) {
 212      return -1.0q;
 213    } else {
 214      return sinq (CONST_PI_Q * x);
 215    }
 216  }
 217  
 218  // @brief PROC (LONG REAL) LONG REAL cospi
 219  
 220  DOUBLE_T cospi_double_real (DOUBLE_T x)
 221  {
 222    x = fmodq (fabsq (x), 2.0q);
 223  // x in [0, 2>.
 224    if (x == 0.5q || x == 1.5q) {
 225      return 0.0q;
 226    } else if (x == 0.0q) {
 227      return 1.0q;
 228    } else if (x == 1.0q) {
 229      return -1.0q;
 230    } else {
 231      return cosq (CONST_PI_Q * x);
 232    }
 233  }
 234  
 235  // @brief PROC (LONG REAL) LONG REAL tanpi
 236  
 237  DOUBLE_T tanpi_double_real (DOUBLE_T x)
 238  {
 239    x = fmodq (x, 1.0q);
 240    if (x <= -0.5q) {
 241      x += 1.0q;
 242    } else if (x > 0.5q) {
 243      x -= 1.0q;
 244    }
 245  // x in <-1/2, 1/2].
 246    A68_OVERFLOW (x == 0.5q);
 247    if (x == -0.25q) {
 248      return -1.0q;
 249    } else if (x == 0) {
 250      return 0.0q;
 251    } else if (x == 0.25q) {
 252      return 1.0q;
 253    } else {
 254      return sinpi_double_real (x) / cospi_double_real (x);
 255    }
 256  }
 257  
 258  // @brief PROC (LONG REAL) LONG REAL cotpi
 259  
 260  DOUBLE_T cotpi_double_real (DOUBLE_T x)
 261  {
 262    x = fmodq (x, 1.0q);
 263    if (x <= -0.5q) {
 264      x += 1.0q;
 265    } else if (x > 0.5q) {
 266      x -= 1.0q;
 267    }
 268  // x in <-1/2, 1/2].
 269    A68_OVERFLOW (x == 0.0q);
 270    if (x == -0.25q) {
 271      return -1.0q;
 272    } else if (x == 0.25q) {
 273      return 1.0q;
 274    } else if (x == 0.5q) {
 275      return 0.0q;
 276    } else {
 277      return cospi_double_real (x) / sinpi_double_real (x);
 278    }
 279  }
 280  
 281  #endif