toolchain/gcc/newlib/libm/math/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  * wrapper jn(int n, double x), yn(int n, double x)
  16  * floating point Bessel's function of the 1st and 2nd kind
  17  * of order n
  18  *
  19  * Special cases:
  20  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  21  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  22  * Note 2. About jn(n,x), yn(n,x)
  23  *      For n=0, j0(x) is called,
  24  *      for n=1, j1(x) is called,
  25  *      for n<x, forward recursion us used starting
  26  *      from values of j0(x) and j1(x).
  27  *      for n>x, a continued fraction approximation to
  28  *      j(n,x)/j(n-1,x) is evaluated and then backward
  29  *      recursion is used starting from a supposed value
  30  *      for j(n,x). The resulting value of j(0,x) is
  31  *      compared with the actual value to correct the
  32  *      supposed value of j(n,x).
  33  *
  34  *      yn(n,x) is similar in all respects, except
  35  *      that forward recursion is used for all
  36  *      values of n>1.
  37  *
  38  */
  39
  40 #include "fdlibm.h"
  41 #include <errno.h>
  42
  43 #ifndef _DOUBLE_IS_32BITS
  44
  45 #ifdef __STDC__
  46         double jn(int n, double x)      /* wrapper jn */
  47 #else
  48         double jn(n,x)                  /* wrapper jn */
  49         double x; int n;
  50 #endif
  51 {
  52 #ifdef _IEEE_LIBM
  53         return __ieee754_jn(n,x);
  54 #else
  55         double z;
  56         struct exception exc;
  57         z = __ieee754_jn(n,x);
  58         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  59         if(fabs(x)>X_TLOSS) {
  60             /* jn(|x|>X_TLOSS) */
  61             exc.type = TLOSS;
  62             exc.name = "jn";
  63             exc.err = 0;
  64             exc.arg1 = n;
  65             exc.arg2 = x;
  66             exc.retval = 0.0;
  67             if (_LIB_VERSION == _POSIX_)
  68                 errno = ERANGE;
  69             else if (!matherr(&exc)) {
  70                errno = ERANGE;
  71             }
  72             if (exc.err != 0)
  73                errno = exc.err;
  74             return exc.retval;
  75         } else
  76             return z;
  77 #endif
  78 }
  79
  80 #ifdef __STDC__
  81         double yn(int n, double x)      /* wrapper yn */
  82 #else
  83         double yn(n,x)                  /* wrapper yn */
  84         double x; int n;
  85 #endif
  86 {
  87 #ifdef _IEEE_LIBM
  88         return __ieee754_yn(n,x);
  89 #else
  90         double z;
  91         struct exception exc;
  92         z = __ieee754_yn(n,x);
  93         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  94         if(x <= 0.0){
  95             /* yn(n,0) = -inf or yn(x<0) = NaN */
  96 #ifndef HUGE_VAL
  97 #define HUGE_VAL inf
  98             double inf = 0.0;
  99
 100             SET_HIGH_WORD(inf,0x7ff00000);      /* set inf to infinite */
 101 #endif
 102             exc.type = DOMAIN;  /* should be SING for IEEE */
 103             exc.name = "yn";
 104             exc.err = 0;
 105             exc.arg1 = n;
 106             exc.arg2 = x;
 107             if (_LIB_VERSION == _SVID_)
 108                 exc.retval = -HUGE;
 109             else
 110                 exc.retval = -HUGE_VAL;
 111             if (_LIB_VERSION == _POSIX_)
 112                 errno = EDOM;
 113             else if (!matherr(&exc)) {
 114                 errno = EDOM;
 115             }
 116             if (exc.err != 0)
 117                errno = exc.err;
 118             return exc.retval;
 119         }
 120         if(x>X_TLOSS) {
 121             /* yn(x>X_TLOSS) */
 122             exc.type = TLOSS;
 123             exc.name = "yn";
 124             exc.err = 0;
 125             exc.arg1 = n;
 126             exc.arg2 = x;
 127             exc.retval = 0.0;
 128             if (_LIB_VERSION == _POSIX_)
 129                 errno = ERANGE;
 130             else if (!matherr(&exc)) {
 131                 errno = ERANGE;
 132             }
 133             if (exc.err != 0)
 134                errno = exc.err;
 135             return exc.retval;
 136         } else
 137             return z;
 138 #endif
 139 }
 140
 141 #endif /* defined(_DOUBLE_IS_32BITS) */