sirtimbly/IronCalc
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

UPDATE: Dump of initial files

Browse Source
This commit is contained in:
Nicolás Hatcher
2023-11-18 21:26:18 +01:00
commit c5b8efd83d
279 changed files with 42654 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
base/src/functions/engineering/transcendental/README.md Normal file
View File

@@ -0,0 +1,29 @@
# Creating tests from transcendental functions
Excel supports a number of transcendental functions like the error functions, gamma nad beta functions.
In this folder we have tests for the Bessel functions.
Some other platform's implementations of those functions are remarkably poor (including Excel), sometimes failing on the third decimal digit. We strive to do better.
To properly test you need to compute some known values with established arbitrary precision arithmetic.
I use for this purpose Arb[1], created by the unrivalled Fredrik Johansson[2]. You might find some python bindings, but I use Julia's Nemo[3]:
```julia
julia> using Nemo
julia> CC = AcbField(200)
julia> besseli(CC("17"), CC("5.6"))
```
If you are new to Julia, just install Julia and in the Julia terminal run:
```
julia> using Pkg; Pkg.add("Nemo")
```
You only need to do that once (like the R philosophy)
Will give you any Bessel function of any order (integer or not) of any value real or complex
[1]: https://arblib.org/
[2]: https://fredrikj.net/
[3]: https://nemocas.github.io/Nemo.jl/latest/

144
base/src/functions/engineering/transcendental/bessel_i.rs Normal file
View File

