ellpk.c
Go to the documentation of this file.
1 /* ellpk.c
2  *
3  * Complete elliptic integral of the first kind
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double m1, y, ellpk();
10  *
11  * y = ellpk( m1 );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Approximates the integral
18  *
19  *
20  *
21  * pi/2
22  * -
23  * | |
24  * | dt
25  * K(m) = | ------------------
26  * | 2
27  * | | sqrt( 1 - m sin t )
28  * -
29  * 0
30  *
31  * where m = 1 - m1, using the approximation
32  *
33  * P(x) - log x Q(x).
34  *
35  * The argument m1 is used internally rather than m so that the logarithmic
36  * singularity at m = 1 will be shifted to the origin; this
37  * preserves maximum accuracy.
38  *
39  * K(0) = pi/2.
40  *
41  * ACCURACY:
42  *
43  * Relative error:
44  * arithmetic domain # trials peak rms
45  * IEEE 0,1 30000 2.5e-16 6.8e-17
46  *
47  * ERROR MESSAGES:
48  *
49  * message condition value returned
50  * ellpk domain x<0, x>1 0.0
51  *
52  */
53 
54 /* ellpk.c */
55 
56 
57 /*
58  * Cephes Math Library, Release 2.0: April, 1987
59  * Copyright 1984, 1987 by Stephen L. Moshier
60  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
61  */
62 
63 #include "mconf.h"
64 
65 static double P[] = {
66  1.37982864606273237150E-4,
67  2.28025724005875567385E-3,
68  7.97404013220415179367E-3,
69  9.85821379021226008714E-3,
70  6.87489687449949877925E-3,
71  6.18901033637687613229E-3,
72  8.79078273952743772254E-3,
73  1.49380448916805252718E-2,
74  3.08851465246711995998E-2,
75  9.65735902811690126535E-2,
76  1.38629436111989062502E0
77 };
78 
79 static double Q[] = {
80  2.94078955048598507511E-5,
81  9.14184723865917226571E-4,
82  5.94058303753167793257E-3,
83  1.54850516649762399335E-2,
84  2.39089602715924892727E-2,
85  3.01204715227604046988E-2,
86  3.73774314173823228969E-2,
87  4.88280347570998239232E-2,
88  7.03124996963957469739E-2,
89  1.24999999999870820058E-1,
90  4.99999999999999999821E-1
91 };
92 
93 static double C1 = 1.3862943611198906188E0; /* log(4) */
94 
95 extern double MACHEP;
96 
97 double ellpk(double x)
98 {
99 
100  if (x < 0.0) {
101  sf_error("ellpk", SF_ERROR_DOMAIN, NULL);
102  return (NAN);
103  }
104 
105  if (x > 1.0) {
106  if (cephes_isinf(x)) {
107  return 0.0;
108  }
109  return ellpk(1/x)/sqrt(x);
110  }
111 
112  if (x > MACHEP) {
113  return (polevl(x, P, 10) - log(x) * polevl(x, Q, 10));
114  }
115  else {
116  if (x == 0.0) {
117  sf_error("ellpk", SF_ERROR_SINGULAR, NULL);
118  return (INFINITY);
119  }
120  else {
121  return (C1 - 0.5 * log(x));
122  }
123  }
124 }
MACHEP
double MACHEP
Definition: const.c:54
P
static double P[]
Definition: ellpk.c:65
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
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
cephes_isinf
#define cephes_isinf(x)
Definition: mconf.h:100
ellpk
double ellpk(double x)
Definition: ellpk.c:97
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
C1
Definition: test_operator_overloading.cpp:97
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
mconf.h
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
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


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:15