struve.c
Go to the documentation of this file.
1 /*
2  * Compute the Struve function.
3  *
4  * Notes
5  * -----
6  *
7  * We use three expansions for the Struve function discussed in [1]:
8  *
9  * - power series
10  * - expansion in Bessel functions
11  * - asymptotic large-z expansion
12  *
13  * Rounding errors are estimated based on the largest terms in the sums.
14  *
15  * ``struve_convergence.py`` plots the convergence regions of the different
16  * expansions.
17  *
18  * (i)
19  *
20  * Looking at the error in the asymptotic expansion, one finds that
21  * it's not worth trying if z ~> 0.7 * v + 12 for v > 0.
22  *
23  * (ii)
24  *
25  * The Bessel function expansion tends to fail for |z| >~ |v| and is not tried
26  * there.
27  *
28  * For Struve H it covers the quadrant v > z where the power series may fail to
29  * produce reasonable results.
30  *
31  * (iii)
32  *
33  * The three expansions together cover for Struve H the region z > 0, v real.
34  *
35  * They also cover Struve L, except that some loss of precision may occur around
36  * the transition region z ~ 0.7 |v|, v < 0, |v| >> 1 where the function changes
37  * rapidly.
38  *
39  * (iv)
40  *
41  * The power series is evaluated in double-double precision. This fixes accuracy
42  * issues in Struve H for |v| << |z| before the asymptotic expansion kicks in.
43  * Moreover, it improves the Struve L behavior for negative v.
44  *
45  *
46  * References
47  * ----------
48  * [1] NIST Digital Library of Mathematical Functions
49  * https://dlmf.nist.gov/11
50  */
51 
52 /*
53  * Copyright (C) 2013 Pauli Virtanen
54  *
55  * Redistribution and use in source and binary forms, with or without
56  * modification, are permitted provided that the following conditions are met:
57  *
58  * a. Redistributions of source code must retain the above copyright notice,
59  * this list of conditions and the following disclaimer.
60  * b. Redistributions in binary form must reproduce the above copyright
61  * notice, this list of conditions and the following disclaimer in the
62  * documentation and/or other materials provided with the distribution.
63  * c. Neither the name of Enthought nor the names of the SciPy Developers
64  * may be used to endorse or promote products derived from this software
65  * without specific prior written permission.
66  *
67  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
68  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
69  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
70  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS
71  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
72  * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
73  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
74  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
75  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
76  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
77  * THE POSSIBILITY OF SUCH DAMAGE.
78  */
79 
80 #include "mconf.h"
81 #include "dd_real.h"
82 
83 // #include "amos_wrappers.h"
84 
85 #define STRUVE_MAXITER 10000
86 #define SUM_EPS 1e-16 /* be sure we are in the tail of the sum */
87 #define SUM_TINY 1e-100
88 #define GOOD_EPS 1e-12
89 #define ACCEPTABLE_EPS 1e-7
90 #define ACCEPTABLE_ATOL 1e-300
91 
92 #define MIN(a, b) ((a) < (b) ? (a) : (b))
93 
94 double struve_power_series(double v, double x, int is_h, double *err);
95 double struve_asymp_large_z(double v, double z, int is_h, double *err);
96 double struve_bessel_series(double v, double z, int is_h, double *err);
97 
98 static double bessel_y(double v, double x);
99 static double bessel_j(double v, double x);
100 static double struve_hl(double v, double x, int is_h);
101 
102 double struve_h(double v, double z)
103 {
104  return struve_hl(v, z, 1);
105 }
106 
107 double struve_l(double v, double z)
108 {
109  return struve_hl(v, z, 0);
110 }
111 
112 static double struve_hl(double v, double z, int is_h)
113 {
114  double value[4], err[4], tmp;
115  int n;
116 
117  if (z < 0) {
118  n = v;
119  if (v == n) {
120  tmp = (n % 2 == 0) ? -1 : 1;
121  return tmp * struve_hl(v, -z, is_h);
122  }
123  else {
124  return NAN;
125  }
126  }
127  else if (z == 0) {
128  if (v < -1) {
129  return gammasgn(v + 1.5) * INFINITY;
130  }
131  else if (v == -1) {
132  return 2 / sqrt(M_PI) / Gamma(0.5);
133  }
134  else {
135  return 0;
136  }
137  }
138 
139  n = -v - 0.5;
140  if (n == -v - 0.5 && n > 0) {
141  if (is_h) {
142  return (n % 2 == 0 ? 1 : -1) * bessel_j(n + 0.5, z);
143  }
144  else {
145  return iv(n + 0.5, z);
146  }
147  }
148 
149  /* Try the asymptotic expansion */
150  if (z >= 0.7*v + 12) {
151  value[0] = struve_asymp_large_z(v, z, is_h, &err[0]);
152  if (err[0] < GOOD_EPS * fabs(value[0])) {
153  return value[0];
154  }
155  }
156  else {
157  err[0] = INFINITY;
158  }
159 
160  /* Try power series */
161  value[1] = struve_power_series(v, z, is_h, &err[1]);
162  if (err[1] < GOOD_EPS * fabs(value[1])) {
163  return value[1];
164  }
165 
166  /* Try bessel series */
167  if (fabs(z) < fabs(v) + 20) {
168  value[2] = struve_bessel_series(v, z, is_h, &err[2]);
169  if (err[2] < GOOD_EPS * fabs(value[2])) {
170  return value[2];
171  }
172  }
173  else {
174  err[2] = INFINITY;
175  }
176 
177  /* Return the best of the three, if it is acceptable */
178  n = 0;
179  if (err[1] < err[n]) n = 1;
180  if (err[2] < err[n]) n = 2;
181  if (err[n] < ACCEPTABLE_EPS * fabs(value[n]) || err[n] < ACCEPTABLE_ATOL) {
182  return value[n];
183  }
184 
185  /* Maybe it really is an overflow? */
186  tmp = -lgam(v + 1.5) + (v + 1)*log(z/2);
187  if (!is_h) {
188  tmp = fabs(tmp);
189  }
190  if (tmp > 700) {
191  sf_error("struve", SF_ERROR_OVERFLOW, NULL);
192  return INFINITY * gammasgn(v + 1.5);
193  }
194 
195  /* Failure */
196  sf_error("struve", SF_ERROR_NO_RESULT, NULL);
197  return NAN;
198 }
199 
200 
201 /*
202  * Power series for Struve H and L
203  * https://dlmf.nist.gov/11.2.1
204  *
205  * Starts to converge roughly at |n| > |z|
206  */
207 double struve_power_series(double v, double z, int is_h, double *err)
208 {
209  int n, sgn;
210  double term, sum, maxterm, scaleexp, tmp;
211  double2 cterm, csum, cdiv, z2, c2v, ctmp;
212 
213  if (is_h) {
214  sgn = -1;
215  }
216  else {
217  sgn = 1;
218  }
219 
220  tmp = -lgam(v + 1.5) + (v + 1)*log(z/2);
221  if (tmp < -600 || tmp > 600) {
222  /* Scale exponent to postpone underflow/overflow */
223  scaleexp = tmp/2;
224  tmp -= scaleexp;
225  }
226  else {
227  scaleexp = 0;
228  }
229 
230  term = 2 / sqrt(M_PI) * exp(tmp) * gammasgn(v + 1.5);
231  sum = term;
232  maxterm = 0;
233 
234  cterm = dd_create_d(term);
235  csum = dd_create_d(sum);
236  z2 = dd_create_d(sgn*z*z);
237  c2v = dd_create_d(2*v);
238 
239  for (n = 0; n < STRUVE_MAXITER; ++n) {
240  /* cdiv = (3 + 2*n) * (3 + 2*n + 2*v)) */
241  cdiv = dd_create_d(3 + 2*n);
242  ctmp = dd_create_d(3 + 2*n);
243  ctmp = dd_add(ctmp, c2v);
244  cdiv = dd_mul(cdiv, ctmp);
245 
246  /* cterm *= z2 / cdiv */
247  cterm = dd_mul(cterm, z2);
248  cterm = dd_div(cterm, cdiv);
249 
250  csum = dd_add(csum, cterm);
251 
252  term = dd_to_double(cterm);
253  sum = dd_to_double(csum);
254 
255  if (fabs(term) > maxterm) {
256  maxterm = fabs(term);
257  }
258  if (fabs(term) < SUM_TINY * fabs(sum) || term == 0 || !isfinite(sum)) {
259  break;
260  }
261  }
262 
263  *err = fabs(term) + fabs(maxterm) * 1e-22;
264 
265  if (scaleexp != 0) {
266  sum *= exp(scaleexp);
267  *err *= exp(scaleexp);
268  }
269 
270  if (sum == 0 && term == 0 && v < 0 && !is_h) {
271  /* Spurious underflow */
272  *err = INFINITY;
273  return NAN;
274  }
275 
276  return sum;
277 }
278 
279 
280 /*
281  * Bessel series
282  * https://dlmf.nist.gov/11.4.19
283  */
284 double struve_bessel_series(double v, double z, int is_h, double *err)
285 {
286  int n;
287  double term, cterm, sum, maxterm;
288 
289  if (is_h && v < 0) {
290  /* Works less reliably in this region */
291  *err = INFINITY;
292  return NAN;
293  }
294 
295  sum = 0;
296  maxterm = 0;
297 
298  cterm = sqrt(z / (2*M_PI));
299 
300  for (n = 0; n < STRUVE_MAXITER; ++n) {
301  if (is_h) {
302  term = cterm * bessel_j(n + v + 0.5, z) / (n + 0.5);
303  cterm *= z/2 / (n + 1);
304  }
305  else {
306  term = cterm * iv(n + v + 0.5, z) / (n + 0.5);
307  cterm *= -z/2 / (n + 1);
308  }
309  sum += term;
310  if (fabs(term) > maxterm) {
311  maxterm = fabs(term);
312  }
313  if (fabs(term) < SUM_EPS * fabs(sum) || term == 0 || !isfinite(sum)) {
314  break;
315  }
316  }
317 
318  *err = fabs(term) + fabs(maxterm) * 1e-16;
319 
320  /* Account for potential underflow of the Bessel functions */
321  *err += 1e-300 * fabs(cterm);
322 
323  return sum;
324 }
325 
326 
327 /*
328  * Large-z expansion for Struve H and L
329  * https://dlmf.nist.gov/11.6.1
330  */
331 double struve_asymp_large_z(double v, double z, int is_h, double *err)
332 {
333  int n, sgn, maxiter;
334  double term, sum, maxterm;
335  double m;
336 
337  if (is_h) {
338  sgn = -1;
339  }
340  else {
341  sgn = 1;
342  }
343 
344  /* Asymptotic expansion divergenge point */
345  m = z/2;
346  if (m <= 0) {
347  maxiter = 0;
348  }
349  else if (m > STRUVE_MAXITER) {
350  maxiter = STRUVE_MAXITER;
351  }
352  else {
353  maxiter = (int)m;
354  }
355  if (maxiter == 0) {
356  *err = INFINITY;
357  return NAN;
358  }
359 
360  if (z < v) {
361  /* Exclude regions where our error estimation fails */
362  *err = INFINITY;
363  return NAN;
364  }
365 
366  /* Evaluate sum */
367  term = -sgn / sqrt(M_PI) * exp(-lgam(v + 0.5) + (v - 1) * log(z/2)) * gammasgn(v + 0.5);
368  sum = term;
369  maxterm = 0;
370 
371  for (n = 0; n < maxiter; ++n) {
372  term *= sgn * (1 + 2*n) * (1 + 2*n - 2*v) / (z*z);
373  sum += term;
374  if (fabs(term) > maxterm) {
375  maxterm = fabs(term);
376  }
377  if (fabs(term) < SUM_EPS * fabs(sum) || term == 0 || !isfinite(sum)) {
378  break;
379  }
380  }
381 
382  if (is_h) {
383  sum += bessel_y(v, z);
384  }
385  else {
386  sum += iv(v, z);
387  }
388 
389  /*
390  * This error estimate is strictly speaking valid only for
391  * n > v - 0.5, but numerical results indicate that it works
392  * reasonably.
393  */
394  *err = fabs(term) + fabs(maxterm) * 1e-16;
395 
396  return sum;
397 }
398 
399 
400 static double bessel_y(double v, double x)
401 {
402  return cbesy_wrap_real(v, x);
403 }
404 
405 static double bessel_j(double v, double x)
406 {
407  return cbesj_wrap_real(v, x);
408 }
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
dd_real.h
bessel_y
static double bessel_y(double v, double x)
Definition: struve.c:400
struve_hl
static double struve_hl(double v, double x, int is_h)
Definition: struve.c:112
double2
Definition: dd_real.h:74
gammasgn
double gammasgn(double x)
Definition: gammasgn.c:3
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
dd_mul
static double2 dd_mul(const double2 a, const double2 b)
Definition: dd_real_idefs.h:404
SF_ERROR_NO_RESULT
@ SF_ERROR_NO_RESULT
Definition: sf_error.h:15
STRUVE_MAXITER
#define STRUVE_MAXITER
Definition: struve.c:85
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
bessel_j
static double bessel_j(double v, double x)
Definition: struve.c:405
Gamma
double Gamma(double x)
Definition: gamma.c:160
SUM_TINY
#define SUM_TINY
Definition: struve.c:87
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
dd_div
static double2 dd_div(const double2 a, const double2 b)
Definition: dd_real_idefs.h:479
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
SUM_EPS
#define SUM_EPS
Definition: struve.c:86
struve_bessel_series
double struve_bessel_series(double v, double z, int is_h, double *err)
Definition: struve.c:284
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
n
int n
Definition: BiCGSTAB_simple.cpp:1
struve_asymp_large_z
double struve_asymp_large_z(double v, double z, int is_h, double *err)
Definition: struve.c:331
GOOD_EPS
#define GOOD_EPS
Definition: struve.c:88
dd_to_double
static double dd_to_double(const double2 a)
Definition: dd_real_idefs.h:90
isfinite
#define isfinite(X)
Definition: main.h:95
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
dd_add
static double2 dd_add(const double2 a, const double2 b)
Definition: dd_real_idefs.h:343
struve_h
double struve_h(double v, double z)
Definition: struve.c:102
ACCEPTABLE_EPS
#define ACCEPTABLE_EPS
Definition: struve.c:89
struve_power_series
double struve_power_series(double v, double x, int is_h, double *err)
Definition: struve.c:207
mconf.h
struve_l
double struve_l(double v, double z)
Definition: struve.c:107
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
dd_create_d
static double2 dd_create_d(double hi)
Definition: dd_real_idefs.h:149
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
iv
double iv(double v, double x)
Definition: scipy_iv.c:80
ACCEPTABLE_ATOL
#define ACCEPTABLE_ATOL
Definition: struve.c:90
NULL
#define NULL
Definition: ccolamd.c:609
M_PI
#define M_PI
Definition: mconf.h:117
test_callbacks.value
value
Definition: test_callbacks.py:160
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
z2
static const Unit3 z2
Definition: testRotateFactor.cpp:43


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:36:27