← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/special/cephes/unity.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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