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__ j0(double x)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__ y0(double x)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 230