Small adjustments.

This commit is contained in:
2016-06-02 17:05:17 +02:00
parent bbd2a2ecaf
commit 410e9612b8

View File

@@ -131,15 +131,15 @@ fn series(m: i32, id: i32) -> f64 {
const EPS : f64 = 1e-17;
let mut ak : f64;
let mut s : f64 = 0.0;
let mut t : f64;
let mut s : f64 = 0.0;
for k in 0..id {
ak = (8 * k + m) as f64;
t = expm((id - k) as f64, ak);
s = s + t / ak;
s = s - (s as i32) as f64;
s += t / ak;
s -= (s as i32) as f64;
}
for k in id..id + 101 {
@@ -150,20 +150,19 @@ fn series(m: i32, id: i32) -> f64 {
break;
}
s = s + t;
s = s - (s as i32) as f64;
s += t;
s -= (s as i32) as f64;
}
s
}
fn expm(p: f64, ak: f64) -> f64 {
const NTP : usize = 25;
if ak == 1.0 {
return 0.0;
}
let mut i : i32 = 0;
let mut p1 : f64 = p;
let mut pt : f64;
let mut r : f64 = 1.0;
const NTP : usize = 25;
lazy_static! {
static ref TP : Vec<f64> = {
@@ -177,32 +176,29 @@ fn expm(p: f64, ak: f64) -> f64 {
};
}
if ak == 1.0 {
return 0.0;
let mut i : usize = 0;
while TP[i] <= p {
i += 1;
}
for n in 0..NTP {
if TP[n] > p {
i = n as i32;
break;
}
}
pt = TP[(i - 1) as usize];
let mut pt = TP[i - 1];
let mut p1 : f64 = p;
let mut r : f64 = 1.0;
for _ in 1..i + 1 {
if p1 >= pt {
r = 16.0 * r;
r = r - ((r / ak) as i32) as f64 * ak;
p1 = p1 - pt;
p1 -= pt;
}
pt = 0.5 * pt;
if pt >= 1.0 {
r = r * r;
r = r - ((r / ak) as i32) as f64 * ak;
r *= r;
r -= ((r / ak) as i32) as f64 * ak;
}
}