1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4
5 #undef NJ
6 #define NJ 9
7
8 #undef NPT
9 #define NPT ((1<<NJ)*10)
10
ibeta(double xcut,double a,double b,double * f1,double * flogx,double * flog1x)11 void ibeta( double xcut , double a , double b ,
12 double *f1 , double *flogx , double *flog1x )
13 {
14 int nn , ii , jj , istep ;
15 double dx , s1,sx,s1x , x, a1=a-1.0,b1=b-1.0 ;
16 double as1[NJ+1] , asx[NJ+1] , as1x[NJ+1] ;
17 double bs1[NJ-1] , bsx[NJ-1] , bs1x[NJ-1] ;
18 double cs1[NJ-3] , csx[NJ-3] , cs1x[NJ-3] ;
19 double xab[NPT] , xabx[NPT] , xab1x[NPT] ;
20
21 if( xcut <= 0.0 || xcut >= 1.0 || a <= 0.0 || b <= 0.0 ) return ;
22
23 if( f1 == NULL && flogx == NULL && flog1x == NULL ) return ;
24
25 dx = xcut / NPT ;
26 for( ii=1 ; ii <= NPT ; ii++ ){
27 x = ii*dx ;
28 /* fprintf(stderr,"ii=%d dx=%g a1=%g b1=%g x=%g\n",ii,dx,a1,b1,x) ; */
29 xab[ii] = pow(x,a1) * pow(1.0-x,b1) ;
30 xabx[ii] = xab[ii] * log(x) ;
31 xab1x[ii] = xab[ii] * log(1.0-x) ;
32
33 /* fprintf(stderr,"ii=%d x=%g xab=%g xabx=%g xab1x=%g a1=%g b1=%g dx=%g\n",
34 ii,x,xab[ii],xabx[ii],xab1x[ii],a1,b1,dx) ; */
35 }
36
37 for( nn=NPT,istep=1,jj=NJ ; jj >= 0 ; jj--,istep*=2,nn/=2 ){
38 s1 = sx = s1x = 0.0 ;
39 for( ii=istep ; ii <= NPT ; ii+=istep ){
40 s1 += xab[ii] ; sx += xabx[ii] ; s1x += xab1x[ii] ;
41 }
42 dx = xcut / nn ;
43 as1[jj] = s1*dx ; asx[jj] = sx*dx ; as1x[jj] = s1x*dx ;
44 }
45
46 #undef AITKEN
47 #define AITKEN(x,y,z) ((x) - ((y)-(x))*((y)-(x)) / (((z)-(y))-((y)-(x))) )
48
49 for( jj=0 ; jj <= NJ-2 ; jj++ ){
50 bs1[jj] = AITKEN( as1[jj] , as1[jj+1] , as1[jj+2] ) ;
51 bsx[jj] = AITKEN( asx[jj] , asx[jj+1] , asx[jj+2] ) ;
52 bs1x[jj] = AITKEN( as1x[jj] , as1x[jj+1] , as1x[jj+2] ) ;
53 }
54
55 for( jj=0 ; jj <= NJ-4 ; jj++ ){
56 cs1[jj] = AITKEN( bs1[jj] , bs1[jj+1] , bs1[jj+2] ) ;
57 csx[jj] = AITKEN( bsx[jj] , bsx[jj+1] , bsx[jj+2] ) ;
58 cs1x[jj] = AITKEN( bs1x[jj] , bs1x[jj+1] , bs1x[jj+2] ) ;
59 }
60
61 #define DMP(tbl,jt) \
62 do{ fprintf(stderr," table " #tbl ":") ; \
63 for( jj=0 ; jj <= jt ; jj++ ) fprintf(stderr," %12.5g",tbl[jj]) ; \
64 fprintf(stderr,"\n") ; } while(0)
65
66 fprintf(stderr,"ibeta(xc=%g,a=%g,b=%g)\n",xcut,a,b) ;
67 DMP(as1,NJ) ; DMP(bs1,NJ-2) ; DMP(cs1,NJ-4) ;
68 DMP(asx,NJ) ; DMP(bsx,NJ-2) ; DMP(csx,NJ-4) ;
69 DMP(as1x,NJ) ; DMP(bs1x,NJ-2) ; DMP(cs1x,NJ-4) ;
70
71 if( f1 != NULL ) *f1 = cs1[NJ-4] ;
72 if( flogx != NULL ) *flogx = csx[NJ-4] ;
73 if( flog1x != NULL ) *flog1x = cs1x[NJ-4] ;
74
75 return ;
76 }
77
main(int argc,char * argv[])78 int main( int argc , char * argv[] )
79 {
80 double xcut , a , b , f1,f2,f3 ;
81
82 if( argc < 4 ) exit(0) ;
83 xcut = strtod(argv[1],NULL) ;
84 a = strtod(argv[2],NULL) ;
85 b = strtod(argv[3],NULL) ;
86
87 ibeta( xcut,a,b , &f1,&f2,&f3 ) ; exit(0) ;
88 }
89