source: sasview/src/sas/models/c_extension/cephes/unity.c @ 79492222

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 79492222 was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.

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