/* Compute bessel function of first kind * Copyright (C) 2006 Bastien ROUCARIES * * COPYRIGHT NOTICE (README FIRST) * ******************************** * * The copyright of this file is quite COMPLICATED: * THIS FILE is put by me bastien ROUCARIES under the * LGPL licence * BUT because I use bessel_I0 and bessel_I1 * from gsl this file becomes de facto GPL * * I would prefere to use a pure LGPL license, but I lack time * to rewrite all the GPL function. * If you think GPL is too strong, * you MUST rewrite the following function and this file will * become LGPL: * * * * COPYRIGHT FOR THIS FILE ALONE * ***************************** * * This library 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 2.1 of the License, or (at your option) any later version. * * This library 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 library; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * * COPYRIGHT FOR EXECUTE AND LINKING THIS FILE (DO SOMETHING USEFUL) * ***************************************************************** * * 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; 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. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * */ /* Please send comments to roucaries.bastien+mathsoft at gmail.com * * ABOUT SPAM and SPAM BOT * If you want to send me spam (you don't and by UE law forbid you to do so) * please send me a mail to my kill list address * address@hidden * END SPAM and SPAM BOT section */ /*!\file cbesseljn.c \brief compute complex bessel J function References: [1] Bessel function of the first kind with complex argument Yousif, Hashim A.; Melka, Richard Computer Physics Communications, vol. 106, Issue 3, pp.199-206 11/1997, ELSEVIER, doi://10.1016/S0010-4655(97)00087-8 [2] Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables Milton Abramowitz and Irene A. Stegun public domain (work of US government) online http://www.math.sfu.ca/~cbm/aands/ [3] Mathematica Manual Bessel, Airy, Struve Functions> BesselJ[nu,z] > General characteristics> Symmetries and periodicities http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/04/02/01/ [4] Mathematica Manual Bessel, Airy, Struve Functions> BesselJ[nu,z] Representations through equivalent functions http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/27/ShowAll.html [5] Algorithms for the evaluation of Bessel functions of complex argument and integer orders G. D. C. Kuiken Applied Mathematics Letters, Volume 2, Issue 4, 1989, Pages 353-356 doi://10.1016/0893-9659(89)90086-4 [6] The numerical computation of Bessel functions of the first and second kind for integer orders and complex arguments du Toit, C.F. Antennas and Propagation, IEEE Transactions on, Volume 38, Issue 9 Sep 1990 pages 1341-1349 doi://10.1109/8.56985 [7] "Bessel Function of the First Kind." Weisstein, Eric W. From MathWorld--A Wolfram Web Resource. online http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html [8] Kahan summation algorithm Wikipedia http://en.wikipedia.org/wiki/Kahan_summation_algorithm [9] Integer sequence A000142 Factorial numbers: n! = 1*2*3*4*...*n http://www.research.att.com/~njas/sequences/table?a=142&fmt=4 [10] N. M. Temme Numerical algorithms for uniform {Airy-type} asymptotic expansions Numerical Algorithms Volume 15, Issue 2 pages 207-225 1997 online http://citeseer.ist.psu.edu/temme97numerical.html */ /* we use special algorithm what do not like optimization */ //#pragma STDC FP_CONTRACT OFF //#pragma STDC CX_LIMITED_RANGE OFF #include #include #include #include #include #include #include #include #include #include #include "debug.h" //#define NDEBUG #define cprint(x) creal(x), cimag(x) #define fmtc "%+4.4e%+4.4e*j" #define forget(a) do {} while(0) // __asm__ __volatile__ ("" :"=m" (a) :"m" (a)) /*!\brief Force a compilation error if condition is true */ #define BUILD_BUG_ON(condition) ((void)sizeof(char[1 - 2*!!(condition)])) #define max(x,y) x < y ? y : x #define sqr(x) ((x) * (x)) #define norm(z) (sqr(creal(z))+sqr(cimag(z))) #define SQR_DBL_EPSILON sqr(DBL_EPSILON) //#define norm(z) (cabs(z)) //#define SQR_DBL_EPSILON (DBL_EPSILON) #define ARRAY_SIZE(x) (sizeof(x)/sizeof(x[0])) #define SMALL_J0_BOUND 1e6 /*!\brief use ascending serie below this magnitude */ #define SMALL_JN_BOUND 5.0 /*!\brief use assymptotic expression above this magnitude */ #define BIG_JN_BOUND 25.0 /*!\brief Tolerance for assert bound checking */ #define TOL_JN_BOUND 1e-2 /*!\brief Arbitrary value for iteration*/ #define MAX_SMALL_ITERATION 1024 /*!\brief Max order for direct approximation */ #define MAX_DIRECT_ORDER 8 /*!\brief Kahan algorithm for summation For better understanding see [8] \todo move to header */ static inline complex kahanaddcomplex(complex R, complex Rk, complex *correction) { complex corrected_next_term; complex new_sum; corrected_next_term = Rk - *(correction); forget(corrected_next_term); new_sum = R + corrected_next_term; forget(new_sum); *correction = (new_sum - R) - corrected_next_term; forget(*correction); return new_sum; } /*\brief The first factorial (\f$n!\f$) 0 to MAXORDER */ static const unsigned long int factorialarray[] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320 }; #ifdef KAHAN_SMALL_ARGUMENT #define declaresmallcorrection(a) complex a = 0.0 #define addsmallRandRk(a, b, c) kahanaddcomplex((a),(b),(c)) #else #define declaresmallcorrection(a) do {} while(0) /*!\brief Return a + b (corrected by c if using Kahan correction) */ #define addsmallRandRk(a, b, c) ((a) + (b)) #endif /*!\brief Implement bessel J function for small arguments according to [5] For small argument we use the following formulae: \f[ J_n(z)=\sum_0^\infty R_k \f] Where: \f{align*} R_{-1}&=1 & R_k &= a_k R_{k-1} \\ a_0&=\frac{\left(\frac{1}{2}z\right)^n}{n!} & a_{+k}&=\frac{-\frac{1}{4} z^2}{k(n+k)} \f} This formula is derivated of [2] eq 9.1.10 p 360: \f[ J_n(z)=(\frac{1}{2}z)^n \sum_{k=0}^\infty \frac{(-\frac{1}{4}z^2)^k}{k!\Gamma(n+k+1)} \f] \note Precision is evaluated to be more than 1e-13 \param[in] n order \param[in] z input argument \note could implement kahan sum if KAHAN_SMALL_ARGUMENT=1 \todo Check if Kahan summation improve precision */ static complex cbesseljn_smallarg (unsigned int n, complex z) { complex ak, Rk; complex R; complex R0; unsigned short int k; /* for 4 * k * (n + k) */ unsigned long diviser; /* For Kahan sum */ declaresmallcorrection(correction); assert(cabs(z) < SMALL_JN_BOUND * ( 1 + TOL_JN_BOUND )); assert(creal(z) >= 0 && cimag(z) >= 0); assert(n <= MAX_DIRECT_ORDER); /* Do not overflow diviser = 4 * k * (k + n) */ BUILD_BUG_ON(4 * MAX_SMALL_ITERATION * (MAX_DIRECT_ORDER + MAX_SMALL_ITERATION) >= ULONG_MAX ); BUILD_BUG_ON(ARRAY_SIZE(factorialarray) != MAX_DIRECT_ORDER + 1); /* a_0 */ errno = 0; ak = cpow (0.5 * z, n); /* underflow NB: errno and exception is propagated from cpow */ if (errno == ERANGE) goto underflow; if (norm(ak) < sqr(DBL_MIN * factorialarray[n])) { /* do the computation if we are lucky we get a subnormal */ ak = ak / factorialarray[n]; goto fullunderflow; } ak = ak / factorialarray[n]; /* R_0 */ R0 = ak * 1.0; Rk = R0; R = R0; DEBUG ("a%-4i: " fmtc " R%-4i: " fmtc " R: " fmtc "\n", 0, cprint (ak), 0, cprint (Rk), cprint (R)); for (k = 1; k < MAX_SMALL_ITERATION; k++) { assert (k < MAX_SMALL_ITERATION - 1); diviser = (4 * k * (n + k)); ak = - sqr(z) / diviser; Rk = ak * Rk; DEBUG("a%-4i: " fmtc " R%-4i: " fmtc " R: " fmtc "\n", k, cprint (ak), k, cprint (Rk), cprint (R)); if (norm(Rk) < norm(R) * SQR_DBL_EPSILON/2) return R; /* R += Rk */ R = addsmallRandRk(R, Rk, &correction); } /* not converged */ errno = ERANGE; (void) feraiseexcept(FE_INVALID); return R; fullunderflow: DEBUG ("Underflow\n"); errno = ERANGE; (void) feraiseexcept(FE_UNDERFLOW); return ak; underflow: DEBUG ("Underflow\n"); return n > 0 ? 0.0 : 1; } #ifdef KAHAN_MEDIUM_ARGUMEMENT #define addmediumsumandak(a, b, c) kahanaddcomplex((a),(b),(c)) #define declaremediumcorrection(a) complex a = 0.0 #else /*!\brief Return a + b (corrected by c if using Kahan correction) */ #define addmediumsumandak(a, b, c) ((a) + (b)) #define declaremediumcorrection(a) do {} while(0) #endif /*!\brief Compute besselj for complex moderate argument and integer even order using poisson expansion According to [5] eq (7) for moderate argument and even order \f[ J_n(z) = \frac{1+(-1)^\frac{n}{2}\cos z}{2m} + \frac{1}{m} \sum_{k=1}^{m-1} \cos (z sin t) cos (nt) \f] Where \f$t\f$ is defined by: \f[ t=\frac{k\pi}{2m} \f] And \f$m\f$ is an integer such as: \f[ m \ge 2 * |z| + \frac{1}{4} (n + |\Im\text{m}\; z|) \f] In order to improve precision tranform \f$1+(-1)^\frac{n}{2}\cos z\f$ to: \f[ \begin{cases} 2\cos^2\frac{z}{2} & \text{if } \frac{n}{2} \text{ even} \\ 2\sin^2\frac{z}{2} & \text{else} \end{cases} \f] Therefore we could transform the evaluated expression to: \f[ \begin{cases} \frac{1}{m}\left[ \cos^2\frac{z}{2} + \sum_{k=1}^{m-1} \cos (z sin t) cos (nt) \right] & \text{if } \frac{n}{2} \text{ even} \\ \frac{1}{m}\left[ \sin^2\frac{z}{2} + \sum_{k=1}^{m-1} \cos (z sin t) cos (nt) \right] & \text{else} \end{cases} \f] \note Precision is evaluated to be better than 1e-13 \note use even m for better precision \note USe 1/2 instead of 1/4 before n (better precision) */ static complex cbesseljn_mediumarg_even (unsigned int n, complex z) { complex sum; complex ak; unsigned int m; unsigned int k; double t; /* For Kahan sum */ declaremediumcorrection(correction); assert(n % 2 == 0); BUILD_BUG_ON(2 * BIG_JN_BOUND + 0.5 * ( MAX_DIRECT_ORDER + BIG_JN_BOUND) >= UINT_MAX); /* use conservative measurement (m even give better error) */ m = ceil(2 * cabs (z) + 0.5 * (n + fabs (cimag (z)))); m = m % 2 == 0 ? m : m + 1; DEBUG("m=%li", (long) m); /* compute (1 + (-1)^n cos(z)) / 2 */ if(n % 4 == 0) /* Use 2 cos^2 (z/2) = (1 + cos z) */ sum = sqr(ccos(z/2)); else /* use 2 sin^2 (z/2) = 1 - cos(z) */ sum = sqr(csin(z/2)); DEBUG("a%-3i=" fmtc " sum=" fmtc "\n", 0, cprint(sum), cprint(sum)); for (k = 1; k <= m - 1; k++) { t = (k * M_PI_2) / m; ak = ccos (z * sin (t)) * cos (n * t); /* R += ak */ sum = addmediumsumandak(sum, ak, &correction); DEBUG("a%-3i=" fmtc " sum=" fmtc "\n", k, cprint(ak), cprint(sum)); } /* first+= second */ return sum / m; } /*!\brief Compute besselj for complex moderate argument and integer odd order using poisson expansion According to [5] eq (7) for moderate argument and odd order \f[ J_n(z) = \frac{(-1)^\frac{n-1}{2}\sin z}{2m} + \frac{1}{m} \sum_{k=1}^{m-1} \sin (z sin t) sin (nt) \f] Where \f$t\f$ is defined by: \f[ t=\frac{k\pi}{2m} \f] And \f$m\f$ is an integer such as: \f[ m \ge 2 |z| + \frac{1}{4} (n + |\Im\text{m}\; z|) \f] \note Precision is evaluated to be better than 1e-13 \note use even m for better precision \note USe 1/2 instead of 1/4 before n (better precision) */ static complex cbesseljn_mediumarg_odd (unsigned int n, complex z) { complex sum; complex ak; unsigned int m; unsigned int k; double t; /* For Kahan sum */ declaremediumcorrection(correction); assert(n % 2 == 1); BUILD_BUG_ON(2 * BIG_JN_BOUND + 0.5 * ( MAX_DIRECT_ORDER + BIG_JN_BOUND) >= UINT_MAX); /* use conservative measurement (m even give better error) */ m = ceil(2 * cabs (z) + 0.5 * (n + fabs (cimag (z)))); m = m % 2 == 0 ? m : m + 1; sum = csin (z) / (2); /* -1^((n-1)/2) */ sum = (n - 1) % 4 == 0 ? sum : - sum; DEBUG("a%-3i=" fmtc " sum=" fmtc "\n", 0, cprint(sum), cprint(sum)); for (k = 1; k <= m - 1; k++) { t = (k * M_PI_2) / m; ak = csin (z * sin (t)) * sin (n * t); /* sum += ak*/ sum = addmediumsumandak(sum, ak, &correction); DEBUG("a%-3i=" fmtc " sum=" fmtc "\n", 0, cprint(sum), cprint(sum)); } return sum / m; } /*!\brief Compute besselj for complex moderate argument and integer order using poisson expansion */ static complex cbesseljn_mediumarg (unsigned int n, complex z) { assert(cabs(z) > SMALL_JN_BOUND * ( 1 - TOL_JN_BOUND )); assert(cabs(z) < BIG_JN_BOUND * ( 1 + TOL_JN_BOUND )); assert(creal(z) >= 0 && cimag(z) >= 0); assert(n <= MAX_DIRECT_ORDER); if (n % 2 == 0) return cbesseljn_mediumarg_even (n, z); else return cbesseljn_mediumarg_odd (n, z); } /*! \brief max iteration avoid overflow*/ #define MAX_LARGE_ITERATION 100 #ifdef KAHAN_LARGE_ARGUMEMENT #define addlargesumandak(a, b, c) kahanaddcomplex((a),(b),(c)) #define declarelargecorrection(a) complex a = 0.0 #else /*!\brief Return a + b (corrected by c if using Kahan correction) */ #define addlargesumandak(a, b, c) ((a) + (b)) #define declarelargecorrection(a) do {} while(0) #endif /*!\brief besselj for large argument According to [6] eq (23), (25), (26) based on [2] eq (9.2.5), (9.2.9), (9.2.10) p 364, for \f$|z|\rightarrow\infty\f$ and \f$\operatorname{arg} z < \pi\f$: \f[ J_n(v)=\sqrt{\frac{2}{\pi z}} \left[ P_n cos \left(z - \frac{n\pi}{2} - \frac{\pi}{4}\right) -Q_n sin \left(z - \frac{n\pi}{2} - \frac{\pi}{4}\right) \right] \f] Where: \f{align*} P_n(z)&= 1 + \sum_{k=1}^\infty (-1)^k \dfrac{ \prod_{r=1}^{2k}\displaylimits \left[4n^2 - (2r - 1) \right]} {(2k)!(8z)^{2k}} \\ Q_n(z)&= \sum_{k=0}^\infty (-1)^k \dfrac{ \prod\limits_{r=1}^{2k+1} \left[4n^2 - (2r - 1) \right]} {(2k+1)!(8z)^{2k+1}} \\ \f} First for n integer we sould note that: \f{align*} \cos \left(z - \frac{n\pi}{2} - \frac{\pi}{4}\right)&= \frac{1}{\sqrt{2}}\cos \left(z - \frac{n\pi}{2}\right) + \frac{1}{\sqrt{2}}\sin \left(z - \frac{n\pi}{2}\right) \\ \sin \left(z - \frac{n\pi}{2} - \frac{\pi}{4}\right)&= \frac{1}{\sqrt{2}}\sin \left(z - \frac{n\pi}{2}\right) - \frac{1}{\sqrt{2}}\cos \left(z - \frac{n\pi}{2}\right) \f} Let \f$l,m\f$ such as: \f{align*} l &= \begin{cases} +1 & \text{if } n\mod 4 \equiv 0,3 \\ -1 & \text{ else} \end{cases}\\ m &= \begin{cases} +1 & \text{if } n\mod 4 \equiv 0,1 \\ -1 & \text{ else} \end{cases} \f} Therefore we find eq (5) of [5]: \f[ \boxed{ J_n(z)= \frac{1}{\sqrt{\pi z}} \left[ (lP_n + mQ_n) \cos z + (mP_n - lQ_n) \sin z \right] } \f] Moreover we could rewrite \f$P_n\f$ using recursivity (see also [5] p 354): \f{align*} P_n(z)&=\sum_{k=0}^\infty P'_k & P'_0&=1 & P'_k&=P'_{k-1} \frac{\left[4n^2 + (4k - 3)^2\right] \left[4n^2 + (4k - 1)^2\right]} {2k(2k-1)} \frac{1}{-64z^2} \f} And for \f$Q_n\f$ (see also [5] p 354) \f{align*} Q_n(z)&=\sum_{k=0}^\infty P'_k & Q'_0&=\frac{4n^2-1}{8z} & Q'_k&=Q'_{k-1} \frac{\left[4n^2 + (4k - 1)^2\right] \left[4n^2 + (4k + 1)^2\right]} {2k(2k+1)} \frac{1}{-64z^2} \f} \note assert \f$|z| > 2 z^2 \f$ \note [5] as a typo, \f$8z^2\f$ should be read \f$(8z)^2\f$ \note valid only if \f$\Re\text{e}\; z > 0 \f$ (asserted) see [6] \note Please note that the serie is not convergent \note According to remark on page 364 afer (9.2.10) of [6] the sum of remainder terms in \f$P\f$ and \f$Q\f$ after rank \f$k\f$ is lesser than respectively \f$Pk\f$ and \f$Qk\f$ if \f$k >\frac{1}{2}n - \frac{1}{4}\f$ for \f$P\f$ and\f$k >\frac{1}{2}n - \frac{3}{4}\f$ for \f$Q\f$. \note According to eq (56) reference [7] this expansion is valid for \f$|z| \gg |m^2 - \frac{1}{4}|\f$. We assert \f$|z| > 10 |m^2 - \frac{1}{4}|\f$ */ static complex cbesseljn_largearg (unsigned int n, complex z) { complex Pk, Pkm1, P, Qk, Qkm1, Q; complex tmp; double denum; double num; complex m8z2; unsigned int k; double l, m; double mu; declarelargecorrection(correctionP); declarelargecorrection(correctionQ); // for debuging using remark of [6] double esterr; bool trust; assert(creal(z) >= 0); assert(cabs(z) >= BIG_JN_BOUND * ( 1 - TOL_JN_BOUND )); assert(norm(z) >= norm(sqr(n) - 1/4) || n <= 8); assert(creal(z) >= 0 && cimag(z) >= 0); /* overflow check */ assert(4.0 * sqr((double) n) < UINT_MAX); mu = 4 * sqr(n); /* P0 & Q0 */ Pk = 1; P = Pk; Pkm1 = Pk; Qk = (mu - 1.0) / (8.0 * z); Q = Qk; Qkm1 = Qk; m8z2 = (-64.0 * sqr (z)); /* P */ for (k = 1; k < MAX_LARGE_ITERATION; k++) { /* Usually converge before overflow */ assert (k <= MAX_LARGE_ITERATION - 1); num = (mu - sqr (4 * k - 3)) * (mu - sqr(4 * k - 1)); denum = 2 * k * (2 * k - 1); Pk = (Pk * num) / (denum * m8z2); DEBUG("num=%e denum=%e" " P=" fmtc " Pk=" fmtc "\n", num, denum , cprint(P), cprint(Pk)); if (norm(Pk) > norm(Pkm1) || norm(Pk) < norm(P) * SQR_DBL_EPSILON) break; Pkm1 = Pk; P = addlargesumandak(P, Pk, &correctionP); } trust = k - 1 > 0.5 * n - 0.25; /* Q */ for (k = 1;k < MAX_LARGE_ITERATION; k++) { /* Usually converge before overflow */ assert (k <= MAX_LARGE_ITERATION - 1); num = (mu - sqr (4 * k - 1)) * (mu - sqr(4 * k + 1)); denum = 2 * k * (2 * k + 1); Qk = (Qk * num) / (denum * m8z2); DEBUG("num=%e denum=%e" " Q=" fmtc " Qk=" fmtc "\n", num, denum , cprint(Q), cprint(Qk)); if (norm(Qk) > norm(Qkm1) || norm(Qk) < norm(Q) * SQR_DBL_EPSILON) break; Qkm1 = Qk; Q = addlargesumandak(Q, Qk, &correctionQ); } trust &= k - 1 > 0.5 * n - 0.75; /* l, m cf [5] eq (5) */ l = (n % 4 == 0) || (n % 4 == 3) ? 1.0 : -1.0; m = (n % 4 == 0) || (n % 4 == 1) ? 1.0 : -1.0; esterr = (cabs(l * Pkm1* ccos(z)) + cabs(m * Qkm1 * ccos (z)) + cabs(m * Pkm1*csin(z)) + cabs(l * Qkm1 * csin(z))) / csqrt (M_PI * z); tmp = (l * P + m * Q) * ccos (z); tmp += (m * P - l * Q) * csin (z); tmp /= csqrt(M_PI *z); DEBUG("max error theoric %e (rel %e) trust it: %s\n", esterr, esterr / cabs(tmp), trust == true ? "true" : "false"); return tmp; } #if 0 /*!\brief Compute \f$\zeta^{3}{2}\f$ for \f$|z| < 1\f$ According to [10] (2.3) p3 and [2] (9.3.38) p368 \f[ \frac{3}{2} \zeta^\frac{3}{2} = \ln \frac{1+\sqrt{1-z^2}}{z} - \sqrt{1-z^2} \f] Therefore: \f[ \zeta^\frac{3}{2} = \frac{2}{3} \left[ \ln \frac{1+\sqrt{1-z^2}}{z} - \sqrt{1-z^2}\right] \f] */ static complex zeta32lesser1(complex z) { complex sqrt1mz2; complex ret; sqrt1mz2 = csqrt(1 - sqr(z)); ret = clog(1 + sqrt1mz2) / z - sqrt1mz2; return 2/3 * ret; } /*!\brief Compute \f$\zeta^\frac{3}{2}\f$ for \f$|z| > 1\f$ According to [10] (2.3) p3 and [2] (9.3.39) p368 \f[ \frac{3}{2} (-\zeta)^\frac{3}{2} = \sqrt(z^2 - 1) - \arccos \frac{1}{z} \f] Therefore: \f[ (\zeta)^\frac{3}{2} = \frac{2i}{3} \left[\sqrt(z^2 - 1) - \arccos \frac{1}{z}\right] \f] */ static complex zeta32greater1(complex z) { complex ret; ret = csqrt(sqr(z) - 1) - cacos(1 / z); return 2/3 * I * ret; } /*!\brief Compute \f$\zeta^\frac{3}{2}\f$ in large order computation Compute \f$\zeta^\frac{3}{2}\f$ where \f$\zeta\f$ is used for computing asymptotic uniform expansion of bessel function as given by (9.3.35) [2, p368] */ static complex zeta32(complex z) { assert(creal(z) > 0); return norm(z) < 1 ? zeta32lesser1(z) : zeta32greater1(z); } /*!\brief Computed \f$lambda_s\f$ \f$lambda_s\f$ is defined by equation (9.3.41) [2,p368] \f[ \lambda_s=\frac{(2s+1)(2s+3)\dotsm(2s+2k+1)\dotsm(6s-1)}{s!(144)^s} \f] \note \f$\lambda_0=1$ see (9.3.40) [2,p368] comment \note Only 11 \f$lambda_s\f$ because we know only 11 \f$u_s\f$ */ static double lambda_s[]= { 1, /* 0 */ 1.0416666666666667e-01, /* 1 */ 8.3550347222222222e-02, /* 2 */ 1.2822657455632716e-01, /* 3 */ 2.9184902646414046e-01, /* 4 */ 8.8162726744375765e-01, /* 5 */ 3.3214082818627675e+00, /* 6 */ 5.6983899350077708e+02, /* 7 */ 2.9990744944402877e+03, /* 8 */ 1.8029158476994044e+04, /* 9 */ 1.2188462345384515e+05, /* 10 */ 9.1528888635321219e+05 /* 11 */ }; /*!\brief Factorized form of \f$u_k\f$ polynom */ struct ukroot { double a_n; /*!< First factor */ complex roots[]; /*!< Root list */ }; /*!\brief \f$u_1\f$ roots and first coefficient \f[ u_1 = 1/8 t - 5/24 t^3 \f] */ static struct ukroot u_1 = { -0.20833333333333333333, /* 5/24 */ { 0, -0.77459666924148338e0, 0.77459666924148338e0 } }; /*!\brief \f$u_2\f$ roots and first coefficient \f[ u_2 &= 9/128 t^2 - 77/192 t^4 + 385/1152 t^6 \f] */ static struct ukroot u_2 = { 0.33420138888888888889, /* 385/1152 */ { 0, 0, -0.9933755698285430168303331391995165308e0, -0.4617412449281712766841302930873661235e0, +0.4617412449281712766841302930873661235e0, +0.9933755698285430168303331391995165308e0, } }; /*!\brief u_3 roots and first coefficient \f[ u_3 = -85085/82944t^9 + 17017/9216 t^7 + 4563/5120 t^5 + 75/1024} t^39/128 t^2 - 77/192 t^4 + 385/1152 t^6 \f] */ static struct ukroot u_3 = { -1.02581259645061728395, /* -85085/82944 */ { 0, 0, 0, -0.138580249705429438745460173524421702e1 - I * 0.101220782843783380746980555686986861e1, -0.138580249705429438745460173524421702e1 + I * 0.101220782843783380746980555686986861e1, -0.907317707191653384781336166975996007e-1, +0.907317707191653384781336166975996007e-1, +0.138580249705429438745460173524421702e1 - I * 0.101220782843783380746980555686986861e1, 0.138580249705429438745460173524421702e1 + I * 0.101220782843783380746980555686986861e1 } }; /*!\brief u_4 roots and first coefficient \f[ u_4 = 37182145/7962624 t^{12} - 7436429/663552 t^{10} + 144001/16384 t^8 - 96833/40960 t^6 + 3675/32768 t^4 \f] */ static struct ukroot u_4 = { 4.66958442342624742798, /* 37182145/7962624 */ { 0, 0, 0, 0, -0.100041841035515907398258771298995e1, -0.941110901446937114426276664565697e0, -0.6736186161701436742651354273089735e0, -0.24435882498737510668615735106017847e0, 0.24435882498737510668615735106017847e0, 0.6736186161701436742651354273089735e0, 0.941110901446937114426276664565697e0, 0.100041841035515907398258771298995e1, } }; /*!\brief \f$u_5\f$ roots and first coefficient \f[ u_5 = -5391411025/191102976 t^{15} + 5391411025/63700992 t^{13} - 108313205/1179648 t^{11} + 250881631/1179648 t^9 - 67608983/9175040 t^7 + 59535/62144 t^5 \f] */ static struct ukroot u_5 = { -28.21207255820024487740, { 0, 0, 0, 0, 0, -0.1670438449895470671100910542996629e1, -0.9270613773884269855139917426769127e0 + I * 0.8782105373619380223726485586305577e0, -0.9270613773884269855139917426769127e0 - I * 0.8782105373619380223726485586305577e0, -0.2052781979771350691712212640237954e0 - I * 0.1597200535972700029703224113919119e0, -0.2052781979771350691712212640237954e0 + I * 0.1597200535972700029703224113919119e0, 0.2052781979771350691712212640237954e0 - I * 0.1597200535972700029703224113919119e0, +0.2052781979771350691712212640237954e0 + I * 0.1597200535972700029703224113919119e0, +0.9270613773884269855139917426769127e0 - I * 0.8782105373619380223726485586305577e0, +0.9270613773884269855139917426769127e0 + I * 0.8782105373619380223726485586305577e0, +0.1670438449895470671100910542996629e1 } }; /*!\brief \f$u_6\f$ roots and first coefficient \f[ u_6 =-6183948445675/27518828544 t^{18} + 415138648925/509607936 t^{16} -4775249765/4194304 t^{14} + 7176153985/9437184 t^{12} -75861726551/314572800 t^{10} + 440748681/4680064 t^8 -2837835/4194304 t^6 \f] */ static struct ukroot u_6 = { -224.71699461288667273874, /* -6183948445675/27518828544 */ { 0, 0, 0, 0, 0, 0, -0.124382324519963642353346130366904e1, -0.104686419493993367114570194049547e1 + I * 0.332780470809332051417451647998361e0, -0.104686419493993367114570194049547e1 - I * 0.332780470809332051417451647998361e0, -0.4886325386868278874127653855249475e0 - I * 0.4342768586722670861475045894935201e0 -0.4886325386868278874127653855249475e0 + I * 0.4342768586722670861475045894935201e0 -0.8554751362758030119554608462468941e-1, +0.8554751362758030119554608462468941e-1, +0.4886325386868278874127653855249475e0 + I * 0.4342768586722670861475045894935201e0, +0.4886325386868278874127653855249475e0 - I * 0.4342768586722670861475045894935201e0, +0.104686419493993367114570194049547e1 + I * 0.332780470809332051417451647998361e0, +0.104686419493993367114570194049547e1 - I * 0.332780470809332051417451647998361e0, 0.124382324519963642353346130366904e1, } }; /*!\brief \f$u_7\f$ roots and first coefficient \f[ u_7 = 1329548915820125/660451885056 t^{21} - 623575990562525/73383542784 t^{19} + 117495020467825/8153726976 t^{17} - 11287669528993/905969664 t^{15} + 4806755210517/838860800 t^{13} - 23169978705569/17616076800 t^{11} + 28375388975/234881024 t^9 - 66891825/33554432 t^7 \f] */ static struct ukroot u_7 = { 2013.08974340710977661787, /* 1329548915820125/660451885056 */ { 0, 0, 0, 0, 0, 0, 0, -0.998039537598447177441282460677148435648e0, -0.99724587573846447076295130893102873649e0, -0.94882272724636817247619094545e0, -0.837374760217824635777173307031e0, -0.6573774609597790399970190455416e0, -0.4196709133247513072532002482099e0, -0.144245206044656036277453505803109e0, +0.144245206044656036277453505803109e0, +0.4196709133247513072532002482099e0, +0.6573774609597790399970190455416e0, +0.837374760217824635777173307031e0, +0.94882272724636817247619094545e0, +0.99724587573846447076295130893102873649e0, +0.998039537598447177441282460677148435648e0, } }; /*!\brief \f$u_8\f$ roots and first coefficient \f[ u_8 = - 2671063771882631125/126806761930752 t^{24} + 59582343274078625/587068342272 t^{22} - 237670164192005545/1174136684544 t^{20} + 42077740772116621/195689447424 t^{18} - 1257093219318079/9663676416 t^{16} + 296916207863309/6710886400 t^{14} - 2904156520889/3758096384 t^{12} + 146540630595/268435456 t^{10} - 14783093325/2147483648 t^8 \f] */ static struct ukroot u_8 = { -21064.04840887960177721512, /* -2671063771882631125/126806761930752 */ { 0, 0, 0, 0, 0, 0, 0, 0, -0.12855209236173717865334911825102e1, -0.11487230378788506995018715201335e1 + I * 0.3055693429633894025779426362946e0, -0.11487230378788506995018715201335e1 - I * 0.3055693429633894025779426362946e0, -0.762452554807904557645736589045814e0 + I * 0.46398205611410723966390073095189e0, -0.762452554807904557645736589045814e0 - I * 0.46398205611410723966390073095189e0, -0.21817886085393558522734608580627e0 - I * 0.251708826111540961305085325918613e0, -0.21817886085393558522734608580627e0 + I * 0.251708826111540961305085325918613e0, -0.112598695910012723571769914377913e0, +0.112598695910012723571769914377913e0, +0.21817886085393558522734608580627e0 - I * 0.251708826111540961305085325918613e0, +0.21817886085393558522734608580627e0 + I * 0.251708826111540961305085325918613e0, +0.762452554807904557645736589045814e0 + I * 0.46398205611410723966390073095189e0, +0.762452554807904557645736589045814e0 - I * 0.46398205611410723966390073095189e0, +0.11487230378788506995018715201335e1 + I * 0.3055693429633894025779426362946e0, +0.11487230378788506995018715201335e1 - I * 0.3055693429633894025779426362946e0, +0.12855209236173717865334911825102e1 } }; /*!\brief \f$u_9\f$ roots and first coefficient \f[ u_9 = 6904699850316601458125/27390260577042432 t^{27} - 461679745093525633675/338151365148672 t^{25} + 29412227933816172445/9393093476352 t^{23} - 37073088786771732695/9393093476352 t^{21} + 6190356955971173843/2087354105856 t^{19} - 519996976064854883/386547056640 t^{17} + 3996930023590772587/11274289152000 t^{15} - 1468251292588911/30064771072 t^{13} + 517406211757245/188978561024 t^{11} - 468131288625/7179869184 t^9 \f] */ static struct ukroot u_9 = { 252085.94970811930830602353, /* 6904699850316601458125/27390260577042432 */ { 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1058716313630550998779067364599140564345e1, -0.1012744489055979013043786254838985367869e1 - I * 0.8079524305020719395243932495366809870049e-1, -0.1012744489055979013043786254838985367869e1 + I * 0.8079524305020719395243932495366809870049e-1, -0.867199766734964955853948753249e0 + I * 0.12151310340817222279307347891e0, -0.867199766734964955853948753249e0 - I * 0.12151310340817222279307347891e0, -0.600543999351638174178690057091e0 + I * 0.107411727631102758807388532378e0, -0.600543999351638174178690057091e0 - I * 0.107411727631102758807388532378e0, -0.216370068593315513557972526683829e0 + I * 0.68924788863076631695122509655846e-1, -0.216370068593315513557972526683829e0 - I * 0.68924788863076631695122509655846e-1, 0.216370068593315513557972526683829e0 - I * 0.68924788863076631695122509655846e-1, 0.216370068593315513557972526683829e0 + I * 0.68924788863076631695122509655846e-1, 0.600543999351638174178690057091e0 + I * 0.107411727631102758807388532378e0, 0.600543999351638174178690057091e0 - I * 0.107411727631102758807388532378e0, 0.867199766734964955853948753249e0 - I * 0.12151310340817222279307347891e0, 0.867199766734964955853948753249e0 + I * 0.12151310340817222279307347891e0, 0.1012744489055979013043786254838985367869e1 - I * 0.8079524305020719395243932495366809870049e-1, 0.1012744489055979013043786254838985367869e1 + I * 0.8079524305020719395243932495366809870049e-1, 0.1058716313630550998779067364599140564345e1 } }; /*!\brief \f$u_{10}\f$ roots and first coefficient \f[ u_{10} = - 4464578923214714502823625/1314732507698036736 t^{30} + 1491741571661309972478475/73040694872113152 t^{28} - 286485860646564818213975/5410421842378752 t^{26} + 17416724745836133903185/225434243432448 t^{24} - 2318810066751674077595/33397665693696 t^{22} + 16487466760916641832507/417470821171200 t^{20} - 45614196450845087013377/3246995275776000 t^{18} + 12736018679229356269/4294967296000 t^{16} - 49042918914307396335/148159191842816 t^{14} + 22818897912239625/1511828488192 t^{12} - 33424574007825/274877906944 t^{10} \f] */ static struct ukroot u_10 = { -3395807.81419312384897239204, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.10001090783214750809440416432541680599188e1, -0.9993294917832275718734757908005237189461e0, -0.9888500788629341914999395888452394835901557e0, -0.9524074579126822822167397126841777753775329e0, -0.882853338709630498587057483665940085863536107e0, -0.7792607086106625162758717324993474593597737525e0, -0.64368000020409505993637643841829812109365971659e0, -0.48077792088388172553910915777979533697628420406e0, -0.2971121517155323652483396314558797798412890267508e0, -0.1005022696406253921602443854248310446528817542823e0, +0.1005022696406253921602443854248310446528817542823e0, +0.2971121517155323652483396314558797798412890267508e0, 0.48077792088388172553910915777979533697628420406e0, 0.64368000020409505993637643841829812109365971659e0, 0.7792607086106625162758717324993474593597737525e0, 0.882853338709630498587057483665940085863536107e0, 0.9524074579126822822167397126841777753775329e0, 0.9888500788629341914999395888452394835901557e0, 0.9993294917832275718734757908005237189461e0, 0.10001090783214750809440416432541680599188e1 } } /*!\brief \f$u_{11}\f$ roots and first coefficient \f[ u_{11} = + 1604407316678887857241980875/31553580184752881664 t^{33} - 392956135061308232223935125/1168651117953810432 t^{31} + 1136426675054854043018733575/1168651117953810432 t^{29} - 209347382482809784715977475/129850124217090048 t^{27} + 4057208263006803008962871/2404631929946112 t^{25} - 7721779642093423553411219/6679533138739200 t^{23} + 24317219180863821793468459/46756731971174400 t^{21} - 185223434015774166216919/1236950581248000 t^{19} + 109880525895631410224729/4233119766937600 t^{17} - 258444671539364377513/107752139522048 t^{15} + 28529686078127770725/314460325543936 t^{13} - 1327867167401775/2199023255552 t^{11} \f] */ /*!\brief Table of \f$u_k\f$ The different u are defined by the following recurence relation: \f{align*} u_0 &= 1 \\ u_{k+1} = \frac{1}{2} t^2 (1-t^2) u'_k + \frac{1}{8} \int_0^t (1-5t^2) u_k(t)\; dt \f} We do not use direct evaluation by monomial formula due to numerical ill behavior \note We compute roots using a arbitrary precision software \f} */ struct ukroot * ukroots[] = { /* u_0 (special case) */ NULL, /* u_1 */ &u_1, /* u_2 */ &u_2, /* u_3 */ &u_3, /* u_4 */ &u_4, /* u_5 */ &u_5, /* u_6 */ &u_6, /* u_7 */ &u_7, /* u_8 */ &u_8, /* u_9 */ &u_9, /* u_10 */ &u_10 }; /*!\brief Bessel function for large order */ complex cbesseljn_largeorder(unsigned int n, complex z) { } #endif /*!\brief Main entry point for besselj function Bessel function of first kind and integer order is The Bessel functions of the first kind \f$J_n(z)\f$ are defined as the solutions to the Bessel differential equation (for \f$n \in \mathbb{N}\f$) \f[ x^2\frac{\text{d}^2y}{\text{d}x^2}+x\frac{\text{d}y}{\text{d}x}+(x^2-n^2)y=0 \f] which are nonsingular at the origin We use two properties of bessel function of first kind in order to compute in the first quadrant: \f[ J_{-v}(z) = (-1)^v J_v(z) \f] Therefore we could compute only whan \f$\Re\text{e}\; z > 0\f$, and apply the following formula: \f{align*} J_{v^*}(z^*) &= (J_v(z))^* & \Re\text{e}\; z &> 0 \f} */ complex cbesseljn (unsigned int n, complex z) { bool retopp = false; bool retconj = false; complex ret; if(isnan(creal(z)) || isnan(cimag(z))) return NAN + I*NAN; /* J_n(-z)=(-1)^n J_n(z) */ if(creal(z) < 0.0) return n % 2 ? cbesseljn (n, -z) : - cbesseljn (n, -z); assert(creal(z) >= 0.0); /* J_n(~z)=~(J_n(z)) */ if(cimag(z) < 0.0) return conj(cbesseljn(n, conj(z))); /* we are in first quadrant */ assert(cimag(z) >= 0.0); /* large argument */ if(norm(z) >= 10 * norm(sqr(n) - 1/4)) return cbesseljn_largearg (n, z); if (norm (z) < norm(SMALL_JN_BOUND)) return cbesseljn_smallarg (n, z); else if (norm (z) > norm(BIG_JN_BOUND)) return cbesseljn_largearg (n, z); else return cbesseljn_mediumarg (n, z); } #ifdef CBESSELJ_TEST #include "test_cbesseljn.c" #endif /* CBESSELJ_TEST */