476 lines
13 KiB
C++
476 lines
13 KiB
C++
/*
|
|
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 <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include "pimathbase.h"
|
|
#include "pitime.h"
|
|
#include <stdlib.h>
|
|
|
|
|
|
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.;
|
|
}
|