URL https://opencores.org/ocsvn/or1k/or1k/trunk

Subversion Repositories or1k

[/] [or1k/] [branches/] [newlib/] [newlib/] [newlib/] [libm/] [math/] [w_jn.c] - Blame information for rev 39

Go to most recent revision | Details | Compare with Previous | View Log


 
/* @(#)w_jn.c 5.1 93/09/24 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
 
/*
 * wrapper jn(int n, double x), yn(int n, double x)
 * floating point Bessel's function of the 1st and 2nd kind
 * of order n
 *
 * Special cases:
 *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
 *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
 * Note 2. About jn(n,x), yn(n,x)
 *      For n=0, j0(x) is called,
 *      for n=1, j1(x) is called,
 *      for n<x, forward recursion us used starting
 *      from values of j0(x) and j1(x).
 *      for n>x, a continued fraction approximation to
 *      j(n,x)/j(n-1,x) is evaluated and then backward
 *      recursion is used starting from a supposed value
 *      for j(n,x). The resulting value of j(0,x) is
 *      compared with the actual value to correct the
 *      supposed value of j(n,x).
 *
 *      yn(n,x) is similar in all respects, except
 *      that forward recursion is used for all
 *      values of n>1.
 *
 */
 
#include "fdlibm.h"
#include <errno.h>
 
#ifndef _DOUBLE_IS_32BITS
 
#ifdef __STDC__
        double jn(int n, double x)      /* wrapper jn */
#else
        double jn(n,x)                  /* wrapper jn */
        double x; int n;
#endif
{
#ifdef _IEEE_LIBM
        return __ieee754_jn(n,x);
#else
        double z;
        struct exception exc;
        z = __ieee754_jn(n,x);
        if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
        if(fabs(x)>X_TLOSS) {
            /* jn(|x|>X_TLOSS) */
            exc.type = TLOSS;
            exc.name = "jn";
            exc.retval = 0.0;
            if (_LIB_VERSION == _POSIX_)
                errno = ERANGE;
            else if (!matherr(&exc)) {
               errno = ERANGE;
            }
            if (exc.err != 0)
               errno = exc.err;
            return exc.retval;
        } else
            return z;
#endif
}
 
#ifdef __STDC__
        double yn(int n, double x)      /* wrapper yn */
#else
        double yn(n,x)                  /* wrapper yn */
        double x; int n;
#endif
{
#ifdef _IEEE_LIBM
        return __ieee754_yn(n,x);
#else
        double z;
        struct exception exc;
        z = __ieee754_yn(n,x);
        if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
        if(x <= 0.0){
            /* yn(n,0) = -inf or yn(x<0) = NaN */
#ifndef HUGE_VAL 
#define HUGE_VAL inf
            double inf = 0.0;
 
            SET_HIGH_WORD(inf,0x7ff00000);      /* set inf to infinite */
#endif
            exc.type = DOMAIN;  /* should be SING for IEEE */
            exc.name = "yn";
            if (_LIB_VERSION == _SVID_)
                exc.retval = -HUGE;
            else
                exc.retval = -HUGE_VAL;
            if (_LIB_VERSION == _POSIX_)
                errno = EDOM;
            else if (!matherr(&exc)) {
                errno = EDOM;
            }
            if (exc.err != 0)
               errno = exc.err;
            return exc.retval;
        }
        if(x>X_TLOSS) {
            /* yn(x>X_TLOSS) */
            exc.type = TLOSS;
            exc.name = "yn";
            exc.retval = 0.0;
            if (_LIB_VERSION == _POSIX_)
                errno = ERANGE;
            else if (!matherr(&exc)) {
                errno = ERANGE;
            }
            if (exc.err != 0)
               errno = exc.err;
            return exc.retval;
        } else
            return z;
#endif
}
 
#endif /* defined(_DOUBLE_IS_32BITS) */

Line No.	Rev	Author	Line
1	39	lampret
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.retval = 0.0;`
64			`if (_LIB_VERSION == _POSIX_)`
65			`errno = ERANGE;`
66			`else if (!matherr(&exc)) {`
67			`errno = ERANGE;`
68			`}`
69			`if (exc.err != 0)`
70			`errno = exc.err;`
71			`return exc.retval;`
72			`} else`
73			`return z;`
74			`#endif`
75			`}`
76
77			`#ifdef __STDC__`
78			`double yn(int n, double x) /* wrapper yn */`
79			`#else`
80			`double yn(n,x) /* wrapper yn */`
81			`double x; int n;`
82			`#endif`
83			`{`
84			`#ifdef _IEEE_LIBM`
85			`return __ieee754_yn(n,x);`
86			`#else`
87			`double z;`
88			`struct exception exc;`
89			`z = __ieee754_yn(n,x);`
90			`if(_LIB_VERSION == _IEEE_ \|\| isnan(x) ) return z;`
91			`if(x <= 0.0){`
92			`/* yn(n,0) = -inf or yn(x<0) = NaN */`
93			`#ifndef HUGE_VAL`
94			`#define HUGE_VAL inf`
95			`double inf = 0.0;`
96
97			`SET_HIGH_WORD(inf,0x7ff00000); /* set inf to infinite */`
98			`#endif`
99			`exc.type = DOMAIN; /* should be SING for IEEE */`
100			`exc.name = "yn";`
101			`if (_LIB_VERSION == _SVID_)`
102			`exc.retval = -HUGE;`
103			`else`
104			`exc.retval = -HUGE_VAL;`
105			`if (_LIB_VERSION == _POSIX_)`
106			`errno = EDOM;`
107			`else if (!matherr(&exc)) {`
108			`errno = EDOM;`
109			`}`
110			`if (exc.err != 0)`
111			`errno = exc.err;`
112			`return exc.retval;`
113			`}`
114			`if(x>X_TLOSS) {`
115			`/* yn(x>X_TLOSS) */`
116			`exc.type = TLOSS;`
117			`exc.name = "yn";`
118			`exc.retval = 0.0;`
119			`if (_LIB_VERSION == _POSIX_)`
120			`errno = ERANGE;`
121			`else if (!matherr(&exc)) {`
122			`errno = ERANGE;`
123			`}`
124			`if (exc.err != 0)`
125			`errno = exc.err;`
126			`return exc.retval;`
127			`} else`
128			`return z;`
129			`#endif`
130			`}`
131
132			`#endif /* defined(_DOUBLE_IS_32BITS) */`

Browse

Tools

Subversion Repositories or1k

[/] [or1k/] [branches/] [newlib/] [newlib/] [newlib/] [libm/] [math/] [w_jn.c] - Blame information for rev 39