1//  Pieced together from Boost C++ and Cephes by
2//  Andreas Kloeckner (C) 2012
3//
4//  Pieces from:
5//
6//  Copyright (c) 2006 Xiaogang Zhang, John Maddock
7//  Use, modification and distribution are subject to the
8//  Boost Software License, Version 1.0. (See
9//  http://www.boost.org/LICENSE_1_0.txt)
10//
11// Cephes Math Library Release 2.8:  June, 2000
12// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
13// What you see here may be used freely, but it comes with no support or
14// guarantee.
15
16
17
18
19
20#pragma once
21
22#include <pyopencl-eval-tbl.cl>
23#include <pyopencl-bessel-j.cl>
24
25typedef double bessel_y_scalar_type;
26
27// {{{ bessel_y0
28
29__constant const bessel_y_scalar_type bessel_y0_P1[] = {
30     1.0723538782003176831e+11,
31    -8.3716255451260504098e+09,
32     2.0422274357376619816e+08,
33    -2.1287548474401797963e+06,
34     1.0102532948020907590e+04,
35    -1.8402381979244993524e+01,
36};
37__constant const bessel_y_scalar_type bessel_y0_Q1[] = {
38     5.8873865738997033405e+11,
39     8.1617187777290363573e+09,
40     5.5662956624278251596e+07,
41     2.3889393209447253406e+05,
42     6.6475986689240190091e+02,
43     1.0,
44};
45__constant const bessel_y_scalar_type bessel_y0_P2[] = {
46    -2.2213976967566192242e+13,
47    -5.5107435206722644429e+11,
48     4.3600098638603061642e+10,
49    -6.9590439394619619534e+08,
50     4.6905288611678631510e+06,
51    -1.4566865832663635920e+04,
52     1.7427031242901594547e+01,
53};
54__constant const bessel_y_scalar_type bessel_y0_Q2[] = {
55     4.3386146580707264428e+14,
56     5.4266824419412347550e+12,
57     3.4015103849971240096e+10,
58     1.3960202770986831075e+08,
59     4.0669982352539552018e+05,
60     8.3030857612070288823e+02,
61     1.0,
62};
63__constant const bessel_y_scalar_type bessel_y0_P3[] = {
64    -8.0728726905150210443e+15,
65     6.7016641869173237784e+14,
66    -1.2829912364088687306e+11,
67    -1.9363051266772083678e+11,
68     2.1958827170518100757e+09,
69    -1.0085539923498211426e+07,
70     2.1363534169313901632e+04,
71    -1.7439661319197499338e+01,
72};
73__constant const bessel_y_scalar_type bessel_y0_Q3[] = {
74     3.4563724628846457519e+17,
75     3.9272425569640309819e+15,
76     2.2598377924042897629e+13,
77     8.6926121104209825246e+10,
78     2.4727219475672302327e+08,
79     5.3924739209768057030e+05,
80     8.7903362168128450017e+02,
81     1.0,
82};
83__constant const bessel_y_scalar_type bessel_y0_PC[] = {
84     2.2779090197304684302e+04,
85     4.1345386639580765797e+04,
86     2.1170523380864944322e+04,
87     3.4806486443249270347e+03,
88     1.5376201909008354296e+02,
89     8.8961548424210455236e-01,
90};
91__constant const bessel_y_scalar_type bessel_y0_QC[] = {
92     2.2779090197304684318e+04,
93     4.1370412495510416640e+04,
94     2.1215350561880115730e+04,
95     3.5028735138235608207e+03,
96     1.5711159858080893649e+02,
97     1.0,
98};
99__constant const bessel_y_scalar_type bessel_y0_PS[] = {
100    -8.9226600200800094098e+01,
101    -1.8591953644342993800e+02,
102    -1.1183429920482737611e+02,
103    -2.2300261666214198472e+01,
104    -1.2441026745835638459e+00,
105    -8.8033303048680751817e-03,
106};
107__constant const bessel_y_scalar_type bessel_y0_QS[] = {
108     5.7105024128512061905e+03,
109     1.1951131543434613647e+04,
110     7.2642780169211018836e+03,
111     1.4887231232283756582e+03,
112     9.0593769594993125859e+01,
113     1.0,
114};
115
116bessel_y_scalar_type bessel_y0(bessel_y_scalar_type x)
117{
118    const bessel_y_scalar_type
119          x1  =  8.9357696627916752158e-01,
120          x2  =  3.9576784193148578684e+00,
121          x3  =  7.0860510603017726976e+00,
122          x11 =  2.280e+02,
123          x12 =  2.9519662791675215849e-03,
124          x21 =  1.0130e+03,
125          x22 =  6.4716931485786837568e-04,
126          x31 =  1.8140e+03,
127          x32 =  1.1356030177269762362e-04;
128
129    bessel_y_scalar_type value, factor, r, rc, rs;
130
131    if (x < 0)
132    {
133       //return policies::raise_domain_error<T>(function,
134       //    "Got x = %1% but x must be non-negative, complex result not supported.", x, pol);
135       return nan((uint)22);
136    }
137    if (x == 0)
138    {
139        return -DBL_MAX;
140    }
141    if (x <= 3)                       // x in (0, 3]
142    {
143        bessel_y_scalar_type y = x * x;
144        bessel_y_scalar_type z = 2 * log(x/x1) * bessel_j0(x) / M_PI;
145        r = boost_evaluate_rational(bessel_y0_P1, bessel_y0_Q1, y);
146        factor = (x + x1) * ((x - x11/256) - x12);
147        value = z + factor * r;
148    }
149    else if (x <= 5.5f)                  // x in (3, 5.5]
150    {
151        bessel_y_scalar_type y = x * x;
152        bessel_y_scalar_type z = 2 * log(x/x2) * bessel_j0(x) / M_PI;
153        r = boost_evaluate_rational(bessel_y0_P2, bessel_y0_Q2, y);
154        factor = (x + x2) * ((x - x21/256) - x22);
155        value = z + factor * r;
156    }
157    else if (x <= 8)                  // x in (5.5, 8]
158    {
159        bessel_y_scalar_type y = x * x;
160        bessel_y_scalar_type z = 2 * log(x/x3) * bessel_j0(x) / M_PI;
161        r = boost_evaluate_rational(bessel_y0_P3, bessel_y0_Q3, y);
162        factor = (x + x3) * ((x - x31/256) - x32);
163        value = z + factor * r;
164    }
165    else                                // x in (8, \infty)
166    {
167        bessel_y_scalar_type y = 8 / x;
168        bessel_y_scalar_type y2 = y * y;
169        bessel_y_scalar_type z = x - 0.25f * M_PI;
170        rc = boost_evaluate_rational(bessel_y0_PC, bessel_y0_QC, y2);
171        rs = boost_evaluate_rational(bessel_y0_PS, bessel_y0_QS, y2);
172        factor = sqrt(2 / (x * M_PI));
173        value = factor * (rc * sin(z) + y * rs * cos(z));
174    }
175
176    return value;
177}
178
179// }}}
180
181// {{{ bessel_y1
182
183__constant const bessel_y_scalar_type bessel_y1_P1[] = {
184     4.0535726612579544093e+13,
185     5.4708611716525426053e+12,
186    -3.7595974497819597599e+11,
187     7.2144548214502560419e+09,
188    -5.9157479997408395984e+07,
189     2.2157953222280260820e+05,
190    -3.1714424660046133456e+02,
191};
192__constant const bessel_y_scalar_type bessel_y1_Q1[] = {
193     3.0737873921079286084e+14,
194     4.1272286200406461981e+12,
195     2.7800352738690585613e+10,
196     1.2250435122182963220e+08,
197     3.8136470753052572164e+05,
198     8.2079908168393867438e+02,
199     1.0,
200};
201__constant const bessel_y_scalar_type bessel_y1_P2[] = {
202     1.1514276357909013326e+19,
203    -5.6808094574724204577e+18,
204    -2.3638408497043134724e+16,
205     4.0686275289804744814e+15,
206    -5.9530713129741981618e+13,
207     3.7453673962438488783e+11,
208    -1.1957961912070617006e+09,
209     1.9153806858264202986e+06,
210    -1.2337180442012953128e+03,
211};
212__constant const bessel_y_scalar_type bessel_y1_Q2[] = {
213     5.3321844313316185697e+20,
214     5.6968198822857178911e+18,
215     3.0837179548112881950e+16,
216     1.1187010065856971027e+14,
217     3.0221766852960403645e+11,
218     6.3550318087088919566e+08,
219     1.0453748201934079734e+06,
220     1.2855164849321609336e+03,
221     1.0,
222};
223__constant const bessel_y_scalar_type bessel_y1_PC[] = {
224    -4.4357578167941278571e+06,
225    -9.9422465050776411957e+06,
226    -6.6033732483649391093e+06,
227    -1.5235293511811373833e+06,
228    -1.0982405543459346727e+05,
229    -1.6116166443246101165e+03,
230     0.0,
231};
232__constant const bessel_y_scalar_type bessel_y1_QC[] = {
233    -4.4357578167941278568e+06,
234    -9.9341243899345856590e+06,
235    -6.5853394797230870728e+06,
236    -1.5118095066341608816e+06,
237    -1.0726385991103820119e+05,
238    -1.4550094401904961825e+03,
239     1.0,
240};
241__constant const bessel_y_scalar_type bessel_y1_PS[] = {
242     3.3220913409857223519e+04,
243     8.5145160675335701966e+04,
244     6.6178836581270835179e+04,
245     1.8494262873223866797e+04,
246     1.7063754290207680021e+03,
247     3.5265133846636032186e+01,
248     0.0,
249};
250__constant const bessel_y_scalar_type bessel_y1_QS[] = {
251     7.0871281941028743574e+05,
252     1.8194580422439972989e+06,
253     1.4194606696037208929e+06,
254     4.0029443582266975117e+05,
255     3.7890229745772202641e+04,
256     8.6383677696049909675e+02,
257     1.0,
258};
259
260bessel_y_scalar_type bessel_y1(bessel_y_scalar_type x)
261{
262    const bessel_y_scalar_type
263      x1  =  2.1971413260310170351e+00,
264      x2  =  5.4296810407941351328e+00,
265      x11 =  5.620e+02,
266      x12 =  1.8288260310170351490e-03,
267      x21 =  1.3900e+03,
268      x22 = -6.4592058648672279948e-06
269    ;
270    bessel_y_scalar_type value, factor, r, rc, rs;
271
272    if (x <= 0)
273    {
274      // domain error
275      return nan((uint)22);
276    }
277    if (x <= 4)                       // x in (0, 4]
278    {
279        bessel_y_scalar_type y = x * x;
280        bessel_y_scalar_type z = 2 * log(x/x1) * bessel_j1(x) / M_PI;
281        r = boost_evaluate_rational(bessel_y1_P1, bessel_y1_Q1, y);
282        factor = (x + x1) * ((x - x11/256) - x12) / x;
283        value = z + factor * r;
284    }
285    else if (x <= 8)                  // x in (4, 8]
286    {
287        bessel_y_scalar_type y = x * x;
288        bessel_y_scalar_type z = 2 * log(x/x2) * bessel_j1(x) / M_PI;
289        r = boost_evaluate_rational(bessel_y1_P2, bessel_y1_Q2, y);
290        factor = (x + x2) * ((x - x21/256) - x22) / x;
291        value = z + factor * r;
292    }
293    else                                // x in (8, \infty)
294    {
295        bessel_y_scalar_type y = 8 / x;
296        bessel_y_scalar_type y2 = y * y;
297        bessel_y_scalar_type z = x - 0.75f * M_PI;
298        rc = boost_evaluate_rational(bessel_y1_PC, bessel_y1_QC, y2);
299        rs = boost_evaluate_rational(bessel_y1_PS, bessel_y1_QS, y2);
300        factor = sqrt(2 / (x * M_PI));
301        value = factor * (rc * sin(z) + y * rs * cos(z));
302    }
303
304    return value;
305}
306
307// }}}
308
309// {{{ bessel_yn
310
311bessel_y_scalar_type bessel_yn_small_z(int n, bessel_y_scalar_type z, bessel_y_scalar_type* scale)
312{
313   //
314   // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
315   //
316   // Note that when called we assume that x < epsilon and n is a positive integer.
317   //
318   // BOOST_ASSERT(n >= 0);
319   // BOOST_ASSERT((z < policies::get_epsilon<T, Policy>()));
320
321   if(n == 0)
322   {
323      return (2 / M_PI) * (log(z / 2) +  M_E);
324   }
325   else if(n == 1)
326   {
327      return (z / M_PI) * log(z / 2)
328         - 2 / (M_PI * z)
329         - (z / (2 * M_PI)) * (1 - 2 * M_E);
330   }
331   else if(n == 2)
332   {
333      return (z * z) / (4 * M_PI) * log(z / 2)
334         - (4 / (M_PI * z * z))
335         - ((z * z) / (8 * M_PI)) * (3./2 - 2 * M_E);
336   }
337   else
338   {
339      bessel_y_scalar_type p = pow(z / 2, (bessel_y_scalar_type) n);
340      bessel_y_scalar_type result = -((tgamma((bessel_y_scalar_type) n) / M_PI));
341      if(p * DBL_MAX < result)
342      {
343         bessel_y_scalar_type div = DBL_MAX / 8;
344         result /= div;
345         *scale /= div;
346         if(p * DBL_MAX < result)
347         {
348            return -DBL_MAX;
349         }
350      }
351      return result / p;
352   }
353}
354
355
356
357
358bessel_y_scalar_type bessel_yn(int n, bessel_y_scalar_type x)
359{
360    //BOOST_MATH_STD_USING
361    bessel_y_scalar_type value, factor, current, prev;
362
363    //using namespace boost::math::tools;
364
365    if ((x == 0) && (n == 0))
366    {
367       return -DBL_MAX;
368    }
369    if (x <= 0)
370    {
371       //return policies::raise_domain_error<T>(function,
372            //"Got x = %1%, but x must be > 0, complex result not supported.", x, pol);
373       return nan((uint)22);
374    }
375
376    //
377    // Reflection comes first:
378    //
379    if (n < 0)
380    {
381        factor = (n & 0x1) ? -1 : 1;  // Y_{-n}(z) = (-1)^n Y_n(z)
382        n = -n;
383    }
384    else
385    {
386        factor = 1;
387    }
388
389    if(x < DBL_EPSILON)
390    {
391       bessel_y_scalar_type scale = 1;
392       value = bessel_yn_small_z(n, x, &scale);
393       if(DBL_MAX * fabs(scale) < fabs(value))
394          return copysign((bessel_y_scalar_type) 1, scale) * copysign((bessel_y_scalar_type) 1, value) * DBL_MAX;
395       value /= scale;
396    }
397    else if (n == 0)
398    {
399        value = bessel_y0(x);
400    }
401    else if (n == 1)
402    {
403        value = factor * bessel_y1(x);
404    }
405    else
406    {
407       prev = bessel_y0(x);
408       current = bessel_y1(x);
409       int k = 1;
410       // BOOST_ASSERT(k < n);
411       do
412       {
413           bessel_y_scalar_type fact = 2 * k / x;
414           if((DBL_MAX - fabs(prev)) / fact < fabs(current))
415           {
416              prev /= current;
417              factor /= current;
418              current = 1;
419           }
420           value = fact * current - prev;
421           prev = current;
422           current = value;
423           ++k;
424       }
425       while(k < n);
426       if(fabs(DBL_MAX * factor) < fabs(value))
427          return sign(value) * sign(value) * DBL_MAX;
428       value /= factor;
429    }
430    return value;
431}
432
433// }}}
434
435// vim: fdm=marker
436