@@ -0,0 +1,144 @@
// This are somewhat lower precision than the BesselJ and BesselY
// needed for BesselK
pub(crate) fn bessel_i0(x: f64) -> f64 {
// Parameters of the polynomial approximation
let p1 = 1.0;
let p2 = 3.5156229;
let p3 = 3.0899424;
let p4 = 1.2067492;
let p5 = 0.2659732;
let p6 = 3.60768e-2;
let p7 = 4.5813e-3;
let q1 = 0.39894228;
let q2 = 1.328592e-2;
let q3 = 2.25319e-3;
let q4 = -1.57565e-3;
let q5 = 9.16281e-3;
let q6 = -2.057706e-2;
let q7 = 2.635537e-2;
let q8 = -1.647633e-2;
let q9 = 3.92377e-3;
let k1 = 3.75;
let ax = x.abs();
if x.is_infinite() {
return 0.0;
}
if ax < k1 {
// let xx = x / k1;
// let y = xx * xx;
// let mut result = 1.0;
// let max_iter = 50;
// let mut term = 1.0;
// for i in 1..max_iter {
// term = term * k1*k1*y /(2.0*i as f64).powi(2);
// result += term;
// };
// result
let xx = x / k1;
let y = xx * xx;
p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))
} else {
let y = k1 / ax;
((ax).exp() / (ax).sqrt())
* (q1
+ y * (q2
+ y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))))
}
}
// needed for BesselK
pub(crate) fn bessel_i1(x: f64) -> f64 {
let p1 = 0.5;
let p2 = 0.87890594;
let p3 = 0.51498869;
let p4 = 0.15084934;
let p5 = 2.658733e-2;
let p6 = 3.01532e-3;
let p7 = 3.2411e-4;
let q1 = 0.39894228;
let q2 = -3.988024e-2;
let q3 = -3.62018e-3;
let q4 = 1.63801e-3;
let q5 = -1.031555e-2;
let q6 = 2.282967e-2;
let q7 = -2.895312e-2;
let q8 = 1.787654e-2;
let q9 = -4.20059e-3;
let k1 = 3.75;
let ax = x.abs();
if ax < k1 {
let xx = x / k1;
let y = xx * xx;
x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))))
} else {
let y = k1 / ax;
let result = ((ax).exp() / (ax).sqrt())
* (q1
+ y * (q2
+ y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
if x < 0.0 {
return -result;
}
result
}
}
pub(crate) fn bessel_i(n: i32, x: f64) -> f64 {
let accuracy = 40;
let large_number = 1e10;
let small_number = 1e-10;
if n < 0 {
return f64::NAN;
}
if n == 0 {
return bessel_i0(x);
}
if x == 0.0 {
// I_n(0) = 0 for all n!= 0
return 0.0;
}
if n == 1 {
return bessel_i1(x);
}
if x.abs() > large_number {
return 0.0;
};
let tox = 2.0 / x.abs();
let mut bip = 0.0;
let mut bi = 1.0;
let mut result = 0.0;
let m = 2 * (((accuracy * n) as f64).sqrt().trunc() as i32 + n);
for j in (1..=m).rev() {
(bip, bi) = (bi, bip + (j as f64) * tox * bi);
// Prevent overflow
if bi.abs() > large_number {
result *= small_number;
bi *= small_number;
bip *= small_number;
}
if j == n {
result = bip;
}
}
result *= bessel_i0(x) / bi;
if (x < 0.0) && (n % 2 == 1) {
result = -result;
}
result
}

402
base/src/functions/engineering/transcendental/bessel_j0_y0.rs Normal file
View File

@@ -0,0 +1,402 @@
/* @(#)e_j0.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* j0(x), y0(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j0(x):
* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
* 2. Reduce x to |x| since j0(x)=j0(-x), and
* for x in (0,2)
* j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
* (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
* for x in (2,inf)
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* as follow:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (cos(x) + sin(x))
* sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse 1.)
*
* 3 Special cases
* j0(nan)= nan
* j0(0) = 1
* j0(inf) = 0
*
* Method -- y0(x):
* 1. For x<2.
* Since
* y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
* therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
* We use the following function to approximate y0,
* y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
* where
* U(z) = u00 + u01*z + ... + u06*z^6
* V(z) = 1 + v01*z + ... + v04*z^4
* with absolute approximation error bounded by 2**-72.
* Note: For tiny x, U/V = u0 and j0(x)~1, hence
* y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
* 2. For x>=2.
* y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* by the method menti1d above.
* 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
*/
use std::f64::consts::FRAC_2_PI;
use super::bessel_util::{high_word, split_words, FRAC_2_SQRT_PI, HUGE};
// R0/S0 on [0, 2.00]
const R02: f64 = 1.562_499_999_999_999_5e-2; // 0x3F8FFFFF, 0xFFFFFFFD
const R03: f64 = -1.899_792_942_388_547_2e-4; // 0xBF28E6A5, 0xB61AC6E9
const R04: f64 = 1.829_540_495_327_006_7e-6; // 0x3EBEB1D1, 0x0C503919
const R05: f64 = -4.618_326_885_321_032e-9; // 0xBE33D5E7, 0x73D63FCE
const S01: f64 = 1.561_910_294_648_900_1e-2; // 0x3F8FFCE8, 0x82C8C2A4
const S02: f64 = 1.169_267_846_633_374_5e-4; // 0x3F1EA6D2, 0xDD57DBF4
const S03: f64 = 5.135_465_502_073_181e-7; // 0x3EA13B54, 0xCE84D5A9
const S04: f64 = 1.166_140_033_337_9e-9; // 0x3E1408BC, 0xF4745D8F
/* The asymptotic expansions of pzero is
* 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
* For x >= 2, We approximate pzero by
* pzero(x) = 1 + (R/S)
* where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
* S = 1 + pS0*s^2 + ... + pS4*s^10
* and
* | pzero(x)-1-R/S | <= 2 ** ( -60.26)
*/
const P_R8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-7.031_249_999_999_004e-2, /* 0xBFB1FFFF, 0xFFFFFD32 */
-8.081_670_412_753_498, /* 0xC02029D0, 0xB44FA779 */
-2.570_631_056_797_048_5e2, /* 0xC0701102, 0x7B19E863 */
-2.485_216_410_094_288e3, /* 0xC0A36A6E, 0xCD4DCAFC */
-5.253_043_804_907_295e3, /* 0xC0B4850B, 0x36CC643D */
];
const P_S8: [f64; 5] = [
1.165_343_646_196_681_8e2, /* 0x405D2233, 0x07A96751 */
3.833_744_753_641_218_3e3, /* 0x40ADF37D, 0x50596938 */
4.059_785_726_484_725_5e4, /* 0x40E3D2BB, 0x6EB6B05F */
1.167_529_725_643_759_2e5, /* 0x40FC810F, 0x8F9FA9BD */
4.762_772_841_467_309_6e4, /* 0x40E74177, 0x4F2C49DC */
];
const P_R5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.141_254_646_918_945e-11, /* 0xBDA918B1, 0x47E495CC */
-7.031_249_408_735_993e-2, /* 0xBFB1FFFF, 0xE69AFBC6 */
-4.159_610_644_705_878, /* 0xC010A370, 0xF90C6BBF */
-6.767_476_522_651_673e1, /* 0xC050EB2F, 0x5A7D1783 */
-3.312_312_996_491_729_7e2, /* 0xC074B3B3, 0x6742CC63 */
-3.464_333_883_656_049e2, /* 0xC075A6EF, 0x28A38BD7 */
];
const P_S5: [f64; 5] = [
6.075_393_826_923_003_4e1, /* 0x404E6081, 0x0C98C5DE */
1.051_252_305_957_045_8e3, /* 0x40906D02, 0x5C7E2864 */
5.978_970_943_338_558e3, /* 0x40B75AF8, 0x8FBE1D60 */
9.625_445_143_577_745e3, /* 0x40C2CCB8, 0xFA76FA38 */
2.406_058_159_229_391e3, /* 0x40A2CC1D, 0xC70BE864 */
];
const P_R3: [f64; 6] = [
/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.547_046_017_719_519e-9, /* 0xBE25E103, 0x6FE1AA86 */
-7.031_196_163_814_817e-2, /* 0xBFB1FFF6, 0xF7C0E24B */
-2.409_032_215_495_296, /* 0xC00345B2, 0xAEA48074 */
-2.196_597_747_348_831e1, /* 0xC035F74A, 0x4CB94E14 */
-5.807_917_047_017_376e1, /* 0xC04D0A22, 0x420A1A45 */
-3.144_794_705_948_885e1, /* 0xC03F72AC, 0xA892D80F */
];
const P_S3: [f64; 5] = [
3.585_603_380_552_097e1, /* 0x4041ED92, 0x84077DD3 */
3.615_139_830_503_038_6e2, /* 0x40769839, 0x464A7C0E */
1.193_607_837_921_115_3e3, /* 0x4092A66E, 0x6D1061D6 */
1.127_996_798_569_074_1e3, /* 0x40919FFC, 0xB8C39B7E */
1.735_809_308_133_357_5e2, /* 0x4065B296, 0xFC379081 */
];
const P_R2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.875_343_330_325_264e-8, /* 0xBE77D316, 0xE927026D */
-7.030_309_954_836_247e-2, /* 0xBFB1FF62, 0x495E1E42 */
-1.450_738_467_809_529_9, /* 0xBFF73639, 0x8A24A843 */
-7.635_696_138_235_278, /* 0xC01E8AF3, 0xEDAFA7F3 */
-1.119_316_688_603_567_5e1, /* 0xC02662E6, 0xC5246303 */
-3.233_645_793_513_353_4, /* 0xC009DE81, 0xAF8FE70F */
];
const P_S2: [f64; 5] = [
2.222_029_975_320_888e1, /* 0x40363865, 0x908B5959 */
1.362_067_942_182_152e2, /* 0x4061069E, 0x0EE8878F */
2.704_702_786_580_835e2, /* 0x4070E786, 0x42EA079B */
1.538_753_942_083_203_3e2, /* 0x40633C03, 0x3AB6FAFF */
1.465_761_769_482_562e1, /* 0x402D50B3, 0x44391809 */
];
// Note: This function is only called for ix>=0x40000000 (see above)
fn pzero(x: f64) -> f64 {
let ix = high_word(x) & 0x7fffffff;
// ix>=0x40000000 && ix<=0x48000000);
let (p, q) = if ix >= 0x40200000 {
(P_R8, P_S8)
} else if ix >= 0x40122E8B {
(P_R5, P_S5)
} else if ix >= 0x4006DB6D {
(P_R3, P_S3)
} else {
(P_R2, P_S2)
};
let z = 1.0 / (x * x);
let r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
let s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
1.0 + r / s
}
/* For x >= 8, the asymptotic expansions of qzero is
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
* We approximate pzero by
* qzero(x) = s*(-1.25 + (R/S))
* where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
* S = 1 + qS0*s^2 + ... + qS5*s^12
* and
* | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
*/
const Q_R8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
7.324_218_749_999_35e-2, /* 0x3FB2BFFF, 0xFFFFFE2C */
1.176_820_646_822_527e1, /* 0x40278952, 0x5BB334D6 */
5.576_733_802_564_019e2, /* 0x40816D63, 0x15301825 */
8.859_197_207_564_686e3, /* 0x40C14D99, 0x3E18F46D */
3.701_462_677_768_878e4, /* 0x40E212D4, 0x0E901566 */
];
const Q_S8: [f64; 6] = [
1.637_760_268_956_898_2e2, /* 0x406478D5, 0x365B39BC */
8.098_344_946_564_498e3, /* 0x40BFA258, 0x4E6B0563 */
1.425_382_914_191_204_8e5, /* 0x41016652, 0x54D38C3F */
8.033_092_571_195_144e5, /* 0x412883DA, 0x83A52B43 */
8.405_015_798_190_605e5, /* 0x4129A66B, 0x28DE0B3D */
-3.438_992_935_378_666e5, /* 0xC114FD6D, 0x2C9530C5 */
];
const Q_R5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
1.840_859_635_945_155_3e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
7.324_217_666_126_848e-2, /* 0x3FB2BFFF, 0xD172B04C */
5.835_635_089_620_569_5, /* 0x401757B0, 0xB9953DD3 */
1.351_115_772_864_498_3e2, /* 0x4060E392, 0x0A8788E9 */
1.027_243_765_961_641e3, /* 0x40900CF9, 0x9DC8C481 */
1.989_977_858_646_053_8e3, /* 0x409F17E9, 0x53C6E3A6 */
];
const Q_S5: [f64; 6] = [
8.277_661_022_365_378e1, /* 0x4054B1B3, 0xFB5E1543 */
2.077_814_164_213_93e3, /* 0x40A03BA0, 0xDA21C0CE */
1.884_728_877_857_181e4, /* 0x40D267D2, 0x7B591E6D */
5.675_111_228_949_473e4, /* 0x40EBB5E3, 0x97E02372 */
3.597_675_384_251_145e4, /* 0x40E19118, 0x1F7A54A0 */
-5.354_342_756_019_448e3, /* 0xC0B4EA57, 0xBEDBC609 */
];
const Q_R3: [f64; 6] = [
/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.377_410_140_897_386e-9, /* 0x3E32CD03, 0x6ADECB82 */
7.324_111_800_429_114e-2, /* 0x3FB2BFEE, 0x0E8D0842 */
3.344_231_375_161_707, /* 0x400AC0FC, 0x61149CF5 */
4.262_184_407_454_126_5e1, /* 0x40454F98, 0x962DAEDD */
1.708_080_913_405_656e2, /* 0x406559DB, 0xE25EFD1F */
1.667_339_486_966_511_7e2, /* 0x4064D77C, 0x81FA21E0 */
];
const Q_S3: [f64; 6] = [
4.875_887_297_245_872e1, /* 0x40486122, 0xBFE343A6 */
7.096_892_210_566_06e2, /* 0x40862D83, 0x86544EB3 */
3.704_148_226_201_113_6e3, /* 0x40ACF04B, 0xE44DFC63 */
6.460_425_167_525_689e3, /* 0x40B93C6C, 0xD7C76A28 */
2.516_333_689_203_689_6e3, /* 0x40A3A8AA, 0xD94FB1C0 */
-1.492_474_518_361_564e2, /* 0xC062A7EB, 0x201CF40F */
];
const Q_R2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.504_444_448_869_832_7e-7, /* 0x3E84313B, 0x54F76BDB */
7.322_342_659_630_793e-2, /* 0x3FB2BEC5, 0x3E883E34 */
1.998_191_740_938_16, /* 0x3FFFF897, 0xE727779C */
1.449_560_293_478_857_4e1, /* 0x402CFDBF, 0xAAF96FE5 */
3.166_623_175_047_815_4e1, /* 0x403FAA8E, 0x29FBDC4A */
1.625_270_757_109_292_7e1, /* 0x403040B1, 0x71814BB4 */
];
const Q_S2: [f64; 6] = [
3.036_558_483_552_192e1, /* 0x403E5D96, 0xF7C07AED */
2.693_481_186_080_498_4e2, /* 0x4070D591, 0xE4D14B40 */
8.447_837_575_953_201e2, /* 0x408A6645, 0x22B3BF22 */
8.829_358_451_124_886e2, /* 0x408B977C, 0x9C5CC214 */
2.126_663_885_117_988_3e2, /* 0x406A9553, 0x0E001365 */
-5.310_954_938_826_669_5, /* 0xC0153E6A, 0xF8B32931 */
];
fn qzero(x: f64) -> f64 {
let ix = high_word(x) & 0x7fffffff;
let (p, q) = if ix >= 0x40200000 {
(Q_R8, Q_S8)
} else if ix >= 0x40122E8B {
(Q_R5, Q_S5)
} else if ix >= 0x4006DB6D {
(Q_R3, Q_S3)
} else {
(Q_R2, Q_S2)
};
let z = 1.0 / (x * x);
let r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
let s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
(-0.125 + r / s) / x
}
const U00: f64 = -7.380_429_510_868_723e-2; /* 0xBFB2E4D6, 0x99CBD01F */
const U01: f64 = 1.766_664_525_091_811_2e-1; /* 0x3FC69D01, 0x9DE9E3FC */
const U02: f64 = -1.381_856_719_455_969e-2; /* 0xBF8C4CE8, 0xB16CFA97 */
const U03: f64 = 3.474_534_320_936_836_5e-4; /* 0x3F36C54D, 0x20B29B6B */
const U04: f64 = -3.814_070_537_243_641_6e-6; /* 0xBECFFEA7, 0x73D25CAD */
const U05: f64 = 1.955_901_370_350_229_2e-8; /* 0x3E550057, 0x3B4EABD4 */
const U06: f64 = -3.982_051_941_321_034e-11; /* 0xBDC5E43D, 0x693FB3C8 */
const V01: f64 = 1.273_048_348_341_237e-2; /* 0x3F8A1270, 0x91C9C71A */
const V02: f64 = 7.600_686_273_503_533e-5; /* 0x3F13ECBB, 0xF578C6C1 */
const V03: f64 = 2.591_508_518_404_578e-7; /* 0x3E91642D, 0x7FF202FD */
const V04: f64 = 4.411_103_113_326_754_7e-10; /* 0x3DFE5018, 0x3BD6D9EF */
pub(crate) fn y0(x: f64) -> f64 {
let (lx, hx) = split_words(x);
let ix = 0x7fffffff & hx;
// Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0
if ix >= 0x7ff00000 {
return 1.0 / (x + x * x);
}
if (ix | lx) == 0 {
return f64::NEG_INFINITY;
}
if hx < 0 {
return f64::NAN;
}
if ix >= 0x40000000 {
// |x| >= 2.0
// y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
// where x0 = x-pi/4
// Better formula:
// cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
// = 1/sqrt(2) * (sin(x) + cos(x))
// sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
// = 1/sqrt(2) * (sin(x) - cos(x))
// To avoid cancellation, use
// sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
// to compute the worse 1.
let s = x.sin();
let c = x.cos();
let mut ss = s - c;
let mut cc = s + c;
// j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
// y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
if ix < 0x7fe00000 {
// make sure x+x not overflow
let z = -(x + x).cos();
if (s * c) < 0.0 {
cc = z / ss;
} else {
ss = z / cc;
}
}
return if ix > 0x48000000 {
FRAC_2_SQRT_PI * ss / x.sqrt()
} else {
let u = pzero(x);
let v = qzero(x);
FRAC_2_SQRT_PI * (u * ss + v * cc) / x.sqrt()
};
}
if ix <= 0x3e400000 {
// x < 2^(-27)
return U00 + FRAC_2_PI * x.ln();
}
let z = x * x;
let u = U00 + z * (U01 + z * (U02 + z * (U03 + z * (U04 + z * (U05 + z * U06)))));
let v = 1.0 + z * (V01 + z * (V02 + z * (V03 + z * V04)));
u / v + FRAC_2_PI * (j0(x) * x.ln())
}
pub(crate) fn j0(x: f64) -> f64 {
let hx = high_word(x);
let ix = hx & 0x7fffffff;
if x.is_nan() {
return x;
} else if x.is_infinite() {
return 0.0;
}
// the function is even
let x = x.abs();
if ix >= 0x40000000 {
// |x| >= 2.0
// let (s, c) = x.sin_cos()
let s = x.sin();
let c = x.cos();
let mut ss = s - c;
let mut cc = s + c;
// makes sure that x+x does not overflow
if ix < 0x7fe00000 {
// |x| < f64::MAX / 2.0
let z = -(x + x).cos();
if s * c < 0.0 {
cc = z / ss;
} else {
ss = z / cc;
}
}
// j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
// y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
return if ix > 0x48000000 {
// x < 5.253807105661922e-287 (2^(-951))
FRAC_2_SQRT_PI * cc / (x.sqrt())
} else {
let u = pzero(x);
let v = qzero(x);
FRAC_2_SQRT_PI * (u * cc - v * ss) / x.sqrt()
};
};
if ix < 0x3f200000 {
// |x| < 2^(-13)
if HUGE + x > 1.0 {
// raise inexact if x != 0
if ix < 0x3e400000 {
return 1.0; // |x|<2^(-27)
} else {
return 1.0 - 0.25 * x * x;
}
}
}
let z = x * x;
let r = z * (R02 + z * (R03 + z * (R04 + z * R05)));
let s = 1.0 + z * (S01 + z * (S02 + z * (S03 + z * S04)));
if ix < 0x3FF00000 {
// |x| < 1.00
1.0 + z * (-0.25 + (r / s))
} else {
let u = 0.5 * x;
(1.0 + u) * (1.0 - u) + z * (r / s)
}
}

