/*							j1l.c
 *
 *	Bessel function of order one
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, j1l();
 *
 * y = j1l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of order one of the argument.
 *
 * The domain is divided into the intervals [0, 9] and
 * (9, infinity). In the first interval the rational approximation
 * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) x P8(x^2) / Q8(x^2),
 * where r, s, t are the first three zeros of the function.
 * In the second interval the expansion is in terms of the
 * modulus M1(x) = sqrt(J1(x)^2 + Y1(x)^2) and phase  P1(x)
 * = atan(Y1(x)/J1(x)).  M1 is approximated by sqrt(1/x)P7(1/x)/Q8(1/x).
 * The approximation to j1 is M1 * cos(x -  3 pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)).
 *
 *
 * ACCURACY:
 *
 *                      Absolute error:
 * arithmetic   domain      # trials      peak         rms
 *    IEEE      0, 30       10000       1.7e-19      5.0e-20
 *
 *
 */
/*							y1l.c
 *
 *	Bessel function of the second kind, order zero
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, y1l();
 *
 * y = y1l( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of the second kind, of order
 * zero, of the argument.
 *
 * The domain is divided into the intervals [0, 4.5>, [4.5,9> and
 * [9, infinity). In the first interval a rational approximation
 * R(x) is employed to compute y0(x)  = R(x) + 2/pi * log(x) * j0(x).
 *
 * In the second interval, the approximation is
 *     (x - p)(x - q)(x - r)(x - s)P9(x)/Q10(x)
 * where p, q, r, s are zeros of y1(x).
 *
 * The third interval uses the same approximations to modulus
 * and phase as j1(x), whence y1(x) = modulus * sin(phase).
 *
 * ACCURACY:
 *
 *  Absolute error, when y0(x) < 1; else relative error:
 *
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       6000        2.4e-19     5.4e-20
 *
 */

/* Copyright 1994 by Stephen L. Moshier (moshier@world.std.com).  */

#include <math.h>

extern long double __polevll ( long double, long double *, int );
extern long double __p1evll ( long double, long double *, int );

/*
j1(x) = (x^2-r0^2)(x^2-r1^2)(x^2-r2^2) x P(x**2)/Q(x**2)
Relative error
n=8, d=8
Peak error =  2e-21
*/
static long double j1n[] = {
-2.634697796221277628971935439379100448941E-4L,
 9.313297622796327912616990544807729349816E-1L,
-1.462801427977939339087861462016763117614E3L,
 1.320001295393312144946720253523751495344E6L,
-7.411832711954540428424579722022697274252E8L,
 2.626500686552841932403375456823890030733E11L,
-5.682630730221834709333595924862737385674E13L,
 6.800062979972634469823465019593889535036E15L,
-3.414700974444745667483143593080298025264E17L
 };
static long double j1d[] = {
/* 1.000000000000000000000000000000000000000E0 */
 2.952679519729437457329616170307633853371E3L,
 4.787239263438296747732226800464625857306E6L,
 5.375447329578075439197694561828270959740E9L,
 4.468662138862678294901755275226730274220E12L,
 2.769597563759616070852381975767135596290E15L,
 1.233678068848311511940047189841753732303E18L,
 3.573258746896955995244858202885402209284E20L,
 5.107790455161415784613110686083256035323E22L
};
/*
sqrt(j0^2(1/x^2) + y0^2(1/x^2)) = z P(z**2)/Q(z**2)
z(x) = 1/sqrt(x)
Relative error
n=7, d=8
Peak error =  1.35e=20
Relative error spread =  9.0e0
*/
static long double modulusn[] = {
-5.041742205078442098873946599818859411841E0L,
 3.918474430130242177354760505592118707545E-1L,
 2.527521168680500659056387541731679323679E0L,
 7.172146812845906480742580751162395880700E0L,
 2.859499532295180940060097685686785747839E0L,
 1.014671139779858141347461961442252531091E0L,
 1.255798064266130869131769135327120168414E-1L,
 1.596507617085714650238013674302102214586E-2L
};
static long double modulusd[] = {
/* 1.000000000000000000000000000000000000000E0 */
-6.233092094568239317498213232200456266608E0L,
-9.214128701852838347001714615803119950789E-1L,
 2.531772200570435289832000302297822691539E0L,
 8.755081357265851765639996579311407253229E0L,
 3.554340386955608261462733259897195898764E0L,
 1.267949948774331531237417342523472412319E0L,
 1.573909467558180942219276563149698450074E-1L,
 2.000925566825407466160190735195649280593E-2L
 };
/*
atan(y1(x)/j1(x))  =  x - 3pi/4 + z P(z**2)/Q(z**2)
z(x) = 1/x
Absolute error
n=5, d=6
Peak error =  4.83e-21
Relative error spread =  1.9e0
*/
static long double phasen[] = {
 2.010456367705144783933016123752706035207E0L,
 1.587378144541918176658121923411891505725E0L,
 2.682837461073751055565177542176613982515E-1L,
 1.472572645054468815027121251529744762234E-2L,
 2.884976126715926258585500789700568171816E-4L,
 1.708502235134706284899170935145399268617E-6L
};
static long double phased[] = {
/* 1.000000000000000000000000000000000000000E0 */
 6.809332495854873089362409327980929501353E0L,
 4.518597941618813112664962212253836700838E0L,
 7.320149039410806471100853405098354700190E-1L,
 3.960155028960712309813713456399687003948E-2L,
 7.713202197319040439861232615648837506269E-4L,
 4.556005960359216767984395285780473266643E-6L
 };

