Back to index

lightning-sunbird  0.9+nobinonly
w_jn.c
Go to the documentation of this file.
00001 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002  *
00003  * ***** BEGIN LICENSE BLOCK *****
00004  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00005  *
00006  * The contents of this file are subject to the Mozilla Public License Version
00007  * 1.1 (the "License"); you may not use this file except in compliance with
00008  * the License. You may obtain a copy of the License at
00009  * http://www.mozilla.org/MPL/
00010  *
00011  * Software distributed under the License is distributed on an "AS IS" basis,
00012  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00013  * for the specific language governing rights and limitations under the
00014  * License.
00015  *
00016  * The Original Code is Mozilla Communicator client code, released
00017  * March 31, 1998.
00018  *
00019  * The Initial Developer of the Original Code is
00020  * Sun Microsystems, Inc.
00021  * Portions created by the Initial Developer are Copyright (C) 1998
00022  * the Initial Developer. All Rights Reserved.
00023  *
00024  * Contributor(s):
00025  *
00026  * Alternatively, the contents of this file may be used under the terms of
00027  * either of the GNU General Public License Version 2 or later (the "GPL"),
00028  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00029  * in which case the provisions of the GPL or the LGPL are applicable instead
00030  * of those above. If you wish to allow use of your version of this file only
00031  * under the terms of either the GPL or the LGPL, and not to allow others to
00032  * use your version of this file under the terms of the MPL, indicate your
00033  * decision by deleting the provisions above and replace them with the notice
00034  * and other provisions required by the GPL or the LGPL. If you do not delete
00035  * the provisions above, a recipient may use your version of this file under
00036  * the terms of any one of the MPL, the GPL or the LGPL.
00037  *
00038  * ***** END LICENSE BLOCK ***** */
00039 
00040 /* @(#)w_jn.c 1.3 95/01/18 */
00041 /*
00042  * ====================================================
00043  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00044  *
00045  * Developed at SunSoft, a Sun Microsystems, Inc. business.
00046  * Permission to use, copy, modify, and distribute this
00047  * software is freely granted, provided that this notice 
00048  * is preserved.
00049  * ====================================================
00050  */
00051 
00052 /*
00053  * wrapper jn(int n, double x), yn(int n, double x)
00054  * floating point Bessel's function of the 1st and 2nd kind
00055  * of order n
00056  *          
00057  * Special cases:
00058  *     y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
00059  *     y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
00060  * Note 2. About jn(n,x), yn(n,x)
00061  *     For n=0, j0(x) is called,
00062  *     for n=1, j1(x) is called,
00063  *     for n<x, forward recursion us used starting
00064  *     from values of j0(x) and j1(x).
00065  *     for n>x, a continued fraction approximation to
00066  *     j(n,x)/j(n-1,x) is evaluated and then backward
00067  *     recursion is used starting from a supposed value
00068  *     for j(n,x). The resulting value of j(0,x) is
00069  *     compared with the actual value to correct the
00070  *     supposed value of j(n,x).
00071  *
00072  *     yn(n,x) is similar in all respects, except
00073  *     that forward recursion is used for all
00074  *     values of n>1.
00075  *     
00076  */
00077 
00078 #include "fdlibm.h"
00079 
00080 #ifdef __STDC__
00081        double fd_jn(int n, double x)      /* wrapper jn */
00082 #else
00083        double fd_jn(n,x)                  /* wrapper jn */
00084        double x; int n;
00085 #endif
00086 {
00087 #ifdef _IEEE_LIBM
00088        return __ieee754_jn(n,x);
00089 #else
00090        double z;
00091        z = __ieee754_jn(n,x);
00092        if(_LIB_VERSION == _IEEE_ || fd_isnan(x) ) return z;
00093        if(fd_fabs(x)>X_TLOSS) {
00094         int err;
00095            return __kernel_standard((double)n,x,38,&err); /* jn(|x|>X_TLOSS,n) */
00096        } else
00097            return z;
00098 #endif
00099 }
00100 
00101 #ifdef __STDC__
00102        double yn(int n, double x)  /* wrapper yn */
00103 #else
00104        double yn(n,x)                     /* wrapper yn */
00105        double x; int n;
00106 #endif
00107 {
00108 #ifdef _IEEE_LIBM
00109        return __ieee754_yn(n,x);
00110 #else
00111        double z;
00112     int err;
00113        z = __ieee754_yn(n,x);
00114        if(_LIB_VERSION == _IEEE_ || fd_isnan(x) ) return z;
00115         if(x <= 0.0){
00116                 if(x==0.0)
00117                     /* d= -one/(x-x); */
00118                     return __kernel_standard((double)n,x,12,&err);
00119                 else
00120                     /* d = zero/(x-x); */
00121                     return __kernel_standard((double)n,x,13,&err);
00122         }
00123        if(x>X_TLOSS) {
00124            return __kernel_standard((double)n,x,39,&err); /* yn(x>X_TLOSS,n) */
00125        } else
00126            return z;
00127 #endif
00128 }