391
base/src/functions/engineering/transcendental/bessel_j1_y1.rs Normal file
View File

@@ -0,0 +1,391 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
* Method -- j1(x):
* 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
* 2. Reduce x to |x| since j1(x)=-j1(-x), and
* for x in (0,2)
* j1(x) = x/2 + x*z*R0/S0, where z = x*x;
* (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 )
* for x in (2,inf)
* j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
* where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
* as follow:
* cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (sin(x) + cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
*
* 3 Special cases
* j1(nan)= nan
* j1(0) = 0
* j1(inf) = 0
*
* Method -- y1(x):
* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
* 2. For x<2.
* Since
* y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
* therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
* We use the following function to approximate y1,
* y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
* where for x in [0,2] (abs err less than 2**-65.89)
* U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
* V(z) = 1 + v0[0]*z + ... + v0[4]*z^5
* Note: For tiny x, 1/x dominate y1 and hence
* y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
* 3. For x>=2.
* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
* where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
* by method mentioned above.
*/
use std::f64::consts::FRAC_2_PI;
use super::bessel_util::{high_word, split_words, FRAC_2_SQRT_PI, HUGE};
// R0/S0 on [0,2]
const R00: f64 = -6.25e-2; // 0xBFB00000, 0x00000000
const R01: f64 = 1.407_056_669_551_897e-3; // 0x3F570D9F, 0x98472C61
const R02: f64 = -1.599_556_310_840_356e-5; // 0xBEF0C5C6, 0xBA169668
const R03: f64 = 4.967_279_996_095_844_5e-8; // 0x3E6AAAFA, 0x46CA0BD9
const S01: f64 = 1.915_375_995_383_634_6e-2; // 0x3F939D0B, 0x12637E53
const S02: f64 = 1.859_467_855_886_309_2e-4; // 0x3F285F56, 0xB9CDF664
const S03: f64 = 1.177_184_640_426_236_8e-6; // 0x3EB3BFF8, 0x333F8498
const S04: f64 = 5.046_362_570_762_170_4e-9; // 0x3E35AC88, 0xC97DFF2C
const S05: f64 = 1.235_422_744_261_379_1e-11; // 0x3DAB2ACF, 0xCFB97ED8
pub(crate) fn j1(x: f64) -> f64 {
let hx = high_word(x);
let ix = hx & 0x7fffffff;
if ix >= 0x7ff00000 {
return 1.0 / x;
}
let y = x.abs();
if ix >= 0x40000000 {
/* |x| >= 2.0 */
let s = y.sin();
let c = y.cos();
let mut ss = -s - c;
let mut cc = s - c;
if ix < 0x7fe00000 {
/* make sure y+y not overflow */
let z = (y + y).cos();
if s * c > 0.0 {
cc = z / ss;
} else {
ss = z / cc;
}
}
// j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
// y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
let z = if ix > 0x48000000 {
FRAC_2_SQRT_PI * cc / y.sqrt()
} else {
let u = pone(y);
let v = qone(y);
FRAC_2_SQRT_PI * (u * cc - v * ss) / y.sqrt()
};
if hx < 0 {
return -z;
} else {
return z;
}
}
if ix < 0x3e400000 {
/* |x|<2**-27 */
if HUGE + x > 1.0 {
return 0.5 * x; /* inexact if x!=0 necessary */
}
}
let z = x * x;
let mut r = z * (R00 + z * (R01 + z * (R02 + z * R03)));
let s = 1.0 + z * (S01 + z * (S02 + z * (S03 + z * (S04 + z * S05))));
r *= x;
x * 0.5 + r / s
}
const U0: [f64; 5] = [
-1.960_570_906_462_389_4e-1, /* 0xBFC91866, 0x143CBC8A */
5.044_387_166_398_113e-2, /* 0x3FA9D3C7, 0x76292CD1 */
-1.912_568_958_757_635_5e-3, /* 0xBF5F55E5, 0x4844F50F */
2.352_526_005_616_105e-5, /* 0x3EF8AB03, 0x8FA6B88E */
-9.190_991_580_398_789e-8, /* 0xBE78AC00, 0x569105B8 */
];
const V0: [f64; 5] = [
1.991_673_182_366_499e-2, /* 0x3F94650D, 0x3F4DA9F0 */
2.025_525_810_251_351_7e-4, /* 0x3F2A8C89, 0x6C257764 */
1.356_088_010_975_162_3e-6, /* 0x3EB6C05A, 0x894E8CA6 */
6.227_414_523_646_215e-9, /* 0x3E3ABF1D, 0x5BA69A86 */
1.665_592_462_079_920_8e-11, /* 0x3DB25039, 0xDACA772A */
];
pub(crate) fn y1(x: f64) -> f64 {
let (lx, hx) = split_words(x);
let ix = 0x7fffffff & hx;
// if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0
if ix >= 0x7ff00000 {
return 1.0 / (x + x * x);
}
if (ix | lx) == 0 {
return f64::NEG_INFINITY;
}
if hx < 0 {
return f64::NAN;
}
if ix >= 0x40000000 {
// |x| >= 2.0
let s = x.sin();
let c = x.cos();
let mut ss = -s - c;
let mut cc = s - c;
if ix < 0x7fe00000 {
// make sure x+x not overflow
let z = (x + x).cos();
if s * c > 0.0 {
cc = z / ss;
} else {
ss = z / cc;
}
}
/* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
* where x0 = x-3pi/4
* Better formula:
* cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (cos(x) + sin(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
return if ix > 0x48000000 {
FRAC_2_SQRT_PI * ss / x.sqrt()
} else {
let u = pone(x);
let v = qone(x);
FRAC_2_SQRT_PI * (u * ss + v * cc) / x.sqrt()
};
}
if ix <= 0x3c900000 {
// x < 2^(-54)
return -FRAC_2_PI / x;
}
let z = x * x;
let u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * U0[4])));
let v = 1.0 + z * (V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * V0[4]))));
x * (u / v) + FRAC_2_PI * (j1(x) * x.ln() - 1.0 / x)
}
/* For x >= 8, the asymptotic expansions of pone is
* 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
* We approximate pone by
* pone(x) = 1 + (R/S)
* where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
* S = 1 + ps0*s^2 + ... + ps4*s^10
* and
* | pone(x)-1-R/S | <= 2 ** ( -60.06)
*/
const PR8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
1.171_874_999_999_886_5e-1, /* 0x3FBDFFFF, 0xFFFFFCCE */
1.323_948_065_930_735_8e1, /* 0x402A7A9D, 0x357F7FCE */
4.120_518_543_073_785_6e2, /* 0x4079C0D4, 0x652EA590 */
3.874_745_389_139_605_3e3, /* 0x40AE457D, 0xA3A532CC */
7.914_479_540_318_917e3, /* 0x40BEEA7A, 0xC32782DD */
];
const PS8: [f64; 5] = [
1.142_073_703_756_784_1e2, /* 0x405C8D45, 0x8E656CAC */
3.650_930_834_208_534_6e3, /* 0x40AC85DC, 0x964D274F */
3.695_620_602_690_334_6e4, /* 0x40E20B86, 0x97C5BB7F */
9.760_279_359_349_508e4, /* 0x40F7D42C, 0xB28F17BB */
3.080_427_206_278_888e4, /* 0x40DE1511, 0x697A0B2D */
];
const PR5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
1.319_905_195_562_435_2e-11, /* 0x3DAD0667, 0xDAE1CA7D */
1.171_874_931_906_141e-1, /* 0x3FBDFFFF, 0xE2C10043 */
6.802_751_278_684_329, /* 0x401B3604, 0x6E6315E3 */
1.083_081_829_901_891_1e2, /* 0x405B13B9, 0x452602ED */
5.176_361_395_331_998e2, /* 0x40802D16, 0xD052D649 */
5.287_152_013_633_375e2, /* 0x408085B8, 0xBB7E0CB7 */
];
const PS5: [f64; 5] = [
5.928_059_872_211_313e1, /* 0x404DA3EA, 0xA8AF633D */
9.914_014_187_336_144e2, /* 0x408EFB36, 0x1B066701 */
5.353_266_952_914_88e3, /* 0x40B4E944, 0x5706B6FB */
7.844_690_317_495_512e3, /* 0x40BEA4B0, 0xB8A5BB15 */
1.504_046_888_103_610_6e3, /* 0x40978030, 0x036F5E51 */
];
const PR3: [f64; 6] = [
3.025_039_161_373_736e-9, /* 0x3E29FC21, 0xA7AD9EDD */
1.171_868_655_672_535_9e-1, /* 0x3FBDFFF5, 0x5B21D17B */
3.932_977_500_333_156_4, /* 0x400F76BC, 0xE85EAD8A */
3.511_940_355_916_369e1, /* 0x40418F48, 0x9DA6D129 */
9.105_501_107_507_813e1, /* 0x4056C385, 0x4D2C1837 */
4.855_906_851_973_649e1, /* 0x4048478F, 0x8EA83EE5 */
];
const PS3: [f64; 5] = [
3.479_130_950_012_515e1, /* 0x40416549, 0xA134069C */
3.367_624_587_478_257_5e2, /* 0x40750C33, 0x07F1A75F */
1.046_871_399_757_751_3e3, /* 0x40905B7C, 0x5037D523 */
8.908_113_463_982_564e2, /* 0x408BD67D, 0xA32E31E9 */
1.037_879_324_396_392_8e2, /* 0x4059F26D, 0x7C2EED53 */
];
const PR2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.077_108_301_068_737_4e-7, /* 0x3E7CE9D4, 0xF65544F4 */
1.171_762_194_626_833_5e-1, /* 0x3FBDFF42, 0xBE760D83 */
2.368_514_966_676_088, /* 0x4002F2B7, 0xF98FAEC0 */
1.224_261_091_482_612_3e1, /* 0x40287C37, 0x7F71A964 */
1.769_397_112_716_877_3e1, /* 0x4031B1A8, 0x177F8EE2 */
5.073_523_125_888_185, /* 0x40144B49, 0xA574C1FE */
];
const PS2: [f64; 5] = [
2.143_648_593_638_214e1, /* 0x40356FBD, 0x8AD5ECDC */
1.252_902_271_684_027_5e2, /* 0x405F5293, 0x14F92CD5 */
2.322_764_690_571_628e2, /* 0x406D08D8, 0xD5A2DBD9 */
1.176_793_732_871_471e2, /* 0x405D6B7A, 0xDA1884A9 */
8.364_638_933_716_183, /* 0x4020BAB1, 0xF44E5192 */
];
/* Note: This function is only called for ix>=0x40000000 (see above) */
fn pone(x: f64) -> f64 {
let ix = high_word(x) & 0x7fffffff;
// ix>=0x40000000 && ix<=0x48000000)
let (p, q) = if ix >= 0x40200000 {
(PR8, PS8)
} else if ix >= 0x40122E8B {
(PR5, PS5)
} else if ix >= 0x4006DB6D {
(PR3, PS3)
} else {
(PR2, PS2)
};
let z = 1.0 / (x * x);
let r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
let s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
1.0 + r / s
}
/* For x >= 8, the asymptotic expansions of qone is
* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
* We approximate pone by
* qone(x) = s*(0.375 + (R/S))
* where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
* S = 1 + qs1*s^2 + ... + qs6*s^12
* and
* | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13)
*/
const QR8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-1.025_390_624_999_927_1e-1, /* 0xBFBA3FFF, 0xFFFFFDF3 */
-1.627_175_345_445_9e1, /* 0xC0304591, 0xA26779F7 */
-7.596_017_225_139_501e2, /* 0xC087BCD0, 0x53E4B576 */
-1.184_980_667_024_295_9e4, /* 0xC0C724E7, 0x40F87415 */
-4.843_851_242_857_503_5e4, /* 0xC0E7A6D0, 0x65D09C6A */
];
const QS8: [f64; 6] = [
1.613_953_697_007_229e2, /* 0x40642CA6, 0xDE5BCDE5 */
7.825_385_999_233_485e3, /* 0x40BE9162, 0xD0D88419 */
1.338_753_362_872_495_8e5, /* 0x4100579A, 0xB0B75E98 */
7.196_577_236_832_409e5, /* 0x4125F653, 0x72869C19 */
6.666_012_326_177_764e5, /* 0x412457D2, 0x7719AD5C */
-2.944_902_643_038_346_4e5, /* 0xC111F969, 0x0EA5AA18 */
];
const QR5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.089_799_311_417_641e-11, /* 0xBDB6FA43, 0x1AA1A098 */
-1.025_390_502_413_754_3e-1, /* 0xBFBA3FFF, 0xCB597FEF */
-8.056_448_281_239_36, /* 0xC0201CE6, 0xCA03AD4B */
-1.836_696_074_748_883_8e2, /* 0xC066F56D, 0x6CA7B9B0 */
-1.373_193_760_655_081_6e3, /* 0xC09574C6, 0x6931734F */
-2.612_444_404_532_156_6e3, /* 0xC0A468E3, 0x88FDA79D */
];
const QS5: [f64; 6] = [
8.127_655_013_843_358e1, /* 0x405451B2, 0xFF5A11B2 */
1.991_798_734_604_859_6e3, /* 0x409F1F31, 0xE77BF839 */
1.746_848_519_249_089e4, /* 0x40D10F1F, 0x0D64CE29 */
4.985_142_709_103_523e4, /* 0x40E8576D, 0xAABAD197 */
2.794_807_516_389_181_2e4, /* 0x40DB4B04, 0xCF7C364B */
-4.719_183_547_951_285e3, /* 0xC0B26F2E, 0xFCFFA004 */
];
const QR3: [f64; 6] = [
-5.078_312_264_617_666e-9, /* 0xBE35CFA9, 0xD38FC84F */
-1.025_378_298_208_370_9e-1, /* 0xBFBA3FEB, 0x51AEED54 */
-4.610_115_811_394_734, /* 0xC01270C2, 0x3302D9FF */
-5.784_722_165_627_836_4e1, /* 0xC04CEC71, 0xC25D16DA */
-2.282_445_407_376_317e2, /* 0xC06C87D3, 0x4718D55F */
-2.192_101_284_789_093_3e2, /* 0xC06B66B9, 0x5F5C1BF6 */
];
const QS3: [f64; 6] = [
4.766_515_503_237_295e1, /* 0x4047D523, 0xCCD367E4 */
6.738_651_126_766_997e2, /* 0x40850EEB, 0xC031EE3E */
3.380_152_866_795_263_4e3, /* 0x40AA684E, 0x448E7C9A */
5.547_729_097_207_228e3, /* 0x40B5ABBA, 0xA61D54A6 */
1.903_119_193_388_108e3, /* 0x409DBC7A, 0x0DD4DF4B */
-1.352_011_914_443_073_4e2, /* 0xC060E670, 0x290A311F */
];
const QR2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.783_817_275_109_588_7e-7, /* 0xBE87F126, 0x44C626D2 */
-1.025_170_426_079_855_5e-1, /* 0xBFBA3E8E, 0x9148B010 */
-2.752_205_682_781_874_6, /* 0xC0060484, 0x69BB4EDA */
-1.966_361_626_437_037_2e1, /* 0xC033A9E2, 0xC168907F */
-4.232_531_333_728_305e1, /* 0xC04529A3, 0xDE104AAA */
-2.137_192_117_037_040_6e1, /* 0xC0355F36, 0x39CF6E52 */
];
const QS2: [f64; 6] = [
2.953_336_290_605_238_5e1, /* 0x403D888A, 0x78AE64FF */
2.529_815_499_821_905_3e2, /* 0x406F9F68, 0xDB821CBA */
7.575_028_348_686_454e2, /* 0x4087AC05, 0xCE49A0F7 */
7.393_932_053_204_672e2, /* 0x40871B25, 0x48D4C029 */
1.559_490_033_366_661_2e2, /* 0x40637E5E, 0x3C3ED8D4 */
-4.959_498_988_226_282, /* 0xC013D686, 0xE71BE86B */
];
// Note: This function is only called for ix>=0x40000000 (see above)
fn qone(x: f64) -> f64 {
let ix = high_word(x) & 0x7fffffff;
// ix>=0x40000000 && ix<=0x48000000
let (p, q) = if ix >= 0x40200000 {
(QR8, QS8)
} else if ix >= 0x40122E8B {
(QR5, QS5)
} else if ix >= 0x4006DB6D {
(QR3, QS3)
} else {
(QR2, QS2)
};
let z = 1.0 / (x * x);
let r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
let s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
(0.375 + r / s) / x
}

