maph/num/
factors.rs

1const FAC_ITER: usize = 100;
2
3pub fn rat_approx_32(f: f32, max: f32) -> (f32, f32) {
4    let frac = f.fract();
5    let int = f - frac;
6    if int > max || int < -max { panic!("Tried to convert an f32 out of bounds for r32, overflowed."); }
7    if frac <= 0.0 {
8        return (int, 1.0);
9    }
10    let dmax = max/int;
11    assert!(0.0 < frac && frac < 1.0);
12    let (mut nl, mut dl) : (f32, f32) = (0.0, 1.0);
13    let (mut nr, mut dr) : (f32, f32) = (1.0, 1.0);
14    let (mut nm, mut dm) : (i64, i64) = (0, 1);
15    let mut loc = None;
16    while (dl <= dmax) && (dr <= dmax) {
17        let ntmp = nl - frac * dl;
18        let dtmp = frac * dr - nr;
19        let mut br = (ntmp/dtmp).floor() as u32;
20        let mut bl = (dtmp/ntmp).floor() as u32;
21        let mut side = 0;
22        if bl == 0 { bl = 1; side = -1; }
23        if br == 0 { br = 1; side = 1; }
24        dm = dl as i64 * bl as i64 + dr as i64 * br as i64;
25        
26        if dm > dmax as i64 {
27            if side == -1 { br = 1.max(((dmax - dl) / dr).floor() as u32); }
28            else if side == 1 { bl = 1.max(((dmax - dr) / dl).floor() as u32); }
29        }
30
31        nm = nl as i64 * bl as i64 + nr as i64 * br as i64;
32        dm = dl as i64 * bl as i64 + dr as i64 * br as i64;
33        if dm > dmax as i64 { break; }
34        let med = (nm as f32)/(dm as f32);
35
36        if frac == med {
37            loc = Some(0);
38            break;
39        } else if frac < med {
40            loc = Some(-1);
41            (nr, dr) = (nm as f32, dm as f32);
42        } else if frac > med {
43            loc = Some(1);
44            (nl, dl) = (nm as f32, dm as f32);
45        }
46    }
47
48    let (mut n, d);
49
50    match loc {
51        Some(0) => {
52            (n, d) = (nm, dm);
53        },
54        _ => {
55            let errl = (frac - (nl as f32/dl as f32)).abs();
56            let errr = (frac - (nr as f32 / dr as f32)).abs();
57            if errl <= errr {
58                (n, d) = (nl as i64, dl as i64);
59            } else {
60                (n, d) = (nr as i64, dr as i64);
61            }
62        },
63    }
64
65
66
67    n += int as i64 * d as i64;
68    (n as f32, d as f32)
69}
70
71pub fn rat_approx_64(f: f64, max: f64) -> (f64, f64) {
72    let frac = f.fract();
73    let int = f - frac;
74    if int > max || int < -max { panic!("Tried to convert an f32 out of bounds for r32, overflowed."); }
75    if frac <= 0.0 {
76        return (int, 1.0);
77    }
78    let dmax = max/int;
79    assert!(0.0 < frac && frac < 1.0);
80    let (mut nl, mut dl) : (f64, f64) = (0.0, 1.0);
81    let (mut nr, mut dr) : (f64, f64) = (1.0, 1.0);
82    let (mut nm, mut dm) : (i128, i128) = (0, 1);
83    let mut loc = None;
84    while (dl <= dmax) && (dr <= dmax) {
85        let ntmp = nl - frac * dl;
86        let dtmp = frac * dr - nr;
87        let mut br = (ntmp/dtmp).floor() as u32;
88        let mut bl = (dtmp/ntmp).floor() as u32;
89        let mut side = 0;
90        if bl == 0 { bl = 1; side = -1; }
91        if br == 0 { br = 1; side = 1; }
92        dm = dl as i128 * bl as i128 + dr as i128 * br as i128;
93        
94        if dm > dmax as i128 {
95            if side == -1 { br = 1.max(((dmax - dl) / dr).floor() as u32); }
96            else if side == 1 { bl = 1.max(((dmax - dr) / dl).floor() as u32); }
97        }
98
99        nm = nl as i128 * bl as i128 + nr as i128 * br as i128;
100        dm = dl as i128 * bl as i128 + dr as i128 * br as i128;
101        if dm > dmax as i128 { break; }
102        let med = (nm as f64)/(dm as f64);
103
104        if frac == med {
105            loc = Some(0);
106            break;
107        } else if frac < med {
108            loc = Some(-1);
109            (nr, dr) = (nm as f64, dm as f64);
110        } else if frac > med {
111            loc = Some(1);
112            (nl, dl) = (nm as f64, dm as f64);
113        }
114    }
115
116    let (mut n, d);
117
118    match loc {
119        Some(0) => {
120            (n, d) = (nm, dm);
121        },
122        _ => {
123            let errl = (frac - (nl as f64/dl as f64)).abs();
124            let errr = (frac - (nr as f64 / dr as f64)).abs();
125            if errl <= errr {
126                (n, d) = (nl as i128, dl as i128);
127            } else {
128                (n, d) = (nr as i128, dr as i128);
129            }
130        },
131    }
132
133
134
135    n += int as i128 * d as i128;
136    (n as f64, d as f64)
137}
138
139///Returns the greatest common denominator of two u16s.
140pub fn gcd16(a: u16, b: u16) -> u16 {
141    let (mut a0, mut b0) = (a, b);
142    while a0 != 0 {
143        let a1 = b0 % a0;
144        b0 = a0;
145        a0 = a1;    
146    }
147    return b0;
148}
149///Returns the greatest common denominator of two u32s.
150pub fn gcd32(a: u32, b: u32) -> u32 {
151    let (mut a0, mut b0) = (a, b);
152    while a0 != 0 {
153        let a1 = b0 % a0;
154        b0 = a0;
155        a0 = a1;    
156    }
157    return b0;
158}
159///Returns the greatest common denominator of two u64s.
160pub fn gcd64(a: u64, b: u64) -> u64 {
161    let (mut a0, mut b0) = (a, b);
162    while a0 != 0 {
163        let a1 = b0 % a0;
164        b0 = a0;
165        a0 = a1;    
166    }
167    return b0;
168}
169///Returns the lowest common multiple of two u32s.
170pub fn lcm16(a: u16, b: u16) -> u16 { (a * b) / gcd16(a, b) }
171///Returns the lowest common multiple of two u32s.
172pub fn lcm32(a: u32, b: u32) -> u32 { (a * b) / gcd32(a, b) }
173///Returns the lowest common multiple of two u64s.
174pub fn lcm64(a: u64, b: u64) -> u64 { (a * b) / gcd64(a, b) }
175///Returns a list of prime factors for a u32.
176pub fn fac32(n: u32) -> Vec<u32> {
177    let mut v = vec!(n);
178    let mut running = true;
179    while running {
180        running = false;
181        let mut v0 = Vec::new();
182        for n0 in &v {
183            match prho32(*n0) {
184                (0, 0) => { v0.push(*n0); },
185                (_, 1) => { v0.push(*n0); },
186                (1, _) => { v0.push(*n0); },
187                (a, b) => { v0.push(a); v0.push(b); running = true; }
188            }
189        }
190        v = v0;
191    }
192    v
193}
194///Pollard Rho algorithm for a u32 - returns a non-trivial factor of a u32 and the original number
195///with the factor removed.
196pub fn prho32(n: u32) -> (u32, u32) {
197    let (mut x, mut y, mut d);
198    let (mut count, mut attempt) = (0, 0);
199    while count < FAC_ITER {
200        count += 1;
201        if n % 2 == 0 { return (2, n/2); }
202        else if n % 3 == 0 { return (3, n/3); }
203        else if n % 5 == 0 { return (5, n/5); }
204        else {
205            x = 2_i32;
206            y = 2_i32;
207            d = 1_i32;
208            while d == 1 {
209                x = (x*x + attempt) % n as i32;
210                y = (y*y + attempt) % n as i32;
211                y = (y*y + attempt) % n as i32;
212                d = gcd32((x - y).abs() as u32, n) as i32;
213            }
214            if d != n as i32 { return (d as u32, n/d as u32); }
215            else { attempt += 1; }
216        }
217    }
218    (0, 0)
219}
220///Returns a list of prime factors for a u64.
221pub fn fac64(n: u64) -> Vec<u64> {
222    let mut v = vec!(n);
223    let mut running = true;
224    while running {
225        running = false;
226        let mut v0 = Vec::new();
227        for n0 in &v {
228            match prho64(*n0) {
229                (0, 0) => { v0.push(*n0); },
230                (_, 1) => { v0.push(*n0); },
231                (1, _) => { v0.push(*n0); },
232                (a, b) => { v0.push(a); v0.push(b); running = true; }
233            }
234        }
235        v = v0;
236    }
237    v
238}
239///Pollard Rho algorithm for a u64 - returns a non-trivial factor of a u64 and the original number
240///with the factor removed.
241pub fn prho64(n: u64) -> (u64, u64) {
242    let (mut x, mut y, mut d);
243    let (mut count, mut attempt) = (0, 0);
244    while count < FAC_ITER {
245        count += 1;
246        if n % 2 == 0 { return (2, n/2); }
247        else if n % 3 == 0 { return (3, n/3); }
248        else if n % 5 == 0 { return (5, n/5); }
249        else {
250            x = 2_i64;
251            y = 2_i64;
252            d = 1_i64;
253            while d == 1 {
254                x = (x*x + attempt) % n as i64;
255                y = (y*y + attempt) % n as i64;
256                y = (y*y + attempt) % n as i64;
257                d = gcd64((x - y).abs() as u64, n) as i64;
258            }
259            if d != n as i64 { return (d as u64, n/d as u64); }
260            else { attempt += 1; }
261        }
262    }
263    (0, 0)
264}
265
266//Returns a vector of square factors for a u32 - each item is a tuple of the unsquared factor and the squared factor.
267pub fn sqfac32(n: u32) -> Vec<(u32, u32)> {
268    let f = fac32(n);
269    let mut v = Vec::new();
270    for fac in f {
271        let mut found = false;
272        for (fac0, n) in &mut v {
273            if fac == *fac0 { *n += 1; found = true; }
274        }
275        if !found { v.push((fac, 1)); }
276    }
277    v.iter().filter(|(_f, n)| n % 2 == 0).map(|(f, n)| (f.pow(*n/2), f.pow(*n))).collect::<Vec<_>>()
278}
279
280//Returns a vector of square factors for a u64 - each item is a tuple of the unsquared factor and the squared factor.
281pub fn sqfac64(n: u64) -> Vec<(u64, u64)> {
282    let f = fac64(n);
283    let mut v = Vec::new();
284    for fac in f {
285        let mut found = false;
286        for (fac0, n) in &mut v {
287            if fac == *fac0 { *n += 1; found = true; }
288        }
289        if !found { v.push((fac, 1)); }
290    }
291    v.iter().filter(|(_f, n)| n % 2 == 0).map(|(f, n)| (f.pow(*n/2), f.pow(*n))).collect::<Vec<_>>()
292}
maph/num/factors.rs

maph/num/
factors.rs
