toolchain/gcc/newlib/libm/math/w_j0.c

   1
   2 /* @(#)w_j0.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 ANSI_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 TRAD_SYNOPSIS
  59 #include <math.h>
  60
  61 double j0(<[x]>)
  62 double <[x]>;
  63 float j0f(<[x]>)
  64 float <[x]>;
  65 double j1(<[x]>)
  66 double <[x]>;
  67 float j1f(<[x]>)
  68 float <[x]>;
  69 double jn(<[n]>, <[x]>)
  70 int <[n]>;
  71 double <[x]>;
  72 float jnf(<[n]>, <[x]>)
  73 int <[n]>;
  74 float <[x]>;
  75
  76 double y0(<[x]>)
  77 double <[x]>;
  78 float y0f(<[x]>)
  79 float <[x]>;
  80 double y1(<[x]>)
  81 double <[x]>;
  82 float y1f(<[x]>)
  83 float <[x]>;
  84 double yn(<[n]>, <[x]>)
  85 int <[n]>;
  86 double <[x]>;
  87 float ynf(<[n]>, <[x]>)
  88 int <[n]>;
  89 float <[x]>;
  90
  91 DESCRIPTION
  92 The Bessel functions are a family of functions that solve the
  93 differential equation
  94 @ifnottex
  95 .  2               2    2
  96 . x  y'' + xy' + (x  - p )y  = 0
  97 @end ifnottex
  98 @tex
  99 $$x^2{d^2y\over dx^2} + x{dy\over dx} + (x^2-p^2)y = 0$$
 100 @end tex
 101 These functions have many applications in engineering and physics.
 102
 103 <<jn>> calculates the Bessel function of the first kind of order
 104 <[n]>.  <<j0>> and <<j1>> are special cases for order 0 and order
 105 1 respectively.
 106
 107 Similarly, <<yn>> calculates the Bessel function of the second kind of
 108 order <[n]>, and <<y0>> and <<y1>> are special cases for order 0 and
 109 1.
 110
 111 <<jnf>>, <<j0f>>, <<j1f>>, <<ynf>>, <<y0f>>, and <<y1f>> perform the
 112 same calculations, but on <<float>> rather than <<double>> values.
 113
 114 RETURNS
 115 The value of each Bessel function at <[x]> is returned.
 116
 117 PORTABILITY
 118 None of the Bessel functions are in ANSI C.
 119 */
 120
 121 /*
 122  * wrapper j0(double x), y0(double x)
 123  */
 124
 125 #include "fdlibm.h"
 126 #include <errno.h>
 127
 128 #ifndef _DOUBLE_IS_32BITS
 129
 130 #ifdef __STDC__
 131         double j0(double x)             /* wrapper j0 */
 132 #else
 133         double j0(x)                    /* wrapper j0 */
 134         double x;
 135 #endif
 136 {
 137 #ifdef _IEEE_LIBM
 138         return __ieee754_j0(x);
 139 #else
 140         struct exception exc;
 141         double z = __ieee754_j0(x);
 142         if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
 143         if(fabs(x)>X_TLOSS) {
 144             /* j0(|x|>X_TLOSS) */
 145             exc.type = TLOSS;
 146             exc.name = "j0";
 147             exc.err = 0;
 148             exc.arg1 = exc.arg2 = x;
 149             exc.retval = 0.0;
 150             if (_LIB_VERSION == _POSIX_)
 151                errno = ERANGE;
 152             else if (!matherr(&exc)) {
 153                errno = ERANGE;
 154             }
 155             if (exc.err != 0)
 156                errno = exc.err;
 157             return exc.retval;
 158         } else
 159             return z;
 160 #endif
 161 }
 162
 163 #ifdef __STDC__
 164         double y0(double x)             /* wrapper y0 */
 165 #else
 166         double y0(x)                    /* wrapper y0 */
 167         double x;
 168 #endif
 169 {
 170 #ifdef _IEEE_LIBM
 171         return __ieee754_y0(x);
 172 #else
 173         double z;
 174         struct exception exc;
 175         z = __ieee754_y0(x);
 176         if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
 177         if(x <= 0.0){
 178 #ifndef HUGE_VAL
 179 #define HUGE_VAL inf
 180             double inf = 0.0;
 181
 182             SET_HIGH_WORD(inf,0x7ff00000);      /* set inf to infinite */
 183 #endif
 184             /* y0(0) = -inf or y0(x<0) = NaN */
 185             exc.type = DOMAIN;  /* should be SING for IEEE y0(0) */
 186             exc.name = "y0";
 187             exc.err = 0;
 188             exc.arg1 = exc.arg2 = x;
 189             if (_LIB_VERSION == _SVID_)
 190                exc.retval = -HUGE;
 191             else
 192                exc.retval = -HUGE_VAL;
 193             if (_LIB_VERSION == _POSIX_)
 194                errno = EDOM;
 195             else if (!matherr(&exc)) {
 196                errno = EDOM;
 197             }
 198             if (exc.err != 0)
 199                errno = exc.err;
 200             return exc.retval;
 201         }
 202         if(x>X_TLOSS) {
 203             /* y0(x>X_TLOSS) */
 204             exc.type = TLOSS;
 205             exc.name = "y0";
 206             exc.err = 0;
 207             exc.arg1 = exc.arg2 = x;
 208             exc.retval = 0.0;
 209             if (_LIB_VERSION == _POSIX_)
 210                errno = ERANGE;
 211             else if (!matherr(&exc)) {
 212                errno = ERANGE;
 213             }
 214             if (exc.err != 0)
 215                errno = exc.err;
 216             return exc.retval;
 217         } else
 218             return z;
 219 #endif
 220 }
 221
 222 #endif /* defined(_DOUBLE_IS_32BITS) */
 223
 224
 225
 226
 227
 228
 229