329
base/src/functions/engineering/transcendental/bessel_jn_yn.rs Normal file
View File

@@ -0,0 +1,329 @@
// https://github.com/JuliaLang/openlibm/blob/master/src/e_jn.c
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
* floating point Bessel's function of the 1st and 2nd kind
* of order n
*
* Special cases:
* y0(0)=y1(0)=yn(n,0) = -inf with division by 0 signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
* Note 2. About jn(n,x), yn(n,x)
* For n=0, j0(x) is called,
* for n=1, j1(x) is called,
* for n<x, forward recursion us used starting
* from values of j0(x) and j1(x).
* for n>x, a continued fraction approximation to
* j(n,x)/j(n-1,x) is evaluated and then backward
* recursion is used starting from a supposed value
* for j(n,x). The resulting value of j(0,x) is
* compared with the actual value to correct the
* supposed value of j(n,x).
*
* yn(n,x) is similar in all respects, except
* that forward recursion is used for all
* values of n>1.
*
*/
use super::{
bessel_j0_y0::{j0, y0},
bessel_j1_y1::{j1, y1},
bessel_util::{split_words, FRAC_2_SQRT_PI},
};
// Special cases are:
//
// $ J_n(n, ±\Infinity) = 0$
// $ J_n(n, NaN} = NaN $
// $ J_n(n, 0) = 0 $
pub(crate) fn jn(n: i32, x: f64) -> f64 {
let (lx, mut hx) = split_words(x);
let ix = 0x7fffffff & hx;
// if J(n,NaN) is NaN
if x.is_nan() {
return x;
}
// if (ix | (/*(u_int32_t)*/(lx | -lx)) >> 31) > 0x7ff00000 {
// return x + x;
// }
let (n, x) = if n < 0 {
// hx ^= 0x80000000;
hx = -hx;
(-n, -x)
} else {
(n, x)
};
if n == 0 {
return j0(x);
}
if n == 1 {
return j1(x);
}
let sign = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
// let sign = if x < 0.0 { -1 } else { 1 };
let x = x.abs();
let b = if (ix | lx) == 0 || ix >= 0x7ff00000 {
// if x is 0 or inf
0.0
} else if n as f64 <= x {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if ix >= 0x52D00000 {
/* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=x.sin(), c=x.cos(),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
let temp = match n & 3 {
0 => x.cos() + x.sin(),
1 => -x.cos() + x.sin(),
2 => -x.cos() - x.sin(),
3 => x.cos() - x.sin(),
_ => {
// Impossible: FIXME!
// panic!("")
0.0
}
};
FRAC_2_SQRT_PI * temp / x.sqrt()
} else {
let mut a = j0(x);
let mut b = j1(x);
for i in 1..n {
let temp = b;
b = b * (((i + i) as f64) / x) - a; /* avoid underflow */
a = temp;
}
b
}
} else {
// x < 2^(-29)
if ix < 0x3e100000 {
// x is tiny, return the first Taylor expansion of J(n,x)
// J(n,x) = 1/n!*(x/2)^n - ...
if n > 33 {
// underflow
0.0
} else {
let temp = x * 0.5;
let mut b = temp;
let mut a = 1;
for i in 2..=n {
a *= i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
b / (a as f64)
}
} else {
/* use backward recurrence */
/* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
* 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
* -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction:
* 1
* = -----------------------
* 1
* w - -----------------
* 1
* w+h - ---------
* w+2h - ...
*
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple
*/
let w = ((n + n) as f64) / x;
let h = 2.0 / x;
let mut q0 = w;
let mut z = w + h;
let mut q1 = w * z - 1.0;
let mut k = 1;
while q1 < 1.0e9 {
k += 1;
z += h;
let tmp = z * q1 - q0;
q0 = q1;
q1 = tmp;
}
let m = n + n;
let mut t = 0.0;
for i in (m..2 * (n + k)).step_by(2).rev() {
t = 1.0 / ((i as f64) / x - t);
}
// for (t=0, i = 2*(n+k); i>=m; i -= 2) t = 1/(i/x-t);
let mut a = t;
let mut b = 1.0;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is
* likely underflow to 0
*/
let mut tmp = n as f64;
let v = 2.0 / x;
tmp = tmp * f64::ln(f64::abs(v * tmp));
if tmp < 7.097_827_128_933_84e2 {
// for(i=n-1, di=(i+i); i>0; i--){
let mut di = 2.0 * ((n - 1) as f64);
for _ in (1..=n - 1).rev() {
let temp = b;
b *= di;
b = b / x - a;
a = temp;
di -= 2.0;
}
} else {
// for(i=n-1, di=(i+i); i>0; i--) {
let mut di = 2.0 * ((n - 1) as f64);
for _ in (1..=n - 1).rev() {
let temp = b;
b *= di;
b = b / x - a;
a = temp;
di -= 2.0;
/* scale b to avoid spurious overflow */
if b > 1e100 {
a /= b;
t /= b;
b = 1.0;
}
}
}
let z = j0(x);
let w = j1(x);
if z.abs() >= w.abs() {
t * z / b
} else {
t * w / a
}
}
};
if sign == 1 {
-b
} else {
b
}
}
// Yn returns the order-n Bessel function of the second kind.
//
// Special cases are:
//
// Y(n, +Inf) = 0
// Y(n ≥ 0, 0) = -Inf
// Y(n < 0, 0) = +Inf if n is odd, -Inf if n is even
// Y(n, x < 0) = NaN
// Y(n, NaN) = NaN
pub(crate) fn yn(n: i32, x: f64) -> f64 {
let (lx, hx) = split_words(x);
let ix = 0x7fffffff & hx;
// if Y(n, NaN) is NaN
if x.is_nan() {
return x;
}
// if (ix | (/*(u_int32_t)*/(lx | -lx)) >> 31) > 0x7ff00000 {
// return x + x;
// }
if (ix | lx) == 0 {
return f64::NEG_INFINITY;
}
if hx < 0 {
return f64::NAN;
}
let (n, sign) = if n < 0 {
(-n, 1 - ((n & 1) << 1))
} else {
(n, 1)
};
if n == 0 {
return y0(x);
}
if n == 1 {
return (sign as f64) * y1(x);
}
if ix == 0x7ff00000 {
return 0.0;
}
let b = if ix >= 0x52D00000 {
// x > 2^302
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=x.sin(), c=x.cos(),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
let temp = match n & 3 {
0 => x.sin() - x.cos(),
1 => -x.sin() - x.cos(),
2 => -x.sin() + x.cos(),
3 => x.sin() + x.cos(),
_ => {
// unreachable
0.0
}
};
FRAC_2_SQRT_PI * temp / x.sqrt()
} else {
let mut a = y0(x);
let mut b = y1(x);
for i in 1..n {
if b.is_infinite() {
break;
}
// if high_word(b) != 0xfff00000 {
// break;
// }
(a, b) = (b, ((2.0 * i as f64) / x) * b - a);
}
b
};
if sign > 0 {
b
} else {
-b
}
}

