123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207 |
- /*
- * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
- * See the copyright notice in the ACK home directory, in the file "Copyright".
- *
- * Author: Ceriel J.H. Jacobs
- */
- /* $Id$ */
- #include <math.h>
- #include <errno.h>
- extern int errno;
- static double
- P1(x)
- double x;
- {
- /* P1(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
- /* Hart & Cheney # 6755 */
- static double p[9] = {
- 0.1000000000000000000000000489e+01,
- 0.5581663300347182292169450071e+01,
- 0.1100186625131173123750501118e+02,
- 0.9727139359130463694593683431e+01,
- 0.4060011483142278994462590992e+01,
- 0.7742832212665311906917358099e+00,
- 0.6021617752811098752098248630e-01,
- 0.1482350677236405118074646993e-02,
- 0.6094215148131061431667573909e-05
- };
- static double q[9] = {
- 0.9999999999999999999999999999e+00,
- 0.5579832245659682292169922224e+01,
- 0.1099168447731617288972771040e+02,
- 0.9707206835125961446797916892e+01,
- 0.4042610016540342097334497865e+01,
- 0.7671965204303836019508430169e+00,
- 0.5893258668794493100786371406e-01,
- 0.1393993644981256852404222530e-02,
- 0.4585597769784750669754696825e-05
- };
- double zsq = 64.0/(x*x);
- return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
- }
- static double
- Q1(x)
- double x;
- {
- /* Q1(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
- /* Hart & Cheney # 7157 */
- /* Probably typerror in Hart & Cheney; it sais:
- Q1(x) = x*P(z*z)/Q(z*z)
- */
- static double p[9] = {
- 0.4687499999999999999999995275e-01,
- 0.3302394516691663879252493748e+00,
- 0.8456888491208195767613862428e+00,
- 0.1008551084218946085420665147e+01,
- 0.5973407972399900690521296181e+00,
- 0.1737697433393258207540273097e+00,
- 0.2303862814819568573893610740e-01,
- 0.1171224207976250587945594946e-02,
- 0.1486418220337492918307904804e-04
- };
- static double q[10] = {
- 0.9999999999999999999999999999e+00,
- 0.7049380763213049609070823421e+01,
- 0.1807129960468949760845562209e+02,
- 0.2159171174362827330505421695e+02,
- 0.1283239297740546866114600499e+02,
- 0.3758349275324260869598403931e+01,
- 0.5055985453754739528620657666e+00,
- 0.2665604326323907148063400439e-01,
- 0.3821140353404633025596424652e-03,
- 0.3206696590241261037875154062e-06
- };
- double zsq = 64.0/(x*x);
- return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
- }
- static double
- smallj1(x)
- double x;
- {
- /* J1(x) = x*P(x*x)/Q(x*x) for x in [0,8] */
- /* Hart & Cheney # 6054 */
- static double p[10] = {
- 0.1921176307760798128049021316e+25,
- -0.2226092031387396254771375773e+24,
- 0.7894463902082476734673226741e+22,
- -0.1269424373753606065436561036e+21,
- 0.1092152214043184787101134641e+19,
- -0.5454629264396819144157448868e+16,
- 0.1634659487571284628830445048e+14,
- -0.2909662785381647825756152444e+11,
- 0.2853433451054763915026471449e+08,
- -0.1197705712815379389149134705e+05
- };
- static double q[10] = {
- 0.3842352615521596256098041912e+25,
- 0.3507567066272028105798868716e+23,
- 0.1611334311633414344007062889e+21,
- 0.4929612313959850319632645381e+18,
- 0.1117536965288162684489793105e+16,
- 0.1969278625584719037168592923e+13,
- 0.2735606122949877990248154504e+10,
- 0.2940957355049651347475558106e+07,
- 0.2274736606126590905134610965e+04,
- 0.1000000000000000000000000000e+01
- };
- double xsq = x*x;
- return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
- }
- double
- j1(x)
- double x;
- {
- /* Use J1(x) = sqrt(2/(pi*x))*(P1(x)*cos(X1)-Q1(x)*sin(X1))
- where X1 = x - 3*pi/4 for |x| > 8.
- Use J1(-x) = -J1(x).
- Use direct approximation of smallj1 for |x| <= 8.
- */
- extern double sqrt(), sin(), cos();
- int negative = x < 0.0;
- if (negative) x = -x;
- if (x > 8.0) {
- double X1 = x - (M_PI - M_PI_4);
- x = sqrt(M_2_PI/x)*(P1(x)*cos(X1) - Q1(x)*sin(X1));
- }
- else x = smallj1(x);
- if (negative) return -x;
- return x;
- }
- static double
- smally1_bar(x)
- double x;
- {
- /* Y1(x) = Y1BAR(x)+(2/pi)*(J1(x)ln(x) - 1/x)
- Approximation of Y1BAR for 0 <= x <= 8:
- Y1BAR(x) = x*P(x*x)/Q(x*x)
- Hart & Cheney # 6449
- */
- static double p[10] = {
- -0.5862655424363443992938931700e+24,
- 0.1570668341992328458208364904e+24,
- -0.7351681299005467428400402479e+22,
- 0.1390658785759080111485190942e+21,
- -0.1339544201526785345938109179e+19,
- 0.7290257386242270629526344379e+16,
- -0.2340575603057015935501295099e+14,
- 0.4411516199185230690878878903e+11,
- -0.4542128738770213026987060358e+08,
- 0.1988612563465350530472715888e+05
- };
- static double q[10] = {
- 0.2990279721605116022908679994e+25,
- 0.2780285010357803058127175655e+23,
- 0.1302687474507355553192845146e+21,
- 0.4071330372239164349602952937e+18,
- 0.9446611865086570116528399283e+15,
- 0.1707657951197456205887347694e+13,
- 0.2440358986882941823431612517e+10,
- 0.2708852767034077697963790196e+07,
- 0.2174361138333330803617969305e+04,
- 0.1000000000000000000000000000e+01
- };
- double xsq = x*x;
- return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
- }
- double
- y1(x)
- double x;
- {
- extern double sqrt(), sin(), cos(), log();
- if (x <= 0.0) {
- errno = EDOM;
- return -HUGE;
- }
- if (x > 8.0) {
- double X1 = x - (M_PI - M_PI_4);
- return sqrt(M_2_PI/x) * (P1(x)*sin(X1)+Q1(x)*cos(X1));
- }
- return smally1_bar(x) + M_2_PI*(j1(x)*log(x) - 1/x);
- }
|