newlib/libm/mathfp/w_jn.c

   1
   2 /* @(#)w_jn.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 /*
  15 FUNCTION
  16 <<jN>>, <<jNf>>, <<yN>>, <<yNf>>---Bessel functions
  17
  18 INDEX
  19 j0
  20 INDEX
  21 j0f
  22 INDEX
  23 j1
  24 INDEX
  25 j1f
  26 INDEX
  27 jn
  28 INDEX
  29 jnf
  30 INDEX
  31 y0
  32 INDEX
  33 y0f
  34 INDEX
  35 y1
  36 INDEX
  37 y1f
  38 INDEX
  39 yn
  40 INDEX
  41 ynf
  42
  43 SYNOPSIS
  44 #include <math.h>
  45 double j0(double <[x]>);
  46 float j0f(float <[x]>);
  47 double j1(double <[x]>);
  48 float j1f(float <[x]>);
  49 double jn(int <[n]>, double <[x]>);
  50 float jnf(int <[n]>, float <[x]>);
  51 double y0(double <[x]>);
  52 float y0f(float <[x]>);
  53 double y1(double <[x]>);
  54 float y1f(float <[x]>);
  55 double yn(int <[n]>, double <[x]>);
  56 float ynf(int <[n]>, float <[x]>);
  57
  58 DESCRIPTION
  59 The Bessel functions are a family of functions that solve the
  60 differential equation
  61 @ifnottex
  62 .  2               2    2
  63 . x  y'' + xy' + (x  - p )y  = 0
  64 @end ifnottex
  65 @tex
  66 $$x^2{d^2y\over dx^2} + x{dy\over dx} + (x^2-p^2)y = 0$$
  67 @end tex
  68 These functions have many applications in engineering and physics.
  69
  70 <<jn>> calculates the Bessel function of the first kind of order
  71 <[n]>.  <<j0>> and <<j1>> are special cases for order 0 and order
  72 1 respectively.
  73
  74 Similarly, <<yn>> calculates the Bessel function of the second kind of
  75 order <[n]>, and <<y0>> and <<y1>> are special cases for order 0 and
  76 1.
  77
  78 <<jnf>>, <<j0f>>, <<j1f>>, <<ynf>>, <<y0f>>, and <<y1f>> perform the
  79 same calculations, but on <<float>> rather than <<double>> values.
  80
  81 RETURNS
  82 The value of each Bessel function at <[x]> is returned.
  83
  84 PORTABILITY
  85 None of the Bessel functions are in ANSI C.
  86 */
  87
  88 /*
  89  * wrapper jn(int n, double x), yn(int n, double x)
  90  * floating point Bessel's function of the 1st and 2nd kind
  91  * of order n
  92  *
  93  * Special cases:
  94  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  95  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  96  * Note 2. About jn(n,x), yn(n,x)
  97  *      For n=0, j0(x) is called,
  98  *      for n=1, j1(x) is called,
  99  *      for n<x, forward recursion us used starting
 100  *      from values of j0(x) and j1(x).
 101  *      for n>x, a continued fraction approximation to
 102  *      j(n,x)/j(n-1,x) is evaluated and then backward
 103  *      recursion is used starting from a supposed value
 104  *      for j(n,x). The resulting value of j(0,x) is
 105  *      compared with the actual value to correct the
 106  *      supposed value of j(n,x).
 107  *
 108  *      yn(n,x) is similar in all respects, except
 109  *      that forward recursion is used for all
 110  *      values of n>1.
 111  *
 112  */
 113
 114 #include "fdlibm.h"
 115 #include <errno.h>
 116
 117 #ifndef _DOUBLE_IS_32BITS
 118
 119 #ifdef __STDC__
 120         double jn(int n, double x)      /* wrapper jn */
 121 #else
 122         double jn(n,x)                  /* wrapper jn */
 123         double x; int n;
 124 #endif
 125 {
 126 #ifdef _IEEE_LIBM
 127         return jn(n,x);
 128 #else
 129         double z;
 130         z = jn(n,x);
 131         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
 132         if(fabs(x)>X_TLOSS) {
 133             /* jn(|x|>X_TLOSS) */
 134             errno = ERANGE;
 135             return 0.0;
 136         } else
 137             return z;
 138 #endif
 139 }
 140
 141 #ifdef __STDC__
 142         double yn(int n, double x)      /* wrapper yn */
 143 #else
 144         double yn(n,x)                  /* wrapper yn */
 145         double x; int n;
 146 #endif
 147 {
 148 #ifdef _IEEE_LIBM
 149         return yn(n,x);
 150 #else
 151         double z;
 152         z = yn(n,x);
 153         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
 154         if(x <= 0.0){
 155             /* yn(n,0) = -inf or yn(x<0) = NaN */
 156 #ifndef HUGE_VAL
 157 #define HUGE_VAL inf
 158             double inf = 0.0;
 159
 160             SET_HIGH_WORD(inf,0x7ff00000);      /* set inf to infinite */
 161 #endif
 162             errno = EDOM;
 163             return -HUGE_VAL;
 164         }
 165         if(x>X_TLOSS) {
 166             /* yn(x>X_TLOSS) */
 167             errno = ERANGE;
 168             return 0.0;
 169         } else
 170             return z;
 171 #endif
 172 }
 173
 174 #endif /* defined(_DOUBLE_IS_32BITS) */