/*
PIP - Platform Independent Primitives
Basic mathematical functions and defines
Ivan Pelipenko peri4ko@yandex.ru
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see .
*/
#include "pimathbase.h"
#include "pitime.h"
#include
double piJ0(const double & v) {
#ifndef PIP_MATH_J0
double x = v;
double xsq;
double nn;
double pzero;
double qzero;
double p1;
double q1;
double result;
if (x < 0) x = -x;
if (x > 8.) {
double xsq_;
double p2;
double q2;
double p3;
double q3;
xsq_ = 64. / (x * x);
p2 = 0.0;
p2 = 2485.271928957404011288128951 + xsq_ * p2;
p2 = 153982.6532623911470917825993 + xsq_ * p2;
p2 = 2016135.283049983642487182349 + xsq_ * p2;
p2 = 8413041.456550439208464315611 + xsq_ * p2;
p2 = 12332384.76817638145232406055 + xsq_ * p2;
p2 = 5393485.083869438325262122897 + xsq_ * p2;
q2 = 1.0;
q2 = 2615.700736920839685159081813 + xsq_ * q2;
q2 = 156001.7276940030940592769933 + xsq_ * q2;
q2 = 2025066.801570134013891035236 + xsq_ * q2;
q2 = 8426449.050629797331554404810 + xsq_ * q2;
q2 = 12338310.22786324960844856182 + xsq_ * q2;
q2 = 5393485.083869438325560444960 + xsq_ * q2;
p3 = -0.0;
p3 = -4.887199395841261531199129300 + xsq_ * p3;
p3 = -226.2630641933704113967255053 + xsq_ * p3;
p3 = -2365.956170779108192723612816 + xsq_ * p3;
p3 = -8239.066313485606568803548860 + xsq_ * p3;
p3 = -10381.41698748464093880530341 + xsq_ * p3;
p3 = -3984.617357595222463506790588 + xsq_ * p3;
q3 = 1.0;
q3 = 408.7714673983499223402830260 + xsq_ * q3;
q3 = 15704.89191515395519392882766 + xsq_ * q3;
q3 = 156021.3206679291652539287109 + xsq_ * q3;
q3 = 533291.3634216897168722255057 + xsq_ * q3;
q3 = 666745.4239319826986004038103 + xsq_ * q3;
q3 = 255015.5108860942382983170882 + xsq_ * q3;
pzero = p2 / q2;
qzero = 8. * p3 / q3 / x;
nn = x - M_PI / 4.;
result = sqrt(2. / M_PI / x) * (pzero * cos(nn) - qzero * sin(nn));
return result;
}
xsq = x * x;
p1 = 26857.86856980014981415848441;
p1 = -40504123.71833132706360663322 + xsq * p1;
p1 = 25071582855.36881945555156435 + xsq * p1;
p1 = -8085222034853.793871199468171 + xsq * p1;
p1 = 1434354939140344.111664316553 + xsq * p1;
p1 = -136762035308817138.6865416609 + xsq * p1;
p1 = 6382059341072356562.289432465 + xsq * p1;
p1 = -117915762910761053603.8440800 + xsq * p1;
p1 = 493378725179413356181.6813446 + xsq * p1;
q1 = 1.;
q1 = 1363.063652328970604442810507 + xsq * q1;
q1 = 1114636.098462985378182402543 + xsq * q1;
q1 = 669998767.2982239671814028660 + xsq * q1;
q1 = 312304311494.1213172572469442 + xsq * q1;
q1 = 112775673967979.8507056031594 + xsq * q1;
q1 = 30246356167094626.98627330784 + xsq * q1;
q1 = 5428918384092285160.200195092 + xsq * q1;
q1 = 493378725179413356211.3278438 + xsq * q1;
return p1 / q1;
#else
return j0(v);
#endif
}
double piJ1(const double & v) {
#ifndef PIP_MATH_J1
double x = v;
double s;
double xsq;
double nn;
double pzero;
double qzero;
double p1;
double q1;
double result;
s = sign(x);
if (x < 0) x = -x;
if (x > 8.) {
double xsq_;
double p2;
double q2;
double p3;
double q3;
xsq_ = 64.0 / (x * x);
p2 = -1611.616644324610116477412898;
p2 = -109824.0554345934672737413139 + xsq_ * p2;
p2 = -1523529.351181137383255105722 + xsq_ * p2;
p2 = -6603373.248364939109255245434 + xsq_ * p2;
p2 = -9942246.505077641195658377899 + xsq_ * p2;
p2 = -4435757.816794127857114720794 + xsq_ * p2;
q2 = 1.0;
q2 = -1455.009440190496182453565068 + xsq_ * q2;
q2 = -107263.8599110382011903063867 + xsq_ * q2;
q2 = -1511809.506634160881644546358 + xsq_ * q2;
q2 = -6585339.479723087072826915069 + xsq_ * q2;
q2 = -9934124.389934585658967556309 + xsq_ * q2;
q2 = -4435757.816794127856828016962 + xsq_ * q2;
p3 = 35.26513384663603218592175580;
p3 = 1706.375429020768002061283546 + xsq_ * p3;
p3 = 18494.26287322386679652009819 + xsq_ * p3;
p3 = 66178.83658127083517939992166 + xsq_ * p3;
p3 = 85145.16067533570196555001171 + xsq_ * p3;
p3 = 33220.91340985722351859704442 + xsq_ * p3;
q3 = 1.0;
q3 = 863.8367769604990967475517183 + xsq_ * q3;
q3 = 37890.22974577220264142952256 + xsq_ * q3;
q3 = 400294.4358226697511708610813 + xsq_ * q3;
q3 = 1419460.669603720892855755253 + xsq_ * q3;
q3 = 1819458.042243997298924553839 + xsq_ * q3;
q3 = 708712.8194102874357377502472 + xsq_ * q3;
pzero = p2 / q2;
qzero = 8 * p3 / q3 / x;
nn = x - 3 * M_PI / 4;
result = sqrt(2 / M_PI / x) * (pzero * cos(nn) - qzero * sin(nn));
if (s < 0) result = -result;
return result;
}
xsq = sqr(x);
p1 = 2701.122710892323414856790990;
p1 = -4695753.530642995859767162166 + xsq * p1;
p1 = 3413234182.301700539091292655 + xsq * p1;
p1 = -1322983480332.126453125473247 + xsq * p1;
p1 = 290879526383477.5409737601689 + xsq * p1;
p1 = -35888175699101060.50743641413 + xsq * p1;
p1 = 2316433580634002297.931815435 + xsq * p1;
p1 = -66721065689249162980.20941484 + xsq * p1;
p1 = 581199354001606143928.050809 + xsq * p1;
q1 = 1.0;
q1 = 1606.931573481487801970916749 + xsq * q1;
q1 = 1501793.594998585505921097578 + xsq * q1;
q1 = 1013863514.358673989967045588 + xsq * q1;
q1 = 524371026216.7649715406728642 + xsq * q1;
q1 = 208166122130760.7351240184229 + xsq * q1;
q1 = 60920613989175217.46105196863 + xsq * q1;
q1 = 11857707121903209998.37113348 + xsq * q1;
q1 = 1162398708003212287858.529400 + xsq * q1;
result = s * x * p1 / q1;
return result;
#else
return j1(v);
#endif
}
double piJn(int n, const double & v) {
#ifndef PIP_MATH_JN
double x = v;
double pkm2;
double pkm1;
double pk;
double xk;
double r;
double ans;
int k;
int sg;
double result;
if (n < 0) {
n = -n;
if (n % 2 == 0)
sg = 1;
else
sg = -1;
} else
sg = 1;
if (x < 0) {
if (n % 2 != 0) sg = -sg;
x = -x;
}
if (n == 0) {
result = sg * piJ0(x);
return result;
}
if (n == 1) {
result = sg * piJ1(x);
return result;
}
if (n == 2) {
if (x == 0)
result = 0;
else
result = sg * (2.0 * piJ1(x) / x - piJ0(x));
return result;
}
if (x < 1E-16) {
result = 0;
return result;
}
k = 53;
pk = 2 * (n + k);
ans = pk;
xk = x * x;
do {
pk = pk - 2.0;
ans = pk - xk / ans;
k = k - 1;
} while (k != 0);
ans = x / ans;
pk = 1.0;
pkm1 = 1.0 / ans;
k = n - 1;
r = 2 * k;
do {
pkm2 = (pkm1 * r - pk * x) / x;
pk = pkm1;
pkm1 = pkm2;
r = r - 2.0;
k = k - 1;
} while (k != 0);
if (fabs(pk) > fabs(pkm1))
ans = piJ1(x) / pk;
else
ans = piJ0(x) / pkm1;
result = sg * ans;
return result;
#else
return jn(n, v);
#endif
}
double piY0(const double & v) {
#ifndef PIP_MATH_Y0
double x = v;
double nn;
double xsq;
double pzero;
double qzero;
double p4;
double q4;
double result;
if (x > 8.) {
double xsq_;
double p2;
double q2;
double p3;
double q3;
xsq_ = 64.0 / (x * x);
p2 = 0.0;
p2 = 2485.271928957404011288128951 + xsq_ * p2;
p2 = 153982.6532623911470917825993 + xsq_ * p2;
p2 = 2016135.283049983642487182349 + xsq_ * p2;
p2 = 8413041.456550439208464315611 + xsq_ * p2;
p2 = 12332384.76817638145232406055 + xsq_ * p2;
p2 = 5393485.083869438325262122897 + xsq_ * p2;
q2 = 1.0;
q2 = 2615.700736920839685159081813 + xsq_ * q2;
q2 = 156001.7276940030940592769933 + xsq_ * q2;
q2 = 2025066.801570134013891035236 + xsq_ * q2;
q2 = 8426449.050629797331554404810 + xsq_ * q2;
q2 = 12338310.22786324960844856182 + xsq_ * q2;
q2 = 5393485.083869438325560444960 + xsq_ * q2;
p3 = -0.0;
p3 = -4.887199395841261531199129300 + xsq_ * p3;
p3 = -226.2630641933704113967255053 + xsq_ * p3;
p3 = -2365.956170779108192723612816 + xsq_ * p3;
p3 = -8239.066313485606568803548860 + xsq_ * p3;
p3 = -10381.41698748464093880530341 + xsq_ * p3;
p3 = -3984.617357595222463506790588 + xsq_ * p3;
q3 = 1.0;
q3 = 408.7714673983499223402830260 + xsq_ * q3;
q3 = 15704.89191515395519392882766 + xsq_ * q3;
q3 = 156021.3206679291652539287109 + xsq_ * q3;
q3 = 533291.3634216897168722255057 + xsq_ * q3;
q3 = 666745.4239319826986004038103 + xsq_ * q3;
q3 = 255015.5108860942382983170882 + xsq_ * q3;
pzero = p2 / q2;
qzero = 8 * p3 / q3 / x;
nn = x - M_PI / 4;
result = sqrt(2 / M_PI / x) * (pzero * sin(nn) + qzero * cos(nn));
return result;
}
xsq = sqr(x);
p4 = -41370.35497933148554125235152;
p4 = 59152134.65686889654273830069 + xsq * p4;
p4 = -34363712229.79040378171030138 + xsq * p4;
p4 = 10255208596863.94284509167421 + xsq * p4;
p4 = -1648605817185729.473122082537 + xsq * p4;
p4 = 137562431639934407.8571335453 + xsq * p4;
p4 = -5247065581112764941.297350814 + xsq * p4;
p4 = 65874732757195549259.99402049 + xsq * p4;
p4 = -27502866786291095837.01933175 + xsq * p4;
q4 = 1.0;
q4 = 1282.452772478993804176329391 + xsq * q4;
q4 = 1001702.641288906265666651753 + xsq * q4;
q4 = 579512264.0700729537480087915 + xsq * q4;
q4 = 261306575504.1081249568482092 + xsq * q4;
q4 = 91620380340751.85262489147968 + xsq * q4;
q4 = 23928830434997818.57439356652 + xsq * q4;
q4 = 4192417043410839973.904769661 + xsq * q4;
q4 = 372645883898616588198.9980 + xsq * q4;
result = p4 / q4 + 2 / M_PI * piJ0(x) * log(x);
return result;
#else
return y0(v);
#endif
}
double piY1(const double & v) {
#ifndef PIP_MATH_Y1
double x = v;
double nn;
double xsq;
double pzero;
double qzero;
double p4;
double q4;
double result;
if (x > 8.) {
double xsq_;
double p2;
double q2;
double p3;
double q3;
xsq_ = 64.0 / (x * x);
p2 = -1611.616644324610116477412898;
p2 = -109824.0554345934672737413139 + xsq_ * p2;
p2 = -1523529.351181137383255105722 + xsq_ * p2;
p2 = -6603373.248364939109255245434 + xsq_ * p2;
p2 = -9942246.505077641195658377899 + xsq_ * p2;
p2 = -4435757.816794127857114720794 + xsq_ * p2;
q2 = 1.0;
q2 = -1455.009440190496182453565068 + xsq_ * q2;
q2 = -107263.8599110382011903063867 + xsq_ * q2;
q2 = -1511809.506634160881644546358 + xsq_ * q2;
q2 = -6585339.479723087072826915069 + xsq_ * q2;
q2 = -9934124.389934585658967556309 + xsq_ * q2;
q2 = -4435757.816794127856828016962 + xsq_ * q2;
p3 = 35.26513384663603218592175580;
p3 = 1706.375429020768002061283546 + xsq_ * p3;
p3 = 18494.26287322386679652009819 + xsq_ * p3;
p3 = 66178.83658127083517939992166 + xsq_ * p3;
p3 = 85145.16067533570196555001171 + xsq_ * p3;
p3 = 33220.91340985722351859704442 + xsq_ * p3;
q3 = 1.0;
q3 = 863.8367769604990967475517183 + xsq_ * q3;
q3 = 37890.22974577220264142952256 + xsq_ * q3;
q3 = 400294.4358226697511708610813 + xsq_ * q3;
q3 = 1419460.669603720892855755253 + xsq_ * q3;
q3 = 1819458.042243997298924553839 + xsq_ * q3;
q3 = 708712.8194102874357377502472 + xsq_ * q3;
pzero = p2 / q2;
qzero = 8 * p3 / q3 / x;
nn = x - 3 * M_PI / 4;
result = sqrt(2 / M_PI / x) * (pzero * sin(nn) + qzero * cos(nn));
return result;
}
xsq = sqr(x);
p4 = -2108847.540133123652824139923;
p4 = 3639488548.124002058278999428 + xsq * p4;
p4 = -2580681702194.450950541426399 + xsq * p4;
p4 = 956993023992168.3481121552788 + xsq * p4;
p4 = -196588746272214065.8820322248 + xsq * p4;
p4 = 21931073399177975921.11427556 + xsq * p4;
p4 = -1212297555414509577913.561535 + xsq * p4;
p4 = 26554738314348543268942.48968 + xsq * p4;
p4 = -99637534243069222259967.44354 + xsq * p4;
q4 = 1.0;
q4 = 1612.361029677000859332072312 + xsq * q4;
q4 = 1563282.754899580604737366452 + xsq * q4;
q4 = 1128686837.169442121732366891 + xsq * q4;
q4 = 646534088126.5275571961681500 + xsq * q4;
q4 = 297663212564727.6729292742282 + xsq * q4;
q4 = 108225825940881955.2553850180 + xsq * q4;
q4 = 29549879358971486742.90758119 + xsq * q4;
q4 = 5435310377188854170800.653097 + xsq * q4;
q4 = 508206736694124324531442.4152 + xsq * q4;
result = x * p4 / q4 + 2 / M_PI * (piJ1(x) * log(x) - 1 / x);
return result;
#else
return y1(v);
#endif
}
double piYn(int n, const double & v) {
#ifndef PIP_MATH_YN
int i;
double x = v;
double a;
double b;
double tmp;
double s;
double result;
s = 1;
if (n < 0) {
n = -n;
if (n % 2 != 0) s = -1;
}
if (n == 0) {
result = piY0(x);
return result;
}
if (n == 1) {
result = s * piY1(x);
return result;
}
a = piY0(x);
b = piY1(x);
for (i = 1; i <= n - 1; i++) {
tmp = b;
b = 2 * i / x * b - a;
a = tmp;
}
result = s * b;
return result;
#else
return yn(n, v);
#endif
}
double randomn(double dv, double sv) {
static bool agen = false;
double s = 2., v0 = 0., v1 = 0.;
if (agen) {
agen = false;
v1 = v1 * sqrt(-2 * log(s) / s);
return v1 * sv + dv;
}
while (s > 1. || s == 0.) {
v0 = randomd();
v1 = randomd();
s = v0 * v0 + v1 * v1;
}
v0 = v0 * sqrt(-2 * log(s) / s);
return v0 * sv + dv;
}
double randomd() {
return (double)rand() / RAND_MAX * 2. - 1.;
}