90
base/src/functions/engineering/transcendental/bessel_k.rs Normal file
View File

@@ -0,0 +1,90 @@
// This are somewhat lower precision than the BesselJ and BesselY
use super::bessel_i::bessel_i0;
use super::bessel_i::bessel_i1;
fn bessel_k0(x: f64) -> f64 {
let p1 = -0.57721566;
let p2 = 0.42278420;
let p3 = 0.23069756;
let p4 = 3.488590e-2;
let p5 = 2.62698e-3;
let p6 = 1.0750e-4;
let p7 = 7.4e-6;
let q1 = 1.25331414;
let q2 = -7.832358e-2;
let q3 = 2.189568e-2;
let q4 = -1.062446e-2;
let q5 = 5.87872e-3;
let q6 = -2.51540e-3;
let q7 = 5.3208e-4;
if x <= 0.0 {
return 0.0;
}
if x <= 2.0 {
let y = x * x / 4.0;
(-(x / 2.0).ln() * bessel_i0(x))
+ (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))))
} else {
let y = 2.0 / x;
((-x).exp() / x.sqrt())
* (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))))
}
}
fn bessel_k1(x: f64) -> f64 {
let p1 = 1.0;
let p2 = 0.15443144;
let p3 = -0.67278579;
let p4 = -0.18156897;
let p5 = -1.919402e-2;
let p6 = -1.10404e-3;
let p7 = -4.686e-5;
let q1 = 1.25331414;
let q2 = 0.23498619;
let q3 = -3.655620e-2;
let q4 = 1.504268e-2;
let q5 = -7.80353e-3;
let q6 = 3.25614e-3;
let q7 = -6.8245e-4;
if x <= 0.0 {
return f64::NAN;
}
if x <= 2.0 {
let y = x * x / 4.0;
((x / 2.0).ln() * bessel_i1(x))
+ (1. / x) * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))))
} else {
let y = 2.0 / x;
((-x).exp() / x.sqrt())
* (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))))
}
}
pub(crate) fn bessel_k(n: i32, x: f64) -> f64 {
if x <= 0.0 || n < 0 {
return f64::NAN;
}
if n == 0 {
return bessel_k0(x);
}
if n == 1 {
return bessel_k1(x);
}
// Perform upward recurrence for all x
let tox = 2.0 / x;
let mut bkm = bessel_k0(x);
let mut bk = bessel_k1(x);
for j in 1..n {
(bkm, bk) = (bk, bkm + (j as f64) * tox * bk);
}
bk
}

