sciforge_lib/maths/complex/
ops.rs

1use super::types::Complex;
2
3pub fn roots_of_unity(n: usize) -> Vec<Complex> {
4    (0..n)
5        .map(|k| {
6            let theta = 2.0 * std::f64::consts::PI * k as f64 / n as f64;
7            Complex::from_polar(1.0, theta)
8        })
9        .collect()
10}
11
12pub fn complex_polynomial_eval(coeffs: &[Complex], z: Complex) -> Complex {
13    let mut result = Complex::zero();
14    let mut zn = Complex::one();
15    for &c in coeffs {
16        result = result + c * zn;
17        zn = zn * z;
18    }
19    result
20}
21
22pub fn mandelbrot_iterate(c: Complex, max_iter: u32) -> u32 {
23    let mut z = Complex::zero();
24    for i in 0..max_iter {
25        if z.norm_sq() > 4.0 {
26            return i;
27        }
28        z = z * z + c;
29    }
30    max_iter
31}
32
33pub fn julia_iterate(z0: Complex, c: Complex, max_iter: u32) -> u32 {
34    let mut z = z0;
35    for i in 0..max_iter {
36        if z.norm_sq() > 4.0 {
37            return i;
38        }
39        z = z * z + c;
40    }
41    max_iter
42}
43
44pub fn newton_fractal_step(z: Complex, coeffs: &[Complex], deriv_coeffs: &[Complex]) -> Complex {
45    let fz = complex_polynomial_eval(coeffs, z);
46    let dfz = complex_polynomial_eval(deriv_coeffs, z);
47    if dfz.norm_sq() < 1e-30 {
48        return z;
49    }
50    z - fz / dfz
51}
52
53pub fn complex_matrix_mul(a: &[Vec<Complex>], b: &[Vec<Complex>]) -> Vec<Vec<Complex>> {
54    let m = a.len();
55    let n = b[0].len();
56    let mut result = vec![vec![Complex::zero(); n]; m];
57    for (i, ri) in result.iter_mut().enumerate() {
58        for (j, rij) in ri.iter_mut().enumerate() {
59            for (k, &aik) in a[i].iter().enumerate() {
60                *rij = *rij + aik * b[k][j];
61            }
62        }
63    }
64    result
65}
66
67pub fn complex_matrix_det(m: &[Vec<Complex>]) -> Complex {
68    let n = m.len();
69    if n == 1 {
70        return m[0][0];
71    }
72    if n == 2 {
73        return m[0][0] * m[1][1] - m[0][1] * m[1][0];
74    }
75    let mut det = Complex::zero();
76    for j in 0..n {
77        let sign = if j % 2 == 0 {
78            Complex::one()
79        } else {
80            -Complex::one()
81        };
82        let minor: Vec<Vec<Complex>> = (1..n)
83            .map(|i| (0..n).filter(|&k| k != j).map(|k| m[i][k]).collect())
84            .collect();
85        det = det + sign * m[0][j] * complex_matrix_det(&minor);
86    }
87    det
88}
89
90pub fn complex_exp(z: Complex) -> Complex {
91    let r = z.re.exp();
92    Complex::new(r * z.im.cos(), r * z.im.sin())
93}
94
95pub fn complex_log(z: Complex) -> Complex {
96    Complex::new(z.norm().ln(), z.im.atan2(z.re))
97}
98
99pub fn complex_sin(z: Complex) -> Complex {
100    Complex::new(z.re.sin() * z.im.cosh(), z.re.cos() * z.im.sinh())
101}
102
103pub fn complex_cos(z: Complex) -> Complex {
104    Complex::new(z.re.cos() * z.im.cosh(), -z.re.sin() * z.im.sinh())
105}
106
107pub fn complex_tan(z: Complex) -> Complex {
108    let s = complex_sin(z);
109    let c = complex_cos(z);
110    if c.norm_sq() < 1e-30 {
111        return Complex::new(f64::INFINITY, 0.0);
112    }
113    s / c
114}
115
116pub fn complex_sinh(z: Complex) -> Complex {
117    Complex::new(z.re.sinh() * z.im.cos(), z.re.cosh() * z.im.sin())
118}
119
120pub fn complex_cosh(z: Complex) -> Complex {
121    Complex::new(z.re.cosh() * z.im.cos(), z.re.sinh() * z.im.sin())
122}
123
124pub fn complex_sqrt(z: Complex) -> Complex {
125    let r = z.norm();
126    let theta = z.im.atan2(z.re);
127    Complex::from_polar(r.sqrt(), theta / 2.0)
128}
129
130pub fn complex_power(z: Complex, w: Complex) -> Complex {
131    if z.norm_sq() < 1e-60 {
132        return Complex::zero();
133    }
134    complex_exp(w * complex_log(z))
135}
136
137pub fn complex_power_real(z: Complex, n: f64) -> Complex {
138    let r = z.norm();
139    let theta = z.im.atan2(z.re);
140    Complex::from_polar(r.powf(n), n * theta)
141}
142
143pub fn mobius_transform(z: Complex, a: Complex, b: Complex, c: Complex, d: Complex) -> Complex {
144    let num = a * z + b;
145    let den = c * z + d;
146    if den.norm_sq() < 1e-60 {
147        return Complex::new(f64::INFINITY, 0.0);
148    }
149    num / den
150}
151
152pub fn bilinear_transform(s: Complex, t_sample: f64) -> Complex {
153    let one = Complex::one();
154    let half_t = Complex::new(t_sample * 0.5, 0.0);
155    let num = one + s * half_t;
156    let den = one - s * half_t;
157    if den.norm_sq() < 1e-60 {
158        return Complex::new(f64::INFINITY, 0.0);
159    }
160    num / den
161}
162
163pub fn complex_contour_integral(f: impl Fn(Complex) -> Complex, path: &[Complex]) -> Complex {
164    let mut result = Complex::zero();
165    for i in 1..path.len() {
166        let dz = path[i] - path[i - 1];
167        let mid = Complex::new(
168            (path[i].re + path[i - 1].re) * 0.5,
169            (path[i].im + path[i - 1].im) * 0.5,
170        );
171        result = result + f(mid) * dz;
172    }
173    result
174}
175
176pub fn complex_conjugate_transpose(m: &[Vec<Complex>]) -> Vec<Vec<Complex>> {
177    let rows = m.len();
178    let cols = if rows > 0 { m[0].len() } else { 0 };
179    let mut result = vec![vec![Complex::zero(); rows]; cols];
180    for (i, mi) in m.iter().enumerate() {
181        for (j, &mij) in mi.iter().enumerate() {
182            result[j][i] = mij.conj();
183        }
184    }
185    result
186}
187
188pub fn complex_matrix_trace(m: &[Vec<Complex>]) -> Complex {
189    let mut trace = Complex::zero();
190    for (i, mi) in m.iter().enumerate() {
191        trace = trace + mi[i];
192    }
193    trace
194}
195
196pub fn complex_matrix_add(a: &[Vec<Complex>], b: &[Vec<Complex>]) -> Vec<Vec<Complex>> {
197    a.iter()
198        .zip(b)
199        .map(|(ra, rb)| ra.iter().zip(rb).map(|(&x, &y)| x + y).collect())
200        .collect()
201}
202
203pub fn complex_matrix_scale(m: &[Vec<Complex>], s: Complex) -> Vec<Vec<Complex>> {
204    m.iter()
205        .map(|row| row.iter().map(|&x| x * s).collect())
206        .collect()
207}
208
209pub fn complex_dft(input: &[Complex]) -> Vec<Complex> {
210    let n = input.len();
211    (0..n)
212        .map(|k| {
213            let mut sum = Complex::zero();
214            for (j, &x) in input.iter().enumerate() {
215                let angle = -2.0 * std::f64::consts::PI * k as f64 * j as f64 / n as f64;
216                sum = sum + x * Complex::from_polar(1.0, angle);
217            }
218            sum
219        })
220        .collect()
221}
222
223pub fn complex_idft(input: &[Complex]) -> Vec<Complex> {
224    let n = input.len();
225    let scale = 1.0 / n as f64;
226    (0..n)
227        .map(|k| {
228            let mut sum = Complex::zero();
229            for (j, &x) in input.iter().enumerate() {
230                let angle = 2.0 * std::f64::consts::PI * k as f64 * j as f64 / n as f64;
231                sum = sum + x * Complex::from_polar(1.0, angle);
232            }
233            sum.scale(scale)
234        })
235        .collect()
236}
sciforge_lib/maths/complex/ops.rs

sciforge_lib/maths/complex/
ops.rs
