scipy_iv.c
Go to the documentation of this file.
1 /* iv.c
2  *
3  * Modified Bessel function of noninteger order
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double v, x, y, iv();
10  *
11  * y = iv( v, x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Returns modified Bessel function of order v of the
18  * argument. If x is negative, v must be integer valued.
19  *
20  */
21 /* iv.c */
22 /* Modified Bessel function of noninteger order */
23 /* If x < 0, then v must be an integer. */
24 
25 
26 /*
27  * Parts of the code are copyright:
28  *
29  * Cephes Math Library Release 2.8: June, 2000
30  * Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
31  *
32  * And other parts:
33  *
34  * Copyright (c) 2006 Xiaogang Zhang
35  * Use, modification and distribution are subject to the
36  * Boost Software License, Version 1.0.
37  *
38  * Boost Software License - Version 1.0 - August 17th, 2003
39  *
40  * Permission is hereby granted, free of charge, to any person or
41  * organization obtaining a copy of the software and accompanying
42  * documentation covered by this license (the "Software") to use, reproduce,
43  * display, distribute, execute, and transmit the Software, and to prepare
44  * derivative works of the Software, and to permit third-parties to whom the
45  * Software is furnished to do so, all subject to the following:
46  *
47  * The copyright notices in the Software and this entire statement,
48  * including the above license grant, this restriction and the following
49  * disclaimer, must be included in all copies of the Software, in whole or
50  * in part, and all derivative works of the Software, unless such copies or
51  * derivative works are solely in the form of machine-executable object code
52  * generated by a source language processor.
53  *
54  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
55  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
56  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
57  * NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE
58  * DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER LIABILITY,
59  * WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
60  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
61  * SOFTWARE.
62  *
63  * And the rest are:
64  *
65  * Copyright (C) 2009 Pauli Virtanen
66  * Distributed under the same license as Scipy.
67  *
68  */
69 
70 #include "mconf.h"
71 #include <float.h>
72 #include <stdlib.h>
73 
74 extern double MACHEP;
75 
76 static double iv_asymptotic(double v, double x);
77 static void ikv_asymptotic_uniform(double v, double x, double *Iv, double *Kv);
78 static void ikv_temme(double v, double x, double *Iv, double *Kv);
79 
80 double iv(double v, double x)
81 {
82  int sign;
83  double t, ax, res;
84 
85  if (isnan(v) || isnan(x)) {
86  return NAN;
87  }
88 
89  /* If v is a negative integer, invoke symmetry */
90  t = floor(v);
91  if (v < 0.0) {
92  if (t == v) {
93  v = -v; /* symmetry */
94  t = -t;
95  }
96  }
97  /* If x is negative, require v to be an integer */
98  sign = 1;
99  if (x < 0.0) {
100  if (t != v) {
101  sf_error("iv", SF_ERROR_DOMAIN, NULL);
102  return (NAN);
103  }
104  if (v != 2.0 * floor(v / 2.0)) {
105  sign = -1;
106  }
107  }
108 
109  /* Avoid logarithm singularity */
110  if (x == 0.0) {
111  if (v == 0.0) {
112  return 1.0;
113  }
114  if (v < 0.0) {
116  return INFINITY;
117  }
118  else
119  return 0.0;
120  }
121 
122  ax = fabs(x);
123  if (fabs(v) > 50) {
124  /*
125  * Uniform asymptotic expansion for large orders.
126  *
127  * This appears to overflow slightly later than the Boost
128  * implementation of Temme's method.
129  */
131  }
132  else {
133  /* Otherwise: Temme's method */
134  ikv_temme(v, ax, &res, NULL);
135  }
136  res *= sign;
137  return res;
138 }
139 
140 
141 /*
142  * Compute Iv from (AMS5 9.7.1), asymptotic expansion for large |z|
143  * Iv ~ exp(x)/sqrt(2 pi x) ( 1 + (4*v*v-1)/8x + (4*v*v-1)(4*v*v-9)/8x/2! + ...)
144  */
145 static double iv_asymptotic(double v, double x)
146 {
147  double mu;
148  double sum, term, prefactor, factor;
149  int k;
150 
151  prefactor = exp(x) / sqrt(2 * M_PI * x);
152 
153  if (prefactor == INFINITY) {
154  return prefactor;
155  }
156 
157  mu = 4 * v * v;
158  sum = 1.0;
159  term = 1.0;
160  k = 1;
161 
162  do {
163  factor = (mu - (2 * k - 1) * (2 * k - 1)) / (8 * x) / k;
164  if (k > 100) {
165  /* didn't converge */
166  sf_error("iv(iv_asymptotic)", SF_ERROR_NO_RESULT, NULL);
167  break;
168  }
169  term *= -factor;
170  sum += term;
171  ++k;
172  } while (fabs(term) > MACHEP * fabs(sum));
173  return sum * prefactor;
174 }
175 
176 
177 /*
178  * Uniform asymptotic expansion factors, (AMS5 9.3.9; AMS5 9.3.10)
179  *
180  * Computed with:
181  * --------------------
182  import numpy as np
183  t = np.poly1d([1,0])
184  def up1(p):
185  return .5*t*t*(1-t*t)*p.deriv() + 1/8. * ((1-5*t*t)*p).integ()
186  us = [np.poly1d([1])]
187  for k in range(10):
188  us.append(up1(us[-1]))
189  n = us[-1].order
190  for p in us:
191  print "{" + ", ".join(["0"]*(n-p.order) + map(repr, p)) + "},"
192  print "N_UFACTORS", len(us)
193  print "N_UFACTOR_TERMS", us[-1].order + 1
194  * --------------------
195  */
196 #define N_UFACTORS 11
197 #define N_UFACTOR_TERMS 31
199  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
200  0, 0, 0, 0, 0, 0, 0, 1},
201  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
202  0, 0, 0, 0, -0.20833333333333334, 0.0, 0.125, 0.0},
203  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
204  0, 0.3342013888888889, 0.0, -0.40104166666666669, 0.0, 0.0703125, 0.0,
205  0.0},
206  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
207  -1.0258125964506173, 0.0, 1.8464626736111112, 0.0,
208  -0.89121093750000002, 0.0, 0.0732421875, 0.0, 0.0, 0.0},
209  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
210  4.6695844234262474, 0.0, -11.207002616222995, 0.0, 8.78912353515625,
211  0.0, -2.3640869140624998, 0.0, 0.112152099609375, 0.0, 0.0, 0.0, 0.0},
212  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -28.212072558200244, 0.0,
213  84.636217674600744, 0.0, -91.818241543240035, 0.0, 42.534998745388457,
214  0.0, -7.3687943594796312, 0.0, 0.22710800170898438, 0.0, 0.0, 0.0,
215  0.0, 0.0},
216  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 212.5701300392171, 0.0,
217  -765.25246814118157, 0.0, 1059.9904525279999, 0.0,
218  -699.57962737613275, 0.0, 218.19051174421159, 0.0,
219  -26.491430486951554, 0.0, 0.57250142097473145, 0.0, 0.0, 0.0, 0.0,
220  0.0, 0.0},
221  {0, 0, 0, 0, 0, 0, 0, 0, 0, -1919.4576623184068, 0.0,
222  8061.7221817373083, 0.0, -13586.550006434136, 0.0, 11655.393336864536,
223  0.0, -5305.6469786134048, 0.0, 1200.9029132163525, 0.0,
224  -108.09091978839464, 0.0, 1.7277275025844574, 0.0, 0.0, 0.0, 0.0, 0.0,
225  0.0, 0.0},
226  {0, 0, 0, 0, 0, 0, 20204.291330966149, 0.0, -96980.598388637503, 0.0,
227  192547.0012325315, 0.0, -203400.17728041555, 0.0, 122200.46498301747,
228  0.0, -41192.654968897557, 0.0, 7109.5143024893641, 0.0,
229  -493.915304773088, 0.0, 6.074042001273483, 0.0, 0.0, 0.0, 0.0, 0.0,
230  0.0, 0.0, 0.0},
231  {0, 0, 0, -242919.18790055133, 0.0, 1311763.6146629769, 0.0,
232  -2998015.9185381061, 0.0, 3763271.2976564039, 0.0,
233  -2813563.2265865342, 0.0, 1268365.2733216248, 0.0,
234  -331645.17248456361, 0.0, 45218.768981362737, 0.0,
235  -2499.8304818112092, 0.0, 24.380529699556064, 0.0, 0.0, 0.0, 0.0, 0.0,
236  0.0, 0.0, 0.0, 0.0},
237  {3284469.8530720375, 0.0, -19706819.11843222, 0.0, 50952602.492664628,
238  0.0, -74105148.211532637, 0.0, 66344512.274729028, 0.0,
239  -37567176.660763353, 0.0, 13288767.166421819, 0.0,
240  -2785618.1280864552, 0.0, 308186.40461266245, 0.0,
241  -13886.089753717039, 0.0, 110.01714026924674, 0.0, 0.0, 0.0, 0.0, 0.0,
242  0.0, 0.0, 0.0, 0.0, 0.0}
243 };
244 
245 
246 /*
247  * Compute Iv, Kv from (AMS5 9.7.7 + 9.7.8), asymptotic expansion for large v
248  */
249 static void ikv_asymptotic_uniform(double v, double x,
250  double *i_value, double *k_value)
251 {
252  double i_prefactor, k_prefactor;
253  double t, t2, eta, z;
254  double i_sum, k_sum, term, divisor;
255  int k, n;
256  int sign = 1;
257 
258  if (v < 0) {
259  /* Negative v; compute I_{-v} and K_{-v} and use (AMS 9.6.2) */
260  sign = -1;
261  v = -v;
262  }
263 
264  z = x / v;
265  t = 1 / sqrt(1 + z * z);
266  t2 = t * t;
267  eta = sqrt(1 + z * z) + log(z / (1 + 1 / t));
268 
269  i_prefactor = sqrt(t / (2 * M_PI * v)) * exp(v * eta);
270  i_sum = 1.0;
271 
272  k_prefactor = sqrt(M_PI * t / (2 * v)) * exp(-v * eta);
273  k_sum = 1.0;
274 
275  divisor = v;
276  for (n = 1; n < N_UFACTORS; ++n) {
277  /*
278  * Evaluate u_k(t) with Horner's scheme;
279  * (using the knowledge about which coefficients are zero)
280  */
281  term = 0;
282  for (k = N_UFACTOR_TERMS - 1 - 3 * n;
283  k < N_UFACTOR_TERMS - n; k += 2) {
284  term *= t2;
285  term += asymptotic_ufactors[n][k];
286  }
287  for (k = 1; k < n; k += 2) {
288  term *= t2;
289  }
290  if (n % 2 == 1) {
291  term *= t;
292  }
293 
294  /* Sum terms */
295  term /= divisor;
296  i_sum += term;
297  k_sum += (n % 2 == 0) ? term : -term;
298 
299  /* Check convergence */
300  if (fabs(term) < MACHEP) {
301  break;
302  }
303 
304  divisor *= v;
305  }
306 
307  if (fabs(term) > 1e-3 * fabs(i_sum)) {
308  /* Didn't converge */
309  sf_error("ikv_asymptotic_uniform", SF_ERROR_NO_RESULT, NULL);
310  }
311  if (fabs(term) > MACHEP * fabs(i_sum)) {
312  /* Some precision lost */
313  sf_error("ikv_asymptotic_uniform", SF_ERROR_LOSS, NULL);
314  }
315 
316  if (k_value != NULL) {
317  /* symmetric in v */
318  *k_value = k_prefactor * k_sum;
319  }
320 
321  if (i_value != NULL) {
322  if (sign == 1) {
323  *i_value = i_prefactor * i_sum;
324  }
325  else {
326  /* (AMS 9.6.2) */
327  *i_value = (i_prefactor * i_sum
328  + (2 / M_PI) * sin(M_PI * v) * k_prefactor * k_sum);
329  }
330  }
331 }
332 
333 
334 /*
335  * The following code originates from the Boost C++ library,
336  * from file `boost/math/special_functions/detail/bessel_ik.hpp`,
337  * converted from C++ to C.
338  */
339 
340 #ifdef DEBUG
341 #define BOOST_ASSERT(a) assert(a)
342 #else
343 #define BOOST_ASSERT(a)
344 #endif
345 
346 /*
347  * Modified Bessel functions of the first and second kind of fractional order
348  *
349  * Calculate K(v, x) and K(v+1, x) by method analogous to
350  * Temme, Journal of Computational Physics, vol 21, 343 (1976)
351  */
352 static int temme_ik_series(double v, double x, double *K, double *K1)
353 {
354  double f, h, p, q, coef, sum, sum1, tolerance;
355  double a, b, c, d, sigma, gamma1, gamma2;
356  unsigned long k;
357  double gp;
358  double gm;
359 
360 
361  /*
362  * |x| <= 2, Temme series converge rapidly
363  * |x| > 2, the larger the |x|, the slower the convergence
364  */
365  BOOST_ASSERT(fabs(x) <= 2);
366  BOOST_ASSERT(fabs(v) <= 0.5f);
367 
368  gp = gamma(v + 1) - 1;
369  gm = gamma(-v + 1) - 1;
370 
371  a = log(x / 2);
372  b = exp(v * a);
373  sigma = -a * v;
374  c = fabs(v) < MACHEP ? 1 : sin(M_PI * v) / (v * M_PI);
375  d = fabs(sigma) < MACHEP ? 1 : sinh(sigma) / sigma;
376  gamma1 = fabs(v) < MACHEP ? -SCIPY_EULER : (0.5f / v) * (gp - gm) * c;
377  gamma2 = (2 + gp + gm) * c / 2;
378 
379  /* initial values */
380  p = (gp + 1) / (2 * b);
381  q = (1 + gm) * b / 2;
382  f = (cosh(sigma) * gamma1 + d * (-a) * gamma2) / c;
383  h = p;
384  coef = 1;
385  sum = coef * f;
386  sum1 = coef * h;
387 
388  /* series summation */
389  tolerance = MACHEP;
390  for (k = 1; k < MAXITER; k++) {
391  f = (k * f + p + q) / (k * k - v * v);
392  p /= k - v;
393  q /= k + v;
394  h = p - k * f;
395  coef *= x * x / (4 * k);
396  sum += coef * f;
397  sum1 += coef * h;
398  if (fabs(coef * f) < fabs(sum) * tolerance) {
399  break;
400  }
401  }
402  if (k == MAXITER) {
403  sf_error("ikv_temme(temme_ik_series)", SF_ERROR_NO_RESULT, NULL);
404  }
405 
406  *K = sum;
407  *K1 = 2 * sum1 / x;
408 
409  return 0;
410 }
411 
412 /* Evaluate continued fraction fv = I_(v+1) / I_v, derived from
413  * Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 */
414 static int CF1_ik(double v, double x, double *fv)
415 {
416  double C, D, f, a, b, delta, tiny, tolerance;
417  unsigned long k;
418 
419 
420  /*
421  * |x| <= |v|, CF1_ik converges rapidly
422  * |x| > |v|, CF1_ik needs O(|x|) iterations to converge
423  */
424 
425  /*
426  * modified Lentz's method, see
427  * Lentz, Applied Optics, vol 15, 668 (1976)
428  */
429  tolerance = 2 * MACHEP;
430  tiny = 1 / sqrt(DBL_MAX);
431  C = f = tiny; /* b0 = 0, replace with tiny */
432  D = 0;
433  for (k = 1; k < MAXITER; k++) {
434  a = 1;
435  b = 2 * (v + k) / x;
436  C = b + a / C;
437  D = b + a * D;
438  if (C == 0) {
439  C = tiny;
440  }
441  if (D == 0) {
442  D = tiny;
443  }
444  D = 1 / D;
445  delta = C * D;
446  f *= delta;
447  if (fabs(delta - 1) <= tolerance) {
448  break;
449  }
450  }
451  if (k == MAXITER) {
452  sf_error("ikv_temme(CF1_ik)", SF_ERROR_NO_RESULT, NULL);
453  }
454 
455  *fv = f;
456 
457  return 0;
458 }
459 
460 /*
461  * Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
462  * z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see
463  * Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
464  */
465 static int CF2_ik(double v, double x, double *Kv, double *Kv1)
466 {
467 
468  double S, C, Q, D, f, a, b, q, delta, tolerance, current, prev;
469  unsigned long k;
470 
471  /*
472  * |x| >= |v|, CF2_ik converges rapidly
473  * |x| -> 0, CF2_ik fails to converge
474  */
475 
476  BOOST_ASSERT(fabs(x) > 1);
477 
478  /*
479  * Steed's algorithm, see Thompson and Barnett,
480  * Journal of Computational Physics, vol 64, 490 (1986)
481  */
482  tolerance = MACHEP;
483  a = v * v - 0.25f;
484  b = 2 * (x + 1); /* b1 */
485  D = 1 / b; /* D1 = 1 / b1 */
486  f = delta = D; /* f1 = delta1 = D1, coincidence */
487  prev = 0; /* q0 */
488  current = 1; /* q1 */
489  Q = C = -a; /* Q1 = C1 because q1 = 1 */
490  S = 1 + Q * delta; /* S1 */
491  for (k = 2; k < MAXITER; k++) { /* starting from 2 */
492  /* continued fraction f = z1 / z0 */
493  a -= 2 * (k - 1);
494  b += 2;
495  D = 1 / (b + a * D);
496  delta *= b * D - 1;
497  f += delta;
498 
499  /* series summation S = 1 + \sum_{n=1}^{\infty} C_n * z_n / z_0 */
500  q = (prev - (b - 2) * current) / a;
501  prev = current;
502  current = q; /* forward recurrence for q */
503  C *= -a / k;
504  Q += C * q;
505  S += Q * delta;
506 
507  /* S converges slower than f */
508  if (fabs(Q * delta) < fabs(S) * tolerance) {
509  break;
510  }
511  }
512  if (k == MAXITER) {
513  sf_error("ikv_temme(CF2_ik)", SF_ERROR_NO_RESULT, NULL);
514  }
515 
516  *Kv = sqrt(M_PI / (2 * x)) * exp(-x) / S;
517  *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
518 
519  return 0;
520 }
521 
522 /* Flags for what to compute */
523 enum {
524  need_i = 0x1,
525  need_k = 0x2
526 };
527 
528 /*
529  * Compute I(v, x) and K(v, x) simultaneously by Temme's method, see
530  * Temme, Journal of Computational Physics, vol 19, 324 (1975)
531  */
532 static void ikv_temme(double v, double x, double *Iv_p, double *Kv_p)
533 {
534  /* Kv1 = K_(v+1), fv = I_(v+1) / I_v */
535  /* Ku1 = K_(u+1), fu = I_(u+1) / I_u */
536  double u, Iv, Kv, Kv1, Ku, Ku1, fv;
537  double W, current, prev, next;
538  int reflect = 0;
539  unsigned n, k;
540  int kind;
541 
542  kind = 0;
543  if (Iv_p != NULL) {
544  kind |= need_i;
545  }
546  if (Kv_p != NULL) {
547  kind |= need_k;
548  }
549 
550  if (v < 0) {
551  reflect = 1;
552  v = -v; /* v is non-negative from here */
553  kind |= need_k;
554  }
555  n = round(v);
556  u = v - n; /* -1/2 <= u < 1/2 */
557 
558  if (x < 0) {
559  if (Iv_p != NULL)
560  *Iv_p = NAN;
561  if (Kv_p != NULL)
562  *Kv_p = NAN;
563  sf_error("ikv_temme", SF_ERROR_DOMAIN, NULL);
564  return;
565  }
566  if (x == 0) {
567  Iv = (v == 0) ? 1 : 0;
568  if (kind & need_k) {
569  sf_error("ikv_temme", SF_ERROR_OVERFLOW, NULL);
570  Kv = INFINITY;
571  }
572  else {
573  Kv = NAN; /* any value will do */
574  }
575 
576  if (reflect && (kind & need_i)) {
577  double z = (u + n % 2);
578 
579  Iv = sin((double)M_PI * z) == 0 ? Iv : INFINITY;
580  if (Iv == INFINITY || Iv == -INFINITY) {
581  sf_error("ikv_temme", SF_ERROR_OVERFLOW, NULL);
582  }
583  }
584 
585  if (Iv_p != NULL) {
586  *Iv_p = Iv;
587  }
588  if (Kv_p != NULL) {
589  *Kv_p = Kv;
590  }
591  return;
592  }
593  /* x is positive until reflection */
594  W = 1 / x; /* Wronskian */
595  if (x <= 2) { /* x in (0, 2] */
596  temme_ik_series(u, x, &Ku, &Ku1); /* Temme series */
597  }
598  else { /* x in (2, \infty) */
599  CF2_ik(u, x, &Ku, &Ku1); /* continued fraction CF2_ik */
600  }
601  prev = Ku;
602  current = Ku1;
603  for (k = 1; k <= n; k++) { /* forward recurrence for K */
604  next = 2 * (u + k) * current / x + prev;
605  prev = current;
606  current = next;
607  }
608  Kv = prev;
609  Kv1 = current;
610  if (kind & need_i) {
611  double lim = (4 * v * v + 10) / (8 * x);
612 
613  lim *= lim;
614  lim *= lim;
615  lim /= 24;
616  if ((lim < MACHEP * 10) && (x > 100)) {
617  /*
618  * x is huge compared to v, CF1 may be very slow
619  * to converge so use asymptotic expansion for large
620  * x case instead. Note that the asymptotic expansion
621  * isn't very accurate - so it's deliberately very hard
622  * to get here - probably we're going to overflow:
623  */
624  Iv = iv_asymptotic(v, x);
625  }
626  else {
627  CF1_ik(v, x, &fv); /* continued fraction CF1_ik */
628  Iv = W / (Kv * fv + Kv1); /* Wronskian relation */
629  }
630  }
631  else {
632  Iv = NAN; /* any value will do */
633  }
634 
635  if (reflect) {
636  double z = (u + n % 2);
637 
638  if (Iv_p != NULL) {
639  *Iv_p = Iv + (2 / M_PI) * sin(M_PI * z) * Kv; /* reflection formula */
640  }
641  if (Kv_p != NULL) {
642  *Kv_p = Kv;
643  }
644  }
645  else {
646  if (Iv_p != NULL) {
647  *Iv_p = Iv;
648  }
649  if (Kv_p != NULL) {
650  *Kv_p = Kv;
651  }
652  }
653  return;
654 }
need_i
@ need_i
Definition: scipy_iv.c:524
D
MatrixXcd D
Definition: EigenSolver_EigenSolver_MatrixType.cpp:14
MAXITER
#define MAXITER
Definition: igam.c:110
BOOST_ASSERT
#define BOOST_ASSERT(a)
Definition: scipy_iv.c:343
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
d
static const double d[K][N]
Definition: igam.h:11
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
mu
double mu
Definition: testBoundingConstraint.cpp:37
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
SF_ERROR_NO_RESULT
@ SF_ERROR_NO_RESULT
Definition: sf_error.h:15
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
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
isnan
#define isnan(X)
Definition: main.h:93
iv_asymptotic
static double iv_asymptotic(double v, double x)
Definition: scipy_iv.c:145
CF1_ik
static int CF1_ik(double v, double x, double *fv)
Definition: scipy_iv.c:414
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
h
const double h
Definition: testSimpleHelicopter.cpp:19
Q
Quaternion Q
Definition: testQuaternion.cpp:27
SCIPY_EULER
#define SCIPY_EULER
Definition: mconf.h:128
sampling::sigma
static const double sigma
Definition: testGaussianBayesNet.cpp:170
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
cosh
const EIGEN_DEVICE_FUNC CoshReturnType cosh() const
Definition: ArrayCwiseUnaryOps.h:353
n
int n
Definition: BiCGSTAB_simple.cpp:1
gtsam.examples.PlanarManipulatorExample.delta
def delta(g0, g1)
Definition: PlanarManipulatorExample.py:45
temme_ik_series
static int temme_ik_series(double v, double x, double *K, double *K1)
Definition: scipy_iv.c:352
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
MACHEP
double MACHEP
Definition: const.c:54
round
double round(double x)
Definition: round.c:38
asymptotic_ufactors
static const double asymptotic_ufactors[N_UFACTORS][N_UFACTOR_TERMS]
Definition: scipy_iv.c:198
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
K1
static Cal3_S2::shared_ptr K1(new Cal3_S2(fov, w, h))
gamma
#define gamma
Definition: mconf.h:85
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
ikv_asymptotic_uniform
static void ikv_asymptotic_uniform(double v, double x, double *Iv, double *Kv)
Definition: scipy_iv.c:249
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
C
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:50
K
#define K
Definition: igam.h:8
gtsam::symbol_shorthand::W
Key W(std::uint64_t j)
Definition: inference/Symbol.h:170
p
float * p
Definition: Tutorial_Map_using.cpp:9
N_UFACTORS
#define N_UFACTORS
Definition: scipy_iv.c:196
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
need_k
@ need_k
Definition: scipy_iv.c:525
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
ikv_temme
static void ikv_temme(double v, double x, double *Iv, double *Kv)
Definition: scipy_iv.c:532
CF2_ik
static int CF2_ik(double v, double x, double *Kv, double *Kv1)
Definition: scipy_iv.c:465
iv
double iv(double v, double x)
Definition: scipy_iv.c:80
NULL
#define NULL
Definition: ccolamd.c:609
N_UFACTOR_TERMS
#define N_UFACTOR_TERMS
Definition: scipy_iv.c:197
M_PI
#define M_PI
Definition: mconf.h:117
sinh
const EIGEN_DEVICE_FUNC SinhReturnType sinh() const
Definition: ArrayCwiseUnaryOps.h:339
align_3::t
Point2 t(10, 10)
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
S
DiscreteKey S(1, 2)
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
gamma2
Vector gamma2
Definition: testSimpleHelicopter.cpp:31
SF_ERROR_LOSS
@ SF_ERROR_LOSS
Definition: sf_error.h:14


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:02