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

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib/] [newlib/] [libm/] [math/] [w_jn.c] - Blame information for rev 56

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

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 56 joel 
            exc.err = 0;
64 
            exc.arg1 = n;
65 
            exc.arg2 = x;
66 39 lampret 
            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 56 joel 
            exc.err = 0;
105 
            exc.arg1 = n;
106 
            exc.arg2 = x;
107 39 lampret 
            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 56 joel 
            exc.err = 0;
125 
            exc.arg1 = n;
126 
            exc.arg2 = x;
127 39 lampret 
            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) */

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.