Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  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) */
  142.