/************************************************************************* Cephes Math Library Release 2.8: June, 2000 Copyright by Stephen L. Moshier Contributors: * Sergey Bochkanov (ALGLIB project). Translation from C to pseudocode. See subroutines comments for additional copyrights. >>> SOURCE LICENSE >>> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (www.fsf.org); either version 2 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 General Public License for more details. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ #ifndef _bessel_h #define _bessel_h #include "ap.h" #include "amp.h" namespace bessel { template amp::ampf besselj0(amp::ampf x); template amp::ampf besselj1(amp::ampf x); template amp::ampf besseljn(int n, amp::ampf x); template amp::ampf bessely0(amp::ampf x); template amp::ampf bessely1(amp::ampf x); template amp::ampf besselyn(int n, amp::ampf x); template amp::ampf besseli0(amp::ampf x); template amp::ampf besseli1(amp::ampf x); template amp::ampf besselk0(amp::ampf x); template amp::ampf besselk1(amp::ampf x); template amp::ampf besselkn(int nn, amp::ampf x); template void besselmfirstcheb(amp::ampf c, amp::ampf& b0, amp::ampf& b1, amp::ampf& b2); template void besselmnextcheb(amp::ampf x, amp::ampf c, amp::ampf& b0, amp::ampf& b1, amp::ampf& b2); template void besselm1firstcheb(amp::ampf c, amp::ampf& b0, amp::ampf& b1, amp::ampf& b2); template void besselm1nextcheb(amp::ampf x, amp::ampf c, amp::ampf& b0, amp::ampf& b1, amp::ampf& b2); template void besselasympt0(amp::ampf x, amp::ampf& pzero, amp::ampf& qzero); template void besselasympt1(amp::ampf x, amp::ampf& pzero, amp::ampf& qzero); /************************************************************************* Bessel function of order zero Returns Bessel function of order zero of the argument. The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval the following rational approximation is used: 2 2 (w - r ) (w - r ) P (w) / Q (w) 1 2 3 8 2 where w = x and the two r's are zeros of the function. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. ACCURACY: Absolute error: arithmetic domain # trials peak rms IEEE 0, 30 60000 4.2e-16 1.1e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf besselj0(amp::ampf x) { amp::ampf result; amp::ampf xsq; amp::ampf nn; amp::ampf pzero; amp::ampf qzero; amp::ampf p1; amp::ampf q1; if( x<0 ) { x = -x; } if( x>amp::ampf("8.0") ) { besselasympt0(x, pzero, qzero); nn = x-amp::pi()/4; result = amp::sqrt(2/amp::pi()/x)*(pzero*amp::cos(nn)-qzero*amp::sin(nn)); return result; } xsq = amp::sqr(x); p1 = amp::ampf("26857.86856980014981415848441"); p1 = -amp::ampf("40504123.71833132706360663322")+xsq*p1; p1 = amp::ampf("25071582855.36881945555156435")+xsq*p1; p1 = -amp::ampf("8085222034853.793871199468171")+xsq*p1; p1 = amp::ampf("1434354939140344.111664316553")+xsq*p1; p1 = -amp::ampf("136762035308817138.6865416609")+xsq*p1; p1 = amp::ampf("6382059341072356562.289432465")+xsq*p1; p1 = -amp::ampf("117915762910761053603.8440800")+xsq*p1; p1 = amp::ampf("493378725179413356181.6813446")+xsq*p1; q1 = amp::ampf("1.0"); q1 = amp::ampf("1363.063652328970604442810507")+xsq*q1; q1 = amp::ampf("1114636.098462985378182402543")+xsq*q1; q1 = amp::ampf("669998767.2982239671814028660")+xsq*q1; q1 = amp::ampf("312304311494.1213172572469442")+xsq*q1; q1 = amp::ampf("112775673967979.8507056031594")+xsq*q1; q1 = amp::ampf("30246356167094626.98627330784")+xsq*q1; q1 = amp::ampf("5428918384092285160.200195092")+xsq*q1; q1 = amp::ampf("493378725179413356211.3278438")+xsq*q1; result = p1/q1; return result; } /************************************************************************* Bessel function of order one Returns Bessel function of order one of the argument. The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5. ACCURACY: Absolute error: arithmetic domain # trials peak rms IEEE 0, 30 30000 2.6e-16 1.1e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf besselj1(amp::ampf x) { amp::ampf result; amp::ampf s; amp::ampf xsq; amp::ampf nn; amp::ampf pzero; amp::ampf qzero; amp::ampf p1; amp::ampf q1; s = amp::sign(x); if( x<0 ) { x = -x; } if( x>amp::ampf("8.0") ) { besselasympt1(x, pzero, qzero); nn = x-3*amp::pi()/4; result = amp::sqrt(2/amp::pi()/x)*(pzero*amp::cos(nn)-qzero*amp::sin(nn)); if( s<0 ) { result = -result; } return result; } xsq = amp::sqr(x); p1 = amp::ampf("2701.122710892323414856790990"); p1 = -amp::ampf("4695753.530642995859767162166")+xsq*p1; p1 = amp::ampf("3413234182.301700539091292655")+xsq*p1; p1 = -amp::ampf("1322983480332.126453125473247")+xsq*p1; p1 = amp::ampf("290879526383477.5409737601689")+xsq*p1; p1 = -amp::ampf("35888175699101060.50743641413")+xsq*p1; p1 = amp::ampf("2316433580634002297.931815435")+xsq*p1; p1 = -amp::ampf("66721065689249162980.20941484")+xsq*p1; p1 = amp::ampf("581199354001606143928.050809")+xsq*p1; q1 = amp::ampf("1.0"); q1 = amp::ampf("1606.931573481487801970916749")+xsq*q1; q1 = amp::ampf("1501793.594998585505921097578")+xsq*q1; q1 = amp::ampf("1013863514.358673989967045588")+xsq*q1; q1 = amp::ampf("524371026216.7649715406728642")+xsq*q1; q1 = amp::ampf("208166122130760.7351240184229")+xsq*q1; q1 = amp::ampf("60920613989175217.46105196863")+xsq*q1; q1 = amp::ampf("11857707121903209998.37113348")+xsq*q1; q1 = amp::ampf("1162398708003212287858.529400")+xsq*q1; result = s*x*p1/q1; return result; } /************************************************************************* Bessel function of integer order Returns Bessel function of order n, where n is a (possibly negative) integer. The ratio of jn(x) to j0(x) is computed by backward recurrence. First the ratio jn/jn-1 is found by a continued fraction expansion. Then the recurrence relating successive orders is applied until j0 or j1 is reached. If n = 0 or 1 the routine for j0 or j1 is called directly. ACCURACY: Absolute error: arithmetic range # trials peak rms IEEE 0, 30 5000 4.4e-16 7.9e-17 Not suitable for large n or x. Use jv() (fractional order) instead. Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf besseljn(int n, amp::ampf x) { amp::ampf result; amp::ampf pkm2; amp::ampf pkm1; amp::ampf pk; amp::ampf xk; amp::ampf r; amp::ampf ans; int k; int sg; 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*besselj0(x); return result; } if( n==1 ) { result = sg*besselj1(x); return result; } if( n==2 ) { if( x==0 ) { result = 0; } else { result = sg*(amp::ampf("2.0")*besselj1(x)/x-besselj0(x)); } return result; } if( x::getAlgoPascalEpsilon() ) { result = 0; return result; } k = 53; pk = 2*(n+k); ans = pk; xk = x*x; do { pk = pk-amp::ampf("2.0"); ans = pk-xk/ans; k = k-1; } while( k!=0 ); ans = x/ans; pk = amp::ampf("1.0"); pkm1 = amp::ampf("1.0")/ans; k = n-1; r = 2*k; do { pkm2 = (pkm1*r-pk*x)/x; pk = pkm1; pkm1 = pkm2; r = r-amp::ampf("2.0"); k = k-1; } while( k!=0 ); if( amp::abs(pk)>amp::abs(pkm1) ) { ans = besselj1(x)/pk; } else { ans = besselj0(x)/pkm1; } result = sg*ans; return result; } /************************************************************************* Bessel function of the second kind, order zero Returns Bessel function of the second kind, of order zero, of the argument. The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval a rational approximation R(x) is employed to compute y0(x) = R(x) + 2 * log(x) * j0(x) / PI. Thus a call to j0() is required. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. ACCURACY: Absolute error, when y0(x) < 1; else relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.3e-15 1.6e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf bessely0(amp::ampf x) { amp::ampf result; amp::ampf nn; amp::ampf xsq; amp::ampf pzero; amp::ampf qzero; amp::ampf p4; amp::ampf q4; if( x>amp::ampf("8.0") ) { besselasympt0(x, pzero, qzero); nn = x-amp::pi()/4; result = amp::sqrt(2/amp::pi()/x)*(pzero*amp::sin(nn)+qzero*amp::cos(nn)); return result; } xsq = amp::sqr(x); p4 = -amp::ampf("41370.35497933148554125235152"); p4 = amp::ampf("59152134.65686889654273830069")+xsq*p4; p4 = -amp::ampf("34363712229.79040378171030138")+xsq*p4; p4 = amp::ampf("10255208596863.94284509167421")+xsq*p4; p4 = -amp::ampf("1648605817185729.473122082537")+xsq*p4; p4 = amp::ampf("137562431639934407.8571335453")+xsq*p4; p4 = -amp::ampf("5247065581112764941.297350814")+xsq*p4; p4 = amp::ampf("65874732757195549259.99402049")+xsq*p4; p4 = -amp::ampf("27502866786291095837.01933175")+xsq*p4; q4 = amp::ampf("1.0"); q4 = amp::ampf("1282.452772478993804176329391")+xsq*q4; q4 = amp::ampf("1001702.641288906265666651753")+xsq*q4; q4 = amp::ampf("579512264.0700729537480087915")+xsq*q4; q4 = amp::ampf("261306575504.1081249568482092")+xsq*q4; q4 = amp::ampf("91620380340751.85262489147968")+xsq*q4; q4 = amp::ampf("23928830434997818.57439356652")+xsq*q4; q4 = amp::ampf("4192417043410839973.904769661")+xsq*q4; q4 = amp::ampf("372645883898616588198.9980")+xsq*q4; result = p4/q4+2/amp::pi()*besselj0(x)*amp::log(x); return result; } /************************************************************************* Bessel function of second kind of order one Returns Bessel function of the second kind of order one of the argument. The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 25 term Chebyshev expansion is used, and a call to j1() is required. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5. ACCURACY: Absolute error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.0e-15 1.3e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf bessely1(amp::ampf x) { amp::ampf result; amp::ampf nn; amp::ampf xsq; amp::ampf pzero; amp::ampf qzero; amp::ampf p4; amp::ampf q4; if( x>amp::ampf("8.0") ) { besselasympt1(x, pzero, qzero); nn = x-3*amp::pi()/4; result = amp::sqrt(2/amp::pi()/x)*(pzero*amp::sin(nn)+qzero*amp::cos(nn)); return result; } xsq = amp::sqr(x); p4 = -amp::ampf("2108847.540133123652824139923"); p4 = amp::ampf("3639488548.124002058278999428")+xsq*p4; p4 = -amp::ampf("2580681702194.450950541426399")+xsq*p4; p4 = amp::ampf("956993023992168.3481121552788")+xsq*p4; p4 = -amp::ampf("196588746272214065.8820322248")+xsq*p4; p4 = amp::ampf("21931073399177975921.11427556")+xsq*p4; p4 = -amp::ampf("1212297555414509577913.561535")+xsq*p4; p4 = amp::ampf("26554738314348543268942.48968")+xsq*p4; p4 = -amp::ampf("99637534243069222259967.44354")+xsq*p4; q4 = amp::ampf("1.0"); q4 = amp::ampf("1612.361029677000859332072312")+xsq*q4; q4 = amp::ampf("1563282.754899580604737366452")+xsq*q4; q4 = amp::ampf("1128686837.169442121732366891")+xsq*q4; q4 = amp::ampf("646534088126.5275571961681500")+xsq*q4; q4 = amp::ampf("297663212564727.6729292742282")+xsq*q4; q4 = amp::ampf("108225825940881955.2553850180")+xsq*q4; q4 = amp::ampf("29549879358971486742.90758119")+xsq*q4; q4 = amp::ampf("5435310377188854170800.653097")+xsq*q4; q4 = amp::ampf("508206736694124324531442.4152")+xsq*q4; result = x*p4/q4+2/amp::pi()*(besselj1(x)*amp::log(x)-1/x); return result; } /************************************************************************* Bessel function of second kind of integer order Returns Bessel function of order n, where n is a (possibly negative) integer. The function is evaluated by forward recurrence on n, starting with values computed by the routines y0() and y1(). If n = 0 or 1 the routine for y0 or y1 is called directly. ACCURACY: Absolute error, except relative when y > 1: arithmetic domain # trials peak rms IEEE 0, 30 30000 3.4e-15 4.3e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf besselyn(int n, amp::ampf x) { amp::ampf result; int i; amp::ampf a; amp::ampf b; amp::ampf tmp; amp::ampf s; s = 1; if( n<0 ) { n = -n; if( n%2!=0 ) { s = -1; } } if( n==0 ) { result = bessely0(x); return result; } if( n==1 ) { result = s*bessely1(x); return result; } a = bessely0(x); b = bessely1(x); for(i=1; i<=n-1; i++) { tmp = b; b = 2*i/x*b-a; a = tmp; } result = s*b; return result; } /************************************************************************* Modified Bessel function of order zero Returns modified Bessel function of order zero of the argument. The function is defined as i0(x) = j0( ix ). The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,30 30000 5.8e-16 1.4e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1984, 1987, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf besseli0(amp::ampf x) { amp::ampf result; amp::ampf y; amp::ampf v; amp::ampf z; amp::ampf b0; amp::ampf b1; amp::ampf b2; if( x<0 ) { x = -x; } if( x<=amp::ampf("8.0") ) { y = x/amp::ampf("2.0")-amp::ampf("2.0"); besselmfirstcheb(-amp::ampf("4.41534164647933937950E-18"), b0, b1, b2); besselmnextcheb(y, amp::ampf("3.33079451882223809783E-17"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("2.43127984654795469359E-16"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.71539128555513303061E-15"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("1.16853328779934516808E-14"), b0, b1, b2); besselmnextcheb(y, amp::ampf("7.67618549860493561688E-14"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("4.85644678311192946090E-13"), b0, b1, b2); besselmnextcheb(y, amp::ampf("2.95505266312963983461E-12"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("1.72682629144155570723E-11"), b0, b1, b2); besselmnextcheb(y, amp::ampf("9.67580903537323691224E-11"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("5.18979560163526290666E-10"), b0, b1, b2); besselmnextcheb(y, amp::ampf("2.65982372468238665035E-9"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("1.30002500998624804212E-8"), b0, b1, b2); besselmnextcheb(y, amp::ampf("6.04699502254191894932E-8"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("2.67079385394061173391E-7"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.11738753912010371815E-6"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("4.41673835845875056359E-6"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.64484480707288970893E-5"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("5.75419501008210370398E-5"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.88502885095841655729E-4"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("5.76375574538582365885E-4"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.63947561694133579842E-3"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("4.32430999505057594430E-3"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.05464603945949983183E-2"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("2.37374148058994688156E-2"), b0, b1, b2); besselmnextcheb(y, amp::ampf("4.93052842396707084878E-2"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("9.49010970480476444210E-2"), b0, b1, b2); besselmnextcheb(y, amp::ampf("1.71620901522208775349E-1"), b0, b1, b2); besselmnextcheb(y, -amp::ampf("3.04682672343198398683E-1"), b0, b1, b2); besselmnextcheb(y, amp::ampf("6.76795274409476084995E-1"), b0, b1, b2); v = amp::ampf("0.5")*(b0-b2); result = amp::exp(x)*v; return result; } z = amp::ampf("32.0")/x-amp::ampf("2.0"); besselmfirstcheb(-amp::ampf("7.23318048787475395456E-18"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("4.83050448594418207126E-18"), b0, b1, b2); besselmnextcheb(z, amp::ampf("4.46562142029675999901E-17"), b0, b1, b2); besselmnextcheb(z, amp::ampf("3.46122286769746109310E-17"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("2.82762398051658348494E-16"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("3.42548561967721913462E-16"), b0, b1, b2); besselmnextcheb(z, amp::ampf("1.77256013305652638360E-15"), b0, b1, b2); besselmnextcheb(z, amp::ampf("3.81168066935262242075E-15"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("9.55484669882830764870E-15"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("4.15056934728722208663E-14"), b0, b1, b2); besselmnextcheb(z, amp::ampf("1.54008621752140982691E-14"), b0, b1, b2); besselmnextcheb(z, amp::ampf("3.85277838274214270114E-13"), b0, b1, b2); besselmnextcheb(z, amp::ampf("7.18012445138366623367E-13"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("1.79417853150680611778E-12"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("1.32158118404477131188E-11"), b0, b1, b2); besselmnextcheb(z, -amp::ampf("3.14991652796324136454E-11"), b0, b1, b2); besselmnextcheb(z, amp::ampf("1.18891471078464383424E-11"), b0, b1, b2); besselmnextcheb(z, amp::ampf("4.94060238822496958910E-10"), b0, b1, b2); besselmnextcheb(z, amp::ampf("3.39623202570838634515E-9"), b0, b1, b2); besselmnextcheb(z, amp::ampf("2.26666899049817806459E-8"), b0, b1, b2); besselmnextcheb(z, amp::ampf("2.04891858946906374183E-7"), b0, b1, b2); besselmnextcheb(z, amp::ampf("2.89137052083475648297E-6"), b0, b1, b2); besselmnextcheb(z, amp::ampf("6.88975834691682398426E-5"), b0, b1, b2); besselmnextcheb(z, amp::ampf("3.36911647825569408990E-3"), b0, b1, b2); besselmnextcheb(z, amp::ampf("8.04490411014108831608E-1"), b0, b1, b2); v = amp::ampf("0.5")*(b0-b2); result = amp::exp(x)*v/amp::sqrt(x); return result; } /************************************************************************* Modified Bessel function of order one Returns modified Bessel function of order one of the argument. The function is defined as i1(x) = -i j1( ix ). The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.9e-15 2.1e-16 Cephes Math Library Release 2.8: June, 2000 Copyright 1985, 1987, 2000 by Stephen L. Moshier *************************************************************************/ template amp::ampf besseli1(amp::ampf x) { amp::ampf result; amp::ampf y; amp::ampf z; amp::ampf v; amp::ampf b0; amp::ampf b1; amp::ampf b2; z = amp::abs(x); if( z<=amp::ampf("8.0") ) { y = z/amp::ampf("2.0")-amp::ampf("2.0"); besselm1firstcheb(amp::ampf("2.77791411276104639959E-18"), b0, b1, b2); besselm1nextcheb(y, -amp::ampf("2.11142121435816608115E-17"), b0, b1, b2); besselm1nextcheb