"Fossies" - the Fresh Open Source Software Archive

Member "ccmath-2.2.1/manual/C09-sfunc" (28 Jan 2000, 23789 Bytes) of package /linux/misc/old/ccmath-2.2.1.tar.gz:


As a special service "Fossies" has tried to format the requested text file into HTML format (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1                                 Chapter 9
    2 
    3                             SPECIAL FUNCTIONS
    4 
    5                                  Summary
    6 
    7                The special functions library segment contains
    8                procedures for the evaluation of many of the
    9                higher transcendental functions encountered in
   10                applications. The functions covered are:
   11 
   12                     o Factorial and Related Functions
   13                     o Elliptic Functions and Integrals
   14                     o Bessel Functions of General Order
   15                     o Spherical Bessel Functions
   16                     o Airy Functions
   17 
   18                The evaluation algorithms employ analytic
   19                approximations in order to cover the widest
   20                possible range of parameters for the functions.
   21                The elliptic integral functions employ the Bartky
   22                parameterization, which significantly simplifies
   23                their use. All the library functions maintain a
   24                high standard of accuracy over the full argument
   25                range.
   26 
   27 -------------------------------------------------------------------------------
   28 
   29  Notes on Contents:
   30 
   31      Functions in this library segment compute the values of some frequently
   32  used higher transcendental functions.  The accuracy standards quoted below
   33  all refer to relative rather than absolute precision.
   34 
   35  o Factorial and Related Functions:
   36 
   37     gaml  -------- evaluate the natural logarithm of the Gamma function.
   38     psi  --------- evaluate the Psi-function for integer argument.
   39     psih  -------- evaluate the Psi-function for half-integer argument.
   40 
   41      The logarithm of the factorial function is computed to an accuracy of 12
   42  or more decimal digits for general positive arguments. Evaluation of the Psi
   43  function is confined to the integer and half-integer arguments required in
   44  the evaluation of other functions in the library.
   45 
   46 
   47  o Elliptic Integrals and Functions:
   48 
   49     amelp  ------- compute the elliptic amplitude function.
   50     felp  -------- compute complete and incomplete elliptic integrals
   51                    of the first and second kind.
   52     gelp  -------- compute a general elliptic integral.
   53     g2elp  ------- compute a general elliptic integral with nonzero
   54                    lower limit.
   55     nome  -------- evaluate the nome q(k) for elliptic functions of
   56                    modulus k.
   57     theta  ------- evaluate the Jacobian theta functions.
   58 
   59      Elliptic integrals are evaluated by using an extension of the Bartky
   60  approach, using arithmetic and geometric means. The accuracy standard for
   61  elliptic integrals is 14 or more decimal digits (ie. nearly full double
   62  precision accuracy).
   63 
   64      Jacobian elliptic functions are evaluated with the aid of theta function
   65  series. The accuracy standard for these evaluations is at least 14 decimal
   66  digits. Rapid convergence of the series yields high execution efficiency.
   67 
   68 
   69  o Bessel Functions of General Order:
   70 
   71     jbes  -------- evaluate the Bessel function of the first kind.
   72     nbes  -------- evaluate the Bessel function of the second kind.
   73     ibes  -------- evaluate the Bessel function of the first kind
   74                    at imaginary argument values.
   75     kbes  -------- evaluate the Bessel function of the third kind
   76                    at imaginary argument values.
   77     drbes  ------- evaluate derivatives of the Bessel functions.
   78     rcbes  ------- compute a sequence of Bessel function values
   79                    via recurrence relations.
   80 
   81      The general Bessel functions maintain a standard of accuracy of at least
   82  12 decimal digits over the full argument range. Series of Bessel functions,
   83  with orders differing by one, are computed using recurrence relations. The
   84  library also contains a function that analytically computes Bessel function
   85  derivatives.
   86 
   87 
   88  o Spherical Bessel Functions:
   89 
   90     jspbes  ------- evaluate the spherical Bessel function of the
   91                     first kind.
   92     yspbes  ------- evaluate the spherical Bessel function of the
   93                     second kind.
   94     kspbes  ------- evaluate spherical Bessel functions of the third
   95                     kind at imaginary argument values.
   96     drspbes  ------ compute the derivatives of spherical Bessel
   97                     functions.
   98     rcspbes  ------ compute a series of spherical Bessel function
   99                     values via recurrence relations.
  100 
  101      The accuracy standard and coverage provided for spherical Bessel
  102  functions matches that used for the Bessel functions.
  103 
  104 
  105  o Airy Functions:
  106 
  107      Airy functions of the first and second kind and their derivatives are
  108  evaluated for all real arguments by a pair of library functions.
  109 
  110     airy  -------- evaluate the Airy function of the first kind or
  111                    its derivative.
  112     biry  -------- evaluate the Airy function of the second kind
  113                    or its derivative.
  114 
  115  The standard of a minimum of 12 accurate decimal digits is maintained for
  116  these functions.
  117 
  118 
  119  -----------------------------------------------------------------------------
  120 
  121  General Technical Comments:
  122 
  123      The primary criteria governing the selection of evaluation algorithms
  124  was the ability to maintain a high standard of accuracy over a comprehensive
  125  parameter range. This dictated the use of analytic approximations rather than
  126  optimized rational approximations. The latter approach can yield highly
  127  efficient execution, but it requires separate approximations for each order
  128  covered.
  129 
  130      The range of orders for which Bessel function evaluation is efficient
  131  extends from zero to one hundred. Efficiency declines as order increases,
  132  however, and the use of asymptotic expansions in order is advised for the
  133  rare application where many functions with order exceeding 100 must be
  134  evaluated.
  135 
  136      General forms of both complete and incomplete elliptic integrals are
  137  evaluated in a uniform manner, that does not require the user to identify
  138  a combination of integrals of the first, second, and third kind. The method
  139  has been extended to the case of singular integrands, by means of an
  140  automatically evoked transformation function (gsing). This approach
  141  represents a significant simplification in the use of this important class
  142  of functions. Its rapid convergence ensures high efficiency.
  143 
  144 -------------------------------------------------------------------------------
  145                              FUNCTION SYNOPSES
  146 -------------------------------------------------------------------------------
  147 
  148      Factorial Functions:
  149 -------------------------------------------------------------------------------
  150 
  151 gaml
  152 
  153      Evaluate the natural logarithm of the Gamma function.
  154 
  155      double gaml(double x)
  156        x = function argument
  157       return value: y = log(Gamma(x))     (logarithm base = e)
  158 
  159           This function is also used in the statistical function library.
  160 
  161      ----------------------------------------------------------------------
  162 
  163 psi
  164 
  165      Evaluate the psi function for integer and half-integer arguments.
  166 
  167                 psi(x) = d {log(Gamma(x))} / dx
  168 
  169      double psi(int m)
  170        m = integer function argument (m>0)
  171       return value: y = psi(m)
  172 
  173 
  174      double psih(double v)
  175        v = half-integer function argument (v=n+.5 for integer n>=0)
  176       return value: y = psi(v)
  177 
  178           The arguments of the psi functions are those needed to
  179           support the other functions in the library.
  180 
  181 -------------------------------------------------------------------------------
  182 
  183      Elliptic Integrals and Functions:
  184 -------------------------------------------------------------------------------
  185 
  186       F(x,k) = Intg(0 to x) 1/sqrt(1 - k^2*[sin(a)]^2) da ,
  187 
  188       E(x,k) = Intg(0 to x) sqrt(1 - k^2*[sin(a)]^2) da , and
  189 
  190       I(x,r,k) = Intg(0 to x) da/{(1-r*[sin(a)]^2)*sqrt(1-k^2*[sin(a)]^2)}
  191 
  192  are the standard definitions of elliptic integrals of the first, second,
  193  and third kind respectively. The parameter k is the modulus of the elliptic
  194  integrals. Complete elliptic integrals are given by
  195 
  196           K = F(pi/2,k)  and  E = E(pi/2,k).
  197 
  198  The Bartky parameterization of elliptic integrals is outlined in the
  199  discussion of gelp below.
  200 -------------------------------------------------------------------------------
  201 
  202 amelp
  203 
  204      Compute the elliptic amplitude function phi = amp(u).
  205 
  206      double amelp(double u,double k)
  207        u = argument of function (u=F(v,k), u>=0)
  208        k = modulus of elliptic functions
  209       return value: v = amp(u)
  210 
  211 
  212           This is the inverse of the elliptic integral of the first kind.
  213           If  u = F(v,k) , then  v = v(u,k) is the amplitude function.
  214           Jacobian elliptic functions satisfy
  215 
  216           sn(u) = sin(v), cn(u) = cos(v),
  217           and dn(u) = sqrt(1-k^2*[sin(v)]^2).
  218           
  219      --------------------------------------------------------------------
  220 
  221 felp
  222 
  223      Compute complete and incomplete elliptic integrals of the first and
  224      second kinds.
  225 
  226      double felp(double an,double k,double *pk,double *pz,double *ph)
  227        an = upper limit of elliptic integrals (an > 0 in radians)
  228        k = modulus of the elliptic integrals
  229        pk = pointer to store for complete integral K
  230              ( output of K, E and E(an,k) can be suppressed
  231                by calling with pk=NULL )
  232        pz = pointer to store for function E(an,k)
  233              ( pz=NULL will suppress output of E(an,k) and E )
  234        ph = pointer to store for complete integral E(k)
  235       return value: f = F(an,k) the elliptic integral of the first kind
  236 
  237 
  238      --------------------------------------------------------------------
  239 
  240      The integrand of a general elliptic integral is parameterized by
  241      Bartky in the following form:
  242 
  243           G(v,k:as,bs,cs) = Intg(0 to v) {M(as,bs,cs,k,y)/R(y,k)} dy  , with
  244 
  245           R(y,k) = sqrt{[cos(y)]^2 + sqrt(1-k^2)*[sin(y)]^2} .
  246 
  247           M = N/D  with,
  248 
  249           N = { as*[cos(y)]^2 + bs*cs[sin(y)]^2 }  and
  250 
  251           D = { as*[cos(y)]^2 + cs*sqrt(1-k^2)*[sin(y)]^2 } .
  252 
  253 
  254      The Bartky parameters B=[as,bs,cs] can be chosen to produce a
  255      general elliptic integral. Integrals of the first, second, and
  256      third kind result from the following initial choices:
  257 
  258           as = 1 and bs = cs = 1  for F(v,k) ,
  259 
  260           as = 1 , bs = 1 - k^2 , and cs = sqrt(1 - k^2)  for E(v,k) , and
  261 
  262           as = 1 , bs = 1/(1-r) , and cs = (1-r)/sqrt(1 - k^2)  for I(v,r,k).
  263 
  264      -------------------------------------------------------------------------
  265 
  266 gelp
  267 
  268      Compute a general elliptic integral G(an,k,B).
  269 
  270      double gelp(double an,double k,double as,double bs,double cs, \
  271                  double *pg,double *pf,double *pk)
  272        an = upper limit of integral (an > 0 in radians)
  273        k = modulus of elliptic integral
  274        as,bs,cs = Bartky parameters of integrand (see note above)
  275        pg = pointer to store for complete integral G=G(pi/2,k,B)
  276              ( output of G, K, and F(an,k) is suppressed by calling
  277                with pg=NULL )
  278        pf,pk = pointers to store for F(an,k) and K respectively
  279                 (output of these integrals of the first kind
  280                  is suppressed when pf=0)
  281       return value: g = G(an,k,B) the incomplete integral
  282 
  283      --------------------------------------------------------------------
  284 
  285 g2elp
  286 
  287      Compute a general elliptic integral with nonzero lower limit
  288      G2(an,bn,k,B).
  289 
  290      double g2elp(double an,double bn,double k,double as,double bs, \
  291                   double cs)
  292        an = lower limit of the integral (radians)
  293        bn = upper limit of the integral (radians)
  294        k = modulus of elliptic integral
  295        as,bs,cs = Bartky parameters of integrand
  296       return value: g2 = G2(an,bn,k,B) = G(bn,k,B)-G(an,k,B)
  297 
  298 
  299          G2(a2,a1,k,B) = Intg(w=a1 to a2) { M(B,k,w)/R(w,k) dw}
  300 
  301      -------------------------------------------------------------------
  302 
  303 
  304      The following two functions are evoked automatically by gelp and g2elp
  305      respectively when a singular integral (cs<0.) is involved. They are not
  306      normally called directly by an application.
  307     
  308 gsng
  309 
  310      Convert a singular elliptic integral (*pc < 0.) to non-singular form.
  311 
  312      double gsng(double *pa,double *pb,double *pc,double b,double an)
  313        pa,pb,pc = pointers to the stores for the integrand's Bartky
  314                   parameters (these parameters correspond to as,bs,cs
  315                   in the note above, with singularity -> cs<0.)
  316        b = modulus parameter (b=sqrt(1-k*k))
  317        an = upper limit of integral (an > 0 in radians)
  318       return value: I = increment to integral after transformation
  319                         (I=HUGE -> limit is coincident with the
  320                          singularity of the integrand)
  321 
  322      -----------------------------------------------------------------------
  323 
  324 gsng2
  325 
  326      Convert an elliptic integral with non-zero lower limit to non-singular
  327      form.
  328 
  329      double gsng2(double *pa,double *pb,double *pc,double b, \
  330                   double an,double bn)
  331      double *pa,*pb,*pc,b,an;
  332        pa,pb,pc = pointers to the stores for the integrand's Bartky
  333                   parameters (these parameters correspond to as,bs,ds
  334                   in the note above, with singularity -> ds<0.)
  335        b = modulus parameter (b=sqrt(1-k*k))
  336        an = lower limit of integral (radians)
  337        bn = upper limit of integral
  338       return value: I2 = increment to integral after transformation
  339                         (I2=HUGE -> limit is coincident with the
  340                          singularity of the integrand)
  341 
  342 
  343           The Bartky parameters *pa,*pb,*pc are altered by these
  344           functions. The value of the integral is given by
  345 
  346             G(v,k:as,bs,cs) = I + G'(v,k:as',bs',cs') ,
  347 
  348             or G2(v2,v1:as,bs,cs) = I2 + G'(v2,v1:as',bs',cs') ,
  349 
  350           with  as', bs', and cs' the transformed parameters
  351           resulting from application of gsng or gsng2. If the
  352           singularity of the integrand lies inside the region of
  353           integration, this result must be interpreted as a Cauchy
  354           principle value integral.
  355 
  356      -----------------------------------------------------------------------
  357 
  358 nome
  359 
  360      Evaluate the nome q(k) for elliptic functions of modulus k.
  361 
  362      double nome(double k,double *pk,double *pkp)
  363        k = modulus of elliptic functions
  364        pk = pointer to store for K
  365        pkp = pointer to store for K1
  366       return value: f = q(k)
  367 
  368 
  369           The nome is defined by
  370 
  371            q(k) = exp( -pi*K1/K) ,  with K = F(pi/2,k),
  372 
  373           and  K1 = F(pi/2,k1) ,  where  k1 = sqrt(1 - k^2) .
  374 
  375      -----------------------------------------------------------
  376 
  377 theta
  378 
  379      Evaluate the Jacobian theta functions.
  380 
  381      double theta(double u,int n)
  382        u = argument of the theta function
  383        n = index of the theta function = 0,1,2,or 3
  384       return value:  theta(u,n)
  385 
  386 
  387      The computation of theta functions must be initialized by
  388      calling the function stheta to initialize the modulus.
  389 
  390 
  391      void stheta(double k)
  392         k = modulus of elliptic functions
  393 
  394 
  395           The theta functions are related to elliptic functions by:
  396 
  397              sn(u) = theta1(u)/(sqrt(k)* theta0(u)) ;
  398 
  399              cn(u) = sqrt(k1/k)* theta2(u)/theta0(u); and
  400 
  401              dn(u) = sqrt(k1)* theta3(u)/theta0(u).
  402 
  403              k is the modulus and k1 = sqrt(1-k^2) .
  404 
  405 
  406 -------------------------------------------------------------------------------
  407 
  408      Bessel Functions:
  409 -------------------------------------------------------------------------------
  410 
  411      The Bessel functions evaluated by the library are defined by
  412 
  413           J(n,x) = {(x/2)^n}/Gamma(n)}* Sum(k=0 to infinity) {
  414                      (-x^2/4)^k / k!*Gamma(n+k) }  ,
  415 
  416           Y(n,x) = N(n,x) = {cos(pi*n)J(n,x) - J(-n,x)}/sin(pi*n) ,
  417 
  418           I(n,x) = (i)^{-n} * J(n,i*x) , and
  419 
  420           K(n,x) = (i*pi/2)(i)^n *{ J(n,i*x) + i*Y(n,i*x) }  .
  421 
  422  These definitions require the use of limit procedures for Y and K when the
  423  order n is an integer.
  424 
  425 -------------------------------------------------------------------------------
  426 
  427 jbes
  428 
  429      Evaluate the Bessel function of the first kind J(n,x).
  430 
  431      double jbes(double v,double x)
  432        v = order of the function (v>=0.)
  433        x = function argument (x>=0.)
  434       return value: f = J(n,x)
  435 
  436      -----------------------------------------------------------
  437 
  438 nbes
  439 
  440      Evaluate the Bessel function of the second kind Y(n,x).
  441      ( N(n,x) is an alternative notation.)
  442 
  443      double nbes(double v,double x)
  444        v = order of the function (v>=0.)
  445        x = function argument (x>=0.)
  446       return value: f = Y(n,x)
  447 
  448      ------------------------------------------------------------
  449 
  450 ibes
  451 
  452      Evaluate the modified Bessel function of the first kind I(n,x).
  453      ( This is a Bessel function with imaginary argument. )
  454 
  455      double ibes(double v,double x)
  456        v = order of the function (v>=0.)
  457        x = function argument (x>=0.)
  458       return value: f = I(n,x)
  459 
  460      -------------------------------------------------------------------
  461 
  462  kbes
  463 
  464      Evaluate the modified Bessel function of the third kind K(n,x).
  465 
  466      double kbes(double v,double x)
  467        v = order of the function (v>=0.)
  468        x = function argument (x>=0.)
  469       return value: f = K(n,x)
  470 
  471      -------------------------------------------------------------------
  472 
  473 drbes
  474 
  475      Evaluate the derivatives of Bessel functions.
  476 
  477      double drbes(double x,double v,int f,double *p)
  478        x = argument of the function (x>=0.)
  479        v = order of the function (v>=0.)
  480        p = pointer to store for supplied function value Z(v,x)
  481             (p=0 -> compute the required function value)
  482        f = function type flag, with:
  483              f='j' -> dJ(v,x)/dx
  484              f='y' -> dY(v,x)/dx
  485              f='i' -> dI(v,x)/dx
  486              f='k' -> dK(v,x)/dx
  487       return value: f = dZ(v,x)/dx , where Z=J,Y,I,or K
  488 
  489      ----------------------------------------------------------
  490 
  491  rcbes
  492 
  493      Recurrence relations for the evaluation of a sequence of Bessel
  494      functions.
  495 
  496      double rcbes()
  497       return value: f = value of the next Bessel function in the
  498                         order of the recurrence (see setrcb below)
  499 
  500 
  501      The recurrence sequence must be initialized by a call of
  502      the function setrcb. This call determines the function
  503      type, recursion direction, the starting order, and the
  504      argument of the sequence.
  505 
  506      void setrcb(double u,double y,int fl,int dr,double *pf, \
  507                  double *ph)
  508        u = starting order of the Bessel functions (u>=0.)
  509        y = argument of the functions (y>=0.)
  510        fl = function type flag, with:
  511              fl='j' -> J(v,x)
  512              fl='y' -> Y(v,x)
  513              fl='i' -> I(v,x)
  514              fl='k' -> K(v,x)
  515        dr = direction flag, with:
  516              dr='u' -> increasing order
  517              dr='d' -> decreasing order
  518        pf,ph = pointers to store for first two function values
  519                 ( *pf=Z(u,y) , *ph=Z(u+1,y) if dr='u' and
  520                   *pf=Z(u,y) , *ph=Z(u-1,y) if dr='d' )
  521 
  522 
  523            The functions J and I are stable for decreasing recurrence,
  524            while Y and K have stable increasing recurrence.
  525 
  526 
  527 -------------------------------------------------------------------------------
  528 
  529      Spherical Bessel Functions:
  530 -------------------------------------------------------------------------------
  531 
  532      The spherical Bessel functions are defined by
  533 
  534           j(n,x) = sqrt(pi/2x)*J(n+1/2,x) ,
  535 
  536           y(n,x) = sqrt(pi/2x)*N(n+1/2,x) , and
  537 
  538           k(n,x) = {exp(-x)/x}* Sum(k=0 to n) {
  539                                 [(n+k)! / (k!(n-k)!]*(2x)^k)] } .
  540 
  541 -------------------------------------------------------------------------------
  542 
  543 jspbes
  544 
  545      Evaluate the spherical Bessel function j(n,x).
  546 
  547      double jspbes(int n,double x)
  548        x = argument of function (x>=0.)
  549        n = order of function (n>=0)
  550       return value: f = j(n,x)
  551 
  552      --------------------------------------------------------
  553 
  554 yspbes
  555 
  556      Evaluate the spherical Bessel function of the second kind y(n,x).
  557 
  558      double yspbes(int n,double x)
  559        x = argument of function (x>=0.)
  560        n = order of function (n>=0)
  561       return value: f = y(n,x)
  562 
  563      -----------------------------------------------------------
  564 
  565 kspbes
  566 
  567      Evaluate the modified spherical Bessel functions k(n,x) and k(n,-x).
  568 
  569      double kspbes(int n,double x)
  570        x = argument of function
  571             ( x>0. -> exponentially decreasing solution, and
  572               x<0. -> exponentially increasing solution)
  573        n = order of function (n>=0)
  574       return value: f = k(n,x)
  575 
  576 
  577           The normalization chosen for k(n,x) yields a leading
  578           asymptotic term = exp(-x)/x at large x.
  579 
  580      -------------------------------------------------------------------
  581 
  582 drspbes
  583 
  584      Evaluate the derivatives of spherical Bessel functions.
  585 
  586      double drspbes(double x,int n,int f,double *p)
  587        x = argument of function (x>=0. for j and y)
  588        n = order of function (n>=0)
  589        p = pointer to store for supplied function value
  590             (p=0 -> compute the required function value)
  591        f = function type flag, with:
  592             f='j' -> dj(n,x)/dx
  593             f='y' -> dy(n,x)/dx
  594             f='k' -> dk(n,x)/dx
  595       return value: d = dz(n,x)/dx for z=j,y,or k
  596 
  597      ----------------------------------------------------------------
  598 
  599 rcspbes
  600 
  601      Recurrence relations for evaluating a sequence of spherical Bessel
  602      functions.
  603 
  604      double rcspbs()
  605       return value: f = value of the next spherical Bessel function
  606                         in the order of the recurrence (see setrcsb
  607                         below)
  608 
  609 
  610      The recurrence sequence must be initialized by a call of
  611      the function setrcsb. This call determines the function
  612      type, recursion direction, the starting order, and the
  613      argument of the sequence.
  614 
  615      setrcsb(int n,double y,int f,int dr,double *pf,double *ph)
  616        n = starting order of the spherical Bessel function
  617            sequence (n>=0)
  618        y = argument of the functions (y>=0.)
  619        f = function type flag, with:
  620              f='j' -> dj(n,x)/dx
  621              f='y' -> dy(n,x)/dx
  622              f='k' -> dk(n,x)/dx
  623        dr = direction flag, with:
  624              dr='u' -> increasing order
  625              dr='d' -> decreasing order
  626        pf,ph = pointers to store for first two function values
  627                 ( *pf=z(n,y) , *ph=z(n+1,y) if dr='u' and
  628                   *pf=z(n,y) , *ph=z(n-1,y) if dr='d' )
  629 
  630 
  631 -------------------------------------------------------------------------------
  632 
  633      Airy Functions:
  634 -------------------------------------------------------------------------------
  635 
  636      The Airy functions can be defined in terms of Bessel functions, with
  637 
  638           Ai(x) = {sqrt(x/3)/pi} K(1/3,u) for x >= 0
  639 
  640                 = {sqrt(-x)/3} { J(1/3,u) + J(-1/3,u) } for x < 0.
  641 
  642           Bi(x) = sqrt(x/3) { I(-1/3,u) + I(1/3,u) } for x >= 0
  643 
  644                 = sqrt(-x/3) { J(-1/3,u) - J1/3,u)) } for x < 0.
  645 
  646           Here u = (2/3)*(|x|)^1.5 .
  647 
  648 -------------------------------------------------------------------------------
  649 
  650 airy
  651 
  652      Evaluate the Airy function Ai(x) or its derivative Ai'(x).
  653 
  654      double airy(double x,int df)
  655        x = argument of the function
  656        df = derivative flag, with
  657              0 -> Ai and 1 -> Ai'
  658       return value: f = Ai(x) if df=0
  659                     f = Ai'(x) if df=1
  660 
  661      -------------------------------------------------------------------
  662 
  663 biry
  664 
  665      Evaluate the Airy function Bi(x) or its derivative Bi'(x).
  666 
  667      double biry(double x,int df)
  668        x = argument of the function
  669        df = derivative flag, with
  670              0 -> Bi and 1 -> Bi'
  671       return value: f = Bi(x) if df=0
  672                     f = Bi'(x) if df=1