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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
/*
* Copyright © 2025 Keith Packard and Bart Massey.
* All Rights Reserved. See the file COPYING in this directory
* for licensing information.
*/
extend namespace Math {
public real acosh(x)
/*
* Returns inverse hyperbolic cosine of 'x'.
*/
{
if (x < 1)
raise invalid_argument("acosh arg < 1", 0, x);
return log(x + sqrt(x*x-1));
}
public real asinh(x)
/*
* Returns inverse hyperbolic sine of 'x'.
*/
{
if (x == 0) return 0;
real sign = 1;
if (x < 0) {
sign = -1;
x = -x;
}
return sign * log(x + sqrt(x*x+1));
}
public real atanh(x)
/*
* Returns inverse hyperbolic tangent of 'x'.
*/
{
if (abs(x) > 1)
raise invalid_argument("atanh arg > 1", 0, x);
return 0.5 * log((1 + x) / (1 - x));
}
public real cosh(x)
/*
* Returns hyperbolic cosine of 'x'.
*/
{
return (exp(x) + exp(-x)) / 2;
}
public real sinh(x)
/*
* Returns hyperbolic sine of 'x'.
*/
{
return (exp(x) - exp(-x)) / 2;
}
public real tanh(x)
/*
* Returns hyperbolic tangent of 'x'.
*/
{
return sinh(x) / cosh(x);
}
real _erf(real x, real off)
{
if (x == off)
return 0;
x = imprecise(x);
int bits = precision(x);
int obits = bits + 512;
real factor = 2 / sqrt(pi_value(obits));
x = imprecise(x, obits);
off = imprecise(off, obits) / factor;
real val = x - off;
for (int n = 1; ; n++) {
int f = 2 * n + 1;
real a = imprecise(((-1)**n * x**f) / (n! * f), obits);
val += a;
if (exponent(val) - exponent(a) > obits)
break;
}
return imprecise(val * factor, bits);
}
public real erf(real x)
/*
* Returns Gauss error function of 'x'.
*/
{
return _erf(x, 0);
}
public real erfc(real x)
/*
* Returns Gauss complementary error function of 'x'.
*/
{
return -_erf(x, 1);
}
public real jn(int n, real x)
/*
* Returns Bessel function of the first kind, order 'n' of 'x'
*/
{
if (n < 0) {
real v = jn(-n, x);
if (n % 2 == 1)
v = -v;
return v;
}
x = imprecise(x);
int bits = precision(x);
int obits = bits + 64;
x = imprecise(x, obits);
real val = imprecise(0, obits);
int s = 1;
int mf = 1;
int mnf = n!;
real x_2 = x/2;
real x_2_pow = x_2 ** n;
real x_2_2 = x_2 ** 2;
for (int m = 0; ;) {
real a = (s / (mf * mnf)) * x_2_pow;
val += a;
if (exponent(val) - exponent(a) > obits)
break;
s = s > 0 ? -1 : 1;
m++;
mf *= m;
mnf *= (m+n);
x_2_pow *= x_2_2;
}
return imprecise(val, bits);
}
public real j0(real x)
/*
* Returns Bessel function of the first kind, order '0' of 'x'
*/
{
return jn(0, x);
}
public real j1(real x)
/*
* Returns Bessel function of the first kind, order '1' of 'x'
*/
{
return jn(1, x);
}
}
|
