oxiproj-core 0.1.2

Foundation types for OxiProj: coordinates, errors, ellipsoids, datums, and units.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

//! Complex polynomial evaluation ported from PROJ `src/zpoly1.cpp`.

use crate::math::complex::Complex;

/// Ported from src/zpoly1.cpp (`pj_zpoly1`).
///
/// Evaluates `z * (c[0] + z*(c[1] + ... + z*c[n]))` — PROJ's complex polynomial
/// (Horner form over `c[0..=n]`, with a final multiply by `z`).
pub fn pj_zpoly1(z: Complex, c: &[Complex], n: usize) -> Complex {
    let mut a = c[n];
    let mut k = n;
    while k > 0 {
        let t = a.r;
        a.r = c[k - 1].r + z.r * t - z.i * a.i;
        a.i = c[k - 1].i + z.r * a.i + z.i * t;
        k -= 1;
    }
    let t = a.r;
    let nr = z.r * t - z.i * a.i;
    let ni = z.r * a.i + z.i * t;
    Complex { r: nr, i: ni }
}

/// Ported from src/zpoly1.cpp (`pj_zpolyd1`).
///
/// Evaluates the same complex polynomial form as [`pj_zpoly1`] and, in the same
/// pass, accumulates its derivative. Returns `(value, derivative)`.
///
/// This mirrors the C routine `COMPLEX pj_zpolyd1(COMPLEX z, const COMPLEX *C,
/// int n, COMPLEX *der)`, which writes the derivative through `der` and returns
/// the value. The `first` flag skips the derivative accumulator update on the
/// first loop iteration, exactly as in the original. The `C` pointer walk
/// (`a = *(C += n)` then `(--C)->r`) is translated to safe indexing over
/// `c: &[Complex]` with an index running from `n` down to `0`.
pub fn pj_zpolyd1(z: Complex, c: &[Complex], n: usize) -> (Complex, Complex) {
    // a = *(C += n);  b = a;  (cursor sits at index `n`)
    let mut a = c[n];
    let mut b = a;
    let mut first = true;
    let mut k = n;
    // while (n-- > 0) { ... }  — executes `n` times.
    while k > 0 {
        if first {
            first = false;
        } else {
            let t = b.r;
            b.r = a.r + z.r * t - z.i * b.i;
            b.i = a.i + z.r * b.i + z.i * t;
        }
        // a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;  — pre-decrement cursor.
        let t = a.r;
        a.r = c[k - 1].r + z.r * t - z.i * a.i;
        a.i = c[k - 1].i + z.r * a.i + z.i * t;
        k -= 1;
    }
    // final derivative update: b.r = a.r + z.r*(t=b.r) - z.i*b.i; ...
    let tb = b.r;
    let der_r = a.r + z.r * tb - z.i * b.i;
    let der_i = a.i + z.r * b.i + z.i * tb;
    let der = Complex { r: der_r, i: der_i };
    // final value: a.r = z.r*(t=a.r) - z.i*a.i; ...
    let ta = a.r;
    let val_r = z.r * ta - z.i * a.i;
    let val_i = z.r * a.i + z.i * ta;
    let val = Complex { r: val_r, i: val_i };
    (val, der)
}

#[cfg(test)]
mod tests {
    use super::*;

    fn cmul(a: Complex, b: Complex) -> Complex {
        Complex::new(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r)
    }

    #[test]
    fn zpoly1_evaluates_z_squared_plus_z_cubed() {
        // c = [0, 1, 1], n = 2.  PROJ's pj_zpoly1 computes
        // z*(c0 + z*(c1 + z*c2)) = z*(0 + z*(1 + z)) = z^2 + z^3.
        let c = [
            Complex::new(0.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
        ];
        let z = Complex::new(0.3, 0.4);
        let result = pj_zpoly1(z, &c, 2);
        // z^2 = (-0.07, 0.24); z^3 = (-0.117, 0.044); sum = (-0.187, 0.284).
        assert!((result.r - (-0.187)).abs() < 1e-12);
        assert!((result.i - 0.284).abs() < 1e-12);
    }

    #[test]
    fn zpoly1_matches_explicit_cmul_form() {
        // Cross-check against the explicit PROJ form z*(c0 + z*(c1 + z*c2)).
        let c = [
            Complex::new(0.0, 0.0),
            Complex::new(0.7, -0.2),
            Complex::new(-0.5, 0.9),
        ];
        let z = Complex::new(0.15, -0.35);
        let result = pj_zpoly1(z, &c, 2);
        let inner1 = Complex::new(c[1].r + cmul(z, c[2]).r, c[1].i + cmul(z, c[2]).i);
        let inner0 = Complex::new(c[0].r + cmul(z, inner1).r, c[0].i + cmul(z, inner1).i);
        let expected = cmul(z, inner0);
        assert!((result.r - expected.r).abs() < 1e-12);
        assert!((result.i - expected.i).abs() < 1e-12);
    }

    #[test]
    fn zpolyd1_value_matches_zpoly1() {
        let c = [
            Complex::new(0.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(0.4, -0.3),
        ];
        let z = Complex::new(0.2, 0.1);
        let (val, _der) = pj_zpolyd1(z, &c, 3);
        let expected = pj_zpoly1(z, &c, 3);
        assert!((val.r - expected.r).abs() < 1e-12);
        assert!((val.i - expected.i).abs() < 1e-12);
    }

    #[test]
    fn zpolyd1_derivative_matches_numeric_difference() {
        let c = [
            Complex::new(0.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(0.4, -0.3),
        ];
        let z = Complex::new(0.2, 0.1);
        let (_val, der) = pj_zpolyd1(z, &c, 3);

        // Numeric derivative along the real axis: (f(z+h) - f(z-h)) / (2h).
        let h = 1e-6;
        let z_plus = Complex::new(z.r + h, z.i);
        let z_minus = Complex::new(z.r - h, z.i);
        let f_plus = pj_zpoly1(z_plus, &c, 3);
        let f_minus = pj_zpoly1(z_minus, &c, 3);
        let der_num_r = (f_plus.r - f_minus.r) / (2.0 * h);
        let der_num_i = (f_plus.i - f_minus.i) / (2.0 * h);

        assert!((der.r - der_num_r).abs() < 1e-6);
        assert!((der.i - der_num_i).abs() < 1e-6);
    }
}