Skip to content

Commit 6d9cf27

Browse files
unageekunageek
authored
Fix derivation and implementation of GCD (#32)
* Fix derivation and implementation of GCD * Workaround for rust-lang/rust#84162 * 💄
1 parent 32fdb3b commit 6d9cf27

File tree

2 files changed

+150
-61
lines changed

2 files changed

+150
-61
lines changed

‎rust-toolchain

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
nightly
1+
nightly-2021-03-25

‎src/interval_set_ops.rs

Lines changed: 149 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ use crate::interval_set::{
44
use gmp_mpfr_sys::mpfr;
55
use inari::{const_dec_interval, const_interval, interval, DecInterval, Decoration, Interval};
66
use rug::Float;
7+
use smallvec::{smallvec, SmallVec};
78
use std::{
89
convert::From,
910
ops::{Add, Mul, Neg, Sub},
@@ -362,36 +363,92 @@ impl TupperIntervalSet {
362363
rs
363364
}
364365

365-
// For x, y ∈ ℚ, gcd(x, y) is defined recursively (the Euclidean algorithm) as:
366+
// For any x, y ∈ ℚ, the GCD (greatest common divisor) of x and y is defined as an extension
367+
// to the integer GCD as:
368+
//
369+
// gcd(x, y) := gcd(q x, q y) / q,
370+
//
371+
// where q is any positive integer that satisfies q x, q y ∈ ℤ (the trivial one is the product
372+
// of the denominators of |x| and |y|). We leave the function undefined for irrational numbers.
373+
// The Euclidean algorithm can be applied to compute the GCD of rational numbers:
366374
//
367375
// gcd(x, y) = | |x| if y = 0,
368376
// | gcd(y, x mod y) otherwise,
369377
//
370-
// assuming the recursion terminates. We leave gcd undefined For irrational numbers.
371-
// Here is an interval extension of it:
378+
// which can be seen from a simple observation:
379+
//
380+
// gcd(q x, q y) / q = | |q x| / q if y = 0,
381+
// | gcd(q y, (q x) mod (q y)) / q otherwise,
382+
// (the Euclidean algorithm for the integer GCD)
383+
// = | |x| if y = 0,
384+
// | gcd(y, ((q x) mod (q y)) / q) otherwise.
385+
// ^^^^^^^^^^^^^^^^^^^^^ = x mod y
386+
//
387+
// We construct an interval extension of the function as follows:
388+
//
389+
// R_0(X, Y) := X,
390+
// R_1(X, Y) := Y,
391+
// R_k(X, Y) := R_{k-2}(X, Y) mod R_{k-1}(X, Y) for any k ∈ ℕ_{≥2},
392+
//
393+
// Z_0(X, Y) := ∅,
394+
// Z_k(X, Y) := | Z_{k-1}(X, Y) ∪ |R_{k-1}(X, Y)| if 0 ∈ R_k(X, Y),
395+
// | Z_{k-1}(X, Y) otherwise,
396+
// for any k ∈ ℕ_{≥1},
397+
//
398+
// gcd(X, Y) := ⋃ Z_k(X, Y).
399+
// k ∈ ℕ_{≥0}
400+
//
401+
// We will denote R_k(X, Y) and Z_k(X, Y) just by R_k and Z_k, respectively.
402+
//
403+
// Proposition. gcd(X, Y) is an interval extension of gcd(x, y).
404+
//
405+
// Proof. Let X and Y be any intervals. There are two possibilities:
406+
//
407+
// (1): X ∩ ℚ = ∅ ∨ Y ∩ ℚ = ∅,
408+
// (2): X ∩ ℚ ≠ ∅ ∧ Y ∩ ℚ ≠ ∅.
409+
//
410+
// Suppose (1). Then, gcd[X, Y] = ∅ ⊆ gcd(X, Y).
411+
// Suppose (2). Let x ∈ X ∩ ℚ, y ∈ Y ∩ ℚ.
412+
// Let r_0 := x, r_1 := y, r_k := r_{k-2} mod r_{k-1} for k ≥ 2.
413+
// Let P(k) :⟺ r_k ∈ R_k. We show that P(k) holds for every k ≥ 0 by induction on k.
414+
// Base cases: From r_0 = x ∈ X = R_0 and r_1 = y ∈ Y = R_1, P(0) and P(1) holds.
415+
// Inductive step: Let k ≥ 0. Suppose P(k), P(k + 1).
416+
// Since X mod Y is an interval extension of x mod y, the following holds:
372417
//
373-
// | |X| if Y = {0} ∨ (X = Y ∧ 0 ∈ X ∧ X mod X ⊆ |X|),
374-
// gcd(X, Y) := | gcd(Y, X mod Y) if 0 ∉ Y,
375-
// | |X| ∪ gcd(Y, X mod Y) otherwise.
418+
// r_{k+2} = r_k mod r_{k+1} ∈ R_k mod R_{k+1} = R_{k+2}.
376419
//
377-
// The second condition of the first case is justified by the following proposition.
378-
// Let X ⊆ ℚ such that 0 ∈ X ∧ X mod X ⊆ |X|. Then
420+
// Thus, P(k + 2) holds. Therefore, P(k) holds for every k ∈ ℕ_{≥0}.
421+
// Since the Euclidean algorithm halts on any input, there exists n ≥ 1 such that
422+
// r_n = 0 ∧ ∀k ∈ {2, …, n-1} : r_k ≠ 0, which leads to |r_{n-1}| = gcd(x, y).
423+
// Let n be such a number. Then from r_n = 0 and r_n ∈ R_n, 0 ∈ R_n. Therefore:
379424
//
380-
// gcd(X, X) ⊆ |X|.
425+
// gcd(x, y) = |r_{n-1}| ∈ |R_{n-1}| ⊆ Z_n ⊆ gcd(X, Y).
381426
//
382-
// Proof: From the definition of gcd,
427+
// Therefore, for any intervals X and Y, gcd[X, Y] ⊆ gcd(X, Y). ■
383428
//
384-
// ∀x, y ∈ X : ∃n ∈ ℕ_{≥1} : ∃z_0, …, z_n ∈ |X| :
385-
// z_0 = |x| ∧ z_1 = |y|
386-
// ∧ ∀k ∈ {2, …, n} : z_k = z_{k-2} mod z_{k-1}
387-
// ∧ z_n = 0 (thus z_{n-1} = gcd(x, y)).
429+
// Proposition. For any intervals X and Y, and any k ≥ 2,
430+
// ∃i ∈ {1, …, k - 1} : R_{i-1} = R_{k-1} ∧ R_i = R_k ∧ Z_{i-1} = Z_{k-1} ⟹ gcd(X, Y) = Z_{k-1}.
431+
// The statement may look a bit awkward, but it makes the implementation easier.
388432
//
389-
// Thus ∀x, y ∈ X : gcd(x, y) ∈ |X|. ■
433+
// Proof. For any j ≥ 1, Z_j can be written in the form:
434+
//
435+
// Z_j = f(Z_{j-1}, R_{j-1}, R_j),
436+
//
437+
// where f is common for every j.
438+
// Suppose ∃i ∈ {1, …, k - 1} : R_{i-1} = R{k-1} ∧ R_i = R_k ∧ Z_{i-1} = Z_{k-1}.
439+
// Let i be such a number. Let n := k - i. Then:
440+
//
441+
// Z_{i+n} = f(Z_{i+n-1}, R_{i+n-1}, R_{i+n})
442+
// = f(Z_{i-1}, R_{i-1}, R_i)
443+
// = Z_i.
444+
//
445+
// By repeating the process, we get ∀m ∈ ℕ_{≥0} : Z_{i+mn} = Z_i.
446+
// Therefore, ∀j ∈ ℕ_{≥i} : Z_j = Z_i.
447+
// Therefore, gcd(X, Y) = Z_i = Z_{k-1}. ■
390448
pub fn gcd(&self, rhs: &Self, site: Option<Site>) -> Self {
391-
const ZERO: Interval = const_interval!(0.0, 0.0);
392449
let mut rs = Self::new();
393-
// gcd(X, Y)
394-
// = gcd(|X|, |Y|)
450+
// {gcd(x, y) | x ∈ X, y ∈ Y}
451+
// = {gcd(x, y) | x ∈ |X|, y ∈ |Y|}
395452
// = {gcd(max(x, y), min(x, y)) | x ∈ |X|, y ∈ |Y|}
396453
// ⊆ {gcd(x, y) | x ∈ max(|X|, |Y|), y ∈ min(|X|, |Y|)}.
397454
let xs = &self.abs();
@@ -406,40 +463,46 @@ impl TupperIntervalSet {
406463
};
407464
let x = DecInterval::set_dec(x.x, dec);
408465
let y = DecInterval::set_dec(y.x, dec);
409-
let mut xs = Self::from(TupperInterval::new(x.max(y), g));
410-
let mut ys = Self::from(TupperInterval::new(x.min(y), g));
411-
let mut prev_xs = xs.clone();
412-
loop {
413-
if ys.iter().any(|y| y.x.contains(0.0)) {
414-
rs.extend(&xs);
415-
416-
if xs == ys || prev_xs == ys {
417-
// Here, in the first iteration, `xs` and `ys` consists of the same,
418-
// single, nonnegative interval that contains zero:
419-
//
420-
// X = Y = |X| = [0, a],
421-
//
422-
// where a ≥ 0. Let x, y ∈ X. Then:
423-
//
424-
// 0 ≤ x mod y < y ≤ a
425-
// ⟹ x mod y ∈ X.
426-
//
427-
// Therefore, 0 ∈ X ∧ X mod X ⊆ |X|.
428-
break;
466+
let mut zs = TupperIntervalSet::new();
467+
let mut zs_prev = zs.clone();
468+
let mut rems: SmallVec<[_; 4]> = smallvec![
469+
Self::from(TupperInterval::new(x.max(y), g)),
470+
Self::from(TupperInterval::new(x.min(y), g))
471+
];
472+
'outer: loop {
473+
// The iteration starts with k = 1.
474+
let xs = &rems[rems.len() - 2];
475+
let ys = &rems[rems.len() - 1];
476+
// R_{k-1} = `xs`, R_k = `ys`, Z_{k-1} = `zs_prev`.
477+
for i_prime in 0..rems.len() - 2 {
478+
// `i_prime` is i with some offset subtracted.
479+
// We have (1): 1 ≤ i < k, (2): Z_{i-1} = Z_i = … = Z_{k-1}.
480+
if &rems[i_prime] == xs && &rems[i_prime + 1] == ys {
481+
// We have R_{i-1} = R_{k-1} ∧ R_i = R_k.
482+
// Therefore, gcd(X, Y) = Z_{k-1}.
483+
break 'outer;
429484
}
430485
}
431486

432-
ys.retain(|y| y.x != ZERO);
433-
if ys.is_empty() {
434-
break;
435-
}
487+
// (used later) R_{k+1} = `rem`.
488+
let mut rem = xs.rem_euclid(&ys, None);
489+
rem.normalize(true);
436490

437-
let xs_rem_ys = xs.rem_euclid(&ys, None);
438-
prev_xs = xs;
439-
xs = ys;
440-
ys = xs_rem_ys;
441-
ys.normalize(true); // To compare it with `xs` in the next iteration.
491+
if ys.iter().any(|y| y.x.contains(0.0)) {
492+
zs.extend(xs);
493+
zs.normalize(true);
494+
// Z_k = `zs`.
495+
if zs != zs_prev {
496+
// Z_k ≠ Z_{k-1}.
497+
// Retain only R_k so that both (1) and (2) will hold
498+
// in subsequent iterations.
499+
rems = rems[rems.len() - 1..].into();
500+
zs_prev = zs.clone();
501+
}
502+
}
503+
rems.push(rem); // […, R_k, R_{k+1}]
442504
}
505+
rs.extend(zs_prev);
443506
}
444507
}
445508
}
@@ -462,16 +525,38 @@ impl TupperIntervalSet {
462525
rs
463526
}
464527

465-
// For x, y ∈ ℚ, lcm(x, y) is defined as:
528+
// For x, y ∈ ℚ, the LCM (least common multiple) of x and y is defined as:
466529
//
467530
// lcm(x, y) = | 0 if x = y = 0,
468531
// | |x y| / gcd(x, y) otherwise.
469532
//
470-
// We leave lcm undefined for irrational numbers.
471-
// Here is an interval extension of it:
533+
// We leave the function undefined for irrational numbers.
534+
// Here is an interval extension of the function:
472535
//
473536
// lcm(X, Y) := | {0} if X = Y = {0},
474537
// | |X Y| / gcd(X, Y) otherwise.
538+
//
539+
// Proposition. lcm(X, Y) is an interval extension of lcm(x, y).
540+
//
541+
// Proof. Let X and Y be any intervals. There are five possibilities:
542+
//
543+
// (1): X ∩ ℚ = ∅ ∨ Y ∩ ℚ = ∅,
544+
// (2): X = Y = {0},
545+
// (3): X = {0} ∧ Y ∩ ℚ\{0} ≠ ∅,
546+
// (4): X ∩ ℚ\{0} ≠ ∅ ∧ Y = {0},
547+
// (5): X ∩ ℚ\{0} ≠ ∅ ∧ Y ∩ ℚ\{0} ≠ ∅.
548+
//
549+
// Suppose (1). Then, lcm[X, Y] = ∅ ⊆ lcm(X, Y).
550+
// Suppose (2). Then, lcm[X, Y] = lcm(X, Y) = {0}.
551+
// Suppose (3). As Y ≠ {0}, lcm(X, Y) = |X Y| / gcd(X, Y).
552+
// Therefore, from 0 ∈ |X Y| and ∃y ∈ Y ∩ ℚ\{0} : |y| ∈ gcd(X, Y), 0 ∈ lcm(X, Y).
553+
// Therefore, lcm[X, Y] = {0} ⊆ lcm(X, Y).
554+
// Suppose (4). In the same manner, lcm[X, Y] ⊆ lcm(X, Y).
555+
// Suppose (5). Let x ∈ X ∩ ℚ\{0}, y ∈ Y ∩ ℚ\{0} ≠ ∅.
556+
// Then, |x y| / gcd(x, y) ∈ lcm(X, Y) = |X Y| / gcd(X, Y).
557+
// Therefore, lcm[X, Y] ⊆ lcm(X, Y).
558+
//
559+
// Hence, the result. ■
475560
pub fn lcm(&self, rhs: &Self, site: Option<Site>) -> Self {
476561
const ZERO: Interval = const_interval!(0.0, 0.0);
477562
let mut rs = TupperIntervalSet::new();
@@ -1551,14 +1636,18 @@ mod tests {
15511636
x.gcd(&y, None)
15521637
}
15531638

1554-
test!(@commut f, @even i!(15.0), @even i!(30.0), (vec![i!(15.0)], Dac));
1555-
test!(@commut f, @even i!(15.0), @even i!(21.0), (vec![i!(3.0)], Dac));
1556-
test!(@commut f, @even i!(15.0), @even i!(17.0), (vec![i!(1.0)], Dac));
1557-
test!(@commut f, @even i!(7.5), @even i!(10.5), (vec![i!(1.5)], Dac));
1558-
test!(@commut f, i!(0.0), @even i!(5.0), (vec![i!(5.0)], Dac));
15591639
test!(f, i!(0.0), i!(0.0), (vec![i!(0.0)], Dac));
1560-
1561-
// This test fails into an infinite loop if you remove `prev_xs == ys`.
1640+
test!(@commut f, i!(0.0), @even i!(5.0), (vec![i!(5.0)], Dac));
1641+
test!(@commut f, @even i!(7.5), @even i!(10.5), (vec![i!(1.5)], Dac));
1642+
test!(@commut f, @even i!(15.0), @even i!(17.0), (vec![i!(1.0)], Dac));
1643+
test!(@commut f, @even i!(15.0), @even i!(21.0), (vec![i!(3.0)], Dac));
1644+
test!(@commut f, @even i!(15.0), @even i!(30.0), (vec![i!(15.0)], Dac));
1645+
test!(
1646+
@commut f,
1647+
@even i!(1348500621.0),
1648+
@even i!(18272779829.0),
1649+
(vec![i!(150991.0)], Dac)
1650+
);
15621651
test!(
15631652
@commut f,
15641653
@even i!(4.0, 6.0),
@@ -1573,10 +1662,10 @@ mod tests {
15731662
x.lcm(&y, None)
15741663
}
15751664

1576-
test!(@commut f, @even i!(3.0), @even i!(5.0), (vec![i!(15.0)], Dac));
1577-
test!(@commut f, @even i!(1.5), @even i!(2.5), (vec![i!(7.5)], Dac));
1578-
test!(@commut f, i!(0.0), @even i!(5.0), (vec![i!(0.0)], Dac));
15791665
test!(f, i!(0.0), i!(0.0), (vec![i!(0.0)], Dac));
1666+
test!(@commut f, i!(0.0), @even i!(5.0), (vec![i!(0.0)], Dac));
1667+
test!(@commut f, @even i!(1.5), @even i!(2.5), (vec![i!(7.5)], Dac));
1668+
test!(@commut f, @even i!(3.0), @even i!(5.0), (vec![i!(15.0)], Dac));
15801669
}
15811670

15821671
#[test]

0 commit comments

Comments
 (0)