19
base/src/functions/engineering/transcendental/bessel_util.rs Normal file
View File

@@ -0,0 +1,19 @@
pub(crate) const HUGE: f64 = 1e300;
pub(crate) const FRAC_2_SQRT_PI: f64 = 5.641_895_835_477_563e-1;
pub(crate) fn high_word(x: f64) -> i32 {
let [_, _, _, _, a1, a2, a3, a4] = x.to_ne_bytes();
// let binding = x.to_ne_bytes();
// let high = <&[u8; 4]>::try_from(&binding[4..8]).expect("");
i32::from_ne_bytes([a1, a2, a3, a4])
}
pub(crate) fn split_words(x: f64) -> (i32, i32) {
let [a1, a2, a3, a4, b1, b2, b3, b4] = x.to_ne_bytes();
// let binding = x.to_ne_bytes();
// let high = <&[u8; 4]>::try_from(&binding[4..8]).expect("");
(
i32::from_ne_bytes([a1, a2, a3, a4]),
i32::from_ne_bytes([b1, b2, b3, b4]),
)
}

14
base/src/functions/engineering/transcendental/create_test.jl Normal file
View File

@@ -0,0 +1,14 @@
# Example file creating testing cases for BesselI
using Nemo
CC = AcbField(100)
values = [1, 2, 3, -2, 5, 30, 2e-8]
for value in values
y_acb = besseli(CC(1), CC(value))
real64 = convert(Float64, real(y_acb))
im64 = convert(Float64, real(y_acb))
println("(", value, ", ", real64, "),")
end

