yn.c
Go to the documentation of this file.
1 /* yn.c
2  *
3  * Bessel function of second kind of integer order
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, yn();
10  * int n;
11  *
12  * y = yn( n, x );
13  *
14  *
15  *
16  * DESCRIPTION:
17  *
18  * Returns Bessel function of order n, where n is a
19  * (possibly negative) integer.
20  *
21  * The function is evaluated by forward recurrence on
22  * n, starting with values computed by the routines
23  * y0() and y1().
24  *
25  * If n = 0 or 1 the routine for y0 or y1 is called
26  * directly.
27  *
28  *
29  *
30  * ACCURACY:
31  *
32  *
33  * Absolute error, except relative
34  * when y > 1:
35  * arithmetic domain # trials peak rms
36  * IEEE 0, 30 30000 3.4e-15 4.3e-16
37  *
38  *
39  * ERROR MESSAGES:
40  *
41  * message condition value returned
42  * yn singularity x = 0 INFINITY
43  * yn overflow INFINITY
44  *
45  * Spot checked against tables for x, n between 0 and 100.
46  *
47  */
48 
49 /*
50  * Cephes Math Library Release 2.8: June, 2000
51  * Copyright 1984, 1987, 2000 by Stephen L. Moshier
52  */
53 
54 #include "mconf.h"
55 extern double MAXLOG;
56 
57 double yn(int n, double x)
58 {
59  double an, anm1, anm2, r;
60  int k, sign;
61 
62  if (n < 0) {
63  n = -n;
64  if ((n & 1) == 0) /* -1**n */
65  sign = 1;
66  else
67  sign = -1;
68  }
69  else
70  sign = 1;
71 
72 
73  if (n == 0)
74  return (sign * y0(x));
75  if (n == 1)
76  return (sign * y1(x));
77 
78  /* test for overflow */
79  if (x == 0.0) {
80  sf_error("yn", SF_ERROR_SINGULAR, NULL);
81  return -INFINITY * sign;
82  }
83  else if (x < 0.0) {
84  sf_error("yn", SF_ERROR_DOMAIN, NULL);
85  return NAN;
86  }
87 
88  /* forward recurrence on n */
89 
90  anm2 = y0(x);
91  anm1 = y1(x);
92  k = 1;
93  r = 2 * k;
94  do {
95  an = r * anm1 / x - anm2;
96  anm2 = anm1;
97  anm1 = an;
98  r += 2.0;
99  ++k;
100  }
101  while (k < n);
102 
103 
104  return (sign * an);
105 }
MAXLOG
double MAXLOG
Definition: const.c:57
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
y0
double y0(double x)
Definition: j0.c:220
n
int n
Definition: BiCGSTAB_simple.cpp:1
y1
double y1(double x)
Definition: j1.c:199
yn
double yn(int n, double x)
Definition: yn.c:57
mconf.h
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16

gtsam
Author(s):
autogenerated on Sun Dec 22 2024 04:18:27