Skip to content

Commit 6cbaa22

Browse files
jackyarndleyjackyarndley
committed
fixed issues with series approximation for now
1 parent 6777e68 commit 6cbaa22

File tree

3 files changed

+73
-72
lines changed

3 files changed

+73
-72
lines changed

‎output.png‎

147 KB
Loading

‎src/math/series_approximation.rs‎

Lines changed: 64 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -46,88 +46,85 @@ impl SeriesApproximation {
4646
}
4747
}
4848

49-
pub fn step(&mut self) -> bool {
50-
self.z.square_mut();
51-
self.z += &self.c;
52-
53-
let z_fixed = to_fixed(&self.z);
54-
55-
if z_fixed.re.abs() < 1e-300 && z_fixed.im.abs() < 1e-300 {
56-
println!("found slow at: {}", self.current_iteration);
57-
}
49+
pub fn run(&mut self) {
50+
self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
51+
self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
52+
self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
53+
self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
5854

59-
self.next_coefficients[0] = to_extended(&self.z);
60-
self.next_coefficients[1] = self.coefficients[0] * self.coefficients[1] * 2.0 + ComplexExtended::new2(1.0, 0.0, 0);
55+
let add_value = ComplexExtended::new2(1.0, 0.0, 0);
6156

62-
// Calculate the new coefficents
63-
for k in 2..=self.order {
64-
let mut sum = self.coefficients[0] * self.coefficients[k];
57+
// Can be changed later into a better loop - this function could also return some more information
58+
while self.current_iteration < self.maximum_iteration {
59+
self.z.square_mut();
60+
self.z += &self.c;
6561

66-
for j in 1..=((k - 1) / 2) {
67-
sum += self.coefficients[j] * self.coefficients[k - j];
68-
}
69-
sum *= 2.0;
62+
let z_fixed = to_fixed(&self.z);
7063

71-
// If even, we include the mid term as well
72-
if k % 2 == 0 {
73-
sum += self.coefficients[k / 2] * self.coefficients[k / 2];
64+
if z_fixed.re.abs() < 1e-300 && z_fixed.im.abs() < 1e-300 {
65+
println!("found slow at: {}", self.current_iteration);
7466
}
7567

76-
sum.reduce();
77-
self.next_coefficients[k] = sum;
78-
}
79-
80-
// this section takes about half of the time
81-
for i in 0..self.original_probes.len() {
82-
// step the probe points using perturbation
83-
self.current_probes[i] = self.current_probes[i] * (self.coefficients[0] * 2.0 + self.current_probes[i]);
84-
self.current_probes[i] += self.original_probes[i];
85-
86-
// TODO this probably does not need to be done every iteration
87-
self.current_probes[i].reduce();
88-
89-
// get the new approximations
90-
let mut series_probe = self.next_coefficients[1] * self.approximation_probes[i][0];
91-
let mut derivative_probe = self.next_coefficients[1] * self.approximation_probes_derivative[i][0];
68+
self.next_coefficients[0] = to_extended(&self.z);
69+
self.next_coefficients[1] = self.coefficients[0] * self.coefficients[1] * 2.0 + add_value;
9270

71+
// Calculate the new coefficents
9372
for k in 2..=self.order {
94-
series_probe += self.next_coefficients[k] * self.approximation_probes[i][k - 1];
95-
derivative_probe += self.next_coefficients[k] * self.approximation_probes_derivative[i][k - 1];
96-
};
73+
let mut sum = self.coefficients[0] * self.coefficients[k];
9774

98-
let relative_error = (self.current_probes[i] - series_probe).norm_square();
99-
let mut derivative = derivative_probe.norm_square();
75+
for j in 1..=((k - 1) / 2) {
76+
sum += self.coefficients[j] * self.coefficients[k - j];
77+
}
78+
sum *= 2.0;
10079

101-
// Check to make sure that the derivative is greater than or equal to 1
102-
if derivative.to_float() < 1.0 {
103-
derivative = FloatExtended::new(1.0, 0);
104-
}
80+
// If even, we include the mid term as well
81+
if k % 2 == 0 {
82+
sum += self.coefficients[k / 2] * self.coefficients[k / 2];
83+
}
10584

106-
// Check that the error over the derivative is less than the pixel spacing
107-
if relative_error / derivative > self.delta_pixel_square {
108-
self.z -= &self.c;
109-
self.z.sqrt_mut();
110-
return false;
85+
sum.reduce();
86+
self.next_coefficients[k] = sum;
11187
}
112-
}
11388

114-
// If the approximation is valid, continue
115-
self.current_iteration += 1;
116-
117-
// Swap the two coefficients buffers
118-
swap(&mut self.coefficients, &mut self.next_coefficients);
119-
self.current_iteration < self.maximum_iteration
120-
}
89+
// this section takes about half of the time
90+
for i in 0..self.original_probes.len() {
91+
// step the probe points using perturbation
92+
self.current_probes[i] = self.current_probes[i] * (self.coefficients[0] * 2.0 + self.current_probes[i]);
93+
self.current_probes[i] += self.original_probes[i];
94+
95+
// TODO this probably does not need to be done every iteration
96+
self.current_probes[i].reduce();
97+
98+
// get the new approximations
99+
let mut series_probe = self.next_coefficients[1] * self.approximation_probes[i][0];
100+
let mut derivative_probe = self.next_coefficients[1] * self.approximation_probes_derivative[i][0];
101+
102+
for k in 2..=self.order {
103+
series_probe += self.next_coefficients[k] * self.approximation_probes[i][k - 1];
104+
derivative_probe += self.next_coefficients[k] * self.approximation_probes_derivative[i][k - 1];
105+
};
106+
107+
let relative_error = (self.current_probes[i] - series_probe).norm_square();
108+
let mut derivative = derivative_probe.norm_square();
109+
110+
// Check to make sure that the derivative is greater than or equal to 1
111+
if derivative.to_float() < 1.0 {
112+
derivative = FloatExtended::new(1.0, 0);
113+
}
114+
115+
// Check that the error over the derivative is less than the pixel spacing
116+
if relative_error / derivative > self.delta_pixel_square {
117+
self.z -= &self.c;
118+
self.z.sqrt_mut();
119+
return;
120+
}
121+
}
121122

122-
pub fn run(&mut self) {
123-
self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
124-
self.add_probe(ComplexExtended::new2(self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
125-
self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
126-
self.add_probe(ComplexExtended::new2(-self.delta_top_left.mantissa.re, -self.delta_top_left.mantissa.im, self.delta_top_left.exponent));
123+
// If the approximation is valid, continue
124+
self.current_iteration += 1;
127125

128-
// Can be changed later into a better loop - this function could also return some more information
129-
while self.step() {
130-
continue;
126+
// Swap the two coefficients buffers
127+
swap(&mut self.coefficients, &mut self.next_coefficients);
131128
}
132129
}
133130

‎src/util/complex_extended.rs‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -10,22 +10,26 @@ pub struct ComplexExtended {
1010
pub exponent: i32,
1111
}
1212

13-
// TODO we can make a struct with shared exponent
1413
impl ComplexExtended {
14+
// maybe set up different versions for these function if the reduce is needed
1515
#[inline]
1616
pub fn new(mantissa: Complex<f64>, exponent: i32) -> Self {
17-
ComplexExtended {
17+
let mut temp = ComplexExtended {
1818
mantissa,
1919
exponent
20-
}
20+
};
21+
temp.reduce();
22+
temp
2123
}
2224

2325
#[inline]
2426
pub fn new2(re: f64, im: f64, exponent: i32) -> Self {
25-
ComplexExtended {
27+
let mut temp = ComplexExtended {
2628
mantissa: Complex::<f64>::new(re, im),
2729
exponent
28-
}
30+
};
31+
temp.reduce();
32+
temp
2933
}
3034

3135
#[inline]

0 commit comments

Comments
 (0)