53
base/src/functions/engineering/transcendental/erf.rs Normal file
View File

@@ -0,0 +1,53 @@
pub(crate) fn erf(x: f64) -> f64 {
let cof = vec![
-1.3026537197817094,
6.419_697_923_564_902e-1,
1.9476473204185836e-2,
-9.561_514_786_808_63e-3,
-9.46595344482036e-4,
3.66839497852761e-4,
4.2523324806907e-5,
-2.0278578112534e-5,
-1.624290004647e-6,
1.303655835580e-6,
1.5626441722e-8,
-8.5238095915e-8,
6.529054439e-9,
5.059343495e-9,
-9.91364156e-10,
-2.27365122e-10,
9.6467911e-11,
2.394038e-12,
-6.886027e-12,
8.94487e-13,
3.13092e-13,
-1.12708e-13,
3.81e-16,
7.106e-15,
-1.523e-15,
-9.4e-17,
1.21e-16,
-2.8e-17,
];
let mut d = 0.0;
let mut dd = 0.0;
let x_abs = x.abs();
let t = 2.0 / (2.0 + x_abs);
let ty = 4.0 * t - 2.0;
for j in (1..=cof.len() - 1).rev() {
let tmp = d;
d = ty * d - dd + cof[j];
dd = tmp;
}
let res = t * f64::exp(-x_abs * x_abs + 0.5 * (cof[0] + ty * d) - dd);
if x < 0.0 {
res - 1.0
} else {
1.0 - res
}
}

