123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106 |
- /* yn.c
- *
- * Bessel function of second kind of integer order
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, yn();
- * int n;
- *
- * y = yn( n, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of order n, where n is a
- * (possibly negative) integer.
- *
- * The function is evaluated by forward recurrence on
- * n, starting with values computed by the routines
- * y0() and y1().
- *
- * If n = 0 or 1 the routine for y0 or y1 is called
- * directly.
- *
- *
- *
- * ACCURACY:
- *
- *
- * Absolute error, except relative
- * when y > 1:
- * arithmetic domain # trials peak rms
- * DEC 0, 30 2200 2.9e-16 5.3e-17
- * IEEE 0, 30 30000 3.4e-15 4.3e-16
- *
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * yn singularity x = 0 MAXNUM
- * yn overflow MAXNUM
- *
- * Spot checked against tables for x, n between 0 and 100.
- *
- */
- /*
- Cephes Math Library Release 2.8: June, 2000
- Copyright 1984, 1987, 2000 by Stephen L. Moshier
- */
- #include "mconf.h"
- #ifdef ANSIPROT
- extern double y0(double);
- extern double y1(double);
- extern double log(double);
- #else
- double y0(), y1(), log();
- #endif
- extern double MAXNUM, MAXLOG;
- double yn(n, x) int n;
- double x;
- {
- double an, anm1, anm2, r;
- int k, sign;
- if (n < 0) {
- n = -n;
- if ((n & 1) == 0) /* -1**n */
- sign = 1;
- else
- sign = -1;
- } else
- sign = 1;
- if (n == 0)
- return (sign * y0(x));
- if (n == 1)
- return (sign * y1(x));
- /* test for overflow */
- if (x <= 0.0) {
- mtherr("yn", SING);
- return (-MAXNUM);
- }
- /* forward recurrence on n */
- anm2 = y0(x);
- anm1 = y1(x);
- k = 1;
- r = 2 * k;
- do {
- an = r * anm1 / x - anm2;
- anm2 = anm1;
- anm1 = an;
- r += 2.0;
- ++k;
- } while (k < n);
- return (sign * an);
- }
|