Le samedi 19 janvier 2008, Mengda Wu a écrit :
> Hi gsler,
>
> Is there a way in GSL to calculate bessel functions with complex numbers
> as input? It seems only real numbers are supported.
Here the bessel file. #if 0 some code that I have not finished...
>
> Thanks,
> Mengda
> _______________________________________________
> Help-gsl mailing list
> [email protected]
> http://lists.gnu.org/mailman/listinfo/help-gsl
--
"ROUCARIES Bastien"
[EMAIL PROTECTED]
-------------------------------------------------------------------------------
DO NOT WRITE TO [EMAIL PROTECTED] OR BE BLACKLISTED
/* 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
* [EMAIL PROTECTED]
* 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 <complex.h>
#include <math.h>
#include <float.h>
#include <fenv.h>
#include <limits.h>
#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#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 */
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
