123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534 |
- #include "slatec-internal.hpp"
- static integer const c__15 = 15;
- static integer const c__5 = 5;
- static integer const c__3 = 3;
- static integer const c__2 = 2;
- static integer const c__1 = 1;
- int dbsknu_(double const *x, double const *fnu, integer const *kode,
- integer const *n, double *y, integer *nz)
- {
-
- static double const x1 = 2.;
- static double const x2 = 17.;
- static double const pi = 3.14159265358979;
- static double const rthpi = 1.2533141373155;
- static double const cc[8] = { .577215664901533,-.0420026350340952,
- -.0421977345555443,.007218943246663,-2.152416741149e-4,
- -2.01348547807e-5,1.133027232e-6,6.116095e-9 };
-
- integer i__1;
-
- double a[160], b[160], f;
- integer i__, j, k;
- double p, q, s, a1, a2, g1, g2, p1, p2, s1, s2 = 0., t1, t2, fc, ak, bk,
- ck = 0., dk, fk;
- integer kk;
- double cx;
- integer nn;
- double ex, tm, pt, st, rx, fhs, fks, dnu, fmu;
- integer inu;
- double sqk, tol, smu, dnu2, coef, elim, flrx;
- integer iflag, koded;
- double etest;
-
- --y;
-
- kk = -i1mach_(15);
- elim = (kk * d1mach_(5) - 3.) * 2.303;
- ak = d1mach_(3);
- tol = max(ak,1e-15);
- if (*x <= 0.) {
- goto L350;
- }
- if (*fnu < 0.) {
- goto L360;
- }
- if (*kode < 1 || *kode > 2) {
- goto L370;
- }
- if (*n < 1) {
- goto L380;
- }
- *nz = 0;
- iflag = 0;
- koded = *kode;
- rx = 2. / *x;
- inu = (integer) (*fnu + .5);
- dnu = *fnu - inu;
- dnu2 = 0.;
- if (abs(dnu) == .5) {
- goto L120;
- }
- if (abs(dnu) < tol) {
- goto L10;
- }
- dnu2 = dnu * dnu;
- L10:
- if (*x > x1) {
- goto L120;
- }
- a1 = 1. - dnu;
- a2 = dnu + 1.;
- t1 = 1. / dgamma_(&a1);
- t2 = 1. / dgamma_(&a2);
- if (abs(dnu) > .1) {
- goto L40;
- }
- s = cc[0];
- ak = 1.;
- for (k = 2; k <= 8; ++k) {
- ak *= dnu2;
- tm = cc[k - 1] * ak;
- s += tm;
- if (abs(tm) < tol) {
- goto L30;
- }
- }
- L30:
- g1 = -s;
- goto L50;
- L40:
- g1 = (t1 - t2) / (dnu + dnu);
- L50:
- g2 = (t1 + t2) * .5;
- smu = 1.;
- fc = 1.;
- flrx = log(rx);
- fmu = dnu * flrx;
- if (dnu == 0.) {
- goto L60;
- }
- fc = dnu * pi;
- fc /= sin(fc);
- if (fmu != 0.) {
- smu = sinh(fmu) / fmu;
- }
- L60:
- f = fc * (g1 * cosh(fmu) + g2 * flrx * smu);
- fc = exp(fmu);
- p = fc * .5 / t2;
- q = .5 / (fc * t1);
- ak = 1.;
- ck = 1.;
- bk = 1.;
- s1 = f;
- s2 = p;
- if (inu > 0 || *n > 1) {
- goto L90;
- }
- if (*x < tol) {
- goto L80;
- }
- cx = *x * *x * .25;
- L70:
- f = (ak * f + p + q) / (bk - dnu2);
- p /= ak - dnu;
- q /= ak + dnu;
- ck = ck * cx / ak;
- t1 = ck * f;
- s1 += t1;
- bk = bk + ak + ak + 1.;
- ak += 1.;
- s = abs(t1) / (abs(s1) + 1.);
- if (s > tol) {
- goto L70;
- }
- L80:
- y[1] = s1;
- if (koded == 1) {
- return 0;
- }
- y[1] = s1 * exp(*x);
- return 0;
- L90:
- if (*x < tol) {
- goto L110;
- }
- cx = *x * *x * .25;
- L100:
- f = (ak * f + p + q) / (bk - dnu2);
- p /= ak - dnu;
- q /= ak + dnu;
- ck = ck * cx / ak;
- t1 = ck * f;
- s1 += t1;
- t2 = ck * (p - ak * f);
- s2 += t2;
- bk = bk + ak + ak + 1.;
- ak += 1.;
- s = abs(t1) / (abs(s1) + 1.) + abs(t2) / (abs(s2) + 1.);
- if (s > tol) {
- goto L100;
- }
- L110:
- s2 *= rx;
- if (koded == 1) {
- goto L170;
- }
- f = exp(*x);
- s1 *= f;
- s2 *= f;
- goto L170;
- L120:
- coef = rthpi / sqrt(*x);
- if (koded == 2) {
- goto L130;
- }
- if (*x > elim) {
- goto L330;
- }
- coef *= exp(-(*x));
- L130:
- if (abs(dnu) == .5) {
- goto L340;
- }
- if (*x > x2) {
- goto L280;
- }
- etest = cos(pi * dnu) / (pi * *x * tol);
- fks = 1.;
- fhs = .25;
- fk = 0.;
- ck = *x + *x + 2.;
- p1 = 0.;
- p2 = 1.;
- k = 0;
- L140:
- ++k;
- fk += 1.;
- ak = (fhs - dnu2) / (fks + fk);
- bk = ck / (fk + 1.);
- pt = p2;
- p2 = bk * p2 - ak * p1;
- p1 = pt;
- a[k - 1] = ak;
- b[k - 1] = bk;
- ck += 2.;
- fks = fks + fk + fk + 1.;
- fhs = fhs + fk + fk;
- if (etest > fk * p1) {
- goto L140;
- }
- kk = k;
- s = 1.;
- p1 = 0.;
- p2 = 1.;
- i__1 = k;
- for (i__ = 1; i__ <= i__1; ++i__) {
- pt = p2;
- p2 = (b[kk - 1] * p2 - p1) / a[kk - 1];
- p1 = pt;
- s += p2;
- --kk;
- }
- s1 = coef * (p2 / s);
- if (inu > 0 || *n > 1) {
- goto L160;
- }
- goto L200;
- L160:
- s2 = s1 * (*x + dnu + .5 - p1 / p2) / *x;
- L170:
- ck = (dnu + dnu + 2.) / *x;
- if (*n == 1) {
- --inu;
- }
- if (inu > 0) {
- goto L180;
- }
- if (*n > 1) {
- goto L200;
- }
- s1 = s2;
- goto L200;
- L180:
- i__1 = inu;
- for (i__ = 1; i__ <= i__1; ++i__) {
- st = s2;
- s2 = ck * s2 + s1;
- s1 = st;
- ck += rx;
- }
- if (*n == 1) {
- s1 = s2;
- }
- L200:
- if (iflag == 1) {
- goto L220;
- }
- y[1] = s1;
- if (*n == 1) {
- return 0;
- }
- y[2] = s2;
- if (*n == 2) {
- return 0;
- }
- i__1 = *n;
- for (i__ = 3; i__ <= i__1; ++i__) {
- y[i__] = ck * y[i__ - 1] + y[i__ - 2];
- ck += rx;
- }
- return 0;
- L220:
- s = -(*x) + log(s1);
- y[1] = 0.;
- *nz = 1;
- if (s < -elim) {
- goto L230;
- }
- y[1] = exp(s);
- *nz = 0;
- L230:
- if (*n == 1) {
- return 0;
- }
- s = -(*x) + log(s2);
- y[2] = 0.;
- ++(*nz);
- if (s < -elim) {
- goto L240;
- }
- --(*nz);
- y[2] = exp(s);
- L240:
- if (*n == 2) {
- return 0;
- }
- kk = 2;
- if (*nz < 2) {
- goto L260;
- }
- i__1 = *n;
- for (i__ = 3; i__ <= i__1; ++i__) {
- kk = i__;
- st = s2;
- s2 = ck * s2 + s1;
- s1 = st;
- ck += rx;
- s = -(*x) + log(s2);
- ++(*nz);
- y[i__] = 0.;
- if (s < -elim) {
- goto L250;
- }
- y[i__] = exp(s);
- --(*nz);
- goto L260;
- L250:
- ;
- }
- return 0;
- L260:
- if (kk == *n) {
- return 0;
- }
- s2 = s2 * ck + s1;
- ck += rx;
- ++kk;
- y[kk] = exp(-(*x) + log(s2));
- if (kk == *n) {
- return 0;
- }
- ++kk;
- i__1 = *n;
- for (i__ = kk; i__ <= i__1; ++i__) {
- y[i__] = ck * y[i__ - 1] + y[i__ - 2];
- ck += rx;
- }
- return 0;
- L280:
- nn = 2;
- if (inu == 0 && *n == 1) {
- nn = 1;
- }
- dnu2 = dnu + dnu;
- fmu = 0.;
- if (abs(dnu2) < tol) {
- goto L290;
- }
- fmu = dnu2 * dnu2;
- L290:
- ex = *x * 8.;
- s2 = 0.;
- i__1 = nn;
- for (k = 1; k <= i__1; ++k) {
- s1 = s2;
- s = 1.;
- ak = 0.;
- ck = 1.;
- sqk = 1.;
- dk = ex;
- for (j = 1; j <= 30; ++j) {
- ck = ck * (fmu - sqk) / dk;
- s += ck;
- dk += ex;
- ak += 8.;
- sqk += ak;
- if (abs(ck) < tol) {
- goto L310;
- }
- }
- L310:
- s2 = s * coef;
- fmu = fmu + dnu * 8. + 4.;
- }
- if (nn > 1) {
- goto L170;
- }
- s1 = s2;
- goto L200;
- L330:
- koded = 2;
- iflag = 1;
- goto L120;
- L340:
- s1 = coef;
- s2 = coef;
- goto L170;
- L350:
- xermsg_("SLATEC", "DBSKNU", "X NOT GREATER THAN ZERO", &c__2, &c__1, (
- ftnlen)6, (ftnlen)6, (ftnlen)23);
- return 0;
- L360:
- xermsg_("SLATEC", "DBSKNU", "FNU NOT ZERO OR POSITIVE", &c__2, &c__1, (
- ftnlen)6, (ftnlen)6, (ftnlen)24);
- return 0;
- L370:
- xermsg_("SLATEC", "DBSKNU", "KODE NOT 1 OR 2", &c__2, &c__1, (ftnlen)6, (
- ftnlen)6, (ftnlen)15);
- return 0;
- L380:
- xermsg_("SLATEC", "DBSKNU", "N NOT GREATER THAN 0", &c__2, &c__1, (ftnlen)
- 6, (ftnlen)6, (ftnlen)20);
- return 0;
- }
|