yggdrasil/libm
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

tests, fixes, format

Browse Source
This commit is contained in:
Andrey Zgarbul 2019-05-07 07:16:24 +03:00
parent b034304632
commit 2c9fd146d3
42 changed files with 4371 additions and 4173 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
build.rs
View File

@ -41,6 +41,19 @@ mod musl_reference_tests {
"rem_pio2.rs",
"rem_pio2_large.rs",
"rem_pio2f.rs",
"remquo.rs",
"remquof.rs",
"lgamma.rs", // lgamma passed, lgamma_r has more than 1 result
"lgammaf.rs", // lgammaf passed, lgammaf_r has more than 1 result
"frexp.rs", // more than 1 result
"frexpf.rs", // more than 1 result
"sincos.rs", // more than 1 result
"sincosf.rs", // more than 1 result
"modf.rs", // more than 1 result
"modff.rs", // more than 1 result
"asinef.rs", // not exists
"jn.rs", // passed, but very slow
"jnf.rs", // passed, but very slow
];
struct Function {
@ -78,12 +91,9 @@ mod musl_reference_tests {
let contents = fs::read_to_string(file).unwrap();
let mut functions = contents.lines().filter(|f| f.starts_with("pub fn"));
let function_to_test = functions.next().unwrap();
if functions.next().is_some() {
panic!("more than one function in");
while let Some(function_to_test) = functions.next() {
math.push(parse(function_to_test));
}
math.push(parse(function_to_test));
}
// Generate a bunch of random inputs for each function. This will
@ -330,7 +340,7 @@ mod musl_reference_tests {
src.push_str(match function.ret {
Ty::F32 => "if _eqf(output, f32::from_bits(*expected as u32)).is_ok() { continue }",
Ty::F64 => "if _eq(output, f64::from_bits(*expected as u64)).is_ok() { continue }",
Ty::I32 => "if output as i64 == expected { continue }",
Ty::I32 => "if output as i64 == *expected { continue }",
Ty::Bool => unreachable!(),
});

0
ci/run-docker.sh Executable file → Normal file
View File

0
ci/run.sh Executable file → Normal file
View File

4
src/math/acosh.rs
View File

@ -11,11 +11,11 @@ pub fn acosh(x: f64) -> f64 {
if e < 0x3ff + 1 {
/* |x| < 2, up to 2ulp error in [1,1.125] */
return log1p(x-1.0+sqrt((x-1.0)*(x-1.0)+2.0*(x-1.0)));
return log1p(x - 1.0 + sqrt((x - 1.0) * (x - 1.0) + 2.0 * (x - 1.0)));
}
if e < 0x3ff + 26 {
/* |x| < 0x1p26 */
return log(2.0*x-1.0/(x+sqrt(x*x-1.0)));
return log(2.0 * x - 1.0 / (x + sqrt(x * x - 1.0)));
}
/* |x| >= 0x1p26 or nan */
return log(x) + LN2;

8
src/math/acoshf.rs
View File

@ -7,14 +7,14 @@ pub fn acoshf(x: f32) -> f32 {
let u = x.to_bits();
let a = u & 0x7fffffff;
if a < 0x3f800000+(1<<23) {
if a < 0x3f800000 + (1 << 23) {
/* |x| < 2, invalid if x < 1 or nan */
/* up to 2ulp error in [1,1.125] */
return log1pf(x-1.0 + sqrtf((x-1.0)*(x-1.0)+2.0*(x-1.0)));
return log1pf(x - 1.0 + sqrtf((x - 1.0) * (x - 1.0) + 2.0 * (x - 1.0)));
}
if a < 0x3f800000+(12<<23) {
if a < 0x3f800000 + (12 << 23) {
/* |x| < 0x1p12 */
return logf(2.0*x - 1.0/(x+sqrtf(x*x-1.0)));
return logf(2.0 * x - 1.0 / (x + sqrtf(x * x - 1.0)));
}
/* x >= 0x1p12 */
return logf(x) + LN2;

8
src/math/asinef.rs
View File

@ -28,9 +28,7 @@ const A: [f32; 2] = [ 0.0, 0.785398163 ];
const B: [f32; 2] = [ 1.570796326, 0.785398163 ];
const Z_ROOTEPS_F: f32 = 1.7263349182589107e-4;
pub fn asinef(x: f32, acosine: usize) -> f32
{
let flag: usize;
pub fn asinef(x: f32, acosine: bool) -> f32 {
let i: usize;
let mut branch: bool = false;
let g: f32;
@ -45,10 +43,10 @@ pub fn asinef(x: f32, acosine: usize) -> f32
}
y = fabsf(x);
flag = acosine;
let flag = acosine;
if y > 0.5 {
i = 1 - flag;
i = (!flag) as usize;
/* Check for range error. */
if y > 1.0 {

4
src/math/asinh.rs
View File

@ -17,10 +17,10 @@ pub fn asinh(mut x: f64) -> f64 {
x = log(x) + LN2;
} else if e >= 0x3ff + 1 {
/* |x| >= 2 */
x = log(2.0*x + 1.0/(sqrt(x*x+1.0)+x));
x = log(2.0 * x + 1.0 / (sqrt(x * x + 1.0) + x));
} else if e >= 0x3ff - 26 {
/* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */
x = log1p(x + x*x/(sqrt(x*x+1.0)+1.0));
x = log1p(x + x * x / (sqrt(x * x + 1.0) + 1.0));
} else {
/* |x| < 0x1p-26, raise inexact if x != 0 */
let x1p120 = f64::from_bits(0x4770000000000000);

12
src/math/asinhf.rs
View File

@ -1,4 +1,4 @@
use super::{logf, log1pf, sqrtf};
use super::{log1pf, logf, sqrtf};
const LN2: f32 = 0.693147180559945309417232121458176568;
@ -11,15 +11,15 @@ pub fn asinhf(mut x: f32) -> f32 {
/* |x| */
x = f32::from_bits(i);
if i >= 0x3f800000 + (12<<23) {
if i >= 0x3f800000 + (12 << 23) {
/* |x| >= 0x1p12 or inf or nan */
x = logf(x) + LN2;
} else if i >= 0x3f800000 + (1<<23) {
} else if i >= 0x3f800000 + (1 << 23) {
/* |x| >= 2 */
x = logf(2.0*x + 1.0/(sqrtf(x*x+1.0)+x));
} else if i >= 0x3f800000 - (12<<23) {
x = logf(2.0 * x + 1.0 / (sqrtf(x * x + 1.0) + x));
} else if i >= 0x3f800000 - (12 << 23) {
/* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
x = log1pf(x + x*x/(sqrtf(x*x+1.0)+1.0));
x = log1pf(x + x * x / (sqrtf(x * x + 1.0) + 1.0));
} else {
/* |x| < 0x1p-12, raise inexact if x!=0 */
let x1p120 = f32::from_bits(0x7b800000);

2
src/math/atan2.rs
View File

@ -53,7 +53,7 @@ pub fn atan2(y: f64, x: f64) -> f64 {
let lx = x.to_bits() as u32;
let mut iy = (y.to_bits() >> 32) as u32;
let ly = y.to_bits() as u32;
if ((ix - 0x3ff00000) | lx) == 0 {
if ((ix.wrapping_sub(0x3ff00000)) | lx) == 0 {
/* x = 1.0 */
return atan(y);
}

19
src/math/atanh.rs
View File

@ -1,33 +1,32 @@
use super::{log1p};
use super::log1p;
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
pub fn atanh(mut x: f64) -> f64 {
let mut u = x.to_bits();
pub fn atanh(x: f64) -> f64 {
let u = x.to_bits();
let e = ((u >> 52) as usize) & 0x7ff;
let sign = (u >> 63) != 0;
/* |x| */
u &= 0x7fffffff;
x = f64::from_bits(u);
let mut y = f64::from_bits(u & 0x7fff_ffff_ffff_ffff);
if e < 0x3ff - 1 {
if e < 0x3ff - 32 {
/* handle underflow */
if e == 0 {
force_eval!(x as f32);
force_eval!(y as f32);
}
} else {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1p(2.0*x + 2.0*x*x/(1.0-x));
y = 0.5 * log1p(2.0 * y + 2.0 * y * y / (1.0 - y));
}
} else {
/* avoid overflow */
x = 0.5*log1p(2.0*(x/(1.0-x)));
y = 0.5 * log1p(2.0 * (y / (1.0 - y)));
}
if sign {
-x
-y
} else {
x
y
}
}

14
src/math/atanhf.rs
View File

@ -1,4 +1,4 @@
use super::{log1pf};
use super::log1pf;
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
pub fn atanhf(mut x: f32) -> f32 {
@ -9,19 +9,19 @@ pub fn atanhf(mut x: f32) -> f32 {
u &= 0x7fffffff;
x = f32::from_bits(u);
if u < 0x3f800000 - (1<<23) {
if u < 0x3f800000 - (32<<23) {
if u < 0x3f800000 - (1 << 23) {
if u < 0x3f800000 - (32 << 23) {
/* handle underflow */
if u < (1<<23) {
force_eval!((x*x) as f32);
if u < (1 << 23) {
force_eval!((x * x) as f32);
}
} else {
/* |x| < 0.5, up to 1.7ulp error */
x = 0.5*log1pf(2.0*x + 2.0*x*x/(1.0-x));
x = 0.5 * log1pf(2.0 * x + 2.0 * x * x / (1.0 - x));
}
} else {
/* avoid overflow */
x = 0.5*log1pf(2.0*(x/(1.0-x)));
x = 0.5 * log1pf(2.0 * (x / (1.0 - x)));
}
if sign {

2
src/math/ceilf.rs
View File

@ -12,7 +12,7 @@ pub fn ceilf(x: f32) -> f32 {
}
}
let mut ui = x.to_bits();
let e = (((ui >> 23) & 0xff) - 0x7f) as i32;
let e = (((ui >> 23) & 0xff).wrapping_sub(0x7f)) as i32;
if e >= 23 {
return x;

2
src/math/copysign.rs
View File

@ -2,6 +2,6 @@ pub fn copysign(x: f64, y: f64) -> f64 {
let mut ux = x.to_bits();
let uy = y.to_bits();
ux &= (!0) >> 1;
ux |= uy & (1<<63);
ux |= uy & (1 << 63);
f64::from_bits(ux)
}

205
src/math/erf.rs
View File

@ -104,73 +104,73 @@ use super::{exp, fabs, get_high_word, with_set_low_word};
* erfc/erf(NaN) is NaN
*/
const ERX: f64 = 8.45062911510467529297e-01; /* 0x3FEB0AC1, 0x60000000 */
const ERX: f64 = 8.45062911510467529297e-01; /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
const EFX8: f64 = 1.02703333676410069053e+00; /* 0x3FF06EBA, 0x8214DB69 */
const PP0: f64 = 1.28379167095512558561e-01; /* 0x3FC06EBA, 0x8214DB68 */
const PP1: f64 = -3.25042107247001499370e-01; /* 0xBFD4CD7D, 0x691CB913 */
const PP2: f64 = -2.84817495755985104766e-02; /* 0xBF9D2A51, 0xDBD7194F */
const PP3: f64 = -5.77027029648944159157e-03; /* 0xBF77A291, 0x236668E4 */
const PP4: f64 = -2.37630166566501626084e-05; /* 0xBEF8EAD6, 0x120016AC */
const QQ1: f64 = 3.97917223959155352819e-01; /* 0x3FD97779, 0xCDDADC09 */
const QQ2: f64 = 6.50222499887672944485e-02; /* 0x3FB0A54C, 0x5536CEBA */
const QQ3: f64 = 5.08130628187576562776e-03; /* 0x3F74D022, 0xC4D36B0F */
const QQ4: f64 = 1.32494738004321644526e-04; /* 0x3F215DC9, 0x221C1A10 */
const QQ5: f64 = -3.96022827877536812320e-06; /* 0xBED09C43, 0x42A26120 */
const EFX8: f64 = 1.02703333676410069053e+00; /* 0x3FF06EBA, 0x8214DB69 */
const PP0: f64 = 1.28379167095512558561e-01; /* 0x3FC06EBA, 0x8214DB68 */
const PP1: f64 = -3.25042107247001499370e-01; /* 0xBFD4CD7D, 0x691CB913 */
const PP2: f64 = -2.84817495755985104766e-02; /* 0xBF9D2A51, 0xDBD7194F */
const PP3: f64 = -5.77027029648944159157e-03; /* 0xBF77A291, 0x236668E4 */
const PP4: f64 = -2.37630166566501626084e-05; /* 0xBEF8EAD6, 0x120016AC */
const QQ1: f64 = 3.97917223959155352819e-01; /* 0x3FD97779, 0xCDDADC09 */
const QQ2: f64 = 6.50222499887672944485e-02; /* 0x3FB0A54C, 0x5536CEBA */
const QQ3: f64 = 5.08130628187576562776e-03; /* 0x3F74D022, 0xC4D36B0F */
const QQ4: f64 = 1.32494738004321644526e-04; /* 0x3F215DC9, 0x221C1A10 */
const QQ5: f64 = -3.96022827877536812320e-06; /* 0xBED09C43, 0x42A26120 */
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
const PA0: f64 = -2.36211856075265944077e-03; /* 0xBF6359B8, 0xBEF77538 */
const PA1: f64 = 4.14856118683748331666e-01; /* 0x3FDA8D00, 0xAD92B34D */
const PA2: f64 = -3.72207876035701323847e-01; /* 0xBFD7D240, 0xFBB8C3F1 */
const PA3: f64 = 3.18346619901161753674e-01; /* 0x3FD45FCA, 0x805120E4 */
const PA4: f64 = -1.10894694282396677476e-01; /* 0xBFBC6398, 0x3D3E28EC */
const PA5: f64 = 3.54783043256182359371e-02; /* 0x3FA22A36, 0x599795EB */
const PA6: f64 = -2.16637559486879084300e-03; /* 0xBF61BF38, 0x0A96073F */
const QA1: f64 = 1.06420880400844228286e-01; /* 0x3FBB3E66, 0x18EEE323 */
const QA2: f64 = 5.40397917702171048937e-01; /* 0x3FE14AF0, 0x92EB6F33 */
const QA3: f64 = 7.18286544141962662868e-02; /* 0x3FB2635C, 0xD99FE9A7 */
const QA4: f64 = 1.26171219808761642112e-01; /* 0x3FC02660, 0xE763351F */
const QA5: f64 = 1.36370839120290507362e-02; /* 0x3F8BEDC2, 0x6B51DD1C */
const QA6: f64 = 1.19844998467991074170e-02; /* 0x3F888B54, 0x5735151D */
const PA0: f64 = -2.36211856075265944077e-03; /* 0xBF6359B8, 0xBEF77538 */
const PA1: f64 = 4.14856118683748331666e-01; /* 0x3FDA8D00, 0xAD92B34D */
const PA2: f64 = -3.72207876035701323847e-01; /* 0xBFD7D240, 0xFBB8C3F1 */
const PA3: f64 = 3.18346619901161753674e-01; /* 0x3FD45FCA, 0x805120E4 */
const PA4: f64 = -1.10894694282396677476e-01; /* 0xBFBC6398, 0x3D3E28EC */
const PA5: f64 = 3.54783043256182359371e-02; /* 0x3FA22A36, 0x599795EB */
const PA6: f64 = -2.16637559486879084300e-03; /* 0xBF61BF38, 0x0A96073F */
const QA1: f64 = 1.06420880400844228286e-01; /* 0x3FBB3E66, 0x18EEE323 */
const QA2: f64 = 5.40397917702171048937e-01; /* 0x3FE14AF0, 0x92EB6F33 */
const QA3: f64 = 7.18286544141962662868e-02; /* 0x3FB2635C, 0xD99FE9A7 */
const QA4: f64 = 1.26171219808761642112e-01; /* 0x3FC02660, 0xE763351F */
const QA5: f64 = 1.36370839120290507362e-02; /* 0x3F8BEDC2, 0x6B51DD1C */
const QA6: f64 = 1.19844998467991074170e-02; /* 0x3F888B54, 0x5735151D */
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
const RA0: f64 = -9.86494403484714822705e-03; /* 0xBF843412, 0x600D6435 */
const RA1: f64 = -6.93858572707181764372e-01; /* 0xBFE63416, 0xE4BA7360 */
const RA2: f64 = -1.05586262253232909814e+01; /* 0xC0251E04, 0x41B0E726 */
const RA3: f64 = -6.23753324503260060396e+01; /* 0xC04F300A, 0xE4CBA38D */
const RA4: f64 = -1.62396669462573470355e+02; /* 0xC0644CB1, 0x84282266 */
const RA5: f64 = -1.84605092906711035994e+02; /* 0xC067135C, 0xEBCCABB2 */
const RA6: f64 = -8.12874355063065934246e+01; /* 0xC0545265, 0x57E4D2F2 */
const RA7: f64 = -9.81432934416914548592e+00; /* 0xC023A0EF, 0xC69AC25C */
const SA1: f64 = 1.96512716674392571292e+01; /* 0x4033A6B9, 0xBD707687 */
const SA2: f64 = 1.37657754143519042600e+02; /* 0x4061350C, 0x526AE721 */
const SA3: f64 = 4.34565877475229228821e+02; /* 0x407B290D, 0xD58A1A71 */
const SA4: f64 = 6.45387271733267880336e+02; /* 0x40842B19, 0x21EC2868 */
const SA5: f64 = 4.29008140027567833386e+02; /* 0x407AD021, 0x57700314 */
const SA6: f64 = 1.08635005541779435134e+02; /* 0x405B28A3, 0xEE48AE2C */
const SA7: f64 = 6.57024977031928170135e+00; /* 0x401A47EF, 0x8E484A93 */
const SA8: f64 = -6.04244152148580987438e-02; /* 0xBFAEEFF2, 0xEE749A62 */
const RA0: f64 = -9.86494403484714822705e-03; /* 0xBF843412, 0x600D6435 */
const RA1: f64 = -6.93858572707181764372e-01; /* 0xBFE63416, 0xE4BA7360 */
const RA2: f64 = -1.05586262253232909814e+01; /* 0xC0251E04, 0x41B0E726 */
const RA3: f64 = -6.23753324503260060396e+01; /* 0xC04F300A, 0xE4CBA38D */
const RA4: f64 = -1.62396669462573470355e+02; /* 0xC0644CB1, 0x84282266 */
const RA5: f64 = -1.84605092906711035994e+02; /* 0xC067135C, 0xEBCCABB2 */
const RA6: f64 = -8.12874355063065934246e+01; /* 0xC0545265, 0x57E4D2F2 */
const RA7: f64 = -9.81432934416914548592e+00; /* 0xC023A0EF, 0xC69AC25C */
const SA1: f64 = 1.96512716674392571292e+01; /* 0x4033A6B9, 0xBD707687 */
const SA2: f64 = 1.37657754143519042600e+02; /* 0x4061350C, 0x526AE721 */
const SA3: f64 = 4.34565877475229228821e+02; /* 0x407B290D, 0xD58A1A71 */
const SA4: f64 = 6.45387271733267880336e+02; /* 0x40842B19, 0x21EC2868 */
const SA5: f64 = 4.29008140027567833386e+02; /* 0x407AD021, 0x57700314 */
const SA6: f64 = 1.08635005541779435134e+02; /* 0x405B28A3, 0xEE48AE2C */
const SA7: f64 = 6.57024977031928170135e+00; /* 0x401A47EF, 0x8E484A93 */
const SA8: f64 = -6.04244152148580987438e-02; /* 0xBFAEEFF2, 0xEE749A62 */
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
const RB0: f64 = -9.86494292470009928597e-03; /* 0xBF843412, 0x39E86F4A */
const RB1: f64 = -7.99283237680523006574e-01; /* 0xBFE993BA, 0x70C285DE */
const RB2: f64 = -1.77579549177547519889e+01; /* 0xC031C209, 0x555F995A */
const RB3: f64 = -1.60636384855821916062e+02; /* 0xC064145D, 0x43C5ED98 */
const RB4: f64 = -6.37566443368389627722e+02; /* 0xC083EC88, 0x1375F228 */
const RB5: f64 = -1.02509513161107724954e+03; /* 0xC0900461, 0x6A2E5992 */
const RB6: f64 = -4.83519191608651397019e+02; /* 0xC07E384E, 0x9BDC383F */
const SB1: f64 = 3.03380607434824582924e+01; /* 0x403E568B, 0x261D5190 */
const SB2: f64 = 3.25792512996573918826e+02; /* 0x40745CAE, 0x221B9F0A */
const SB3: f64 = 1.53672958608443695994e+03; /* 0x409802EB, 0x189D5118 */
const SB4: f64 = 3.19985821950859553908e+03; /* 0x40A8FFB7, 0x688C246A */
const SB5: f64 = 2.55305040643316442583e+03; /* 0x40A3F219, 0xCEDF3BE6 */
const SB6: f64 = 4.74528541206955367215e+02; /* 0x407DA874, 0xE79FE763 */
const SB7: f64 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
const RB0: f64 = -9.86494292470009928597e-03; /* 0xBF843412, 0x39E86F4A */
const RB1: f64 = -7.99283237680523006574e-01; /* 0xBFE993BA, 0x70C285DE */
const RB2: f64 = -1.77579549177547519889e+01; /* 0xC031C209, 0x555F995A */
const RB3: f64 = -1.60636384855821916062e+02; /* 0xC064145D, 0x43C5ED98 */
const RB4: f64 = -6.37566443368389627722e+02; /* 0xC083EC88, 0x1375F228 */
const RB5: f64 = -1.02509513161107724954e+03; /* 0xC0900461, 0x6A2E5992 */
const RB6: f64 = -4.83519191608651397019e+02; /* 0xC07E384E, 0x9BDC383F */
const SB1: f64 = 3.03380607434824582924e+01; /* 0x403E568B, 0x261D5190 */
const SB2: f64 = 3.25792512996573918826e+02; /* 0x40745CAE, 0x221B9F0A */
const SB3: f64 = 1.53672958608443695994e+03; /* 0x409802EB, 0x189D5118 */
const SB4: f64 = 3.19985821950859553908e+03; /* 0x40A8FFB7, 0x688C246A */
const SB5: f64 = 2.55305040643316442583e+03; /* 0x40A3F219, 0xCEDF3BE6 */
const SB6: f64 = 4.74528541206955367215e+02; /* 0x407DA874, 0xE79FE763 */
const SB7: f64 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
fn erfc1(x: f64) -> f64 {
let s: f64;
@ -178,10 +178,10 @@ fn erfc1(x: f64) -> f64 {
let q: f64;
s = fabs(x) - 1.0;
p = PA0+s*(PA1+s*(PA2+s*(PA3+s*(PA4+s*(PA5+s*PA6)))));
q = 1.0+s*(QA1+s*(QA2+s*(QA3+s*(QA4+s*(QA5+s*QA6)))));
p = PA0 + s * (PA1 + s * (PA2 + s * (PA3 + s * (PA4 + s * (PA5 + s * PA6)))));
q = 1.0 + s * (QA1 + s * (QA2 + s * (QA3 + s * (QA4 + s * (QA5 + s * QA6)))));
1.0 - ERX - p/q
1.0 - ERX - p / q
}
fn erfc2(ix: u32, mut x: f64) -> f64 {
@ -190,26 +190,28 @@ fn erfc2(ix: u32, mut x: f64) -> f64 {
let big_s: f64;
let z: f64;
if ix < 0x3ff40000 { /* |x| < 1.25 */
if ix < 0x3ff40000 {
/* |x| < 1.25 */
return erfc1(x);
}
x = fabs(x);
s = 1.0/(x*x);
if ix < 0x4006db6d { /* |x| < 1/.35 ~ 2.85714 */
r = RA0+s*(RA1+s*(RA2+s*(RA3+s*(RA4+s*(
RA5+s*(RA6+s*RA7))))));
big_s = 1.0+s*(SA1+s*(SA2+s*(SA3+s*(SA4+s*(
SA5+s*(SA6+s*(SA7+s*SA8)))))));
} else { /* |x| > 1/.35 */
r = RB0+s*(RB1+s*(RB2+s*(RB3+s*(RB4+s*(
RB5+s*RB6)))));
big_s = 1.0+s*(SB1+s*(SB2+s*(SB3+s*(SB4+s*(
SB5+s*(SB6+s*SB7))))));
s = 1.0 / (x * x);
if ix < 0x4006db6d {
/* |x| < 1/.35 ~ 2.85714 */
r = RA0 + s * (RA1 + s * (RA2 + s * (RA3 + s * (RA4 + s * (RA5 + s * (RA6 + s * RA7))))));
big_s = 1.0
+ s * (SA1
+ s * (SA2 + s * (SA3 + s * (SA4 + s * (SA5 + s * (SA6 + s * (SA7 + s * SA8)))))));
} else {
/* |x| > 1/.35 */
r = RB0 + s * (RB1 + s * (RB2 + s * (RB3 + s * (RB4 + s * (RB5 + s * RB6)))));
big_s =
1.0 + s * (SB1 + s * (SB2 + s * (SB3 + s * (SB4 + s * (SB5 + s * (SB6 + s * SB7))))));
}
z = with_set_low_word(x, 0);
exp(-z*z-0.5625)*exp((z-x)*(z+x)+r/big_s)/x
exp(-z * z - 0.5625) * exp((z - x) * (z + x) + r / big_s) / x
}
pub fn erf(x: f64) -> f64 {
@ -221,25 +223,28 @@ pub fn erf(x: f64) -> f64 {
let sign: usize;
ix = get_high_word(x);
sign = (ix>>31) as usize;
sign = (ix >> 31) as usize;
ix &= 0x7fffffff;
if ix >= 0x7ff00000 {
/* erf(nan)=nan, erf(+-inf)=+-1 */
return 1.0-2.0*(sign as f64) + 1.0/x;
return 1.0 - 2.0 * (sign as f64) + 1.0 / x;
}
if ix < 0x3feb0000 { /* |x| < 0.84375 */
if ix < 0x3e300000 { /* |x| < 2**-28 */
if ix < 0x3feb0000 {
/* |x| < 0.84375 */
if ix < 0x3e300000 {
/* |x| < 2**-28 */
/* avoid underflow */
return 0.125*(8.0*x + EFX8*x);
return 0.125 * (8.0 * x + EFX8 * x);
}
z = x*x;
r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
y = r/s;
return x + x*y;
z = x * x;
r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4)));
s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5))));
y = r / s;
return x + x * y;
}
if ix < 0x40180000 { /* 0.84375 <= |x| < 6 */
y = 1.0 - erfc2(ix,x);
if ix < 0x40180000 {
/* 0.84375 <= |x| < 6 */
y = 1.0 - erfc2(ix, x);
} else {
let x1p_1022 = f64::from_bits(0x0010000000000000);
y = 1.0 - x1p_1022;
@ -261,30 +266,34 @@ pub fn erfc(x: f64) -> f64 {
let sign: usize;
ix = get_high_word(x);
sign = (ix>>31) as usize;
sign = (ix >> 31) as usize;
ix &= 0x7fffffff;
if ix >= 0x7ff00000 {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
return 2.0*(sign as f64) + 1.0/x;
return 2.0 * (sign as f64) + 1.0 / x;
}
if ix < 0x3feb0000 { /* |x| < 0.84375 */
if ix < 0x3c700000 { /* |x| < 2**-56 */
if ix < 0x3feb0000 {
/* |x| < 0.84375 */
if ix < 0x3c700000 {
/* |x| < 2**-56 */
return 1.0 - x;
}
z = x*x;
r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
y = r/s;
if sign != 0 || ix < 0x3fd00000 { /* x < 1/4 */
return 1.0 - (x+x*y);
z = x * x;
r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4)));
s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5))));
y = r / s;
if sign != 0 || ix < 0x3fd00000 {
/* x < 1/4 */
return 1.0 - (x + x * y);
}
return 0.5 - (x - 0.5 + x*y);
return 0.5 - (x - 0.5 + x * y);
}
if ix < 0x403c0000 { /* 0.84375 <= |x| < 28 */
if ix < 0x403c0000 {
/* 0.84375 <= |x| < 28 */
if sign != 0 {
return 2.0 - erfc2(ix,x);
return 2.0 - erfc2(ix, x);
} else {
return erfc2(ix,x);
return erfc2(ix, x);
}
}
@ -292,6 +301,6 @@ pub fn erfc(x: f64) -> f64 {
if sign != 0 {
2.0 - x1p_1022
} else {
x1p_1022*x1p_1022
x1p_1022 * x1p_1022
}
}

206
src/math/erff.rs
View File

@ -15,73 +15,73 @@
use super::{expf, fabsf};
const ERX: f32 = 8.4506291151e-01; /* 0x3f58560b */
const ERX: f32 = 8.4506291151e-01; /* 0x3f58560b */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
const EFX8: f32 = 1.0270333290e+00; /* 0x3f8375d4 */
const PP0: f32 = 1.2837916613e-01; /* 0x3e0375d4 */
const PP1: f32 = -3.2504209876e-01; /* 0xbea66beb */
const PP2: f32 = -2.8481749818e-02; /* 0xbce9528f */
const PP3: f32 = -5.7702702470e-03; /* 0xbbbd1489 */
const PP4: f32 = -2.3763017452e-05; /* 0xb7c756b1 */
const QQ1: f32 = 3.9791721106e-01; /* 0x3ecbbbce */
const QQ2: f32 = 6.5022252500e-02; /* 0x3d852a63 */
const QQ3: f32 = 5.0813062117e-03; /* 0x3ba68116 */
const QQ4: f32 = 1.3249473704e-04; /* 0x390aee49 */
const QQ5: f32 = -3.9602282413e-06; /* 0xb684e21a */
const EFX8: f32 = 1.0270333290e+00; /* 0x3f8375d4 */
const PP0: f32 = 1.2837916613e-01; /* 0x3e0375d4 */
const PP1: f32 = -3.2504209876e-01; /* 0xbea66beb */
const PP2: f32 = -2.8481749818e-02; /* 0xbce9528f */
const PP3: f32 = -5.7702702470e-03; /* 0xbbbd1489 */
const PP4: f32 = -2.3763017452e-05; /* 0xb7c756b1 */
const QQ1: f32 = 3.9791721106e-01; /* 0x3ecbbbce */
const QQ2: f32 = 6.5022252500e-02; /* 0x3d852a63 */
const QQ3: f32 = 5.0813062117e-03; /* 0x3ba68116 */
const QQ4: f32 = 1.3249473704e-04; /* 0x390aee49 */
const QQ5: f32 = -3.9602282413e-06; /* 0xb684e21a */
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
const PA0: f32 = -2.3621185683e-03; /* 0xbb1acdc6 */
const PA1: f32 = 4.1485610604e-01; /* 0x3ed46805 */
const PA2: f32 = -3.7220788002e-01; /* 0xbebe9208 */
const PA3: f32 = 3.1834661961e-01; /* 0x3ea2fe54 */
const PA4: f32 = -1.1089469492e-01; /* 0xbde31cc2 */
const PA5: f32 = 3.5478305072e-02; /* 0x3d1151b3 */
const PA6: f32 = -2.1663755178e-03; /* 0xbb0df9c0 */
const QA1: f32 = 1.0642088205e-01; /* 0x3dd9f331 */
const QA2: f32 = 5.4039794207e-01; /* 0x3f0a5785 */
const QA3: f32 = 7.1828655899e-02; /* 0x3d931ae7 */
const QA4: f32 = 1.2617121637e-01; /* 0x3e013307 */
const QA5: f32 = 1.3637083583e-02; /* 0x3c5f6e13 */
const QA6: f32 = 1.1984500103e-02; /* 0x3c445aa3 */
const PA0: f32 = -2.3621185683e-03; /* 0xbb1acdc6 */
const PA1: f32 = 4.1485610604e-01; /* 0x3ed46805 */
const PA2: f32 = -3.7220788002e-01; /* 0xbebe9208 */
const PA3: f32 = 3.1834661961e-01; /* 0x3ea2fe54 */
const PA4: f32 = -1.1089469492e-01; /* 0xbde31cc2 */
const PA5: f32 = 3.5478305072e-02; /* 0x3d1151b3 */
const PA6: f32 = -2.1663755178e-03; /* 0xbb0df9c0 */
const QA1: f32 = 1.0642088205e-01; /* 0x3dd9f331 */
const QA2: f32 = 5.4039794207e-01; /* 0x3f0a5785 */
const QA3: f32 = 7.1828655899e-02; /* 0x3d931ae7 */
const QA4: f32 = 1.2617121637e-01; /* 0x3e013307 */
const QA5: f32 = 1.3637083583e-02; /* 0x3c5f6e13 */
const QA6: f32 = 1.1984500103e-02; /* 0x3c445aa3 */
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
const RA0: f32 = -9.8649440333e-03; /* 0xbc21a093 */
const RA1: f32 = -6.9385856390e-01; /* 0xbf31a0b7 */
const RA2: f32 = -1.0558626175e+01; /* 0xc128f022 */
const RA3: f32 = -6.2375331879e+01; /* 0xc2798057 */
const RA4: f32 = -1.6239666748e+02; /* 0xc322658c */
const RA5: f32 = -1.8460508728e+02; /* 0xc3389ae7 */
const RA6: f32 = -8.1287437439e+01; /* 0xc2a2932b */
const RA7: f32 = -9.8143291473e+00; /* 0xc11d077e */
const SA1: f32 = 1.9651271820e+01; /* 0x419d35ce */
const SA2: f32 = 1.3765776062e+02; /* 0x4309a863 */
const SA3: f32 = 4.3456588745e+02; /* 0x43d9486f */
const SA4: f32 = 6.4538726807e+02; /* 0x442158c9 */
const SA5: f32 = 4.2900814819e+02; /* 0x43d6810b */
const SA6: f32 = 1.0863500214e+02; /* 0x42d9451f */
const SA7: f32 = 6.5702495575e+00; /* 0x40d23f7c */
const SA8: f32 = -6.0424413532e-02; /* 0xbd777f97 */
const RA0: f32 = -9.8649440333e-03; /* 0xbc21a093 */
const RA1: f32 = -6.9385856390e-01; /* 0xbf31a0b7 */
const RA2: f32 = -1.0558626175e+01; /* 0xc128f022 */
const RA3: f32 = -6.2375331879e+01; /* 0xc2798057 */
const RA4: f32 = -1.6239666748e+02; /* 0xc322658c */
const RA5: f32 = -1.8460508728e+02; /* 0xc3389ae7 */
const RA6: f32 = -8.1287437439e+01; /* 0xc2a2932b */
const RA7: f32 = -9.8143291473e+00; /* 0xc11d077e */
const SA1: f32 = 1.9651271820e+01; /* 0x419d35ce */
const SA2: f32 = 1.3765776062e+02; /* 0x4309a863 */
const SA3: f32 = 4.3456588745e+02; /* 0x43d9486f */
const SA4: f32 = 6.4538726807e+02; /* 0x442158c9 */
const SA5: f32 = 4.2900814819e+02; /* 0x43d6810b */
const SA6: f32 = 1.0863500214e+02; /* 0x42d9451f */
const SA7: f32 = 6.5702495575e+00; /* 0x40d23f7c */
const SA8: f32 = -6.0424413532e-02; /* 0xbd777f97 */
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
const RB0: f32 = -9.8649431020e-03; /* 0xbc21a092 */
const RB1: f32 = -7.9928326607e-01; /* 0xbf4c9dd4 */
const RB2: f32 = -1.7757955551e+01; /* 0xc18e104b */
const RB3: f32 = -1.6063638306e+02; /* 0xc320a2ea */
const RB4: f32 = -6.3756646729e+02; /* 0xc41f6441 */
const RB5: f32 = -1.0250950928e+03; /* 0xc480230b */
const RB6: f32 = -4.8351919556e+02; /* 0xc3f1c275 */
const SB1: f32 = 3.0338060379e+01; /* 0x41f2b459 */
const SB2: f32 = 3.2579251099e+02; /* 0x43a2e571 */
const SB3: f32 = 1.5367296143e+03; /* 0x44c01759 */
const SB4: f32 = 3.1998581543e+03; /* 0x4547fdbb */
const SB5: f32 = 2.5530502930e+03; /* 0x451f90ce */
const SB6: f32 = 4.7452853394e+02; /* 0x43ed43a7 */
const SB7: f32 = -2.2440952301e+01; /* 0xc1b38712 */
const RB0: f32 = -9.8649431020e-03; /* 0xbc21a092 */
const RB1: f32 = -7.9928326607e-01; /* 0xbf4c9dd4 */
const RB2: f32 = -1.7757955551e+01; /* 0xc18e104b */
const RB3: f32 = -1.6063638306e+02; /* 0xc320a2ea */
const RB4: f32 = -6.3756646729e+02; /* 0xc41f6441 */
const RB5: f32 = -1.0250950928e+03; /* 0xc480230b */
const RB6: f32 = -4.8351919556e+02; /* 0xc3f1c275 */
const SB1: f32 = 3.0338060379e+01; /* 0x41f2b459 */
const SB2: f32 = 3.2579251099e+02; /* 0x43a2e571 */
const SB3: f32 = 1.5367296143e+03; /* 0x44c01759 */
const SB4: f32 = 3.1998581543e+03; /* 0x4547fdbb */
const SB5: f32 = 2.5530502930e+03; /* 0x451f90ce */
const SB6: f32 = 4.7452853394e+02; /* 0x43ed43a7 */
const SB7: f32 = -2.2440952301e+01; /* 0xc1b38712 */
fn erfc1(x: f32) -> f32 {
let s: f32;
@ -89,9 +89,9 @@ fn erfc1(x: f32) -> f32 {
let q: f32;
s = fabsf(x) - 1.0;
p = PA0+s*(PA1+s*(PA2+s*(PA3+s*(PA4+s*(PA5+s*PA6)))));
q = 1.0+s*(QA1+s*(QA2+s*(QA3+s*(QA4+s*(QA5+s*QA6)))));
return 1.0 - ERX - p/q;
p = PA0 + s * (PA1 + s * (PA2 + s * (PA3 + s * (PA4 + s * (PA5 + s * PA6)))));
q = 1.0 + s * (QA1 + s * (QA2 + s * (QA3 + s * (QA4 + s * (QA5 + s * QA6)))));
return 1.0 - ERX - p / q;
}
fn erfc2(mut ix: u32, mut x: f32) -> f32 {
@ -100,31 +100,32 @@ fn erfc2(mut ix: u32, mut x: f32) -> f32 {
let big_s: f32;
let z: f32;
if ix < 0x3fa00000 { /* |x| < 1.25 */
if ix < 0x3fa00000 {
/* |x| < 1.25 */
return erfc1(x);
}
x = fabsf(x);
s = 1.0/(x*x);
if ix < 0x4036db6d { /* |x| < 1/0.35 */
r = RA0+s*(RA1+s*(RA2+s*(RA3+s*(RA4+s*(
RA5+s*(RA6+s*RA7))))));
big_s = 1.0+s*(SA1+s*(SA2+s*(SA3+s*(SA4+s*(
SA5+s*(SA6+s*(SA7+s*SA8)))))));
} else { /* |x| >= 1/0.35 */
r = RB0+s*(RB1+s*(RB2+s*(RB3+s*(RB4+s*(
RB5+s*RB6)))));
big_s = 1.0+s*(SB1+s*(SB2+s*(SB3+s*(SB4+s*(
SB5+s*(SB6+s*SB7))))));
s = 1.0 / (x * x);
if ix < 0x4036db6d {
/* |x| < 1/0.35 */
r = RA0 + s * (RA1 + s * (RA2 + s * (RA3 + s * (RA4 + s * (RA5 + s * (RA6 + s * RA7))))));
big_s = 1.0
+ s * (SA1
+ s * (SA2 + s * (SA3 + s * (SA4 + s * (SA5 + s * (SA6 + s * (SA7 + s * SA8)))))));
} else {
/* |x| >= 1/0.35 */
r = RB0 + s * (RB1 + s * (RB2 + s * (RB3 + s * (RB4 + s * (RB5 + s * RB6)))));
big_s =
1.0 + s * (SB1 + s * (SB2 + s * (SB3 + s * (SB4 + s * (SB5 + s * (SB6 + s * SB7))))));
}
ix = x.to_bits();
z = f32::from_bits(ix&0xffffe000);
z = f32::from_bits(ix & 0xffffe000);
expf(-z*z - 0.5625) * expf((z-x)*(z+x) + r/big_s)/x
expf(-z * z - 0.5625) * expf((z - x) * (z + x) + r / big_s) / x
}
pub fn erff(x: f32) -> f32
{
pub fn erff(x: f32) -> f32 {
let r: f32;
let s: f32;
let z: f32;
@ -133,25 +134,28 @@ pub fn erff(x: f32) -> f32
let sign: usize;
ix = x.to_bits();
sign = (ix>>31) as usize;
sign = (ix >> 31) as usize;
ix &= 0x7fffffff;
if ix >= 0x7f800000 {
/* erf(nan)=nan, erf(+-inf)=+-1 */
return 1.0-2.0*(sign as f32) + 1.0/x;
return 1.0 - 2.0 * (sign as f32) + 1.0 / x;
}
if ix < 0x3f580000 { /* |x| < 0.84375 */
if ix < 0x31800000 { /* |x| < 2**-28 */
if ix < 0x3f580000 {
/* |x| < 0.84375 */
if ix < 0x31800000 {
/* |x| < 2**-28 */
/*avoid underflow */
return 0.125*(8.0*x + EFX8*x);
return 0.125 * (8.0 * x + EFX8 * x);
}
z = x*x;
r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
y = r/s;
return x + x*y;
z = x * x;
r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4)));
s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5))));
y = r / s;
return x + x * y;
}
if ix < 0x40c00000 { /* |x| < 6 */
y = 1.0 - erfc2(ix,x);
if ix < 0x40c00000 {
/* |x| < 6 */
y = 1.0 - erfc2(ix, x);
} else {
let x1p_120 = f32::from_bits(0x03800000);
y = 1.0 - x1p_120;
@ -173,27 +177,31 @@ pub fn erfcf(x: f32) -> f32 {
let sign: usize;
ix = x.to_bits();
sign = (ix>>31) as usize;
sign = (ix >> 31) as usize;
ix &= 0x7fffffff;
if ix >= 0x7f800000 {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
return 2.0*(sign as f32) + 1.0/x;
return 2.0 * (sign as f32) + 1.0 / x;
}
if ix < 0x3f580000 { /* |x| < 0.84375 */
if ix < 0x23800000 { /* |x| < 2**-56 */
if ix < 0x3f580000 {
/* |x| < 0.84375 */
if ix < 0x23800000 {
/* |x| < 2**-56 */
return 1.0 - x;
}
z = x*x;
r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
y = r/s;
if sign != 0 || ix < 0x3e800000 { /* x < 1/4 */
return 1.0 - (x+x*y);
z = x * x;
r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4)));
s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5))));
y = r / s;
if sign != 0 || ix < 0x3e800000 {
/* x < 1/4 */
return 1.0 - (x + x * y);
}
return 0.5 - (x - 0.5 + x*y);
return 0.5 - (x - 0.5 + x * y);
}
if ix < 0x41e00000 { /* |x| < 28 */
if ix < 0x41e00000 {
/* |x| < 28 */
if sign != 0 {
return 2.0 - erfc2(ix, x);
} else {
@ -205,6 +213,6 @@ pub fn erfcf(x: f32) -> f32 {
if sign != 0 {
2.0 - x1p_120
} else {
x1p_120*x1p_120
x1p_120 * x1p_120
}
}

11
src/math/exp10.rs
View File

@ -2,18 +2,15 @@ use super::{exp2, modf, pow};
const LN10: f64 = 3.32192809488736234787031942948939;
const P10: &[f64] = &[
1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15
1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15,
];
pub fn exp10(x: f64) -> f64
{
pub fn exp10(x: f64) -> f64 {
let (mut y, n) = modf(x);
let u: u64 = n.to_bits();
/* fabs(n) < 16 without raising invalid on nan */
if (u>>52 & 0x7ff) < 0x3ff+4 {
if (u >> 52 & 0x7ff) < 0x3ff + 4 {
if y == 0.0 {
return P10[((n as isize) + 15) as usize];
}

7
src/math/exp10f.rs
View File

@ -3,17 +3,16 @@ use super::{exp2, exp2f, modff};
const LN10_F32: f32 = 3.32192809488736234787031942948939;
const LN10_F64: f64 = 3.32192809488736234787031942948939;
const P10: &[f32] = &[
1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,
];
pub fn exp10f(x: f32) -> f32 {
let (mut y, n) = modff(x);
let u = n.to_bits();
/* fabsf(n) < 8 without raising invalid on nan */
if (u>>23 & 0xff) < 0x7f+3 {
if (u >> 23 & 0xff) < 0x7f + 3 {
if y == 0.0 {
return P10[((n as isize) + 7) as usize]
return P10[((n as isize) + 7) as usize];
}
y = exp2f(LN10_F32 * y);
return y * P10[((n as isize) + 7) as usize];

4
src/math/fma.rs
View File

@ -43,7 +43,7 @@ fn mul(x: u64, y: u64) -> (u64, u64) {
t1 = xlo * ylo;
t2 = xlo * yhi + xhi * ylo;
t3 = xhi * yhi;
let lo = t1 + (t2 << 32);
let lo = t1.wrapping_add(t2 << 32);
let hi = t3 + (t2 >> 32) + (t1 > lo) as u64;
(hi, lo)
}
@ -116,7 +116,7 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 {
let mut nonzero: i32 = 1;
if samesign {
/* r += z */
rlo += zlo;
rlo = rlo.wrapping_add(zlo);
rhi += zhi + (rlo < zlo) as u64;
} else {
/* r -= z */

4
src/math/fmod.rs
View File

@ -46,7 +46,7 @@ pub fn fmod(x: f64, y: f64) -> f64 {
/* x mod y */
while ex > ey {
i = uxi - uyi;
i = uxi.wrapping_sub(uyi);
if i >> 63 == 0 {
if i == 0 {
return 0.0 * x;
@ -56,7 +56,7 @@ pub fn fmod(x: f64, y: f64) -> f64 {
uxi <<= 1;
ex -= 1;
}
i = uxi - uyi;
i = uxi.wrapping_sub(uyi);
if i >> 63 == 0 {
if i == 0 {
return 0.0 * x;

4
src/math/fmodf.rs
View File

@ -52,7 +52,7 @@ pub fn fmodf(x: f32, y: f32) -> f32 {
/* x mod y */
while ex > ey {
i = uxi - uyi;
i = uxi.wrapping_sub(uyi);
if i >> 31 == 0 {
if i == 0 {
return 0.0 * x;
@ -64,7 +64,7 @@ pub fn fmodf(x: f32, y: f32) -> f32 {
ex -= 1;
}
i = uxi - uyi;
i = uxi.wrapping_sub(uyi);
if i >> 31 == 0 {
if i == 0 {
return 0.0 * x;

6
src/math/frexp.rs
View File

@ -1,11 +1,11 @@
pub fn frexp(x: f64) -> (f64, isize) {
pub fn frexp(x: f64) -> (f64, i32) {
let mut y = x.to_bits();
let ee = ((y>>52) & 0x7ff) as isize;
let ee = ((y >> 52) & 0x7ff) as i32;
if ee == 0 {
if x != 0.0 {
let x1p64 = f64::from_bits(0x43f0000000000000);
let (x, e) = frexp(x*x1p64);
let (x, e) = frexp(x * x1p64);
return (x, e - 64);
}
return (x, 0);

8
src/math/frexpf.rs
View File

@ -1,11 +1,11 @@
pub fn frexpf(x: f32) -> (f32, isize) {
pub fn frexpf(x: f32) -> (f32, i32) {
let mut y = x.to_bits();
let ee: isize = ((y>>23) & 0xff) as isize;
let ee: i32 = ((y >> 23) & 0xff) as i32;
if ee == 0 {
if x != 0.0 {
let x1p64 = f32::from_bits(0x5f800000);
let (x, e) = frexpf(x*x1p64);
let (x, e) = frexpf(x * x1p64);
return (x, e - 64);
} else {
return (x, 0);
@ -17,5 +17,5 @@ pub fn frexpf(x: f32) -> (f32, isize) {
let e = ee - 0x7e;
y &= 0x807fffff;
y |= 0x3f000000;
return (f32::from_bits(y), e);
(f32::from_bits(y), e)
}

18
src/math/ilogb.rs
View File

@ -1,30 +1,30 @@
const FP_ILOGBNAN: isize = -1 - (((!0) >> 1));
const FP_ILOGB0: isize = FP_ILOGBNAN;
const FP_ILOGBNAN: i32 = -1 - ((!0) >> 1);
const FP_ILOGB0: i32 = FP_ILOGBNAN;
pub fn ilogb(x: f64) -> isize {
pub fn ilogb(x: f64) -> i32 {
let mut i: u64 = x.to_bits();
let e = ((i>>52) & 0x7ff) as isize;
let e = ((i >> 52) & 0x7ff) as i32;
if e == 0 {
i <<= 12;
if i == 0 {
force_eval!(0.0/0.0);
force_eval!(0.0 / 0.0);
return FP_ILOGB0;
}
/* subnormal x */
let mut e = -0x3ff;
while (i>>63) == 0 {
while (i >> 63) == 0 {
e -= 1;
i <<= 1;
}
return e;
}
if e == 0x7ff {
force_eval!(0.0/0.0);
if (i<<12) != 0 {
force_eval!(0.0 / 0.0);
if (i << 12) != 0 {
return FP_ILOGBNAN;
} else {
return isize::max_value();
return i32::max_value();
}
}
return e - 0x3ff;

18
src/math/ilogbf.rs
View File

@ -1,30 +1,30 @@
const FP_ILOGBNAN: isize = -1 - (((!0) >> 1));
const FP_ILOGB0: isize = FP_ILOGBNAN;
const FP_ILOGBNAN: i32 = -1 - ((!0) >> 1);
const FP_ILOGB0: i32 = FP_ILOGBNAN;
pub fn ilogbf(x: f32) -> isize {
pub fn ilogbf(x: f32) -> i32 {
let mut i = x.to_bits();
let e = ((i>>23) & 0xff) as isize;
let e = ((i >> 23) & 0xff) as i32;
if e == 0 {
i <<= 9;
if i == 0 {
force_eval!(0.0/0.0);
force_eval!(0.0 / 0.0);
return FP_ILOGB0;
}
/* subnormal x */
let mut e = -0x7f;
while (i>>31) == 0 {
while (i >> 31) == 0 {
e -= 1;
i <<= 1;
}
return e;
}
if e == 0xff {
force_eval!(0.0/0.0);
if (i<<9) != 0 {
force_eval!(0.0 / 0.0);
if (i << 9) != 0 {
return FP_ILOGBNAN;
} else {
return isize::max_value();
return i32::max_value();
}
}
return e - 0x7f;

380
src/math/j0.rs
View File

@ -54,9 +54,9 @@
* 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
*/
use super::{cos, get_low_word, get_high_word, fabs, log, sin, sqrt};
use super::{cos, fabs, get_high_word, get_low_word, log, sin, sqrt};
const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
const TPI: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
const TPI: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
/* common method when |x|>=2 */
fn common(ix: u32, x: f64, y0: bool) -> f64 {
@ -79,38 +79,37 @@ fn common(ix: u32, x: f64, y0: bool) -> f64 {
if y0 {
c = -c;
}
cc = s+c;
cc = s + c;
/* avoid overflow in 2*x, big ulp error when x>=0x1p1023 */
if ix < 0x7fe00000 {
ss = s-c;
z = -cos(2.0*x);
if s*c < 0.0 {
cc = z/ss;
ss = s - c;
z = -cos(2.0 * x);
if s * c < 0.0 {
cc = z / ss;
} else {
ss = z/cc;
ss = z / cc;
}
if ix < 0x48000000 {
if y0 {
ss = -ss;
}
cc = pzero(x)*cc-qzero(x)*ss;
cc = pzero(x) * cc - qzero(x) * ss;
}
}
return INVSQRTPI*cc/sqrt(x);
return INVSQRTPI * cc / sqrt(x);
}
/* R0/S0 on [0, 2.00] */
const R02: f64 = 1.56249999999999947958e-02; /* 0x3F8FFFFF, 0xFFFFFFFD */
const R02: f64 = 1.56249999999999947958e-02; /* 0x3F8FFFFF, 0xFFFFFFFD */
const R03: f64 = -1.89979294238854721751e-04; /* 0xBF28E6A5, 0xB61AC6E9 */
const R04: f64 = 1.82954049532700665670e-06; /* 0x3EBEB1D1, 0x0C503919 */
const R04: f64 = 1.82954049532700665670e-06; /* 0x3EBEB1D1, 0x0C503919 */
const R05: f64 = -4.61832688532103189199e-09; /* 0xBE33D5E7, 0x73D63FCE */
const S01: f64 = 1.56191029464890010492e-02; /* 0x3F8FFCE8, 0x82C8C2A4 */
const S02: f64 = 1.16926784663337450260e-04; /* 0x3F1EA6D2, 0xDD57DBF4 */
const S03: f64 = 5.13546550207318111446e-07; /* 0x3EA13B54, 0xCE84D5A9 */
const S04: f64 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
const S01: f64 = 1.56191029464890010492e-02; /* 0x3F8FFCE8, 0x82C8C2A4 */
const S02: f64 = 1.16926784663337450260e-04; /* 0x3F1EA6D2, 0xDD57DBF4 */
const S03: f64 = 5.13546550207318111446e-07; /* 0x3EA13B54, 0xCE84D5A9 */
const S04: f64 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
pub fn j0(mut x: f64) -> f64
{
pub fn j0(mut x: f64) -> f64 {
let z: f64;
let r: f64;
let s: f64;
@ -121,47 +120,49 @@ pub fn j0(mut x: f64) -> f64
/* j0(+-inf)=0, j0(nan)=nan */
if ix >= 0x7ff00000 {
return 1.0/(x*x);
return 1.0 / (x * x);
}
x = fabs(x);
if ix >= 0x40000000 { /* |x| >= 2 */
if ix >= 0x40000000 {
/* |x| >= 2 */
/* large ulp error near zeros: 2.4, 5.52, 8.6537,.. */
return common(ix,x,false);
return common(ix, x, false);
}
/* 1 - x*x/4 + x*x*R(x^2)/S(x^2) */
if ix >= 0x3f200000 { /* |x| >= 2**-13 */
if ix >= 0x3f200000 {
/* |x| >= 2**-13 */
/* up to 4ulp error close to 2 */
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
return (1.0+x/2.0)*(1.0-x/2.0) + z*(r/s);
z = x * x;
r = z * (R02 + z * (R03 + z * (R04 + z * R05)));
s = 1.0 + z * (S01 + z * (S02 + z * (S03 + z * S04)));
return (1.0 + x / 2.0) * (1.0 - x / 2.0) + z * (r / s);
}
/* 1 - x*x/4 */
/* prevent underflow */
/* inexact should be raised when x!=0, this is not done correctly */
if ix >= 0x38000000 { /* |x| >= 2**-127 */
x = 0.25*x*x;
if ix >= 0x38000000 {
/* |x| >= 2**-127 */
x = 0.25 * x * x;
}
return 1.0 - x;
}
const U00: f64 = -7.38042951086872317523e-02; /* 0xBFB2E4D6, 0x99CBD01F */
const U01: f64 = 1.76666452509181115538e-01; /* 0x3FC69D01, 0x9DE9E3FC */
const U02: f64 = -1.38185671945596898896e-02; /* 0xBF8C4CE8, 0xB16CFA97 */
const U03: f64 = 3.47453432093683650238e-04; /* 0x3F36C54D, 0x20B29B6B */
const U04: f64 = -3.81407053724364161125e-06; /* 0xBECFFEA7, 0x73D25CAD */
const U05: f64 = 1.95590137035022920206e-08; /* 0x3E550057, 0x3B4EABD4 */
const U06: f64 = -3.98205194132103398453e-11; /* 0xBDC5E43D, 0x693FB3C8 */
const V01: f64 = 1.27304834834123699328e-02; /* 0x3F8A1270, 0x91C9C71A */
const V02: f64 = 7.60068627350353253702e-05; /* 0x3F13ECBB, 0xF578C6C1 */
const V03: f64 = 2.59150851840457805467e-07; /* 0x3E91642D, 0x7FF202FD */
const V04: f64 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
const U00: f64 = -7.38042951086872317523e-02; /* 0xBFB2E4D6, 0x99CBD01F */
const U01: f64 = 1.76666452509181115538e-01; /* 0x3FC69D01, 0x9DE9E3FC */
const U02: f64 = -1.38185671945596898896e-02; /* 0xBF8C4CE8, 0xB16CFA97 */
const U03: f64 = 3.47453432093683650238e-04; /* 0x3F36C54D, 0x20B29B6B */
const U04: f64 = -3.81407053724364161125e-06; /* 0xBECFFEA7, 0x73D25CAD */
const U05: f64 = 1.95590137035022920206e-08; /* 0x3E550057, 0x3B4EABD4 */
const U06: f64 = -3.98205194132103398453e-11; /* 0xBDC5E43D, 0x693FB3C8 */
const V01: f64 = 1.27304834834123699328e-02; /* 0x3F8A1270, 0x91C9C71A */
const V02: f64 = 7.60068627350353253702e-05; /* 0x3F13ECBB, 0xF578C6C1 */
const V03: f64 = 2.59150851840457805467e-07; /* 0x3E91642D, 0x7FF202FD */
const V04: f64 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
pub fn y0(x: f64) -> f64
{
pub fn y0(x: f64) -> f64 {
let z: f64;
let u: f64;
let v: f64;
@ -172,30 +173,32 @@ pub fn y0(x: f64) -> f64
lx = get_low_word(x);
/* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */
if ((ix<<1) | lx) == 0 {
return -1.0/0.0;
if ((ix << 1) | lx) == 0 {
return -1.0 / 0.0;
}
if (ix>>31) != 0 {
return 0.0/0.0;
if (ix >> 31) != 0 {
return 0.0 / 0.0;
}
if ix >= 0x7ff00000 {
return 1.0/x;
return 1.0 / x;
}
if ix >= 0x40000000 { /* x >= 2 */
if ix >= 0x40000000 {
/* x >= 2 */
/* large ulp errors near zeros: 3.958, 7.086,.. */
return common(ix,x,true);
return common(ix, x, true);
}
/* U(x^2)/V(x^2) + (2/pi)*j0(x)*log(x) */
if ix >= 0x3e400000 { /* x >= 2**-27 */
if ix >= 0x3e400000 {
/* x >= 2**-27 */
/* large ulp error near the first zero, x ~= 0.89 */
z = x*x;
u = U00+z*(U01+z*(U02+z*(U03+z*(U04+z*(U05+z*U06)))));
v = 1.0+z*(V01+z*(V02+z*(V03+z*V04)));
return u/v + TPI*(j0(x)*log(x));
z = x * x;
u = U00 + z * (U01 + z * (U02 + z * (U03 + z * (U04 + z * (U05 + z * U06)))));
v = 1.0 + z * (V01 + z * (V02 + z * (V03 + z * V04)));
return u / v + TPI * (j0(x) * log(x));
}
return U00 + TPI*log(x);
return U00 + TPI * log(x);
}
/* The asymptotic expansions of pzero is
@ -207,72 +210,75 @@ pub fn y0(x: f64) -> f64
* and
* | pzero(x)-1-R/S | <= 2 ** ( -60.26)
*/
const PR8: [f64; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */
-8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */
-2.57063105679704847262e+02, /* 0xC0701102, 0x7B19E863 */
-2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */
-5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */
const PR8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */
-8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */
-2.57063105679704847262e+02, /* 0xC0701102, 0x7B19E863 */
-2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */
-5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */
];
const PS8: [f64; 5] = [
1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */
3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */
4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */
1.16752972564375915681e+05, /* 0x40FC810F, 0x8F9FA9BD */
4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */
1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */
3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */
4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */
1.16752972564375915681e+05, /* 0x40FC810F, 0x8F9FA9BD */
4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */
];
const PR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */
-7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */
-4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */
-6.76747652265167261021e+01, /* 0xC050EB2F, 0x5A7D1783 */
-3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */
-3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */
const PR5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */
-7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */
-4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */
-6.76747652265167261021e+01, /* 0xC050EB2F, 0x5A7D1783 */
-3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */
-3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */
];
const PS5: [f64; 5] = [
6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */
1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */
5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */
9.62544514357774460223e+03, /* 0x40C2CCB8, 0xFA76FA38 */
2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */
6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */
1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */
5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */
9.62544514357774460223e+03, /* 0x40C2CCB8, 0xFA76FA38 */
2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */
];
const PR3: [f64; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */
-7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */
-2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */
-2.19659774734883086467e+01, /* 0xC035F74A, 0x4CB94E14 */
-5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */
-3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */
const PR3: [f64; 6] = [
/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */
-7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */
-2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */
-2.19659774734883086467e+01, /* 0xC035F74A, 0x4CB94E14 */
-5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */
-3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */
];
const PS3: [f64; 5] = [
3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */
3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */
1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */
1.12799679856907414432e+03, /* 0x40919FFC, 0xB8C39B7E */
1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */
3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */
3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */
1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */
1.12799679856907414432e+03, /* 0x40919FFC, 0xB8C39B7E */
1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */
];
const PR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */
-7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */
-1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */
-7.63569613823527770791e+00, /* 0xC01E8AF3, 0xEDAFA7F3 */
-1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */
-3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */
const PR2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */
-7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */
-1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */
-7.63569613823527770791e+00, /* 0xC01E8AF3, 0xEDAFA7F3 */
-1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */
-3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */
];
const PS2: [f64; 5] = [
2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
];
fn pzero(x: f64) -> f64
{
fn pzero(x: f64) -> f64 {
let p: &[f64; 6];
let q: &[f64; 5];
let z: f64;
@ -282,17 +288,27 @@ fn pzero(x: f64) -> f64
ix = get_high_word(x);
ix &= 0x7fffffff;
if ix >= 0x40200000 {p = &PR8; q = &PS8;}
else if ix >= 0x40122E8B {p = &PR5; q = &PS5;}
else if ix >= 0x4006DB6D {p = &PR3; q = &PS3;}
else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0 + r/s;
if ix >= 0x40200000 {
p = &PR8;
q = &PS8;
} else if ix >= 0x40122E8B {
p = &PR5;
q = &PS5;
} else if ix >= 0x4006DB6D {
p = &PR3;
q = &PS3;
} else
/*ix >= 0x40000000*/
{
p = &PR2;
q = &PS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
return 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
@ -302,76 +318,79 @@ fn pzero(x: f64) -> f64
* and
* | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
*/
const QR8: [f64; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */
1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */
5.57673380256401856059e+02, /* 0x40816D63, 0x15301825 */
8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */
3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */
const QR8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */
1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */
5.57673380256401856059e+02, /* 0x40816D63, 0x15301825 */
8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */
3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */
];
const QS8: [f64; 6] = [
1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */
8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */
1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */
8.03309257119514397345e+05, /* 0x412883DA, 0x83A52B43 */
8.40501579819060512818e+05, /* 0x4129A66B, 0x28DE0B3D */
-3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */
1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */
8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */
1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */
8.03309257119514397345e+05, /* 0x412883DA, 0x83A52B43 */
8.40501579819060512818e+05, /* 0x4129A66B, 0x28DE0B3D */
-3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */
];
const QR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */
5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */
1.35111577286449829671e+02, /* 0x4060E392, 0x0A8788E9 */
1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */
1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */
const QR5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */
5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */
1.35111577286449829671e+02, /* 0x4060E392, 0x0A8788E9 */
1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */
1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */
];
const QS5: [f64; 6] = [
8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */
2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */
1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */
5.67511122894947329769e+04, /* 0x40EBB5E3, 0x97E02372 */
3.59767538425114471465e+04, /* 0x40E19118, 0x1F7A54A0 */
-5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */
8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */
2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */
1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */
5.67511122894947329769e+04, /* 0x40EBB5E3, 0x97E02372 */
3.59767538425114471465e+04, /* 0x40E19118, 0x1F7A54A0 */
-5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */
];
const QR3: [f64; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */
7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */
3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */
4.26218440745412650017e+01, /* 0x40454F98, 0x962DAEDD */
1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */
1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */
const QR3: [f64; 6] = [
/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */
7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */
3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */
4.26218440745412650017e+01, /* 0x40454F98, 0x962DAEDD */
1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */
1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */
];
const QS3: [f64; 6] = [
4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */
7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */
3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */
6.46042516752568917582e+03, /* 0x40B93C6C, 0xD7C76A28 */
2.51633368920368957333e+03, /* 0x40A3A8AA, 0xD94FB1C0 */
-1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */
4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */
7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */
3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */
6.46042516752568917582e+03, /* 0x40B93C6C, 0xD7C76A28 */
2.51633368920368957333e+03, /* 0x40A3A8AA, 0xD94FB1C0 */
-1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */
];
const QR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */
7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */
1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */
1.44956029347885735348e+01, /* 0x402CFDBF, 0xAAF96FE5 */
3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */
1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */
const QR2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */
7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */
1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */
1.44956029347885735348e+01, /* 0x402CFDBF, 0xAAF96FE5 */
3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */
1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */
];
const QS2: [f64; 6] = [
3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
-5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
-5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
];
fn qzero(x: f64) -> f64
{
fn qzero(x: f64) -> f64 {
let p: &[f64; 6];
let q: &[f64; 6];
let s: f64;
@ -381,12 +400,23 @@ fn qzero(x: f64) -> f64
ix = get_high_word(x);
ix &= 0x7fffffff;
if ix >= 0x40200000 {p = &QR8; q = &QS8;}
else if ix >= 0x40122E8B {p = &QR5; q = &QS5;}
else if ix >= 0x4006DB6D {p = &QR3; q = &QS3;}
else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-0.125 + r/s)/x;
if ix >= 0x40200000 {
p = &QR8;
q = &QS8;
} else if ix >= 0x40122E8B {
p = &QR5;
q = &QS5;
} else if ix >= 0x4006DB6D {
p = &QR3;
q = &QS3;
} else
/*ix >= 0x40000000*/
{
p = &QR2;
q = &QS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
return (-0.125 + r / s) / x;
}

377
src/math/j0f.rs
View File

@ -16,10 +16,9 @@
use super::{cosf, fabsf, logf, sinf, sqrtf};
const INVSQRTPI: f32 = 5.6418961287e-01; /* 0x3f106ebb */
const TPI: f32 = 6.3661974669e-01; /* 0x3f22f983 */
const TPI: f32 = 6.3661974669e-01; /* 0x3f22f983 */
fn common(ix: u32, x: f32, y0: bool) -> f32
{
fn common(ix: u32, x: f32, y0: bool) -> f32 {
let z: f32;
let s: f32;
let mut c: f32;
@ -34,37 +33,36 @@ fn common(ix: u32, x: f32, y0: bool) -> f32
if y0 {
c = -c;
}
cc = s+c;
cc = s + c;
if ix < 0x7f000000 {
ss = s-c;
z = -cosf(2.0*x);
if s*c < 0.0 {
cc = z/ss;
ss = s - c;
z = -cosf(2.0 * x);
if s * c < 0.0 {
cc = z / ss;
} else {
ss = z/cc;
ss = z / cc;
}
if ix < 0x58800000 {
if y0 {
ss = -ss;
}
cc = pzerof(x)*cc-qzerof(x)*ss;
cc = pzerof(x) * cc - qzerof(x) * ss;
}
}
return INVSQRTPI*cc/sqrtf(x);
return INVSQRTPI * cc / sqrtf(x);
}
/* R0/S0 on [0, 2.00] */
const R02: f32 = 1.5625000000e-02; /* 0x3c800000 */
const R02: f32 = 1.5625000000e-02; /* 0x3c800000 */
const R03: f32 = -1.8997929874e-04; /* 0xb947352e */
const R04: f32 = 1.8295404516e-06; /* 0x35f58e88 */
const R04: f32 = 1.8295404516e-06; /* 0x35f58e88 */
const R05: f32 = -4.6183270541e-09; /* 0xb19eaf3c */
const S01: f32 = 1.5619102865e-02; /* 0x3c7fe744 */
const S02: f32 = 1.1692678527e-04; /* 0x38f53697 */
const S03: f32 = 5.1354652442e-07; /* 0x3509daa6 */
const S04: f32 = 1.1661400734e-09; /* 0x30a045e8 */
const S01: f32 = 1.5619102865e-02; /* 0x3c7fe744 */
const S02: f32 = 1.1692678527e-04; /* 0x38f53697 */
const S03: f32 = 5.1354652442e-07; /* 0x3509daa6 */
const S04: f32 = 1.1661400734e-09; /* 0x30a045e8 */
pub fn j0f(mut x: f32) -> f32
{
pub fn j0f(mut x: f32) -> f32 {
let z: f32;
let r: f32;
let s: f32;
@ -73,41 +71,43 @@ pub fn j0f(mut x: f32) -> f32
ix = x.to_bits();
ix &= 0x7fffffff;
if ix >= 0x7f800000 {
return 1.0/(x*x);
return 1.0 / (x * x);
}
x = fabsf(x);
if ix >= 0x40000000 { /* |x| >= 2 */
if ix >= 0x40000000 {
/* |x| >= 2 */
/* large ulp error near zeros */
return common(ix, x, false);
}
if ix >= 0x3a000000 { /* |x| >= 2**-11 */
if ix >= 0x3a000000 {
/* |x| >= 2**-11 */
/* up to 4ulp error near 2 */
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
return (1.0+x/2.0)*(1.0-x/2.0) + z*(r/s);
z = x * x;
r = z * (R02 + z * (R03 + z * (R04 + z * R05)));
s = 1.0 + z * (S01 + z * (S02 + z * (S03 + z * S04)));
return (1.0 + x / 2.0) * (1.0 - x / 2.0) + z * (r / s);
}
if ix >= 0x21800000 { /* |x| >= 2**-60 */
x = 0.25*x*x;
if ix >= 0x21800000 {
/* |x| >= 2**-60 */
x = 0.25 * x * x;
}
return 1.0 - x;
}
const U00: f32 = -7.3804296553e-02; /* 0xbd9726b5 */
const U01: f32 = 1.7666645348e-01; /* 0x3e34e80d */
const U02: f32 = -1.3818567619e-02; /* 0xbc626746 */
const U03: f32 = 3.4745343146e-04; /* 0x39b62a69 */
const U04: f32 = -3.8140706238e-06; /* 0xb67ff53c */
const U05: f32 = 1.9559013964e-08; /* 0x32a802ba */
const U06: f32 = -3.9820518410e-11; /* 0xae2f21eb */
const V01: f32 = 1.2730483897e-02; /* 0x3c509385 */
const V02: f32 = 7.6006865129e-05; /* 0x389f65e0 */
const V03: f32 = 2.5915085189e-07; /* 0x348b216c */
const V04: f32 = 4.4111031494e-10; /* 0x2ff280c2 */
const U00: f32 = -7.3804296553e-02; /* 0xbd9726b5 */
const U01: f32 = 1.7666645348e-01; /* 0x3e34e80d */
const U02: f32 = -1.3818567619e-02; /* 0xbc626746 */
const U03: f32 = 3.4745343146e-04; /* 0x39b62a69 */
const U04: f32 = -3.8140706238e-06; /* 0xb67ff53c */
const U05: f32 = 1.9559013964e-08; /* 0x32a802ba */
const U06: f32 = -3.9820518410e-11; /* 0xae2f21eb */
const V01: f32 = 1.2730483897e-02; /* 0x3c509385 */
const V02: f32 = 7.6006865129e-05; /* 0x389f65e0 */
const V03: f32 = 2.5915085189e-07; /* 0x348b216c */
const V04: f32 = 4.4111031494e-10; /* 0x2ff280c2 */
pub fn y0f(x: f32) -> f32
{
pub fn y0f(x: f32) -> f32 {
let z: f32;
let u: f32;
let v: f32;
@ -115,26 +115,28 @@ pub fn y0f(x: f32) -> f32
ix = x.to_bits();
if (ix & 0x7fffffff) == 0 {
return -1.0/0.0;
return -1.0 / 0.0;
}
if (ix>>31) !=0 {
return 0.0/0.0;
if (ix >> 31) != 0 {
return 0.0 / 0.0;
}
if ix >= 0x7f800000 {
return 1.0/x;
return 1.0 / x;
}
if ix >= 0x40000000 { /* |x| >= 2.0 */
if ix >= 0x40000000 {
/* |x| >= 2.0 */
/* large ulp error near zeros */
return common(ix,x,true);
return common(ix, x, true);
}
if ix >= 0x39000000 { /* x >= 2**-13 */
if ix >= 0x39000000 {
/* x >= 2**-13 */
/* large ulp error at x ~= 0.89 */
z = x*x;
u = U00+z*(U01+z*(U02+z*(U03+z*(U04+z*(U05+z*U06)))));
v = 1.0+z*(V01+z*(V02+z*(V03+z*V04)));
return u/v + TPI*(j0f(x)*logf(x));
z = x * x;
u = U00 + z * (U01 + z * (U02 + z * (U03 + z * (U04 + z * (U05 + z * U06)))));
v = 1.0 + z * (V01 + z * (V02 + z * (V03 + z * V04)));
return u / v + TPI * (j0f(x) * logf(x));
}
return U00 + TPI*logf(x);
return U00 + TPI * logf(x);
}
/* The asymptotic expansions of pzero is
@ -146,71 +148,74 @@ pub fn y0f(x: f32) -> f32
* and
* | pzero(x)-1-R/S | <= 2 ** ( -60.26)
*/
const PR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
-7.0312500000e-02, /* 0xbd900000 */
-8.0816707611e+00, /* 0xc1014e86 */
-2.5706311035e+02, /* 0xc3808814 */
-2.4852163086e+03, /* 0xc51b5376 */
-5.2530439453e+03, /* 0xc5a4285a */
const PR8: [f32; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
-7.0312500000e-02, /* 0xbd900000 */
-8.0816707611e+00, /* 0xc1014e86 */
-2.5706311035e+02, /* 0xc3808814 */
-2.4852163086e+03, /* 0xc51b5376 */
-5.2530439453e+03, /* 0xc5a4285a */
];
const PS8: [f32; 5] = [
1.1653436279e+02, /* 0x42e91198 */
3.8337448730e+03, /* 0x456f9beb */
4.0597855469e+04, /* 0x471e95db */
1.1675296875e+05, /* 0x47e4087c */
4.7627726562e+04, /* 0x473a0bba */
1.1653436279e+02, /* 0x42e91198 */
3.8337448730e+03, /* 0x456f9beb */
4.0597855469e+04, /* 0x471e95db */
1.1675296875e+05, /* 0x47e4087c */
4.7627726562e+04, /* 0x473a0bba */
];
const PR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.1412546255e-11, /* 0xad48c58a */
-7.0312492549e-02, /* 0xbd8fffff */
-4.1596107483e+00, /* 0xc0851b88 */
-6.7674766541e+01, /* 0xc287597b */
-3.3123129272e+02, /* 0xc3a59d9b */
-3.4643338013e+02, /* 0xc3ad3779 */
const PR5: [f32; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
-1.1412546255e-11, /* 0xad48c58a */
-7.0312492549e-02, /* 0xbd8fffff */
-4.1596107483e+00, /* 0xc0851b88 */
-6.7674766541e+01, /* 0xc287597b */
-3.3123129272e+02, /* 0xc3a59d9b */
-3.4643338013e+02, /* 0xc3ad3779 */
];
const PS5: [f32; 5] = [
6.0753936768e+01, /* 0x42730408 */
1.0512523193e+03, /* 0x44836813 */
5.9789707031e+03, /* 0x45bad7c4 */
9.6254453125e+03, /* 0x461665c8 */
2.4060581055e+03, /* 0x451660ee */
6.0753936768e+01, /* 0x42730408 */
1.0512523193e+03, /* 0x44836813 */
5.9789707031e+03, /* 0x45bad7c4 */
9.6254453125e+03, /* 0x461665c8 */
2.4060581055e+03, /* 0x451660ee */
];
const PR3: [f32; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.5470459075e-09, /* 0xb12f081b */
-7.0311963558e-02, /* 0xbd8fffb8 */
-2.4090321064e+00, /* 0xc01a2d95 */
-2.1965976715e+01, /* 0xc1afba52 */
-5.8079170227e+01, /* 0xc2685112 */
-3.1447946548e+01, /* 0xc1fb9565 */
const PR3: [f32; 6] = [
/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-2.5470459075e-09, /* 0xb12f081b */
-7.0311963558e-02, /* 0xbd8fffb8 */
-2.4090321064e+00, /* 0xc01a2d95 */
-2.1965976715e+01, /* 0xc1afba52 */
-5.8079170227e+01, /* 0xc2685112 */
-3.1447946548e+01, /* 0xc1fb9565 */
];
const PS3: [f32; 5] = [
3.5856033325e+01, /* 0x420f6c94 */
3.6151397705e+02, /* 0x43b4c1ca */
1.1936077881e+03, /* 0x44953373 */
1.1279968262e+03, /* 0x448cffe6 */
1.7358093262e+02, /* 0x432d94b8 */
3.5856033325e+01, /* 0x420f6c94 */
3.6151397705e+02, /* 0x43b4c1ca */
1.1936077881e+03, /* 0x44953373 */
1.1279968262e+03, /* 0x448cffe6 */
1.7358093262e+02, /* 0x432d94b8 */
];
const PR2: [f32; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.8753431271e-08, /* 0xb3be98b7 */
-7.0303097367e-02, /* 0xbd8ffb12 */
-1.4507384300e+00, /* 0xbfb9b1cc */
-7.6356959343e+00, /* 0xc0f4579f */
-1.1193166733e+01, /* 0xc1331736 */
-3.2336456776e+00, /* 0xc04ef40d */
const PR2: [f32; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
-8.8753431271e-08, /* 0xb3be98b7 */
-7.0303097367e-02, /* 0xbd8ffb12 */
-1.4507384300e+00, /* 0xbfb9b1cc */
-7.6356959343e+00, /* 0xc0f4579f */
-1.1193166733e+01, /* 0xc1331736 */
-3.2336456776e+00, /* 0xc04ef40d */
];
const PS2: [f32; 5] = [
2.2220300674e+01, /* 0x41b1c32d */
1.3620678711e+02, /* 0x430834f0 */
2.7047027588e+02, /* 0x43873c32 */
1.5387539673e+02, /* 0x4319e01a */
1.4657617569e+01, /* 0x416a859a */
2.2220300674e+01, /* 0x41b1c32d */
1.3620678711e+02, /* 0x430834f0 */
2.7047027588e+02, /* 0x43873c32 */
1.5387539673e+02, /* 0x4319e01a */
1.4657617569e+01, /* 0x416a859a */
];
fn pzerof(x: f32) -> f32
{
fn pzerof(x: f32) -> f32 {
let p: &[f32; 6];
let q: &[f32; 5];
let z: f32;
@ -220,17 +225,27 @@ fn pzerof(x: f32) -> f32
ix = x.to_bits();
ix &= 0x7fffffff;
if ix >= 0x41000000 {p = &PR8; q = &PS8;}
else if ix >= 0x409173eb {p = &PR5; q = &PS5;}
else if ix >= 0x4036d917 {p = &PR3; q = &PS3;}
else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0 + r/s;
if ix >= 0x41000000 {
p = &PR8;
q = &PS8;
} else if ix >= 0x409173eb {
p = &PR5;
q = &PS5;
} else if ix >= 0x4036d917 {
p = &PR3;
q = &PS3;
} else
/*ix >= 0x40000000*/
{
p = &PR2;
q = &PS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
return 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
@ -240,76 +255,79 @@ fn pzerof(x: f32) -> f32
* and
* | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
*/
const QR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
7.3242187500e-02, /* 0x3d960000 */
1.1768206596e+01, /* 0x413c4a93 */
5.5767340088e+02, /* 0x440b6b19 */
8.8591972656e+03, /* 0x460a6cca */
3.7014625000e+04, /* 0x471096a0 */
const QR8: [f32; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
7.3242187500e-02, /* 0x3d960000 */
1.1768206596e+01, /* 0x413c4a93 */
5.5767340088e+02, /* 0x440b6b19 */
8.8591972656e+03, /* 0x460a6cca */
3.7014625000e+04, /* 0x471096a0 */
];
const QS8: [f32; 6] = [
1.6377603149e+02, /* 0x4323c6aa */
8.0983447266e+03, /* 0x45fd12c2 */
1.4253829688e+05, /* 0x480b3293 */
8.0330925000e+05, /* 0x49441ed4 */
8.4050156250e+05, /* 0x494d3359 */
-3.4389928125e+05, /* 0xc8a7eb69 */
1.6377603149e+02, /* 0x4323c6aa */
8.0983447266e+03, /* 0x45fd12c2 */
1.4253829688e+05, /* 0x480b3293 */
8.0330925000e+05, /* 0x49441ed4 */
8.4050156250e+05, /* 0x494d3359 */
-3.4389928125e+05, /* 0xc8a7eb69 */
];
const QR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.8408595828e-11, /* 0x2da1ec79 */
7.3242180049e-02, /* 0x3d95ffff */
5.8356351852e+00, /* 0x40babd86 */
1.3511157227e+02, /* 0x43071c90 */
1.0272437744e+03, /* 0x448067cd */
1.9899779053e+03, /* 0x44f8bf4b */
const QR5: [f32; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
1.8408595828e-11, /* 0x2da1ec79 */
7.3242180049e-02, /* 0x3d95ffff */
5.8356351852e+00, /* 0x40babd86 */
1.3511157227e+02, /* 0x43071c90 */
1.0272437744e+03, /* 0x448067cd */
1.9899779053e+03, /* 0x44f8bf4b */
];
const QS5: [f32; 6] = [
8.2776611328e+01, /* 0x42a58da0 */
2.0778142090e+03, /* 0x4501dd07 */
1.8847289062e+04, /* 0x46933e94 */
5.6751113281e+04, /* 0x475daf1d */
3.5976753906e+04, /* 0x470c88c1 */
-5.3543427734e+03, /* 0xc5a752be */
8.2776611328e+01, /* 0x42a58da0 */
2.0778142090e+03, /* 0x4501dd07 */
1.8847289062e+04, /* 0x46933e94 */
5.6751113281e+04, /* 0x475daf1d */
3.5976753906e+04, /* 0x470c88c1 */
-5.3543427734e+03, /* 0xc5a752be */
];
const QR3: [f32; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.3774099900e-09, /* 0x3196681b */
7.3241114616e-02, /* 0x3d95ff70 */
3.3442313671e+00, /* 0x405607e3 */
4.2621845245e+01, /* 0x422a7cc5 */
1.7080809021e+02, /* 0x432acedf */
1.6673394775e+02, /* 0x4326bbe4 */
const QR3: [f32; 6] = [
/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
4.3774099900e-09, /* 0x3196681b */
7.3241114616e-02, /* 0x3d95ff70 */
3.3442313671e+00, /* 0x405607e3 */
4.2621845245e+01, /* 0x422a7cc5 */
1.7080809021e+02, /* 0x432acedf */
1.6673394775e+02, /* 0x4326bbe4 */
];
const QS3: [f32; 6] = [
4.8758872986e+01, /* 0x42430916 */
7.0968920898e+02, /* 0x44316c1c */
3.7041481934e+03, /* 0x4567825f */
6.4604252930e+03, /* 0x45c9e367 */
2.5163337402e+03, /* 0x451d4557 */
-1.4924745178e+02, /* 0xc3153f59 */
4.8758872986e+01, /* 0x42430916 */
7.0968920898e+02, /* 0x44316c1c */
3.7041481934e+03, /* 0x4567825f */
6.4604252930e+03, /* 0x45c9e367 */
2.5163337402e+03, /* 0x451d4557 */
-1.4924745178e+02, /* 0xc3153f59 */
];
const QR2: [f32; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.5044444979e-07, /* 0x342189db */
7.3223426938e-02, /* 0x3d95f62a */
1.9981917143e+00, /* 0x3fffc4bf */
1.4495602608e+01, /* 0x4167edfd */
3.1666231155e+01, /* 0x41fd5471 */
1.6252708435e+01, /* 0x4182058c */
const QR2: [f32; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.5044444979e-07, /* 0x342189db */
7.3223426938e-02, /* 0x3d95f62a */
1.9981917143e+00, /* 0x3fffc4bf */
1.4495602608e+01, /* 0x4167edfd */
3.1666231155e+01, /* 0x41fd5471 */
1.6252708435e+01, /* 0x4182058c */
];
const QS2: [f32; 6] = [
3.0365585327e+01, /* 0x41f2ecb8 */
2.6934811401e+02, /* 0x4386ac8f */
8.4478375244e+02, /* 0x44533229 */
8.8293585205e+02, /* 0x445cbbe5 */
2.1266638184e+02, /* 0x4354aa98 */
-5.3109550476e+00, /* 0xc0a9f358 */
3.0365585327e+01, /* 0x41f2ecb8 */
2.6934811401e+02, /* 0x4386ac8f */
8.4478375244e+02, /* 0x44533229 */
8.8293585205e+02, /* 0x445cbbe5 */
2.1266638184e+02, /* 0x4354aa98 */
-5.3109550476e+00, /* 0xc0a9f358 */
];
fn qzerof(x: f32) -> f32
{
fn qzerof(x: f32) -> f32 {
let p: &[f32; 6];
let q: &[f32; 6];
let s: f32;
@ -319,12 +337,23 @@ fn qzerof(x: f32) -> f32
ix = x.to_bits();
ix &= 0x7fffffff;
if ix >= 0x41000000 {p = &QR8; q = &QS8;}
else if ix >= 0x409173eb {p = &QR5; q = &QS5;}
else if ix >= 0x4036d917 {p = &QR3; q = &QS3;}
else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-0.125 + r/s)/x;
if ix >= 0x41000000 {
p = &QR8;
q = &QS8;
} else if ix >= 0x409173eb {
p = &QR5;
q = &QS5;
} else if ix >= 0x4036d917 {
p = &QR3;
q = &QS3;
} else
/*ix >= 0x40000000*/
{
p = &QR2;
q = &QS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
return (-0.125 + r / s) / x;
}

371
src/math/j1.rs
View File

@ -54,13 +54,12 @@
* by method mentioned above.
*/
use super::{cos, get_high_word, get_low_word, fabs, log, sin, sqrt};
use super::{cos, fabs, get_high_word, get_low_word, log, sin, sqrt};
const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
const TPI: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
const TPI: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
fn common(ix: u32, x: f64, y1: bool, sign: bool) -> f64
{
fn common(ix: u32, x: f64, y1: bool, sign: bool) -> f64 {
let z: f64;
let mut s: f64;
let c: f64;
@ -80,42 +79,41 @@ fn common(ix: u32, x: f64, y1: bool, sign: bool) -> f64
s = -s;
}
c = cos(x);
cc = s-c;
cc = s - c;
if ix < 0x7fe00000 {
/* avoid overflow in 2*x */
ss = -s-c;
z = cos(2.0*x);
if s*c > 0.0 {
cc = z/ss;
ss = -s - c;
z = cos(2.0 * x);
if s * c > 0.0 {
cc = z / ss;
} else {
ss = z/cc;
ss = z / cc;
}
if ix < 0x48000000 {
if y1 {
ss = -ss;
}
cc = pone(x)*cc-qone(x)*ss;
cc = pone(x) * cc - qone(x) * ss;
}
}
if sign {
cc = -cc;
}
return INVSQRTPI*cc/sqrt(x);
return INVSQRTPI * cc / sqrt(x);
}
/* R0/S0 on [0,2] */
const R00: f64 = -6.25000000000000000000e-02; /* 0xBFB00000, 0x00000000 */
const R01: f64 = 1.40705666955189706048e-03; /* 0x3F570D9F, 0x98472C61 */
const R01: f64 = 1.40705666955189706048e-03; /* 0x3F570D9F, 0x98472C61 */
const R02: f64 = -1.59955631084035597520e-05; /* 0xBEF0C5C6, 0xBA169668 */
const R03: f64 = 4.96727999609584448412e-08; /* 0x3E6AAAFA, 0x46CA0BD9 */
const S01: f64 = 1.91537599538363460805e-02; /* 0x3F939D0B, 0x12637E53 */
const S02: f64 = 1.85946785588630915560e-04; /* 0x3F285F56, 0xB9CDF664 */
const S03: f64 = 1.17718464042623683263e-06; /* 0x3EB3BFF8, 0x333F8498 */
const S04: f64 = 5.04636257076217042715e-09; /* 0x3E35AC88, 0xC97DFF2C */
const S05: f64 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
const R03: f64 = 4.96727999609584448412e-08; /* 0x3E6AAAFA, 0x46CA0BD9 */
const S01: f64 = 1.91537599538363460805e-02; /* 0x3F939D0B, 0x12637E53 */
const S02: f64 = 1.85946785588630915560e-04; /* 0x3F285F56, 0xB9CDF664 */
const S03: f64 = 1.17718464042623683263e-06; /* 0x3EB3BFF8, 0x333F8498 */
const S04: f64 = 5.04636257076217042715e-09; /* 0x3E35AC88, 0xC97DFF2C */
const S05: f64 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
pub fn j1(x: f64) -> f64
{
pub fn j1(x: f64) -> f64 {
let mut z: f64;
let r: f64;
let s: f64;
@ -123,43 +121,44 @@ pub fn j1(x: f64) -> f64
let sign: bool;
ix = get_high_word(x);
sign = (ix>>31) != 0;
sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
if ix >= 0x7ff00000 {
return 1.0/(x*x);
return 1.0 / (x * x);
}
if ix >= 0x40000000 { /* |x| >= 2 */
if ix >= 0x40000000 {
/* |x| >= 2 */
return common(ix, fabs(x), false, sign);
}
if ix >= 0x38000000 { /* |x| >= 2**-127 */
z = x*x;
r = z*(R00+z*(R01+z*(R02+z*R03)));
s = 1.0+z*(S01+z*(S02+z*(S03+z*(S04+z*S05))));
z = r/s;
if ix >= 0x38000000 {
/* |x| >= 2**-127 */
z = x * x;
r = z * (R00 + z * (R01 + z * (R02 + z * R03)));
s = 1.0 + z * (S01 + z * (S02 + z * (S03 + z * (S04 + z * S05))));
z = r / s;
} else {
/* avoid underflow, raise inexact if x!=0 */
z = x;
}
return (0.5 + z)*x;
return (0.5 + z) * x;
}
const U0: [f64; 5] = [
-1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
-1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
-9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
-1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
-1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
-9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
];
const V0: [f64; 5] = [
1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */
1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */
1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
];
pub fn y1(x: f64) -> f64
{
pub fn y1(x: f64) -> f64 {
let z: f64;
let u: f64;
let v: f64;
@ -170,26 +169,28 @@ pub fn y1(x: f64) -> f64
lx = get_low_word(x);
/* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */
if (ix<<1 | lx) == 0 {
return -1.0/0.0;
if (ix << 1 | lx) == 0 {
return -1.0 / 0.0;
}
if (ix>>31) != 0 {
return 0.0/0.0;
if (ix >> 31) != 0 {
return 0.0 / 0.0;
}
if ix >= 0x7ff00000 {
return 1.0/x;
return 1.0 / x;
}
if ix >= 0x40000000 { /* x >= 2 */
if ix >= 0x40000000 {
/* x >= 2 */
return common(ix, x, true, false);
}
if ix < 0x3c900000 { /* x < 2**-54 */
return -TPI/x;
if ix < 0x3c900000 {
/* x < 2**-54 */
return -TPI / x;
}
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + TPI*(j1(x)*log(x)-1.0/x);
z = x * x;
u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * U0[4])));
v = 1.0 + z * (V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * V0[4]))));
return x * (u / v) + TPI * (j1(x) * log(x) - 1.0 / x);
}
/* For x >= 8, the asymptotic expansions of pone is
@ -202,72 +203,74 @@ pub fn y1(x: f64) -> f64
* | 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.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */
3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
const PR8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */
3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
];
const PS8: [f64; 5] = [
1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */
3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */
3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
];
const PR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */
5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
const PR5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */
5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
];
const PS5: [f64; 5] = [
5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */
1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */
1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
];
const PR3: [f64; 6] = [
3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */
9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */
9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
];
const PS3: [f64; 5] = [
3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */
1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */
1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
];
const PR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */
1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
const PR2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */
1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
];
const PS2: [f64; 5] = [
2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
];
fn pone(x: f64) -> f64
{
fn pone(x: f64) -> f64 {
let p: &[f64; 6];
let q: &[f64; 5];
let z: f64;
@ -277,14 +280,25 @@ fn pone(x: f64) -> f64
ix = get_high_word(x);
ix &= 0x7fffffff;
if ix >= 0x40200000 {p = &PR8; q = &PS8;}
else if ix >= 0x40122E8B {p = &PR5; q = &PS5;}
else if ix >= 0x4006DB6D {p = &PR3; q = &PS3;}
else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0+ r/s;
if ix >= 0x40200000 {
p = &PR8;
q = &PS8;
} else if ix >= 0x40122E8B {
p = &PR5;
q = &PS5;
} else if ix >= 0x4006DB6D {
p = &PR3;
q = &PS3;
} else
/*ix >= 0x40000000*/
{
p = &PR2;
q = &PS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
return 1.0 + r / s;
}
/* For x >= 8, the asymptotic expansions of qone is
@ -297,76 +311,78 @@ fn pone(x: f64) -> f64
* | 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.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
-1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
-7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */
-1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
-4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
const QR8: [f64; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
-1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
-7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */
-1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
-4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
];
const QS8: [f64; 6] = [
1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */
6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */
-2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */
6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */
-2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
];
const QR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
-1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
-8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
-1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */
-1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
-2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
const QR5: [f64; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
-1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
-8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
-1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */
-1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
-2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
];
const QS5: [f64; 6] = [
8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */
2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */
-4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */
2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */
-4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
];
const QR3: [f64; 6] = [
-5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
-1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
-4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
-5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */
-2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
-2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
-5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
-1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
-4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
-5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */
-2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
-2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
];
const QS3: [f64; 6] = [
4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */
1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */
-1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */
1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */
-1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
];
const QR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
-1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
-2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
-1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */
-4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
-2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
const QR2: [f64; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
-1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
-2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
-1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */
-4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
-2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
];
const QS2: [f64; 6] = [
2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
-4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
-4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
];
fn qone(x: f64) -> f64
{
fn qone(x: f64) -> f64 {
let p: &[f64; 6];
let q: &[f64; 6];
let s: f64;
@ -376,12 +392,23 @@ fn qone(x: f64) -> f64
ix = get_high_word(x);
ix &= 0x7fffffff;
if ix >= 0x40200000 {p = &QR8; q = &QS8;}
else if ix >= 0x40122E8B {p = &QR5; q = &QS5;}
else if ix >= 0x4006DB6D {p = &QR3; q = &QS3;}
else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (0.375 + r/s)/x;
if ix >= 0x40200000 {
p = &QR8;
q = &QS8;
} else if ix >= 0x40122E8B {
p = &QR5;
q = &QS5;
} else if ix >= 0x4006DB6D {
p = &QR3;
q = &QS3;
} else
/*ix >= 0x40000000*/
{
p = &QR2;
q = &QS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
return (0.375 + r / s) / x;
}

369
src/math/j1f.rs
View File

@ -16,10 +16,9 @@
use super::{cosf, fabsf, logf, sinf, sqrtf};
const INVSQRTPI: f32 = 5.6418961287e-01; /* 0x3f106ebb */
const TPI: f32 = 6.3661974669e-01; /* 0x3f22f983 */
const TPI: f32 = 6.3661974669e-01; /* 0x3f22f983 */
fn common(ix: u32, x: f32, y1: bool, sign: bool) -> f32
{
fn common(ix: u32, x: f32, y1: bool, sign: bool) -> f32 {
let z: f64;
let mut s: f64;
let c: f64;
@ -31,41 +30,40 @@ fn common(ix: u32, x: f32, y1: bool, sign: bool) -> f32
s = -s;
}
c = cosf(x) as f64;
cc = s-c;
cc = s - c;
if ix < 0x7f000000 {
ss = -s-c;
z = cosf(2.0*x) as f64;
if s*c > 0.0 {
cc = z/ss;
ss = -s - c;
z = cosf(2.0 * x) as f64;
if s * c > 0.0 {
cc = z / ss;
} else {
ss = z/cc;
ss = z / cc;
}
if ix < 0x58800000 {
if y1 {
ss = -ss;
}
cc = (ponef(x) as f64)*cc-(qonef(x) as f64)*ss;
cc = (ponef(x) as f64) * cc - (qonef(x) as f64) * ss;
}
}
if sign {
cc = -cc;
}
return INVSQRTPI*(cc as f32)/sqrtf(x);
return INVSQRTPI * (cc as f32) / sqrtf(x);
}
/* R0/S0 on [0,2] */
const R00: f32 = -6.2500000000e-02; /* 0xbd800000 */
const R01: f32 = 1.4070566976e-03; /* 0x3ab86cfd */
const R01: f32 = 1.4070566976e-03; /* 0x3ab86cfd */
const R02: f32 = -1.5995563444e-05; /* 0xb7862e36 */
const R03: f32 = 4.9672799207e-08; /* 0x335557d2 */
const S01: f32 = 1.9153760746e-02; /* 0x3c9ce859 */
const S02: f32 = 1.8594678841e-04; /* 0x3942fab6 */
const S03: f32 = 1.1771846857e-06; /* 0x359dffc2 */
const S04: f32 = 5.0463624390e-09; /* 0x31ad6446 */
const S05: f32 = 1.2354227016e-11; /* 0x2d59567e */
const R03: f32 = 4.9672799207e-08; /* 0x335557d2 */
const S01: f32 = 1.9153760746e-02; /* 0x3c9ce859 */
const S02: f32 = 1.8594678841e-04; /* 0x3942fab6 */
const S03: f32 = 1.1771846857e-06; /* 0x359dffc2 */
const S04: f32 = 5.0463624390e-09; /* 0x31ad6446 */
const S05: f32 = 1.2354227016e-11; /* 0x2d59567e */
pub fn j1f(x: f32) -> f32
{
pub fn j1f(x: f32) -> f32 {
let mut z: f32;
let r: f32;
let s: f32;
@ -73,42 +71,43 @@ pub fn j1f(x: f32) -> f32
let sign: bool;
ix = x.to_bits();
sign = (ix>>31) != 0;
sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
if ix >= 0x7f800000 {
return 1.0/(x*x);
return 1.0 / (x * x);
}
if ix >= 0x40000000 { /* |x| >= 2 */
if ix >= 0x40000000 {
/* |x| >= 2 */
return common(ix, fabsf(x), false, sign);
}
if ix >= 0x39000000 { /* |x| >= 2**-13 */
z = x*x;
r = z*(R00+z*(R01+z*(R02+z*R03)));
s = 1.0+z*(S01+z*(S02+z*(S03+z*(S04+z*S05))));
z = 0.5 + r/s;
if ix >= 0x39000000 {
/* |x| >= 2**-13 */
z = x * x;
r = z * (R00 + z * (R01 + z * (R02 + z * R03)));
s = 1.0 + z * (S01 + z * (S02 + z * (S03 + z * (S04 + z * S05))));
z = 0.5 + r / s;
} else {
z = 0.5;
}
return z*x;
return z * x;
}
const U0: [f32; 5] = [
-1.9605709612e-01, /* 0xbe48c331 */
5.0443872809e-02, /* 0x3d4e9e3c */
-1.9125689287e-03, /* 0xbafaaf2a */
2.3525259166e-05, /* 0x37c5581c */
-9.1909917899e-08, /* 0xb3c56003 */
-1.9605709612e-01, /* 0xbe48c331 */
5.0443872809e-02, /* 0x3d4e9e3c */
-1.9125689287e-03, /* 0xbafaaf2a */
2.3525259166e-05, /* 0x37c5581c */
-9.1909917899e-08, /* 0xb3c56003 */
];
const V0: [f32; 5] = [
1.9916731864e-02, /* 0x3ca3286a */
2.0255257550e-04, /* 0x3954644b */
1.3560879779e-06, /* 0x35b602d4 */
6.2274145840e-09, /* 0x31d5f8eb */
1.6655924903e-11, /* 0x2d9281cf */
1.9916731864e-02, /* 0x3ca3286a */
2.0255257550e-04, /* 0x3954644b */
1.3560879779e-06, /* 0x35b602d4 */
6.2274145840e-09, /* 0x31d5f8eb */
1.6655924903e-11, /* 0x2d9281cf */
];
pub fn y1f(x: f32) -> f32
{
pub fn y1f(x: f32) -> f32 {
let z: f32;
let u: f32;
let v: f32;
@ -116,24 +115,26 @@ pub fn y1f(x: f32) -> f32
ix = x.to_bits();
if (ix & 0x7fffffff) == 0 {
return -1.0/0.0;
return -1.0 / 0.0;
}
if (ix>>31) != 0{
return 0.0/0.0;
if (ix >> 31) != 0 {
return 0.0 / 0.0;
}
if ix >= 0x7f800000 {
return 1.0/x;
return 1.0 / x;
}
if ix >= 0x40000000 { /* |x| >= 2.0 */
return common(ix,x,true,false);
if ix >= 0x40000000 {
/* |x| >= 2.0 */
return common(ix, x, true, false);
}
if ix < 0x33000000 { /* x < 2**-25 */
return -TPI/x;
if ix < 0x33000000 {
/* x < 2**-25 */
return -TPI / x;
}
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + TPI*(j1f(x)*logf(x)-1.0/x);
z = x * x;
u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * U0[4])));
v = 1.0 + z * (V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * V0[4]))));
return x * (u / v) + TPI * (j1f(x) * logf(x) - 1.0 / x);
}
/* For x >= 8, the asymptotic expansions of pone is
@ -146,72 +147,74 @@ pub fn y1f(x: f32) -> f32
* | pone(x)-1-R/S | <= 2 ** ( -60.06)
*/
const PR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
1.1718750000e-01, /* 0x3df00000 */
1.3239480972e+01, /* 0x4153d4ea */
4.1205184937e+02, /* 0x43ce06a3 */
3.8747453613e+03, /* 0x45722bed */
7.9144794922e+03, /* 0x45f753d6 */
const PR8: [f32; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
1.1718750000e-01, /* 0x3df00000 */
1.3239480972e+01, /* 0x4153d4ea */
4.1205184937e+02, /* 0x43ce06a3 */
3.8747453613e+03, /* 0x45722bed */
7.9144794922e+03, /* 0x45f753d6 */
];
const PS8: [f32; 5] = [
1.1420736694e+02, /* 0x42e46a2c */
3.6509309082e+03, /* 0x45642ee5 */
3.6956207031e+04, /* 0x47105c35 */
9.7602796875e+04, /* 0x47bea166 */
3.0804271484e+04, /* 0x46f0a88b */
1.1420736694e+02, /* 0x42e46a2c */
3.6509309082e+03, /* 0x45642ee5 */
3.6956207031e+04, /* 0x47105c35 */
9.7602796875e+04, /* 0x47bea166 */
3.0804271484e+04, /* 0x46f0a88b */
];
const PR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
1.3199052094e-11, /* 0x2d68333f */
1.1718749255e-01, /* 0x3defffff */
6.8027510643e+00, /* 0x40d9b023 */
1.0830818176e+02, /* 0x42d89dca */
5.1763616943e+02, /* 0x440168b7 */
5.2871520996e+02, /* 0x44042dc6 */
const PR5: [f32; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
1.3199052094e-11, /* 0x2d68333f */
1.1718749255e-01, /* 0x3defffff */
6.8027510643e+00, /* 0x40d9b023 */
1.0830818176e+02, /* 0x42d89dca */
5.1763616943e+02, /* 0x440168b7 */
5.2871520996e+02, /* 0x44042dc6 */
];
const PS5: [f32; 5] = [
5.9280597687e+01, /* 0x426d1f55 */
9.9140142822e+02, /* 0x4477d9b1 */
5.3532670898e+03, /* 0x45a74a23 */
7.8446904297e+03, /* 0x45f52586 */
1.5040468750e+03, /* 0x44bc0180 */
5.9280597687e+01, /* 0x426d1f55 */
9.9140142822e+02, /* 0x4477d9b1 */
5.3532670898e+03, /* 0x45a74a23 */
7.8446904297e+03, /* 0x45f52586 */
1.5040468750e+03, /* 0x44bc0180 */
];
const PR3: [f32; 6] = [
3.0250391081e-09, /* 0x314fe10d */
1.1718686670e-01, /* 0x3defffab */
3.9329774380e+00, /* 0x407bb5e7 */
3.5119403839e+01, /* 0x420c7a45 */
9.1055007935e+01, /* 0x42b61c2a */
4.8559066772e+01, /* 0x42423c7c */
3.0250391081e-09, /* 0x314fe10d */
1.1718686670e-01, /* 0x3defffab */
3.9329774380e+00, /* 0x407bb5e7 */
3.5119403839e+01, /* 0x420c7a45 */
9.1055007935e+01, /* 0x42b61c2a */
4.8559066772e+01, /* 0x42423c7c */
];
const PS3: [f32; 5] = [
3.4791309357e+01, /* 0x420b2a4d */
3.3676245117e+02, /* 0x43a86198 */
1.0468714600e+03, /* 0x4482dbe3 */
8.9081134033e+02, /* 0x445eb3ed */
1.0378793335e+02, /* 0x42cf936c */
3.4791309357e+01, /* 0x420b2a4d */
3.3676245117e+02, /* 0x43a86198 */
1.0468714600e+03, /* 0x4482dbe3 */
8.9081134033e+02, /* 0x445eb3ed */
1.0378793335e+02, /* 0x42cf936c */
];
const PR2: [f32; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.0771083225e-07, /* 0x33e74ea8 */
1.1717621982e-01, /* 0x3deffa16 */
2.3685150146e+00, /* 0x401795c0 */
1.2242610931e+01, /* 0x4143e1bc */
1.7693971634e+01, /* 0x418d8d41 */
5.0735230446e+00, /* 0x40a25a4d */
const PR2: [f32; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
1.0771083225e-07, /* 0x33e74ea8 */
1.1717621982e-01, /* 0x3deffa16 */
2.3685150146e+00, /* 0x401795c0 */
1.2242610931e+01, /* 0x4143e1bc */
1.7693971634e+01, /* 0x418d8d41 */
5.0735230446e+00, /* 0x40a25a4d */
];
const PS2: [f32; 5] = [
2.1436485291e+01, /* 0x41ab7dec */
1.2529022980e+02, /* 0x42fa9499 */
2.3227647400e+02, /* 0x436846c7 */
1.1767937469e+02, /* 0x42eb5bd7 */
8.3646392822e+00, /* 0x4105d590 */
2.1436485291e+01, /* 0x41ab7dec */
1.2529022980e+02, /* 0x42fa9499 */
2.3227647400e+02, /* 0x436846c7 */
1.1767937469e+02, /* 0x42eb5bd7 */
8.3646392822e+00, /* 0x4105d590 */
];
fn ponef(x: f32) -> f32
{
fn ponef(x: f32) -> f32 {
let p: &[f32; 6];
let q: &[f32; 5];
let z: f32;
@ -221,14 +224,25 @@ fn ponef(x: f32) -> f32
ix = x.to_bits();
ix &= 0x7fffffff;
if ix >= 0x41000000 {p = &PR8; q = &PS8;}
else if ix >= 0x409173eb {p = &PR5; q = &PS5;}
else if ix >= 0x4036d917 {p = &PR3; q = &PS3;}
else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0 + r/s;
if ix >= 0x41000000 {
p = &PR8;
q = &PS8;
} else if ix >= 0x409173eb {
p = &PR5;
q = &PS5;
} else if ix >= 0x4036d917 {
p = &PR3;
q = &PS3;
} else
/*ix >= 0x40000000*/
{
p = &PR2;
q = &PS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4]))));
return 1.0 + r / s;
}
/* For x >= 8, the asymptotic expansions of qone is
@ -241,76 +255,78 @@ fn ponef(x: f32) -> f32
* | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13)
*/
const QR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
-1.0253906250e-01, /* 0xbdd20000 */
-1.6271753311e+01, /* 0xc1822c8d */
-7.5960174561e+02, /* 0xc43de683 */
-1.1849806641e+04, /* 0xc639273a */
-4.8438511719e+04, /* 0xc73d3683 */
const QR8: [f32; 6] = [
/* for x in [inf, 8]=1/[0,0.125] */
0.0000000000e+00, /* 0x00000000 */
-1.0253906250e-01, /* 0xbdd20000 */
-1.6271753311e+01, /* 0xc1822c8d */
-7.5960174561e+02, /* 0xc43de683 */
-1.1849806641e+04, /* 0xc639273a */
-4.8438511719e+04, /* 0xc73d3683 */
];
const QS8: [f32; 6] = [
1.6139537048e+02, /* 0x43216537 */
7.8253862305e+03, /* 0x45f48b17 */
1.3387534375e+05, /* 0x4802bcd6 */
7.1965775000e+05, /* 0x492fb29c */
6.6660125000e+05, /* 0x4922be94 */
-2.9449025000e+05, /* 0xc88fcb48 */
1.6139537048e+02, /* 0x43216537 */
7.8253862305e+03, /* 0x45f48b17 */
1.3387534375e+05, /* 0x4802bcd6 */
7.1965775000e+05, /* 0x492fb29c */
6.6660125000e+05, /* 0x4922be94 */
-2.9449025000e+05, /* 0xc88fcb48 */
];
const QR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.0897993405e-11, /* 0xadb7d219 */
-1.0253904760e-01, /* 0xbdd1fffe */
-8.0564479828e+00, /* 0xc100e736 */
-1.8366960144e+02, /* 0xc337ab6b */
-1.3731937256e+03, /* 0xc4aba633 */
-2.6124443359e+03, /* 0xc523471c */
const QR5: [f32; 6] = [
/* for x in [8,4.5454]=1/[0.125,0.22001] */
-2.0897993405e-11, /* 0xadb7d219 */
-1.0253904760e-01, /* 0xbdd1fffe */
-8.0564479828e+00, /* 0xc100e736 */
-1.8366960144e+02, /* 0xc337ab6b */
-1.3731937256e+03, /* 0xc4aba633 */
-2.6124443359e+03, /* 0xc523471c */
];
const QS5: [f32; 6] = [
8.1276550293e+01, /* 0x42a28d98 */
1.9917987061e+03, /* 0x44f8f98f */
1.7468484375e+04, /* 0x468878f8 */
4.9851425781e+04, /* 0x4742bb6d */
2.7948074219e+04, /* 0x46da5826 */
-4.7191835938e+03, /* 0xc5937978 */
8.1276550293e+01, /* 0x42a28d98 */
1.9917987061e+03, /* 0x44f8f98f */
1.7468484375e+04, /* 0x468878f8 */
4.9851425781e+04, /* 0x4742bb6d */
2.7948074219e+04, /* 0x46da5826 */
-4.7191835938e+03, /* 0xc5937978 */
];
const QR3: [f32; 6] = [
-5.0783124372e-09, /* 0xb1ae7d4f */
-1.0253783315e-01, /* 0xbdd1ff5b */
-4.6101160049e+00, /* 0xc0938612 */
-5.7847221375e+01, /* 0xc267638e */
-2.2824453735e+02, /* 0xc3643e9a */
-2.1921012878e+02, /* 0xc35b35cb */
-5.0783124372e-09, /* 0xb1ae7d4f */
-1.0253783315e-01, /* 0xbdd1ff5b */
-4.6101160049e+00, /* 0xc0938612 */
-5.7847221375e+01, /* 0xc267638e */
-2.2824453735e+02, /* 0xc3643e9a */
-2.1921012878e+02, /* 0xc35b35cb */
];
const QS3: [f32; 6] = [
4.7665153503e+01, /* 0x423ea91e */
6.7386511230e+02, /* 0x4428775e */
3.3801528320e+03, /* 0x45534272 */
5.5477290039e+03, /* 0x45ad5dd5 */
1.9031191406e+03, /* 0x44ede3d0 */
-1.3520118713e+02, /* 0xc3073381 */
4.7665153503e+01, /* 0x423ea91e */
6.7386511230e+02, /* 0x4428775e */
3.3801528320e+03, /* 0x45534272 */
5.5477290039e+03, /* 0x45ad5dd5 */
1.9031191406e+03, /* 0x44ede3d0 */
-1.3520118713e+02, /* 0xc3073381 */
];
const QR2: [f32; 6] = [ /* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.7838172539e-07, /* 0xb43f8932 */
-1.0251704603e-01, /* 0xbdd1f475 */
-2.7522056103e+00, /* 0xc0302423 */
-1.9663616180e+01, /* 0xc19d4f16 */
-4.2325313568e+01, /* 0xc2294d1f */
-2.1371921539e+01, /* 0xc1aaf9b2 */
const QR2: [f32; 6] = [
/* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.7838172539e-07, /* 0xb43f8932 */
-1.0251704603e-01, /* 0xbdd1f475 */
-2.7522056103e+00, /* 0xc0302423 */
-1.9663616180e+01, /* 0xc19d4f16 */
-4.2325313568e+01, /* 0xc2294d1f */
-2.1371921539e+01, /* 0xc1aaf9b2 */
];
const QS2: [f32; 6] = [
2.9533363342e+01, /* 0x41ec4454 */
2.5298155212e+02, /* 0x437cfb47 */
7.5750280762e+02, /* 0x443d602e */
7.3939318848e+02, /* 0x4438d92a */
1.5594900513e+02, /* 0x431bf2f2 */
-4.9594988823e+00, /* 0xc09eb437 */
2.9533363342e+01, /* 0x41ec4454 */
2.5298155212e+02, /* 0x437cfb47 */
7.5750280762e+02, /* 0x443d602e */
7.3939318848e+02, /* 0x4438d92a */
1.5594900513e+02, /* 0x431bf2f2 */
-4.9594988823e+00, /* 0xc09eb437 */
];
fn qonef(x: f32) -> f32
{
fn qonef(x: f32) -> f32 {
let p: &[f32; 6];
let q: &[f32; 6];
let s: f32;
@ -320,12 +336,23 @@ fn qonef(x: f32) -> f32
ix = x.to_bits();
ix &= 0x7fffffff;
if ix >= 0x41000000 {p = &QR8; q = &QS8;}
else if ix >= 0x409173eb {p = &QR5; q = &QS5;}
else if ix >= 0x4036d917 {p = &QR3; q = &QS3;}
else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (0.375 + r/s)/x;
if ix >= 0x41000000 {
p = &QR8;
q = &QS8;
} else if ix >= 0x409173eb {
p = &QR5;
q = &QS5;
} else if ix >= 0x4036d917 {
p = &QR3;
q = &QS3;
} else
/*ix >= 0x40000000*/
{
p = &QR2;
q = &QS2;
}
z = 1.0 / (x * x);
r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5]))));
s = 1.0 + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5])))));
return (0.375 + r / s) / x;
}

119
src/math/jn.rs
View File

@ -34,16 +34,15 @@
* values of n>1.
*/
use super::{cos, fabs, get_high_word, get_low_word, log, j0, j1, sin, sqrt, y0, y1};
use super::{cos, fabs, get_high_word, get_low_word, j0, j1, log, sin, sqrt, y0, y1};
const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
pub fn jn(n: isize, mut x: f64) -> f64
{
pub fn jn(n: i32, mut x: f64) -> f64 {
let mut ix: u32;
let lx: u32;
let nm1: isize;
let mut i: isize;
let nm1: i32;
let mut i: i32;
let mut sign: bool;
let mut a: f64;
let mut b: f64;
@ -51,11 +50,12 @@ pub fn jn(n: isize, mut x: f64) -> f64
ix = get_high_word(x);
lx = get_low_word(x);
sign = (ix>>31) != 0;
sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
// -lx == !lx + 1
if (ix | (lx|(!lx+1))>>31) > 0x7ff00000 { /* nan */
if (ix | (lx | (!lx + 1)) >> 31) > 0x7ff00000 {
/* nan */
return x;
}
@ -67,23 +67,25 @@ pub fn jn(n: isize, mut x: f64) -> f64
return j0(x);
}
if n < 0 {
nm1 = -(n+1);
nm1 = -(n + 1);
x = -x;
sign = !sign;
} else {
nm1 = n-1;
nm1 = n - 1;
}
if nm1 == 0 {
return j1(x);
}
sign &= (n & 1) != 0; /* even n: 0, odd n: signbit(x) */
sign &= (n & 1) != 0; /* even n: 0, odd n: signbit(x) */
x = fabs(x);
if (ix|lx) == 0 || ix == 0x7ff00000 { /* if x is 0 or inf */
if (ix | lx) == 0 || ix == 0x7ff00000 {
/* if x is 0 or inf */
b = 0.0;
} else if (nm1 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 */
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)
@ -97,13 +99,13 @@ pub fn jn(n: isize, mut x: f64) -> f64
* 2 -s+c -c-s
* 3 s+c c-s
*/
temp = match nm1&3 {
0 => -cos(x)+sin(x),
1 => -cos(x)-sin(x),
2 => cos(x)-sin(x),
3 | _ => cos(x)+sin(x),
temp = match nm1 & 3 {
0 => -cos(x) + sin(x),
1 => -cos(x) - sin(x),
2 => cos(x) - sin(x),
3 | _ => cos(x) + sin(x),
};
b = INVSQRTPI*temp/sqrt(x);
b = INVSQRTPI * temp / sqrt(x);
} else {
a = j0(x);
b = j1(x);
@ -111,28 +113,30 @@ pub fn jn(n: isize, mut x: f64) -> f64
while i < nm1 {
i += 1;
temp = b;
b = b*(2.0*(i as f64)/x) - a; /* avoid underflow */
b = b * (2.0 * (i as f64) / x) - a; /* avoid underflow */
a = temp;
}
}
} else {
if ix < 0x3e100000 { /* x < 2**-29 */
if ix < 0x3e100000 {
/* x < 2**-29 */
/* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if nm1 > 32 { /* underflow */
if nm1 > 32 {
/* underflow */
b = 0.0;
} else {
temp = x*0.5;
temp = x * 0.5;
b = temp;
a = 1.0;
i = 2;
while i <= nm1 + 1 {
a *= i as f64; /* a = n! */
b *= temp; /* b = (x/2)^n */
b *= temp; /* b = (x/2)^n */
i += 1;
}
b = b/a;
b = b / a;
}
} else {
/* use backward recurrence */
@ -173,26 +177,26 @@ pub fn jn(n: isize, mut x: f64) -> f64
let mut tmp: f64;
let nf: f64;
let mut k: isize;
let mut k: i32;
nf = (nm1 as f64) + 1.0;
w = 2.0*nf/x;
h = 2.0/x;
z = w+h;
w = 2.0 * nf / x;
h = 2.0 / x;
z = w + h;
q0 = w;
q1 = w*z - 1.0;
q1 = w * z - 1.0;
k = 1;
while q1 < 1.0e9 {
k += 1;
z += h;
tmp = z*q1 - q0;
tmp = z * q1 - q0;
q0 = q1;
q1 = tmp;
}
t = 0.0;
i = k;
while i >= 0 {
t = 1.0/(2.0*((i as f64)+nf)/x - t);
t = 1.0 / (2.0 * ((i as f64) + nf) / x - t);
i -= 1;
}
a = t;
@ -205,12 +209,12 @@ pub fn jn(n: isize, mut x: f64) -> f64
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = nf*log(fabs(w));
tmp = nf * log(fabs(w));
if tmp < 7.09782712893383973096e+02 {
i = nm1;
while i > 0 {
temp = b;
b = b*(2.0*(i as f64))/x - a;
b = b * (2.0 * (i as f64)) / x - a;
a = temp;
i -= 1;
}
@ -218,14 +222,14 @@ pub fn jn(n: isize, mut x: f64) -> f64
i = nm1;
while i > 0 {
temp = b;
b = b*(2.0*(i as f64))/x - a;
b = b * (2.0 * (i as f64)) / x - a;
a = temp;
/* scale b to avoid spurious overflow */
let x1p500 = f64::from_bits(0x5f30000000000000); // 0x1p500 == 2^500
if b > x1p500 {
a /= b;
t /= b;
b = 1.0;
b = 1.0;
}
i -= 1;
}
@ -233,9 +237,9 @@ pub fn jn(n: isize, mut x: f64) -> f64
z = j0(x);
w = j1(x);
if fabs(z) >= fabs(w) {
b = t*z/b;
b = t * z / b;
} else {
b = t*w/a;
b = t * w / a;
}
}
}
@ -247,30 +251,30 @@ pub fn jn(n: isize, mut x: f64) -> f64
}
}
pub fn yn(n: isize, x: f64) -> f64
{
pub fn yn(n: i32, x: f64) -> f64 {
let mut ix: u32;
let lx: u32;
let mut ib: u32;
let nm1: isize;
let nm1: i32;
let mut sign: bool;
let mut i: isize;
let mut i: i32;
let mut a: f64;
let mut b: f64;
let mut temp: f64;
ix = get_high_word(x);
lx = get_low_word(x);
sign = (ix>>31) != 0;
sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
// -lx == !lx + 1
if (ix | (lx|(!lx+1))>>31) > 0x7ff00000 { /* nan */
if (ix | (lx | (!lx + 1)) >> 31) > 0x7ff00000 {
/* nan */
return x;
}
if sign && (ix|lx) != 0 { /* x < 0 */
return 0.0/0.0;
if sign && (ix | lx) != 0 {
/* x < 0 */
return 0.0 / 0.0;
}
if ix == 0x7ff00000 {
return 0.0;
@ -280,10 +284,10 @@ pub fn yn(n: isize, x: f64) -> f64
return y0(x);
}
if n < 0 {
nm1 = -(n+1);
sign = (n&1) != 0;
nm1 = -(n + 1);
sign = (n & 1) != 0;
} else {
nm1 = n-1;
nm1 = n - 1;
sign = false;
}
if nm1 == 0 {
@ -294,7 +298,8 @@ pub fn yn(n: isize, x: f64) -> f64
}
}
if ix >= 0x52d00000 { /* x > 2**302 */
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)
@ -308,13 +313,13 @@ pub fn yn(n: isize, x: f64) -> f64
* 2 -s+c -c-s
* 3 s+c c-s
*/
temp = match nm1&3 {
0 => -sin(x)-cos(x),
1 => -sin(x)+cos(x),
2 => sin(x)+cos(x),
3 | _ => sin(x)-cos(x),
temp = match nm1 & 3 {
0 => -sin(x) - cos(x),
1 => -sin(x) + cos(x),
2 => sin(x) + cos(x),
3 | _ => sin(x) - cos(x),
};
b = INVSQRTPI*temp/sqrt(x);
b = INVSQRTPI * temp / sqrt(x);
} else {
a = y0(x);
b = y1(x);
@ -324,7 +329,7 @@ pub fn yn(n: isize, x: f64) -> f64
while i < nm1 && ib != 0xfff00000 {
i += 1;
temp = b;
b = (2.0*(i as f64)/x)*b - a;
b = (2.0 * (i as f64) / x) * b - a;
ib = get_high_word(b);
a = temp;
}

86
src/math/jnf.rs
View File

@ -15,20 +15,20 @@
use super::{fabsf, j0f, j1f, logf, y0f, y1f};
pub fn jnf(n: isize, mut x: f32) -> f32
{
pub fn jnf(n: i32, mut x: f32) -> f32 {
let mut ix: u32;
let mut nm1: isize;
let mut nm1: i32;
let mut sign: bool;
let mut i: isize;
let mut i: i32;
let mut a: f32;
let mut b: f32;
let mut temp: f32;
ix = x.to_bits();
sign = (ix>>31) != 0;
sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
if ix > 0x7f800000 { /* nan */
if ix > 0x7f800000 {
/* nan */
return x;
}
@ -37,19 +37,20 @@ pub fn jnf(n: isize, mut x: f32) -> f32
return j0f(x);
}
if n < 0 {
nm1 = -(n+1);
nm1 = -(n + 1);
x = -x;
sign = !sign;
} else {
nm1 = n-1;
nm1 = n - 1;
}
if nm1 == 0 {
return j1f(x);
}
sign &= (n&1) != 0; /* even n: 0, odd n: signbit(x) */
sign &= (n & 1) != 0; /* even n: 0, odd n: signbit(x) */
x = fabsf(x);
if ix == 0 || ix == 0x7f800000 { /* if x is 0 or inf */
if ix == 0 || ix == 0x7f800000 {
/* if x is 0 or inf */
b = 0.0;
} else if (nm1 as f32) < x {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
@ -59,15 +60,17 @@ pub fn jnf(n: isize, mut x: f32) -> f32
while i < nm1 {
i += 1;
temp = b;
b = b*(2.0*(i as f32)/x) - a;
b = b * (2.0 * (i as f32) / x) - a;
a = temp;
}
} else {
if ix < 0x35800000 { /* x < 2**-20 */
if ix < 0x35800000 {
/* x < 2**-20 */
/* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if nm1 > 8 { /* underflow */
if nm1 > 8 {
/* underflow */
nm1 = 8;
}
temp = 0.5 * x;
@ -75,11 +78,11 @@ pub fn jnf(n: isize, mut x: f32) -> f32
a = 1.0;
i = 2;
while i <= nm1 + 1 {
a *= i as f32; /* a = n! */
b *= temp; /* b = (x/2)^n */
a *= i as f32; /* a = n! */
b *= temp; /* b = (x/2)^n */
i += 1;
}
b = b/a;
b = b / a;
} else {
/* use backward recurrence */
/* x x^2 x^2
@ -118,26 +121,26 @@ pub fn jnf(n: isize, mut x: f32) -> f32
let mut z: f32;
let mut tmp: f32;
let nf: f32;
let mut k: isize;
let mut k: i32;
nf = (nm1 as f32)+1.0;
w = 2.0*(nf as f32)/x;
h = 2.0/x;
z = w+h;
nf = (nm1 as f32) + 1.0;
w = 2.0 * (nf as f32) / x;
h = 2.0 / x;
z = w + h;
q0 = w;
q1 = w*z - 1.0;
q1 = w * z - 1.0;
k = 1;
while q1 < 1.0e4 {
k += 1;
z += h;
tmp = z*q1 - q0;
tmp = z * q1 - q0;
q0 = q1;
q1 = tmp;
}
t = 0.0;
i = k;
while i >= 0 {
t = 1.0/(2.0*((i as f32)+nf)/x-t);
t = 1.0 / (2.0 * ((i as f32) + nf) / x - t);
i -= 1;
}
a = t;
@ -150,12 +153,12 @@ pub fn jnf(n: isize, mut x: f32) -> f32
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = nf*logf(fabsf(w));
tmp = nf * logf(fabsf(w));
if tmp < 88.721679688 {
i = nm1;
while i > 0 {
temp = b;
b = 2.0*(i as f32)*b/x - a;
b = 2.0 * (i as f32) * b / x - a;
a = temp;
i -= 1;
}
@ -163,7 +166,7 @@ pub fn jnf(n: isize, mut x: f32) -> f32
i = nm1;
while i > 0 {
temp = b;
b = 2.0*(i as f32)*b/x - a;
b = 2.0 * (i as f32) * b / x - a;
a = temp;
/* scale b to avoid spurious overflow */
let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
@ -178,9 +181,9 @@ pub fn jnf(n: isize, mut x: f32) -> f32
z = j0f(x);
w = j1f(x);
if fabsf(z) >= fabsf(w) {
b = t*z/b;
b = t * z / b;
} else {
b = t*w/a;
b = t * w / a;
}
}
}
@ -192,25 +195,26 @@ pub fn jnf(n: isize, mut x: f32) -> f32
}
}
pub fn ynf(n: isize, x: f32) -> f32
{
pub fn ynf(n: i32, x: f32) -> f32 {
let mut ix: u32;
let mut ib: u32;
let nm1: isize;
let nm1: i32;
let mut sign: bool;
let mut i: isize;
let mut i: i32;
let mut a: f32;
let mut b: f32;
let mut temp: f32;
ix = x.to_bits();
sign = (ix>>31) != 0;
sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
if ix > 0x7f800000 { /* nan */
if ix > 0x7f800000 {
/* nan */
return x;
}
if sign && ix != 0 { /* x < 0 */
return 0.0/0.0;
if sign && ix != 0 {
/* x < 0 */
return 0.0 / 0.0;
}
if ix == 0x7f800000 {
return 0.0;
@ -220,10 +224,10 @@ pub fn ynf(n: isize, x: f32) -> f32
return y0f(x);
}
if n < 0 {
nm1 = -(n+1);
sign = (n&1) != 0;
nm1 = -(n + 1);
sign = (n & 1) != 0;
} else {
nm1 = n-1;
nm1 = n - 1;
sign = false;
}
if nm1 == 0 {
@ -242,7 +246,7 @@ pub fn ynf(n: isize, x: f32) -> f32
while i < nm1 && ib != 0xff800000 {
i += 1;
temp = b;
b = (2.0*(i as f32)/x)*b - a;
b = (2.0 * (i as f32) / x) * b - a;
ib = b.to_bits();
a = temp;
}

248
src/math/lgamma.rs
View File

@ -80,88 +80,87 @@
use super::{floor, k_cos, k_sin, log};
const PI: f64 = 3.14159265358979311600e+00; /* 0x400921FB, 0x54442D18 */
const A0: f64 = 7.72156649015328655494e-02; /* 0x3FB3C467, 0xE37DB0C8 */
const A1: f64 = 3.22467033424113591611e-01; /* 0x3FD4A34C, 0xC4A60FAD */
const A2: f64 = 6.73523010531292681824e-02; /* 0x3FB13E00, 0x1A5562A7 */
const A3: f64 = 2.05808084325167332806e-02; /* 0x3F951322, 0xAC92547B */
const A4: f64 = 7.38555086081402883957e-03; /* 0x3F7E404F, 0xB68FEFE8 */
const A5: f64 = 2.89051383673415629091e-03; /* 0x3F67ADD8, 0xCCB7926B */
const A6: f64 = 1.19270763183362067845e-03; /* 0x3F538A94, 0x116F3F5D */
const A7: f64 = 5.10069792153511336608e-04; /* 0x3F40B6C6, 0x89B99C00 */
const A8: f64 = 2.20862790713908385557e-04; /* 0x3F2CF2EC, 0xED10E54D */
const A9: f64 = 1.08011567247583939954e-04; /* 0x3F1C5088, 0x987DFB07 */
const A10: f64 = 2.52144565451257326939e-05; /* 0x3EFA7074, 0x428CFA52 */
const A11: f64 = 4.48640949618915160150e-05; /* 0x3F07858E, 0x90A45837 */
const TC: f64 = 1.46163214496836224576e+00; /* 0x3FF762D8, 0x6356BE3F */
const TF: f64 = -1.21486290535849611461e-01; /* 0xBFBF19B9, 0xBCC38A42 */
const PI: f64 = 3.14159265358979311600e+00; /* 0x400921FB, 0x54442D18 */
const A0: f64 = 7.72156649015328655494e-02; /* 0x3FB3C467, 0xE37DB0C8 */
const A1: f64 = 3.22467033424113591611e-01; /* 0x3FD4A34C, 0xC4A60FAD */
const A2: f64 = 6.73523010531292681824e-02; /* 0x3FB13E00, 0x1A5562A7 */
const A3: f64 = 2.05808084325167332806e-02; /* 0x3F951322, 0xAC92547B */
const A4: f64 = 7.38555086081402883957e-03; /* 0x3F7E404F, 0xB68FEFE8 */
const A5: f64 = 2.89051383673415629091e-03; /* 0x3F67ADD8, 0xCCB7926B */
const A6: f64 = 1.19270763183362067845e-03; /* 0x3F538A94, 0x116F3F5D */
const A7: f64 = 5.10069792153511336608e-04; /* 0x3F40B6C6, 0x89B99C00 */
const A8: f64 = 2.20862790713908385557e-04; /* 0x3F2CF2EC, 0xED10E54D */
const A9: f64 = 1.08011567247583939954e-04; /* 0x3F1C5088, 0x987DFB07 */
const A10: f64 = 2.52144565451257326939e-05; /* 0x3EFA7074, 0x428CFA52 */
const A11: f64 = 4.48640949618915160150e-05; /* 0x3F07858E, 0x90A45837 */
const TC: f64 = 1.46163214496836224576e+00; /* 0x3FF762D8, 0x6356BE3F */
const TF: f64 = -1.21486290535849611461e-01; /* 0xBFBF19B9, 0xBCC38A42 */
/* tt = -(tail of TF) */
const TT: f64 = -3.63867699703950536541e-18; /* 0xBC50C7CA, 0xA48A971F */
const T0: f64 = 4.83836122723810047042e-01; /* 0x3FDEF72B, 0xC8EE38A2 */
const T1: f64 = -1.47587722994593911752e-01; /* 0xBFC2E427, 0x8DC6C509 */
const T2: f64 = 6.46249402391333854778e-02; /* 0x3FB08B42, 0x94D5419B */
const T3: f64 = -3.27885410759859649565e-02; /* 0xBFA0C9A8, 0xDF35B713 */
const T4: f64 = 1.79706750811820387126e-02; /* 0x3F9266E7, 0x970AF9EC */
const T5: f64 = -1.03142241298341437450e-02; /* 0xBF851F9F, 0xBA91EC6A */
const T6: f64 = 6.10053870246291332635e-03; /* 0x3F78FCE0, 0xE370E344 */
const T7: f64 = -3.68452016781138256760e-03; /* 0xBF6E2EFF, 0xB3E914D7 */
const T8: f64 = 2.25964780900612472250e-03; /* 0x3F6282D3, 0x2E15C915 */
const T9: f64 = -1.40346469989232843813e-03; /* 0xBF56FE8E, 0xBF2D1AF1 */
const T10: f64 = 8.81081882437654011382e-04; /* 0x3F4CDF0C, 0xEF61A8E9 */
const TT: f64 = -3.63867699703950536541e-18; /* 0xBC50C7CA, 0xA48A971F */
const T0: f64 = 4.83836122723810047042e-01; /* 0x3FDEF72B, 0xC8EE38A2 */
const T1: f64 = -1.47587722994593911752e-01; /* 0xBFC2E427, 0x8DC6C509 */
const T2: f64 = 6.46249402391333854778e-02; /* 0x3FB08B42, 0x94D5419B */
const T3: f64 = -3.27885410759859649565e-02; /* 0xBFA0C9A8, 0xDF35B713 */
const T4: f64 = 1.79706750811820387126e-02; /* 0x3F9266E7, 0x970AF9EC */
const T5: f64 = -1.03142241298341437450e-02; /* 0xBF851F9F, 0xBA91EC6A */
const T6: f64 = 6.10053870246291332635e-03; /* 0x3F78FCE0, 0xE370E344 */
const T7: f64 = -3.68452016781138256760e-03; /* 0xBF6E2EFF, 0xB3E914D7 */
const T8: f64 = 2.25964780900612472250e-03; /* 0x3F6282D3, 0x2E15C915 */
const T9: f64 = -1.40346469989232843813e-03; /* 0xBF56FE8E, 0xBF2D1AF1 */
const T10: f64 = 8.81081882437654011382e-04; /* 0x3F4CDF0C, 0xEF61A8E9 */
const T11: f64 = -5.38595305356740546715e-04; /* 0xBF41A610, 0x9C73E0EC */
const T12: f64 = 3.15632070903625950361e-04; /* 0x3F34AF6D, 0x6C0EBBF7 */
const T12: f64 = 3.15632070903625950361e-04; /* 0x3F34AF6D, 0x6C0EBBF7 */
const T13: f64 = -3.12754168375120860518e-04; /* 0xBF347F24, 0xECC38C38 */
const T14: f64 = 3.35529192635519073543e-04; /* 0x3F35FD3E, 0xE8C2D3F4 */
const U0: f64 = -7.72156649015328655494e-02; /* 0xBFB3C467, 0xE37DB0C8 */
const U1: f64 = 6.32827064025093366517e-01; /* 0x3FE4401E, 0x8B005DFF */
const U2: f64 = 1.45492250137234768737e+00; /* 0x3FF7475C, 0xD119BD6F */
const U3: f64 = 9.77717527963372745603e-01; /* 0x3FEF4976, 0x44EA8450 */
const U4: f64 = 2.28963728064692451092e-01; /* 0x3FCD4EAE, 0xF6010924 */
const U5: f64 = 1.33810918536787660377e-02; /* 0x3F8B678B, 0xBF2BAB09 */
const V1: f64 = 2.45597793713041134822e+00; /* 0x4003A5D7, 0xC2BD619C */
const V2: f64 = 2.12848976379893395361e+00; /* 0x40010725, 0xA42B18F5 */
const V3: f64 = 7.69285150456672783825e-01; /* 0x3FE89DFB, 0xE45050AF */
const V4: f64 = 1.04222645593369134254e-01; /* 0x3FBAAE55, 0xD6537C88 */
const V5: f64 = 3.21709242282423911810e-03; /* 0x3F6A5ABB, 0x57D0CF61 */
const S0: f64 = -7.72156649015328655494e-02; /* 0xBFB3C467, 0xE37DB0C8 */
const S1: f64 = 2.14982415960608852501e-01; /* 0x3FCB848B, 0x36E20878 */
const S2: f64 = 3.25778796408930981787e-01; /* 0x3FD4D98F, 0x4F139F59 */
const S3: f64 = 1.46350472652464452805e-01; /* 0x3FC2BB9C, 0xBEE5F2F7 */
const S4: f64 = 2.66422703033638609560e-02; /* 0x3F9B481C, 0x7E939961 */
const S5: f64 = 1.84028451407337715652e-03; /* 0x3F5E26B6, 0x7368F239 */
const S6: f64 = 3.19475326584100867617e-05; /* 0x3F00BFEC, 0xDD17E945 */
const R1: f64 = 1.39200533467621045958e+00; /* 0x3FF645A7, 0x62C4AB74 */
const R2: f64 = 7.21935547567138069525e-01; /* 0x3FE71A18, 0x93D3DCDC */
const R3: f64 = 1.71933865632803078993e-01; /* 0x3FC601ED, 0xCCFBDF27 */
const R4: f64 = 1.86459191715652901344e-02; /* 0x3F9317EA, 0x742ED475 */
const R5: f64 = 7.77942496381893596434e-04; /* 0x3F497DDA, 0xCA41A95B */
const R6: f64 = 7.32668430744625636189e-06; /* 0x3EDEBAF7, 0xA5B38140 */
const W0: f64 = 4.18938533204672725052e-01; /* 0x3FDACFE3, 0x90C97D69 */
const W1: f64 = 8.33333333333329678849e-02; /* 0x3FB55555, 0x5555553B */
const W2: f64 = -2.77777777728775536470e-03; /* 0xBF66C16C, 0x16B02E5C */
const W3: f64 = 7.93650558643019558500e-04; /* 0x3F4A019F, 0x98CF38B6 */
const W4: f64 = -5.95187557450339963135e-04; /* 0xBF4380CB, 0x8C0FE741 */
const W5: f64 = 8.36339918996282139126e-04; /* 0x3F4B67BA, 0x4CDAD5D1 */
const W6: f64 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
const T14: f64 = 3.35529192635519073543e-04; /* 0x3F35FD3E, 0xE8C2D3F4 */
const U0: f64 = -7.72156649015328655494e-02; /* 0xBFB3C467, 0xE37DB0C8 */
const U1: f64 = 6.32827064025093366517e-01; /* 0x3FE4401E, 0x8B005DFF */
const U2: f64 = 1.45492250137234768737e+00; /* 0x3FF7475C, 0xD119BD6F */
const U3: f64 = 9.77717527963372745603e-01; /* 0x3FEF4976, 0x44EA8450 */
const U4: f64 = 2.28963728064692451092e-01; /* 0x3FCD4EAE, 0xF6010924 */
const U5: f64 = 1.33810918536787660377e-02; /* 0x3F8B678B, 0xBF2BAB09 */
const V1: f64 = 2.45597793713041134822e+00; /* 0x4003A5D7, 0xC2BD619C */
const V2: f64 = 2.12848976379893395361e+00; /* 0x40010725, 0xA42B18F5 */
const V3: f64 = 7.69285150456672783825e-01; /* 0x3FE89DFB, 0xE45050AF */
const V4: f64 = 1.04222645593369134254e-01; /* 0x3FBAAE55, 0xD6537C88 */
const V5: f64 = 3.21709242282423911810e-03; /* 0x3F6A5ABB, 0x57D0CF61 */
const S0: f64 = -7.72156649015328655494e-02; /* 0xBFB3C467, 0xE37DB0C8 */
const S1: f64 = 2.14982415960608852501e-01; /* 0x3FCB848B, 0x36E20878 */
const S2: f64 = 3.25778796408930981787e-01; /* 0x3FD4D98F, 0x4F139F59 */
const S3: f64 = 1.46350472652464452805e-01; /* 0x3FC2BB9C, 0xBEE5F2F7 */
const S4: f64 = 2.66422703033638609560e-02; /* 0x3F9B481C, 0x7E939961 */
const S5: f64 = 1.84028451407337715652e-03; /* 0x3F5E26B6, 0x7368F239 */
const S6: f64 = 3.19475326584100867617e-05; /* 0x3F00BFEC, 0xDD17E945 */
const R1: f64 = 1.39200533467621045958e+00; /* 0x3FF645A7, 0x62C4AB74 */
const R2: f64 = 7.21935547567138069525e-01; /* 0x3FE71A18, 0x93D3DCDC */
const R3: f64 = 1.71933865632803078993e-01; /* 0x3FC601ED, 0xCCFBDF27 */
const R4: f64 = 1.86459191715652901344e-02; /* 0x3F9317EA, 0x742ED475 */
const R5: f64 = 7.77942496381893596434e-04; /* 0x3F497DDA, 0xCA41A95B */
const R6: f64 = 7.32668430744625636189e-06; /* 0x3EDEBAF7, 0xA5B38140 */
const W0: f64 = 4.18938533204672725052e-01; /* 0x3FDACFE3, 0x90C97D69 */
const W1: f64 = 8.33333333333329678849e-02; /* 0x3FB55555, 0x5555553B */
const W2: f64 = -2.77777777728775536470e-03; /* 0xBF66C16C, 0x16B02E5C */
const W3: f64 = 7.93650558643019558500e-04; /* 0x3F4A019F, 0x98CF38B6 */
const W4: f64 = -5.95187557450339963135e-04; /* 0xBF4380CB, 0x8C0FE741 */
const W5: f64 = 8.36339918996282139126e-04; /* 0x3F4B67BA, 0x4CDAD5D1 */
const W6: f64 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
/* sin(PI*x) assuming x > 2^-100, if sin(PI*x)==0 the sign is arbitrary */
fn sin_pi(mut x: f64) -> f64
{
let mut n: isize;
fn sin_pi(mut x: f64) -> f64 {
let mut n: i32;
/* spurious inexact if odd int */
x = 2.0*(x*0.5 - floor(x*0.5)); /* x mod 2.0 */
x = 2.0 * (x * 0.5 - floor(x * 0.5)); /* x mod 2.0 */
n = (x*4.0) as isize;
n = (n+1)/2;
x -= (n as f64)*0.5;
n = (x * 4.0) as i32;
n = (n + 1) / 2;
x -= (n as f64) * 0.5;
x *= PI;
match n {
1 => k_cos(x, 0.0),
2 => k_sin(-x, 0.0, 0),
3 => -k_cos(x, 0.0),
0|_ => k_sin(x, 0.0, 0),
1 => k_cos(x, 0.0),
2 => k_sin(-x, 0.0, 0),
3 => -k_cos(x, 0.0),
0 | _ => k_sin(x, 0.0, 0),
}
}
@ -169,8 +168,7 @@ pub fn lgamma(x: f64) -> f64 {
lgamma_r(x).0
}
pub fn lgamma_r(mut x: f64) -> (f64, isize)
{
pub fn lgamma_r(mut x: f64) -> (f64, i32) {
let u: u64 = x.to_bits();
let mut t: f64;
let y: f64;
@ -185,17 +183,18 @@ pub fn lgamma_r(mut x: f64) -> (f64, isize)
let w: f64;
let ix: u32;
let sign: bool;
let i: isize;
let mut signgam: isize;
let i: i32;
let mut signgam: i32;
/* purge off +-inf, NaN, +-0, tiny and negative arguments */
signgam = 1;
sign = (u>>63) != 0;
ix = ((u>>32) as u32) & 0x7fffffff;
sign = (u >> 63) != 0;
ix = ((u >> 32) as u32) & 0x7fffffff;
if ix >= 0x7ff00000 {
return (x*x, signgam);
return (x * x, signgam);
}
if ix < (0x3ff-70)<<20 { /* |x|<2**-70, return -log(|x|) */
if ix < (0x3ff - 70) << 20 {
/* |x|<2**-70, return -log(|x|) */
if sign {
x = -x;
signgam = -1;
@ -205,15 +204,16 @@ pub fn lgamma_r(mut x: f64) -> (f64, isize)
if sign {
x = -x;
t = sin_pi(x);
if t == 0.0 { /* -integer */
return (1.0/(x-x), signgam);
if t == 0.0 {
/* -integer */
return (1.0 / (x - x), signgam);
}
if t > 0.0 {
signgam = -1;
} else {
t = -t;
}
nadj = log(PI/(t*x));
nadj = log(PI / (t * x));
} else {
nadj = 0.0;
}
@ -224,13 +224,14 @@ pub fn lgamma_r(mut x: f64) -> (f64, isize)
}
/* for x < 2.0 */
else if ix < 0x40000000 {
if ix <= 0x3feccccc { /* lgamma(x) = lgamma(x+1)-log(x) */
if ix <= 0x3feccccc {
/* lgamma(x) = lgamma(x+1)-log(x) */
r = -log(x);
if ix >= 0x3FE76944 {
y = 1.0 - x;
i = 0;
} else if ix >= 0x3FCDA661 {
y = x - (TC-1.0);
y = x - (TC - 1.0);
i = 1;
} else {
y = x;
@ -238,10 +239,12 @@ pub fn lgamma_r(mut x: f64) -> (f64, isize)
}
} else {
r = 0.0;
if ix >= 0x3FFBB4C3 { /* [1.7316,2] */
if ix >= 0x3FFBB4C3 {
/* [1.7316,2] */
y = 2.0 - x;
i = 0;
} else if ix >= 0x3FF3B4C4 { /* [1.23,1.73] */
} else if ix >= 0x3FF3B4C4 {
/* [1.23,1.73] */
y = x - TC;
i = 1;
} else {
@ -251,56 +254,67 @@ pub fn lgamma_r(mut x: f64) -> (f64, isize)
}
match i {
0 => {
z = y*y;
p1 = A0+z*(A2+z*(A4+z*(A6+z*(A8+z*A10))));
p2 = z*(A1+z*(A3+z*(A5+z*(A7+z*(A9+z*A11)))));
p = y*p1+p2;
r += p-0.5*y;
z = y * y;
p1 = A0 + z * (A2 + z * (A4 + z * (A6 + z * (A8 + z * A10))));
p2 = z * (A1 + z * (A3 + z * (A5 + z * (A7 + z * (A9 + z * A11)))));
p = y * p1 + p2;
r += p - 0.5 * y;
}
1 => {
z = y*y;
w = z*y;
p1 = T0+w*(T3+w*(T6+w*(T9 +w*T12))); /* parallel comp */
p2 = T1+w*(T4+w*(T7+w*(T10+w*T13)));
p3 = T2+w*(T5+w*(T8+w*(T11+w*T14)));
p = z*p1-(TT-w*(p2+y*p3));
z = y * y;
w = z * y;
p1 = T0 + w * (T3 + w * (T6 + w * (T9 + w * T12))); /* parallel comp */
p2 = T1 + w * (T4 + w * (T7 + w * (T10 + w * T13)));
p3 = T2 + w * (T5 + w * (T8 + w * (T11 + w * T14)));
p = z * p1 - (TT - w * (p2 + y * p3));
r += TF + p;
}
2 => {
p1 = y*(U0+y*(U1+y*(U2+y*(U3+y*(U4+y*U5)))));
p2 = 1.0+y*(V1+y*(V2+y*(V3+y*(V4+y*V5))));
r += -0.5*y + p1/p2;
p1 = y * (U0 + y * (U1 + y * (U2 + y * (U3 + y * (U4 + y * U5)))));
p2 = 1.0 + y * (V1 + y * (V2 + y * (V3 + y * (V4 + y * V5))));
r += -0.5 * y + p1 / p2;
}
#[cfg(feature = "checked")]
_ => unreachable!(),
#[cfg(not(feature = "checked"))]
_ => {}
}
} else if ix < 0x40200000 { /* x < 8.0 */
i = x as isize;
} else if ix < 0x40200000 {
/* x < 8.0 */
i = x as i32;
y = x - (i as f64);
p = y*(S0+y*(S1+y*(S2+y*(S3+y*(S4+y*(S5+y*S6))))));
q = 1.0+y*(R1+y*(R2+y*(R3+y*(R4+y*(R5+y*R6)))));
r = 0.5*y+p/q;
z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
p = y * (S0 + y * (S1 + y * (S2 + y * (S3 + y * (S4 + y * (S5 + y * S6))))));
q = 1.0 + y * (R1 + y * (R2 + y * (R3 + y * (R4 + y * (R5 + y * R6)))));
r = 0.5 * y + p / q;
z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
// TODO: In C, this was implemented using switch jumps with fallthrough.
// Does this implementation have performance problems?
if i >= 7 { z *= y + 6.0; }
if i >= 6 { z *= y + 5.0; }
if i >= 5 { z *= y + 4.0; }
if i >= 4 { z *= y + 3.0; }
if i >= 7 {
z *= y + 6.0;
}
if i >= 6 {
z *= y + 5.0;
}
if i >= 5 {
z *= y + 4.0;
}
if i >= 4 {
z *= y + 3.0;
}
if i >= 3 {
z *= y + 2.0;
r += log(z);
}
} else if ix < 0x43900000 { /* 8.0 <= x < 2**58 */
} else if ix < 0x43900000 {
/* 8.0 <= x < 2**58 */
t = log(x);
z = 1.0/x;
y = z*z;
w = W0+z*(W1+y*(W2+y*(W3+y*(W4+y*(W5+y*W6)))));
r = (x-0.5)*(t-1.0)+w;
} else { /* 2**58 <= x <= inf */
r = x*(log(x)-1.0);
z = 1.0 / x;
y = z * z;
w = W0 + z * (W1 + y * (W2 + y * (W3 + y * (W4 + y * (W5 + y * W6)))));
r = (x - 0.5) * (t - 1.0) + w;
} else {
/* 2**58 <= x <= inf */
r = x * (log(x) - 1.0);
}
if sign {
r = nadj - r;

238
src/math/lgammaf.rs
View File

@ -15,88 +15,87 @@
use super::{floorf, k_cosf, k_sinf, logf};
const PI: f32 = 3.1415927410e+00; /* 0x40490fdb */
const A0: f32 = 7.7215664089e-02; /* 0x3d9e233f */
const A1: f32 = 3.2246702909e-01; /* 0x3ea51a66 */
const A2: f32 = 6.7352302372e-02; /* 0x3d89f001 */
const A3: f32 = 2.0580807701e-02; /* 0x3ca89915 */
const A4: f32 = 7.3855509982e-03; /* 0x3bf2027e */
const A5: f32 = 2.8905137442e-03; /* 0x3b3d6ec6 */
const A6: f32 = 1.1927076848e-03; /* 0x3a9c54a1 */
const A7: f32 = 5.1006977446e-04; /* 0x3a05b634 */
const A8: f32 = 2.2086278477e-04; /* 0x39679767 */
const A9: f32 = 1.0801156895e-04; /* 0x38e28445 */
const A10: f32 = 2.5214456400e-05; /* 0x37d383a2 */
const A11: f32 = 4.4864096708e-05; /* 0x383c2c75 */
const TC: f32 = 1.4616321325e+00; /* 0x3fbb16c3 */
const TF: f32 = -1.2148628384e-01; /* 0xbdf8cdcd */
const PI: f32 = 3.1415927410e+00; /* 0x40490fdb */
const A0: f32 = 7.7215664089e-02; /* 0x3d9e233f */
const A1: f32 = 3.2246702909e-01; /* 0x3ea51a66 */
const A2: f32 = 6.7352302372e-02; /* 0x3d89f001 */
const A3: f32 = 2.0580807701e-02; /* 0x3ca89915 */
const A4: f32 = 7.3855509982e-03; /* 0x3bf2027e */
const A5: f32 = 2.8905137442e-03; /* 0x3b3d6ec6 */
const A6: f32 = 1.1927076848e-03; /* 0x3a9c54a1 */
const A7: f32 = 5.1006977446e-04; /* 0x3a05b634 */
const A8: f32 = 2.2086278477e-04; /* 0x39679767 */
const A9: f32 = 1.0801156895e-04; /* 0x38e28445 */
const A10: f32 = 2.5214456400e-05; /* 0x37d383a2 */
const A11: f32 = 4.4864096708e-05; /* 0x383c2c75 */
const TC: f32 = 1.4616321325e+00; /* 0x3fbb16c3 */
const TF: f32 = -1.2148628384e-01; /* 0xbdf8cdcd */
/* TT = -(tail of TF) */
const TT: f32 = 6.6971006518e-09; /* 0x31e61c52 */
const T0: f32 = 4.8383611441e-01; /* 0x3ef7b95e */
const T1: f32 = -1.4758771658e-01; /* 0xbe17213c */
const T2: f32 = 6.4624942839e-02; /* 0x3d845a15 */
const T3: f32 = -3.2788541168e-02; /* 0xbd064d47 */
const T4: f32 = 1.7970675603e-02; /* 0x3c93373d */
const T5: f32 = -1.0314224288e-02; /* 0xbc28fcfe */
const T6: f32 = 6.1005386524e-03; /* 0x3bc7e707 */
const T7: f32 = -3.6845202558e-03; /* 0xbb7177fe */
const T8: f32 = 2.2596477065e-03; /* 0x3b141699 */
const T9: f32 = -1.4034647029e-03; /* 0xbab7f476 */
const T10: f32 = 8.8108185446e-04; /* 0x3a66f867 */
const TT: f32 = 6.6971006518e-09; /* 0x31e61c52 */
const T0: f32 = 4.8383611441e-01; /* 0x3ef7b95e */
const T1: f32 = -1.4758771658e-01; /* 0xbe17213c */
const T2: f32 = 6.4624942839e-02; /* 0x3d845a15 */
const T3: f32 = -3.2788541168e-02; /* 0xbd064d47 */
const T4: f32 = 1.7970675603e-02; /* 0x3c93373d */
const T5: f32 = -1.0314224288e-02; /* 0xbc28fcfe */
const T6: f32 = 6.1005386524e-03; /* 0x3bc7e707 */
const T7: f32 = -3.6845202558e-03; /* 0xbb7177fe */
const T8: f32 = 2.2596477065e-03; /* 0x3b141699 */
const T9: f32 = -1.4034647029e-03; /* 0xbab7f476 */
const T10: f32 = 8.8108185446e-04; /* 0x3a66f867 */
const T11: f32 = -5.3859531181e-04; /* 0xba0d3085 */
const T12: f32 = 3.1563205994e-04; /* 0x39a57b6b */
const T12: f32 = 3.1563205994e-04; /* 0x39a57b6b */
const T13: f32 = -3.1275415677e-04; /* 0xb9a3f927 */
const T14: f32 = 3.3552918467e-04; /* 0x39afe9f7 */
const U0: f32 = -7.7215664089e-02; /* 0xbd9e233f */
const U1: f32 = 6.3282704353e-01; /* 0x3f2200f4 */
const U2: f32 = 1.4549225569e+00; /* 0x3fba3ae7 */
const U3: f32 = 9.7771751881e-01; /* 0x3f7a4bb2 */
const U4: f32 = 2.2896373272e-01; /* 0x3e6a7578 */
const U5: f32 = 1.3381091878e-02; /* 0x3c5b3c5e */
const V1: f32 = 2.4559779167e+00; /* 0x401d2ebe */
const V2: f32 = 2.1284897327e+00; /* 0x4008392d */
const V3: f32 = 7.6928514242e-01; /* 0x3f44efdf */
const V4: f32 = 1.0422264785e-01; /* 0x3dd572af */
const V5: f32 = 3.2170924824e-03; /* 0x3b52d5db */
const S0: f32 = -7.7215664089e-02; /* 0xbd9e233f */
const S1: f32 = 2.1498242021e-01; /* 0x3e5c245a */
const S2: f32 = 3.2577878237e-01; /* 0x3ea6cc7a */
const S3: f32 = 1.4635047317e-01; /* 0x3e15dce6 */
const S4: f32 = 2.6642270386e-02; /* 0x3cda40e4 */
const S5: f32 = 1.8402845599e-03; /* 0x3af135b4 */
const S6: f32 = 3.1947532989e-05; /* 0x3805ff67 */
const R1: f32 = 1.3920053244e+00; /* 0x3fb22d3b */
const R2: f32 = 7.2193557024e-01; /* 0x3f38d0c5 */
const R3: f32 = 1.7193385959e-01; /* 0x3e300f6e */
const R4: f32 = 1.8645919859e-02; /* 0x3c98bf54 */
const R5: f32 = 7.7794247773e-04; /* 0x3a4beed6 */
const R6: f32 = 7.3266842264e-06; /* 0x36f5d7bd */
const W0: f32 = 4.1893854737e-01; /* 0x3ed67f1d */
const W1: f32 = 8.3333335817e-02; /* 0x3daaaaab */
const W2: f32 = -2.7777778450e-03; /* 0xbb360b61 */
const W3: f32 = 7.9365057172e-04; /* 0x3a500cfd */
const W4: f32 = -5.9518753551e-04; /* 0xba1c065c */
const W5: f32 = 8.3633989561e-04; /* 0x3a5b3dd2 */
const W6: f32 = -1.6309292987e-03; /* 0xbad5c4e8 */
const T14: f32 = 3.3552918467e-04; /* 0x39afe9f7 */
const U0: f32 = -7.7215664089e-02; /* 0xbd9e233f */
const U1: f32 = 6.3282704353e-01; /* 0x3f2200f4 */
const U2: f32 = 1.4549225569e+00; /* 0x3fba3ae7 */
const U3: f32 = 9.7771751881e-01; /* 0x3f7a4bb2 */
const U4: f32 = 2.2896373272e-01; /* 0x3e6a7578 */
const U5: f32 = 1.3381091878e-02; /* 0x3c5b3c5e */
const V1: f32 = 2.4559779167e+00; /* 0x401d2ebe */
const V2: f32 = 2.1284897327e+00; /* 0x4008392d */
const V3: f32 = 7.6928514242e-01; /* 0x3f44efdf */
const V4: f32 = 1.0422264785e-01; /* 0x3dd572af */
const V5: f32 = 3.2170924824e-03; /* 0x3b52d5db */
const S0: f32 = -7.7215664089e-02; /* 0xbd9e233f */
const S1: f32 = 2.1498242021e-01; /* 0x3e5c245a */
const S2: f32 = 3.2577878237e-01; /* 0x3ea6cc7a */
const S3: f32 = 1.4635047317e-01; /* 0x3e15dce6 */
const S4: f32 = 2.6642270386e-02; /* 0x3cda40e4 */
const S5: f32 = 1.8402845599e-03; /* 0x3af135b4 */
const S6: f32 = 3.1947532989e-05; /* 0x3805ff67 */
const R1: f32 = 1.3920053244e+00; /* 0x3fb22d3b */
const R2: f32 = 7.2193557024e-01; /* 0x3f38d0c5 */
const R3: f32 = 1.7193385959e-01; /* 0x3e300f6e */
const R4: f32 = 1.8645919859e-02; /* 0x3c98bf54 */
const R5: f32 = 7.7794247773e-04; /* 0x3a4beed6 */
const R6: f32 = 7.3266842264e-06; /* 0x36f5d7bd */
const W0: f32 = 4.1893854737e-01; /* 0x3ed67f1d */
const W1: f32 = 8.3333335817e-02; /* 0x3daaaaab */
const W2: f32 = -2.7777778450e-03; /* 0xbb360b61 */
const W3: f32 = 7.9365057172e-04; /* 0x3a500cfd */
const W4: f32 = -5.9518753551e-04; /* 0xba1c065c */
const W5: f32 = 8.3633989561e-04; /* 0x3a5b3dd2 */
const W6: f32 = -1.6309292987e-03; /* 0xbad5c4e8 */
/* sin(PI*x) assuming x > 2^-100, if sin(PI*x)==0 the sign is arbitrary */
fn sin_pi(mut x: f32) -> f32
{
fn sin_pi(mut x: f32) -> f32 {
let mut y: f64;
let mut n: isize;
/* spurious inexact if odd int */
x = 2.0*(x*0.5 - floorf(x*0.5)); /* x mod 2.0 */
x = 2.0 * (x * 0.5 - floorf(x * 0.5)); /* x mod 2.0 */
n = (x*4.0) as isize;
n = (n+1)/2;
y = (x as f64) - (n as f64)*0.5;
n = (x * 4.0) as isize;
n = (n + 1) / 2;
y = (x as f64) - (n as f64) * 0.5;
y *= 3.14159265358979323846;
match n {
1 => k_cosf(y),
2 => k_sinf(-y),
3 => -k_cosf(y),
0|_ => k_sinf(y),
1 => k_cosf(y),
2 => k_sinf(-y),
3 => -k_cosf(y),
0 | _ => k_sinf(y),
}
}
@ -104,8 +103,7 @@ pub fn lgammaf(x: f32) -> f32 {
lgammaf_r(x).0
}
pub fn lgammaf_r(mut x: f32) -> (f32, isize)
{
pub fn lgammaf_r(mut x: f32) -> (f32, isize) {
let u = x.to_bits();
let mut t: f32;
let y: f32;
@ -125,12 +123,13 @@ pub fn lgammaf_r(mut x: f32) -> (f32, isize)
/* purge off +-inf, NaN, +-0, tiny and negative arguments */
signgam = 1;
sign = (u>>31) != 0;
sign = (u >> 31) != 0;
ix = u & 0x7fffffff;
if ix >= 0x7f800000 {
return (x*x, signgam);
return (x * x, signgam);
}
if ix < 0x35000000 { /* |x| < 2**-21, return -log(|x|) */
if ix < 0x35000000 {
/* |x| < 2**-21, return -log(|x|) */
if sign {
signgam = -1;
x = -x;
@ -140,15 +139,16 @@ pub fn lgammaf_r(mut x: f32) -> (f32, isize)
if sign {
x = -x;
t = sin_pi(x);
if t == 0.0 { /* -integer */
return (1.0/(x-x), signgam);
if t == 0.0 {
/* -integer */
return (1.0 / (x - x), signgam);
}
if t > 0.0 {
signgam = -1;
} else {
t = -t;
}
nadj = logf(PI/(t*x));
nadj = logf(PI / (t * x));
} else {
nadj = 0.0;
}
@ -159,13 +159,14 @@ pub fn lgammaf_r(mut x: f32) -> (f32, isize)
}
/* for x < 2.0 */
else if ix < 0x40000000 {
if ix <= 0x3f666666 { /* lgamma(x) = lgamma(x+1)-log(x) */
if ix <= 0x3f666666 {
/* lgamma(x) = lgamma(x+1)-log(x) */
r = -logf(x);
if ix >= 0x3f3b4a20 {
y = 1.0 - x;
i = 0;
} else if ix >= 0x3e6d3308 {
y = x - (TC-1.0);
y = x - (TC - 1.0);
i = 1;
} else {
y = x;
@ -173,10 +174,12 @@ pub fn lgammaf_r(mut x: f32) -> (f32, isize)
}
} else {
r = 0.0;
if ix >= 0x3fdda618 { /* [1.7316,2] */
if ix >= 0x3fdda618 {
/* [1.7316,2] */
y = 2.0 - x;
i = 0;
} else if ix >= 0x3F9da620 { /* [1.23,1.73] */
} else if ix >= 0x3F9da620 {
/* [1.23,1.73] */
y = x - TC;
i = 1;
} else {
@ -186,56 +189,67 @@ pub fn lgammaf_r(mut x: f32) -> (f32, isize)
}
match i {
0 => {
z = y*y;
p1 = A0+z*(A2+z*(A4+z*(A6+z*(A8+z*A10))));
p2 = z*(A1+z*(A3+z*(A5+z*(A7+z*(A9+z*A11)))));
p = y*p1+p2;
r += p - 0.5*y;
z = y * y;
p1 = A0 + z * (A2 + z * (A4 + z * (A6 + z * (A8 + z * A10))));
p2 = z * (A1 + z * (A3 + z * (A5 + z * (A7 + z * (A9 + z * A11)))));
p = y * p1 + p2;
r += p - 0.5 * y;
}
1 => {
z = y*y;
w = z*y;
p1 = T0+w*(T3+w*(T6+w*(T9 +w*T12))); /* parallel comp */
p2 = T1+w*(T4+w*(T7+w*(T10+w*T13)));
p3 = T2+w*(T5+w*(T8+w*(T11+w*T14)));
p = z*p1-(TT-w*(p2+y*p3));
z = y * y;
w = z * y;
p1 = T0 + w * (T3 + w * (T6 + w * (T9 + w * T12))); /* parallel comp */
p2 = T1 + w * (T4 + w * (T7 + w * (T10 + w * T13)));
p3 = T2 + w * (T5 + w * (T8 + w * (T11 + w * T14)));
p = z * p1 - (TT - w * (p2 + y * p3));
r += TF + p;
}
2 => {
p1 = y*(U0+y*(U1+y*(U2+y*(U3+y*(U4+y*U5)))));
p2 = 1.0+y*(V1+y*(V2+y*(V3+y*(V4+y*V5))));
r += -0.5*y + p1/p2;
p1 = y * (U0 + y * (U1 + y * (U2 + y * (U3 + y * (U4 + y * U5)))));
p2 = 1.0 + y * (V1 + y * (V2 + y * (V3 + y * (V4 + y * V5))));
r += -0.5 * y + p1 / p2;
}
#[cfg(feature = "checked")]
_ => unreachable!(),
#[cfg(not(feature = "checked"))]
_ => {}
}
} else if ix < 0x41000000 { /* x < 8.0 */
} else if ix < 0x41000000 {
/* x < 8.0 */
i = x as isize;
y = x - (i as f32);
p = y*(S0+y*(S1+y*(S2+y*(S3+y*(S4+y*(S5+y*S6))))));
q = 1.0+y*(R1+y*(R2+y*(R3+y*(R4+y*(R5+y*R6)))));
r = 0.5*y+p/q;
z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
p = y * (S0 + y * (S1 + y * (S2 + y * (S3 + y * (S4 + y * (S5 + y * S6))))));
q = 1.0 + y * (R1 + y * (R2 + y * (R3 + y * (R4 + y * (R5 + y * R6)))));
r = 0.5 * y + p / q;
z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
// TODO: In C, this was implemented using switch jumps with fallthrough.
// Does this implementation have performance problems?
if i >= 7 { z *= y + 6.0; }
if i >= 6 { z *= y + 5.0; }
if i >= 5 { z *= y + 4.0; }
if i >= 4 { z *= y + 3.0; }
if i >= 7 {
z *= y + 6.0;
}
if i >= 6 {
z *= y + 5.0;
}
if i >= 5 {
z *= y + 4.0;
}
if i >= 4 {
z *= y + 3.0;
}
if i >= 3 {
z *= y + 2.0;
r += logf(z);
}
} else if ix < 0x5c800000 { /* 8.0 <= x < 2**58 */
} else if ix < 0x5c800000 {
/* 8.0 <= x < 2**58 */
t = logf(x);
z = 1.0/x;
y = z*z;
w = W0+z*(W1+y*(W2+y*(W3+y*(W4+y*(W5+y*W6)))));
r = (x-0.5)*(t-1.0)+w;
} else { /* 2**58 <= x <= inf */
r = x*(logf(x)-1.0);
z = 1.0 / x;
y = z * z;
w = W0 + z * (W1 + y * (W2 + y * (W3 + y * (W4 + y * (W5 + y * W6)))));
r = (x - 0.5) * (t - 1.0) + w;
} else {
/* 2**58 <= x <= inf */
r = x * (logf(x) - 1.0);
}
if sign {
r = nadj - r;

10
src/math/mod.rs
View File

@ -118,8 +118,6 @@ mod frexp;
mod frexpf;
mod hypot;
mod hypotf;
mod ldexp;
mod ldexpf;
mod ilogb;
mod ilogbf;
mod j0;
@ -128,6 +126,8 @@ mod j1;
mod j1f;
mod jn;
mod jnf;
mod ldexp;
mod ldexpf;
mod lgamma;
mod lgammaf;
mod log;
@ -192,8 +192,8 @@ pub use self::cosh::cosh;
pub use self::coshf::coshf;
pub use self::erf::erf;
pub use self::erf::erfc;
pub use self::erff::erff;
pub use self::erff::erfcf;
pub use self::erff::erff;
pub use self::exp::exp;
pub use self::exp10::exp10;
pub use self::exp10f::exp10f;
@ -216,8 +216,6 @@ pub use self::frexp::frexp;
pub use self::frexpf::frexpf;
pub use self::hypot::hypot;
pub use self::hypotf::hypotf;
pub use self::ldexp::ldexp;
pub use self::ldexpf::ldexpf;
pub use self::ilogb::ilogb;
pub use self::ilogbf::ilogbf;
pub use self::j0::j0;
@ -232,6 +230,8 @@ pub use self::jn::jn;
pub use self::jn::yn;
pub use self::jnf::jnf;
pub use self::jnf::ynf;
pub use self::ldexp::ldexp;
pub use self::ldexpf::ldexpf;
pub use self::lgamma::lgamma;
pub use self::lgamma::lgamma_r;
pub use self::lgammaf::lgammaf;

13
src/math/modf.rs
View File

@ -2,29 +2,30 @@ pub fn modf(x: f64) -> (f64, f64) {
let rv2: f64;
let mut u = x.to_bits();
let mask: u64;
let e = ((u>>52 & 0x7ff) as isize) - 0x3ff;
let e = ((u >> 52 & 0x7ff) as isize) - 0x3ff;
/* no fractional part */
if e >= 52 {
rv2 = x;
if e == 0x400 && (u<<12) != 0 { /* nan */
if e == 0x400 && (u << 12) != 0 {
/* nan */
return (x, rv2);
}
u &= 1<<63;
u &= 1 << 63;
return (f64::from_bits(u), rv2);
}
/* no integral part*/
if e < 0 {
u &= 1<<63;
u &= 1 << 63;
rv2 = f64::from_bits(u);
return (x, rv2);
}
mask = ((!0)>>12)>>e;
mask = ((!0) >> 12) >> e;
if (u & mask) == 0 {
rv2 = x;
u &= 1<<63;
u &= 1 << 63;
return (f64::from_bits(u), rv2);
}
u &= !mask;

7
src/math/modff.rs
View File

@ -2,12 +2,13 @@ pub fn modff(x: f32) -> (f32, f32) {
let rv2: f32;
let mut u: u32 = x.to_bits();
let mask: u32;
let e = ((u>>23 & 0xff) as isize) - 0x7f;
let e = ((u >> 23 & 0xff) as isize) - 0x7f;
/* no fractional part */
if e >= 23 {
rv2 = x;
if e == 0x80 && (u<<9) != 0 { /* nan */
if e == 0x80 && (u << 9) != 0 {
/* nan */
return (x, rv2);
}
u &= 0x80000000;
@ -20,7 +21,7 @@ pub fn modff(x: f32) -> (f32, f32) {
return (x, rv2);
}
mask = 0x007fffff>>e;
mask = 0x007fffff >> e;
if (u & mask) == 0 {
rv2 = x;
u &= 0x80000000;

33
src/math/remquo.rs
View File

@ -1,26 +1,25 @@
pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize)
{
pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize) {
let ux: u64 = x.to_bits();
let mut uy: u64 = y.to_bits();
let mut ex = ((ux>>52) & 0x7ff) as isize;
let mut ey = ((uy>>52) & 0x7ff) as isize;
let sx = (ux>>63) != 0;
let sy = (uy>>63) != 0;
let mut ex = ((ux >> 52) & 0x7ff) as isize;
let mut ey = ((uy >> 52) & 0x7ff) as isize;
let sx = (ux >> 63) != 0;
let sy = (uy >> 63) != 0;
let mut q: u32;
let mut i: u64;
let mut uxi: u64 = ux;
if (uy<<1) == 0 || y.is_nan() || ex == 0x7ff {
return ((x*y)/(x*y), 0);
if (uy << 1) == 0 || y.is_nan() || ex == 0x7ff {
return ((x * y) / (x * y), 0);
}
if (ux<<1) == 0 {
if (ux << 1) == 0 {
return (x, 0);
}
/* normalize x and y */
if ex == 0 {
i = uxi << 12;
while (i>>63) == 0 {
while (i >> 63) == 0 {
ex -= 1;
i <<= 1;
}
@ -30,8 +29,8 @@ pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize)
uxi |= 1 << 52;
}
if ey == 0 {
i = uy<<12;
while (i>>63) == 0 {
i = uy << 12;
while (i >> 63) == 0 {
ey -= 1;
i <<= 1;
}
@ -43,14 +42,14 @@ pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize)
q = 0;
if ex+1 != ey {
if ex + 1 != ey {
if ex < ey {
return (x, 0);
}
/* x mod y */
while ex > ey {
i = uxi - uy;
if (i>>63) == 0 {
if (i >> 63) == 0 {
uxi = i;
q += 1;
}
@ -59,14 +58,14 @@ pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize)
ex -= 1;
}
i = uxi - uy;
if (i>>63) == 0 {
if (i >> 63) == 0 {
uxi = i;
q += 1;
}
if uxi == 0 {
ex = -60;
} else {
while (uxi>>52) == 0 {
while (uxi >> 52) == 0 {
uxi <<= 1;
ex -= 1;
}
@ -84,7 +83,7 @@ pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize)
if sy {
y = -y;
}
if ex == ey || (ex+1 == ey && (2.0*x > y || (2.0*x == y && (q%2) != 0))) {
if ex == ey || (ex + 1 == ey && (2.0 * x > y || (2.0 * x == y && (q % 2) != 0))) {
x -= y;
q += 1;
}

37
src/math/remquof.rs
View File

@ -1,26 +1,25 @@
pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize)
{
pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize) {
let ux: u32 = x.to_bits();
let mut uy: u32 = y.to_bits();
let mut ex = ((ux>>23) & 0xff) as isize;
let mut ey = ((uy>>23) & 0xff) as isize;
let sx = (ux>>31) != 0;
let sy = (uy>>31) != 0;
let mut ex = ((ux >> 23) & 0xff) as isize;
let mut ey = ((uy >> 23) & 0xff) as isize;
let sx = (ux >> 31) != 0;
let sy = (uy >> 31) != 0;
let mut q: u32;
let mut i: u32;
let mut uxi: u32 = ux;
if (uy<<1) == 0 || y.is_nan() || ex == 0xff {
return ((x*y)/(x*y), 0);
if (uy << 1) == 0 || y.is_nan() || ex == 0xff {
return ((x * y) / (x * y), 0);
}
if (ux<<1) == 0 {
if (ux << 1) == 0 {
return (x, 0);
}
/* normalize x and y */
if ex == 0 {
i = uxi<<9;
while (i>>31) == 0 {
i = uxi << 9;
while (i >> 31) == 0 {
ex -= 1;
i <<= 1;
}
@ -30,8 +29,8 @@ pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize)
uxi |= 1 << 23;
}
if ey == 0 {
i = uy<<9;
while (i>>31) == 0 {
i = uy << 9;
while (i >> 31) == 0 {
ey -= 1;
i <<= 1;
}
@ -42,14 +41,14 @@ pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize)
}
q = 0;
if ex+1 != ey {
if ex + 1 != ey {
if ex < ey {
return (x, 0);
}
/* x mod y */
while ex > ey {
i = uxi - uy;
if (i>>31) == 0 {
if (i >> 31) == 0 {
uxi = i;
q += 1;
}
@ -58,14 +57,14 @@ pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize)
ex -= 1;
}
i = uxi - uy;
if (i>>31) == 0 {
if (i >> 31) == 0 {
uxi = i;
q += 1;
}
if uxi == 0 {
ex = -30;
} else {
while (uxi>>23) == 0 {
while (uxi >> 23) == 0 {
uxi <<= 1;
ex -= 1;
}
@ -83,12 +82,12 @@ pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize)
if sy {
y = -y;
}
if ex == ey || (ex+1 == ey && (2.0*x > y || (2.0*x == y && (q%2) != 0))) {
if ex == ey || (ex + 1 == ey && (2.0 * x > y || (2.0 * x == y && (q % 2) != 0))) {
x -= y;
q += 1;
}
q &= 0x7fffffff;
let quo = if sx^sy { -(q as isize) } else { q as isize };
let quo = if sx ^ sy { -(q as isize) } else { q as isize };
if sx {
(-x, quo)
} else {

9
src/math/sincos.rs
View File

@ -12,8 +12,7 @@
use super::{get_high_word, k_cos, k_sin, rem_pio2};
pub fn sincos(x: f64) -> (f64, f64)
{
pub fn sincos(x: f64) -> (f64, f64) {
let s: f64;
let c: f64;
let mut ix: u32;
@ -28,9 +27,9 @@ pub fn sincos(x: f64) -> (f64, f64)
/* raise inexact if x!=0 and underflow if subnormal */
let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120 == 2^120
if ix < 0x00100000 {
force_eval!(x/x1p120);
force_eval!(x / x1p120);
} else {
force_eval!(x+x1p120);
force_eval!(x + x1p120);
}
return (x, 1.0);
}
@ -47,7 +46,7 @@ pub fn sincos(x: f64) -> (f64, f64)
let (n, y0, y1) = rem_pio2(x);
s = k_sin(y0, y1, 1);
c = k_cos(y0, y1);
match n&3 {
match n & 3 {
0 => (s, c),
1 => (c, -s),
2 => (-s, -c),

23
src/math/sincosf.rs
View File

@ -18,13 +18,12 @@ use super::{k_cosf, k_sinf, rem_pio2f};
/* Small multiples of pi/2 rounded to double precision. */
const PI_2: f32 = 0.5 * 3.1415926535897931160E+00;
const S1PIO2: f32 = 1.0*PI_2; /* 0x3FF921FB, 0x54442D18 */
const S2PIO2: f32 = 2.0*PI_2; /* 0x400921FB, 0x54442D18 */
const S3PIO2: f32 = 3.0*PI_2; /* 0x4012D97C, 0x7F3321D2 */
const S4PIO2: f32 = 4.0*PI_2; /* 0x401921FB, 0x54442D18 */
const S1PIO2: f32 = 1.0 * PI_2; /* 0x3FF921FB, 0x54442D18 */
const S2PIO2: f32 = 2.0 * PI_2; /* 0x400921FB, 0x54442D18 */
const S3PIO2: f32 = 3.0 * PI_2; /* 0x4012D97C, 0x7F3321D2 */
const S4PIO2: f32 = 4.0 * PI_2; /* 0x401921FB, 0x54442D18 */
pub fn sincosf(x: f32) -> (f32, f32)
{
pub fn sincosf(x: f32) -> (f32, f32) {
let s: f32;
let c: f32;
let mut ix: u32;
@ -42,9 +41,9 @@ pub fn sincosf(x: f32) -> (f32, f32)
let x1p120 = f32::from_bits(0x7b800000); // 0x1p120 == 2^120
if ix < 0x00100000 {
force_eval!(x/x1p120);
force_eval!(x / x1p120);
} else {
force_eval!(x+x1p120);
force_eval!(x + x1p120);
}
return (x, 1.0);
}
@ -53,7 +52,8 @@ pub fn sincosf(x: f32) -> (f32, f32)
/* |x| ~<= 5*pi/4 */
if ix <= 0x407b53d1 {
if ix <= 0x4016cbe3 { /* |x| ~<= 3pi/4 */
if ix <= 0x4016cbe3 {
/* |x| ~<= 3pi/4 */
if sign {
s = -k_cosf((x + S1PIO2) as f64);
c = k_sinf((x + S1PIO2) as f64);
@ -78,7 +78,8 @@ pub fn sincosf(x: f32) -> (f32, f32)
/* |x| ~<= 9*pi/4 */
if ix <= 0x40e231d5 {
if ix <= 0x40afeddf { /* |x| ~<= 7*pi/4 */
if ix <= 0x40afeddf {
/* |x| ~<= 7*pi/4 */
if sign {
s = k_cosf((x + S3PIO2) as f64);
c = -k_sinf((x + S3PIO2) as f64);
@ -109,7 +110,7 @@ pub fn sincosf(x: f32) -> (f32, f32)
let (n, y) = rem_pio2f(x);
s = k_sinf(y);
c = k_cosf(y);
match n&3 {
match n & 3 {
0 => (s, c),
1 => (c, -s),
2 => (-s, -c),

84
src/math/tgamma.rs
View File

@ -28,8 +28,7 @@ use super::{exp, floor, k_cos, k_sin, pow};
const PI: f64 = 3.141592653589793238462643383279502884;
/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
fn sinpi(mut x: f64) -> f64
{
fn sinpi(mut x: f64) -> f64 {
let mut n: isize;
/* argument reduction: x = |x| mod 2 */
@ -39,22 +38,22 @@ fn sinpi(mut x: f64) -> f64
/* reduce x into [-.25,.25] */
n = (4.0 * x) as isize;
n = (n+1)/2;
n = (n + 1) / 2;
x -= (n as f64) * 0.5;
x *= PI;
match n {
1 => k_cos(x, 0.0),
2 => k_sin(-x, 0.0, 0),
3 => -k_cos(x, 0.0),
0|_ => k_sin(x, 0.0, 0),
1 => k_cos(x, 0.0),
2 => k_sin(-x, 0.0, 0),
3 => -k_cos(x, 0.0),
0 | _ => k_sin(x, 0.0, 0),
}
}
const N: usize = 12;
//static const double g = 6.024680040776729583740234375;
const GMHALF: f64 = 5.524680040776729583740234375;
const SNUM: [f64; N+1] = [
const SNUM: [f64; N + 1] = [
23531376880.410759688572007674451636754734846804940,
42919803642.649098768957899047001988850926355848959,
35711959237.355668049440185451547166705960488635843,
@ -69,21 +68,50 @@ const SNUM: [f64; N+1] = [
210.82427775157934587250973392071336271166969580291,
2.5066282746310002701649081771338373386264310793408,
];
const SDEN: [f64; N+1] = [
0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0,
45995730.0, 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0,
const SDEN: [f64; N + 1] = [
0.0,
39916800.0,
120543840.0,
150917976.0,
105258076.0,
45995730.0,
13339535.0,
2637558.0,
357423.0,
32670.0,
1925.0,
66.0,
1.0,
];
/* n! for small integer n */
const FACT: [f64; 23] = [
1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, 3628800.0,
39916800.0, 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0,
20922789888000.0, 355687428096000.0, 6402373705728000.0, 121645100408832000.0,
2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
1.0,
1.0,
2.0,
6.0,
24.0,
120.0,
720.0,
5040.0,
40320.0,
362880.0,
3628800.0,
39916800.0,
479001600.0,
6227020800.0,
87178291200.0,
1307674368000.0,
20922789888000.0,
355687428096000.0,
6402373705728000.0,
121645100408832000.0,
2432902008176640000.0,
51090942171709440000.0,
1124000727777607680000.0,
];
/* S(x) rational function for positive x */
fn s(x: f64) -> f64
{
fn s(x: f64) -> f64 {
let mut num: f64 = 0.0;
let mut den: f64 = 0.0;
@ -99,11 +127,10 @@ fn s(x: f64) -> f64
den = den / x + SDEN[i];
}
}
return num/den;
return num / den;
}
pub fn tgamma(mut x: f64) -> f64
{
pub fn tgamma(mut x: f64) -> f64 {
let u: u64 = x.to_bits();
let absx: f64;
let mut y: f64;
@ -111,23 +138,23 @@ pub fn tgamma(mut x: f64) -> f64
let mut z: f64;
let mut r: f64;
let ix: u32 = ((u >> 32) as u32) & 0x7fffffff;
let sign: bool = (u>>64) != 0;
let sign: bool = (u >> 63) != 0;
/* special cases */
if ix >= 0x7ff00000 {
/* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
return x + core::f64::INFINITY;
}
if ix < ((0x3ff-54)<<20) {
if ix < ((0x3ff - 54) << 20) {
/* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
return 1.0/x;
return 1.0 / x;
}
/* integer arguments */
/* raise inexact when non-integer */
if x == floor(x) {
if sign {
return 0.0/0.0;
return 0.0 / 0.0;
}
if x <= FACT.len() as f64 {
return FACT[(x as usize) - 1];
@ -136,10 +163,11 @@ pub fn tgamma(mut x: f64) -> f64
/* x >= 172: tgamma(x)=inf with overflow */
/* x =< -184: tgamma(x)=+-0 with underflow */
if ix >= 0x40670000 { /* |x| >= 184 */
if ix >= 0x40670000 {
/* |x| >= 184 */
if sign {
let x1p_126 = f64::from_bits(0x3810000000000000); // 0x1p-126 == 2^-126
force_eval!((x1p_126/x) as f32);
force_eval!((x1p_126 / x) as f32);
if floor(x) * 0.5 == floor(x * 0.5) {
return 0.0;
} else {
@ -172,8 +200,8 @@ pub fn tgamma(mut x: f64) -> f64
dy = -dy;
z = -z;
}
r += dy * (GMHALF+0.5) * r / y;
z = pow(y, 0.5*z);
r += dy * (GMHALF + 0.5) * r / y;
z = pow(y, 0.5 * z);
y = r * z * z;
return y;
}

2
src/math/tgammaf.rs
View File

@ -1,4 +1,4 @@
use super::{tgamma};
use super::tgamma;
pub fn tgammaf(x: f32) -> f32 {
tgamma(x as f64) as f32