16
base/src/functions/engineering/transcendental/mod.rs Normal file
View File

@@ -0,0 +1,16 @@
mod bessel_i;
mod bessel_j0_y0;
mod bessel_j1_y1;
mod bessel_jn_yn;
mod bessel_k;
mod bessel_util;
mod erf;
#[cfg(test)]
mod test_bessel;
pub(crate) use bessel_i::bessel_i;
pub(crate) use bessel_jn_yn::jn as bessel_j;
pub(crate) use bessel_jn_yn::yn as bessel_y;
pub(crate) use bessel_k::bessel_k;
pub(crate) use erf::erf;

183
base/src/functions/engineering/transcendental/test_bessel.rs Normal file
View File

@@ -0,0 +1,183 @@
use crate::functions::engineering::transcendental::bessel_k;
use super::{
bessel_i::bessel_i,
bessel_j0_y0::{j0, y0},
bessel_j1_y1::j1,
bessel_jn_yn::{jn, yn},
};
const EPS: f64 = 1e-13;
const EPS_LOW: f64 = 1e-6;
// Known values computed with Arb via Nemo.jl in Julia
// You can also use Mathematica
/// But please do not use Excel or any other software without arbitrary precision
fn numbers_are_close(a: f64, b: f64) -> bool {
if a == b {
// avoid underflow if a = b = 0.0
return true;
}
(a - b).abs() / ((a * a + b * b).sqrt()) < EPS
}
fn numbers_are_somewhat_close(a: f64, b: f64) -> bool {
if a == b {
// avoid underflow if a = b = 0.0
return true;
}
(a - b).abs() / ((a * a + b * b).sqrt()) < EPS_LOW
}
#[test]
fn bessel_j0_known_values() {
let cases = [
(2.4, 0.002507683297243813),
(0.5, 0.9384698072408129),
(1.0, 0.7651976865579666),
(1.12345, 0.7084999488947348),
(27.0, 0.07274191800588709),
(33.0, 0.09727067223550946),
(2e-4, 0.9999999900000001),
(0.0, 1.0),
(1e10, 2.175591750246892e-6),
];
for (value, known) in cases {
let f = j0(value);
assert!(
numbers_are_close(f, known),
"Got: {f}, expected: {known} for j0({value})"
);
}
}
#[test]
fn bessel_y0_known_values() {
let cases = [
(2.4, 0.5104147486657438),
(0.5, -0.4445187335067065),
(1.0, 0.08825696421567692),
(1.12345, 0.1783162909790613),
(27.0, 0.1352149762078722),
(33.0, 0.0991348255208796),
(2e-4, -5.496017824512429),
(1e10, -7.676508175792937e-6),
(1e-300, -439.8351636227653),
];
for (value, known) in cases {
let f = y0(value);
assert!(
numbers_are_close(f, known),
"Got: {f}, expected: {known} for y0({value})"
);
}
assert!(y0(0.0).is_infinite());
}
#[test]
fn bessel_j1_known_values() {
// Values computed with Maxima, the computer algebra system
// TODO: Recompute
let cases = [
(2.4, 0.5201852681819311),
(0.5, 0.2422684576748738),
(1.0, 0.4400505857449335),
(1.17232, 0.4910665691824317),
(27.5, 0.1521418932046569),
(42.0, -0.04599388822188721),
(3e-5, 1.499999999831249E-5),
(350.0, -0.02040531295214455),
(0.0, 0.0),
(1e12, -7.913802683850441e-7),
];
for (value, known) in cases {
let f = j1(value);
assert!(
numbers_are_close(f, known),
"Got: {f}, expected: {known} for j1({value})"
);
}
}
#[test]
fn bessel_jn_known_values() {
// Values computed with Maxima, the computer algebra system
// TODO: Recompute
let cases = [
(3, 0.5, 0.002_563_729_994_587_244),
(4, 0.5, 0.000_160_736_476_364_287_6),
(-3, 0.5, -0.002_563_729_994_587_244),
(-4, 0.5, 0.000_160_736_476_364_287_6),
(3, 30.0, 0.129211228759725),
(-3, 30.0, -0.129211228759725),
(4, 30.0, -0.052609000321320355),
(20, 30.0, 0.0048310199934040645),
(7, 0.0, 0.0),
];
for (n, value, known) in cases {
let f = jn(n, value);
assert!(
numbers_are_close(f, known),
"Got: {f}, expected: {known} for jn({n}, {value})"
);
}
}
#[test]
fn bessel_yn_known_values() {
let cases = [
(3, 0.5, -42.059494304723883),
(4, 0.5, -499.272_560_819_512_3),
(-3, 0.5, 42.059494304723883),
(-4, 0.5, -499.272_560_819_512_3),
(3, 35.0, -0.13191405300596323),
(-12, 12.2, -0.310438011314211),
(7, 1e12, 1.016_712_505_197_956_3e-7),
(35, 3.0, -6.895_879_073_343_495e31),
];
for (n, value, known) in cases {
let f = yn(n, value);
assert!(
numbers_are_close(f, known),
"Got: {f}, expected: {known} for yn({n}, {value})"
);
}
}
#[test]
fn bessel_in_known_values() {
let cases = [
(1, 0.5, 0.2578943053908963),
(3, 0.5, 0.002645111968990286),
(7, 0.2, 1.986608521182497e-11),
(7, 0.0, 0.0),
(0, -0.5, 1.0634833707413236),
// worse case scenario
(0, 3.7499, 9.118167894541882),
(0, 3.7501, 9.119723897590003),
];
for (n, value, known) in cases {
let f = bessel_i(n, value);
assert!(
numbers_are_somewhat_close(f, known),
"Got: {f}, expected: {known} for in({n}, {value})"
);
}
}
#[test]
fn bessel_kn_known_values() {
let cases = [
(1, 0.5, 1.656441120003301),
(0, 0.5, 0.9244190712276659),
(3, 0.5, 62.05790952993026),
];
for (n, value, known) in cases {
let f = bessel_k(n, value);
assert!(
numbers_are_somewhat_close(f, known),
"Got: {f}, expected: {known} for kn({n}, {value})"
);
}
}