00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <math.h>
00022 #include <errno.h>
00023 #include <float.h>
00024 #include "moremath.h"
00025
00026 #ifdef ACC
00027 #define ACCOLD ACC
00028 #define ACC 40
00029 #else
00030 #undef ACCOLD
00031 #define ACC 40
00032 #endif
00033
00034 #ifdef BIGNO
00035 #define BIGNOOLD BIGNO
00036 #define BIGNO 1.0e10
00037 #else
00038 #undef BIGNOOLD
00039 #define BIGNO 1.0e10
00040 #endif
00041
00042 #ifdef BIGNI
00043 #define BIGNIOLD BIGNI
00044 #define BIGNI 1.0e-10
00045 #else
00046 #undef BIGNIOLD
00047 #define BIGNI 1.0e-10
00048 #endif
00049
00050 #ifdef ITMAX
00051 #define ITMAXOLD ITMAX
00052 #define ITMAX 100
00053 #else
00054 #undef ITMAXOLD
00055 #define ITMAX 100
00056 #endif
00057
00058 #ifdef EPS
00059 #define EPSOLD EPS
00060 #define EPS 3.0e-7
00061 #else
00062 #undef EPSOLD
00063 #define EPS 3.0e-7
00064 #endif
00065
00066 #ifdef FPMIN
00067 #define FPMINOLD FPMIN
00068 #define FPMIN 1.0e-30
00069 #else
00070 #undef FPMINOLD
00071 #define FPMIN 1.0e-30
00072 #endif
00073
00074
00075 double ansi_asinh(double arg)
00076 {
00077 arg = log(arg+sqrt(arg*arg+1));
00078 return arg;
00079 }
00080
00081
00082 double ansi_acosh(double arg)
00083 {
00084 if (arg < 1)
00085 {
00086 errno = EDOM;
00087 return 0;
00088 };
00089 arg = log(arg+sqrt(arg*arg-1));
00090 return arg;
00091 }
00092
00093
00094 double ansi_atanh(double arg)
00095 {
00096 if ((arg <= -1) || (arg >= 1))
00097 {
00098 errno = EDOM;
00099 return 0;
00100 };
00101 arg = 0.5*log(1+arg)-0.5*log(1-arg);
00102 return arg;
00103 }
00104
00105
00106 double ansi_lgamma(double arg)
00107 {
00108 double x,y,ser,tmp;
00109 static double coeff[6] = {76.18009172947146,-86.50532032941677,24.01409824083091,
00110 -1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
00111 int j;
00112
00113 y = x = arg;
00114 tmp = x+5.5;
00115 tmp -= (x+0.5)*log(tmp);
00116 ser = 1.000000000190015;
00117 for (j = 0; j <= 5;j++) ser += coeff[j]/++y;
00118 return -tmp+log(2.5066282746310005*ser/x);
00119 }
00120
00121
00122
00123 double ansi_j0(double arg)
00124 {
00125 double absarg,z,x,y,ans,ans1,ans2;
00126
00127 if ( (absarg = fabs(arg)) < 8 ){
00128 y = arg*arg;
00129 ans1 = 57568490574.0 + y*(-13362590354.0 + y*(651619640.7 + y*(-11214424.18 + y*(77392.33017 - y*184.9052456))));
00130 ans2 = 57568490411.0 + y*(1029532985.0 + y*(9494680.718 + y*(59272.64853 + y*(267.8532712+y))));
00131 ans = ans1/ans2;
00132 }
00133 else {
00134 z = 8/absarg;
00135 y = z*z;
00136 x = absarg - 0.785398164;
00137 ans1 = 1.0 + y*(-0.1098628627e-2 + y*(0.2734510407e-4 + y*(-0.2073370639e-5 + y*0.2093887211e-6)));
00138 ans2 = -0.1562499995e-1 + y*(0.1430488765e-3 + y*(-0.6911147651e-5 + y*(0.7621095161e-6 - y*0.934935152e-7)));
00139 ans = sqrt( 0.636619772 / absarg ) * (cos(x)*ans1 - z*sin(x)*ans2);
00140 }
00141 return ans;
00142 }
00143
00144
00145
00146 double ansi_j1(double arg)
00147 {
00148 double absarg,z,x,y,ans,ans1,ans2;
00149
00150 if ( (absarg = fabs(arg)) < 8 ){
00151 y = arg*arg;
00152 ans1 = arg*(72362614232.0 + y*(-7895059235.0 + y*(242396853.1 + y*(-2972611.439 + y*(15704.48260 + y*(-30.16036606))))));
00153 ans2 = 144725228442.0 + y*(2300535178.0 + y*(18583304.74 + y*(99447.43394 + y*(376.9991397 + y))));
00154 ans = ans1/ans2;
00155 }
00156 else {
00157 z = 8/absarg;
00158 y = z*z;
00159 x = absarg - 2.356194491;
00160 ans1 = 1.0 + y*(0.183105e-2 + y*(-0.3516396496e-4 + y*(0.2457520174e-5 + y*(-0.240337019e-6))));
00161 ans2 = 0.04687499995 + y*(-0.2002690873e-3 + y*(0.8449199096e-5 + y*(-0.88228987e-6 + y*0.105787412e-6)));
00162 ans = sqrt( 0.636619772 / absarg ) * (cos(x)*ans1 - z*sin(x)*ans2);
00163 if (x < 0) ans = -ans;
00164 }
00165 return ans;
00166 }
00167
00168
00169
00170 double ansi_jn(int n, double arg)
00171 {
00172 int i,m,sumsign;
00173 double absarg,toarg,bj1,bj2,bj3,sum,ans;
00174
00175 if (n < 0) {
00176 errno = EDOM;
00177 return 0;
00178 }
00179 if (n == 0) return ansi_j0(arg);
00180 if (n == 1) return ansi_j1(arg);
00181 absarg = fabs(arg);
00182 if (absarg == 0) return 0;
00183 else {
00184 toarg = 2/absarg;
00185 if ( absarg > n ) {
00186 bj1 = ansi_j0(arg);
00187 bj2 = ansi_j1(arg);
00188 for (i = 1; i < n; i++) {
00189 bj3 = i*toarg*bj2-bj1;
00190 bj1 = bj2;
00191 bj2 = bj3;
00192 }
00193 ans = bj3;
00194 }
00195 else {
00196 m = 2*((n + (int) sqrt(double(ACC*n)))/2 );
00197 sumsign = 0;
00198 bj3 = ans = sum = 0;
00199 bj2 = 1;
00200 for(i = m; i >0; i--) {
00201 bj1 = i*toarg*bj2 - bj3;
00202 bj3 = bj2;
00203 bj2 = bj1;
00204 if (fabs(bj2) > BIGNO ) {
00205 bj2 *= BIGNI;
00206 bj3 *= BIGNI;
00207 ans *= BIGNI;
00208 sum *= BIGNI;
00209 }
00210 if (sumsign) sum += bj2;
00211 sumsign = !sumsign;
00212 if (i == n) ans = bj3;
00213 }
00214 sum = 2*sum - bj2;
00215 ans /= sum;
00216 }
00217 }
00218 return ((arg < 0) && (n & 1)) ? -ans : ans;
00219 }
00220
00221
00222
00223 double ansi_gamma_series(double a, double x)
00224 {
00225 int i;
00226 double sum,del,ap;
00227
00228 if (x <= 0) {
00229 if (x < 0) errno = EDOM;
00230 return 0;
00231 }
00232 else {
00233 ap = a;
00234 del = sum = 1/a;
00235 for (i = 1; i <= ITMAX; i++) {
00236 ap++;
00237 del *= x/ap;
00238 sum += del;
00239 if (fabs(del) < fabs(sum)*EPS) return ( sum*exp(-x + a*log(x) - ansi_lgamma(a) ) );
00240 }
00241 errno = EDOM;
00242 return 0;
00243 }
00244 }
00245
00246
00247
00248 double ansi_gamma_cf(double a, double x)
00249 {
00250 int i;
00251 double an,b,c,d,del,h;
00252
00253 b = x+1-a;
00254 c = 1/FPMIN;
00255 h = d = 1/b;
00256 for (i = 1; i<=ITMAX; i++) {
00257 an = -i*(i-a);
00258 b += 2;
00259 d = an*d + b;
00260 if (fabs(d) < FPMIN ) d = FPMIN;
00261 c = b + an/c;
00262 if (fabs(c) < FPMIN ) c = FPMIN;
00263 d = 1/d;
00264 del = d*c;
00265 h *= del;
00266 if (fabs(del-1) < EPS) break;
00267 }
00268 if (i > ITMAX) errno = EDOM;
00269 return ( exp( -x + a*log(x) - ansi_lgamma(a) ) * h );
00270 }
00271
00272
00273
00274 double ansi_gammap(double a, double x)
00275 {
00276 if ((x < 0) || (a <= 0)) {
00277 errno = EDOM;
00278 return 0;
00279 }
00280 if ( x < a+1 ) return ansi_gamma_series(a,x);
00281 else return (1 - ansi_gamma_cf(a,x));
00282 }
00283
00284
00285
00286 double ansi_erf(double arg)
00287 {
00288 return (arg<0) ? -ansi_gammap(0.5,arg*arg) : ansi_gammap(0.5,arg*arg);
00289 }
00290
00291
00292
00293 double ansi_erfc(double arg)
00294 {
00295 return (arg<0) ? 1+ansi_gammap(0.5,arg*arg) : ansi_gammap(0.5,arg*arg);
00296 }
00297
00298 #ifdef ACCOLD
00299 #define ACC ACCOLD
00300 #endif
00301
00302 #ifdef BIGNOOLD
00303 #define BIGNO BIGNOOLD
00304 #endif
00305
00306 #ifdef BIGNIOLD
00307 #define BIGNI BIGNIOLD
00308 #endif
00309
00310 #ifdef ITMAXOLD
00311 #define ITMAX ITMAXOLD
00312 #endif
00313
00314 #ifdef EPSOLD
00315 #define EPS EPSOLD
00316 #endif
00317
00318 #ifdef FPMINOLD
00319 #define FPMIN FPMINOLD
00320 #endif
00321
00322
00323 double factorial(double arg)
00324 {
00325 double Ans,Temp;
00326 static double a[13] = {1,1,2,6,24,120,720,5040,40320,362880,3628800,
00327 39916800,479001600};
00328 if ((arg < 0) || (floor(arg) != arg))
00329 {
00330 errno = EDOM;
00331 return 0;
00332 };
00333 if (arg > 170)
00334 {
00335 errno = ERANGE;
00336 return HUGE_VAL;
00337 };
00338 if (arg <= 12) return a[(int)arg];
00339 Ans = 1;
00340 for (Temp = 1;Temp <= arg;Temp++)
00341 {
00342 Ans = Ans*Temp;
00343 };
00344 return Ans;
00345 }
00346
00347
00348 double combination(double n,double r)
00349 {
00350 double ans;
00351 if ((n < 0) || (r < 0) || (n < r))
00352 {
00353 errno = ERANGE;
00354 return 0;
00355 };
00356 if (n > 170)
00357 {
00358 errno = EDOM;
00359 return HUGE_VAL;
00360 };
00361 ans = factorial(n);
00362 ans = ans/factorial(r);
00363 ans = ans/factorial(n-r);
00364 return ans;
00365 }
00366
00367
00368 double permutation(double n,double r)
00369 {
00370 double ans;
00371 if ((n < 0) || (r < 0) || (n < r))
00372 {
00373 errno = ERANGE;
00374 return 0;
00375 };
00376 if (n > 170)
00377 {
00378 errno = EDOM;
00379 return HUGE_VAL;
00380 };
00381 ans = factorial(n);
00382 ans = ans/factorial(r);
00383 return ans;
00384 }
00385
00386
00387 double integ(double arg)
00388 {
00389 double ans;
00390 modf(arg,&ans);
00391 return ans;
00392 }
00393
00394
00395 double frac(double arg)
00396 {
00397 double ans;
00398 return modf(arg,&ans);
00399 }
00400
00401
00402 double step(double arg)
00403 {
00404 return (arg<0) ? 0 : 1;
00405 }
00406
00407
00408 double sec(double arg)
00409 {
00410 return 1/cos(arg);
00411 }
00412
00413
00414 double cosec(double arg)
00415 {
00416 return 1/sin(arg);
00417 }
00418
00419
00420 double cot(double arg)
00421 {
00422 return 1/tan(arg);
00423 }
00424