1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
/* ellpef.c
*
* Complete elliptic integral of the second kind
*
*
*
* SYNOPSIS:
*
* float m1, y, ellpef();
*
* y = ellpef( m1 );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
* pi/2
* -
* | | 2
* E(m) = | sqrt( 1 - m sin t ) dt
* | |
* -
* 0
*
* Where m = 1 - m1, using the approximation
*
* P(x) - x log x Q(x).
*
* Though there are no singularities, the argument m1 is used
* rather than m for compatibility with ellpk().
*
* E(1) = 1; E(0) = pi/2.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 1 30000 1.1e-7 3.9e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* ellpef domain x<0, x>1 0.0
*
*/
/* ellpe.c */
/* Elliptic integral of second kind */
/*
Cephes Math Library, Release 2.1: February, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <math.h>
static float P[] = {
1.53552577301013293365E-4,
2.50888492163602060990E-3,
8.68786816565889628429E-3,
1.07350949056076193403E-2,
7.77395492516787092951E-3,
7.58395289413514708519E-3,
1.15688436810574127319E-2,
2.18317996015557253103E-2,
5.68051945617860553470E-2,
4.43147180560990850618E-1,
1.00000000000000000299E0
};
static float Q[] = {
3.27954898576485872656E-5,
1.00962792679356715133E-3,
6.50609489976927491433E-3,
1.68862163993311317300E-2,
2.61769742454493659583E-2,
3.34833904888224918614E-2,
4.27180926518931511717E-2,
5.85936634471101055642E-2,
9.37499997197644278445E-2,
2.49999999999888314361E-1
};
float polevlf(float, float *, int), logf(float);
float ellpef( float xx)
{
float x;
x = xx;
if( (x <= 0.0) || (x > 1.0) )
{
if( x == 0.0 )
return( 1.0 );
mtherr( "ellpef", DOMAIN );
return( 0.0 );
}
return( polevlf(x,P,10) - logf(x) * (x * polevlf(x,Q,9)) );
}
|
