"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/cephes/unity.c" (2 Oct 2015, 2265 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "unity.c" see the Fossies "Dox" file reference documentation.

    1 /*                          unity.c
    2  *
    3  * Relative error approximations for function arguments near
    4  * unity.
    5  *
    6  *    log1p(x) = log(1+x)
    7  *    expm1(x) = exp(x) - 1
    8  *    cosm1(x) = cos(x) - 1
    9  *
   10  */
   11 
   12 #include "mconf.h"
   13 
   14 /* log1p(x) = log(1 + x)  */
   15 
   16 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
   17  * 1/sqrt(2) <= x < sqrt(2)
   18  * Theoretical peak relative error = 2.32e-20
   19  */
   20 
   21 static double LP[] = {
   22  4.5270000862445199635215E-5,
   23  4.9854102823193375972212E-1,
   24  6.5787325942061044846969E0,
   25  2.9911919328553073277375E1,
   26  6.0949667980987787057556E1,
   27  5.7112963590585538103336E1,
   28  2.0039553499201281259648E1,
   29 };
   30 
   31 static double LQ[] = {
   32 /* 1.0000000000000000000000E0,*/
   33  1.5062909083469192043167E1,
   34  8.3047565967967209469434E1,
   35  2.2176239823732856465394E2,
   36  3.0909872225312059774938E2,
   37  2.1642788614495947685003E2,
   38  6.0118660497603843919306E1,
   39 };
   40 
   41 #define SQRTH 0.70710678118654752440
   42 #define SQRT2 1.41421356237309504880
   43 
   44 double cephes_log (double x)
   45 {
   46     double z = 1.0 + x;
   47 
   48     if (z < SQRTH || z > SQRT2) {
   49     return log(z);
   50     }
   51 
   52     z = x*x;
   53     z = -0.5 * z + x * (z * polevl(x, LP, 6) / p1evl(x, LQ, 6));
   54 
   55     return x + z;
   56 }
   57 
   58 /* expm1(x) = exp(x) - 1  */
   59 
   60 /*  e^x =  1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
   61  * -0.5 <= x <= 0.5
   62  */
   63 
   64 static double EP[3] = {
   65  1.2617719307481059087798E-4,
   66  3.0299440770744196129956E-2,
   67  9.9999999999999999991025E-1,
   68 };
   69 
   70 static double EQ[4] = {
   71  3.0019850513866445504159E-6,
   72  2.5244834034968410419224E-3,
   73  2.2726554820815502876593E-1,
   74  2.0000000000000000000897E0,
   75 };
   76 
   77 double cephes_exp (double x)
   78 {
   79     double r, xx;
   80 
   81     if (isnan(x)) {
   82     return x;
   83     }
   84 
   85     if (!isfinite(x)) {
   86     return (x < 0)? -1.0 : x;
   87     }
   88 
   89     if (x < -0.5 || x > 0.5) {
   90     return exp(x) - 1.0;
   91     }
   92 
   93     xx = x * x;
   94     r = x * polevl(xx, EP, 2);
   95     r = r/(polevl(xx, EQ, 3) - r);
   96 
   97     return r + r;
   98 }
   99 
  100 /* cosm1(x) = cos(x) - 1  */
  101 
  102 static double coscof[7] = {
  103  4.7377507964246204691685E-14,
  104 -1.1470284843425359765671E-11,
  105  2.0876754287081521758361E-9,
  106 -2.7557319214999787979814E-7,
  107  2.4801587301570552304991E-5,
  108 -1.3888888888888872993737E-3,
  109  4.1666666666666666609054E-2,
  110 };
  111 
  112 double cosm1 (double x)
  113 {
  114     double xx;
  115 
  116     if (x < -PIO4 || x > PIO4) {
  117     return cos(x) - 1.0;
  118     }
  119 
  120     xx = x * x;
  121     xx = -0.5*xx + xx * xx * polevl(xx, coscof, 6);
  122 
  123     return xx;
  124 }