SourceBrowser

experimental.c /size: 6348 b last modification: 2025-02-21 11:03
1/******************************************************************************/
2/*  Experimental code                                                         */
3/******************************************************************************/
4
5/*
6    Compute w_of_z via Fourier integration using Ooura-Mori transform.
7    Agreement with Johnson's code usually < 1E-15, so far always < 1E-13.
8    Todo:
9    - sign for negative x or y
10    - determine application limits
11    - more systematical comparison with Johnson's code
12    - comparison with Abrarov&Quine
13 */
14
15#define max_iter_int 10
16#define num_range 5
17#define PI           3.14159265358979323846L  /* pi */
18#define SQR(x) ((x)*(x))
19#include <errno.h>
20
21double cerf_experimental_integration( int kind, double x, double y )
22// kind: 0 cos, 1 sin transform (precomputing arrays[2] depend on this)
23{
24    // unused parameters
25    static int mu = 0;
26    int intgr_debug = 0;
27    static double intgr_delta=2.2e-16, intgr_eps=5.5e-20;
28
29    if( x<0 || y<0 ) {
30        fprintf( stderr, "negative arguments not yet implemented\n" );
31        exit( EDOM );
32    }
33
34    double w = sqrt(2)*x;
35    double gamma = sqrt(2)*y;
36
37    int iter;
38    int kaux;
39    int isig;
40    int N;
41    int j;               // range
42    long double S=0;     // trapezoid sum
43    long double S_last;  // - in last iteration
44    long double s;       // term contributing to S
45    long double T;       // sum of abs(s)
46    // precomputed coefficients
47    static int firstCall=1;
48    static int iterDone[2][num_range]; // Nm,Np,ak,bk are precomputed up to this
49    static int Nm[num_range][max_iter_int];
50    static int Np[num_range][max_iter_int];
51    static long double *ak[2][num_range][max_iter_int];
52    static long double *bk[2][num_range][max_iter_int];
53    // auxiliary for computing ak and bk
54    long double u;
55    long double e;
56    long double tk;
57    long double chi;
58    long double dchi;
59    long double h;
60    long double k;
61    long double f;
62    long double ahk;
63    long double chk;
64    long double dhk;
65    double p;
66    double q;
67    const double Smin=2e-20; // to assess worst truncation error
68
69    // dynamic initialization upon first call
70    if ( firstCall ) {
71        for ( j=0; j<num_range; ++ j ) {
72            iterDone[0][j] = -1;
73            iterDone[1][j] = -1;
74        }
75        firstCall = 0;
76    }
77
78    // determine range, set p,q
79    j=1; p=1.4; q=0.6;
80
81    // iterative integration
82    if( intgr_debug & 4 )
83        N = 100;
84    else
85        N = 40;
86    for ( iter=0; iter<max_iter_int; ++iter ) {
87        // static initialisation of Nm, Np, ak, bk for given 'iter'
88        if ( iter>iterDone[kind][j] ) {
89            if ( N>1e6 )
90                return -3; // integral limits overflow
91            Nm[j][iter] = N;
92            Np[j][iter] = N;
93            if ( !( ak[kind][j][iter]=malloc((sizeof(long double))*
94                                         (Nm[j][iter]+1+Np[j][iter])) ) ||
95                !( bk[kind][j][iter]=malloc((sizeof(long double))*
96                                         (Nm[j][iter]+1+Np[j][iter])) ) ) {
97                fprintf( stderr, "Workspace allocation failed\n" );
98                exit( ENOMEM );
99            }
100            h = logl( logl( 42*N/intgr_delta/Smin ) / p ) / N; // 42=(pi+1)*10
101            isig=1-2*(Nm[j][iter]&1);
102            for ( kaux=-Nm[j][iter]; kaux<=Np[j][iter]; ++kaux ) {
103                k = kaux;
104                if( !kind )
105                    k -= 0.5;
106                u = k*h;
107                chi  = 2*p*sinhl(u) + 2*q*u;
108                dchi = 2*p*coshl(u) + 2*q;
109                if ( u==0 ) {
110                    if ( k!=0 )
111                        return -4; // integration variable underflow
112                    // special treatment to bridge singularity at u=0
113                    ahk = PI/h/dchi;
114                    dhk = 0.5;
115                    chk = sin( ahk );
116                } else {
117                    if ( -chi>DBL_MAX_EXP/2 )
118                        return -5; // integral transformation overflow
119                    e = expl( -chi );
120                    ahk = PI/h * u/(1-e);
121                    dhk = 1/(1-e) - u*e*dchi/SQR(1-e);
122                    chk = e>1 ?
123                        ( kind ? sinl( PI*k/(1-e) ) : cosl( PI*k/(1-e) ) ) :
124                        isig * sinl( PI*k*e/(1-e) );
125                }
126                ak[kind][j][iter][kaux+Nm[j][iter]] = ahk;
127                bk[kind][j][iter][kaux+Nm[j][iter]] = dhk * chk;
128                isig = -isig;
129            }
130            iterDone[kind][j] = iter;
131        }
132        // integrate according to trapezoidal rule
133        S_last = S;
134        S = 0;
135        T = 0;
136        for ( kaux=-Nm[j][iter]; kaux<=Np[j][iter]; ++kaux ) {
137            tk = ak[kind][j][iter][kaux+Nm[j][iter]] / w;
138            f = expl(-tk*gamma-SQR(tk)/2); // Fourier kernel
139            if ( mu )
140                f /= tk; // TODO
141            s = bk[kind][j][iter][kaux+Nm[j][iter]] * f;
142            S += s;
143            T += fabsl(s);
144            if( intgr_debug & 2 )
145                printf( "%2i %6i %12.4Lg %12.4Lg"
146                        " %12.4Lg %12.4Lg %12.4Lg %12.4Lg\n",
147                        iter, kaux, ak[kind][j][iter][kaux+Nm[j][iter]],
148                        bk[kind][j][iter][kaux+Nm[j][iter]], f, s, S, T );
149        }
150        if( intgr_debug & 1 )
151            printf( "%23.17Le  %23.17Le\n", S, T );
152        // intgr_num_of_terms += Np[j][iter]-(-Nm[j][iter])+1;
153        // termination criteria
154        if      ( intgr_debug & 4 )
155            return -1; // we want to inspect just one sum
156        else if ( S < 0 )
157            return -6; // cancelling terms lead to negative S
158        else if ( intgr_eps*T > intgr_delta*fabs(S) )
159            return -2; // cancellation
160        else if ( iter && fabs(S-S_last) + intgr_eps*T < intgr_delta*fabs(S) )
161            return S*sqrt(2*PI)/w; // success
162            // factor 2 from int_-infty^+infty = 2 * int_0^+infty
163            // factor pi/w from formula 48 in kww paper
164            // factor 1/sqrt(2*pi) from Gaussian
165        N *= 2; // retry with more points
166    }
167    return -9; // not converged
168}
169
170double cerf_experimental_imw( double x, double y )
171{
172    return cerf_experimental_integration( 1, x, y );
173}
174
175double cerf_experimental_rew( double x, double y )
176{
177    return cerf_experimental_integration( 0, x, y );
178}
179
Source Browser ?