#define JZ1 1.4681970642123893257219777768630106989681587E1L
#define JZ2 4.9218456321694603670267082846388138932425412E1L
#define JZ3 1.0349945389513658033222363253561305574983502E2L

#define THPIO4L  2.35619449019234492885L

long double j1l(long double x)
{
long double xx, y, z, modulus, phase;

xx = x * x;
if( xx < 81.0L )
  {
    y = (xx - JZ1) * (xx - JZ2) * (xx - JZ3);
    y *= x * __polevll( xx, j1n, 8 ) / __p1evll( xx, j1d, 8 );
    return y;
  }

y = fabsl(x);
xx = 1.0/xx;
phase = __polevll( xx, phasen, 5 ) / __p1evll( xx, phased, 6 );

z = 1.0/y;
modulus = __polevll( z, modulusn, 7 ) / __p1evll( z, modulusd, 8 );

y = modulus * cosl( y -  THPIO4L + z*phase) / sqrtl(y);
if( x < 0 )
  y = -y;
return y;
}

static long double y1n[] = {
-1.288901054372751879530799223223505257692E5L,
 9.914315981558815369372336678918231933103E7L,
-2.906793378120403577273915213389829871402E10L,
 3.954354656937677136265837177266774549541E12L,
-2.445982226888344140153571265516582079837E14L,
 5.685362960165615942886131910189352872528E15L,
-2.158855258453711703120399853735329573110E16L
};
static long double y1d[] = {
/* 1.000000000000000000000000000000000000000E0 */
 8.926354644853231136072991636781412786125E2L,
 4.679841933793707979658529184983101722128E5L,
 1.775133253792677466650594382943712081101E8L,
 5.089532584184822833415809008343104067480E10L,
 1.076474894829072923244442357703372832195E13L,
 1.525917240904692387994338842633536676896E15L,
 1.101136026928555260167699409793359571793E17L
};
static long double y159n[] = {
-6.806634906054210550895796232243478627879E-1L,
 4.306669585790359450531735079332477032239E1L,
-9.230477746767243316013573089219930917796E2L,
 6.171186628598134035236603420969585848071E3L,
 2.096869860275353982829414221776143235491E4L,
-1.238961670382216747944146253091813185559E5L,
-1.781314136808997406108546510088125743955E6L,
-1.803400156074242435454306040864758160815E6L,
-1.155761550219364178626921880155041761746E6L,
 3.112221202330688509818373491140459519285E5L
};
static long double y159d[] = {
/* 1.000000000000000000000000000000000000000E0 */
-6.181482377814679766978324825945117512812E1L,
 2.238187927382180589098927722621122192486E3L,
-5.225317824142187494325968315136322493497E4L,
 9.217235006983512475118388815131670860571E5L,
-1.183757638771741974520896241943158149865E7L,
 1.208072488974110742912058706805732230596E8L,
-8.193431077523942651172530209520018834293E8L,
 4.282669747880013349980886241568066930618E9L,
-1.171523459555524458807562674565375794753E9L,
 1.078445545755236785691969932840418530054E8L
};

#define MAXNUML 1.189731495357231765021263853E4932L
#define TWOOPI 6.36619772367581343075535E-1L
#define THPIO4 2.35619449019234492885L
#define Y1Z1 2.1971413260310170351490335626989662730530183E0L
#define Y1Z2 5.4296810407941351327720051908525841965837575E0L
#define Y1Z3 8.5960058683311689264296061801639678511029216E0L
#define Y1Z4 1.1749154830839881243399421939922350714301166E1L

long double y1l(long double x)
{
long double xx, y, z, modulus, phase;

if( x < 0.0 )
  {
    return (-MAXNUML);
  }
z = 1.0/x;
xx = x * x;
if( xx < 81.0L )
  {
    if( xx < 20.25L )
      {
	y = TWOOPI * (logl(x) * j1l(x) - z);
	y += x * __polevll( xx, y1n, 6 ) / __p1evll( xx, y1d, 7 );
      }
    else
      {
	y = (x - Y1Z1)*(x - Y1Z2)*(x - Y1Z3)*(x - Y1Z4);
	y *= __polevll( x, y159n, 9 ) / __p1evll( x, y159d, 10 );
      }
    return y;
  }

xx = 1.0/xx;
phase = __polevll( xx, phasen, 5 ) / __p1evll( xx, phased, 6 );

modulus = __polevll( z, modulusn, 7 ) / __p1evll( z, modulusd, 8 );

z = modulus * sinl( x -  THPIO4L + z*phase) / sqrtl(x);
return z;
}