nettrash-a 0.0.27

A library of algorithms in Rust
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

//! The portable kernels of the greatest common divisor, for Miri and for the targets without an
//! assembly kernel: a product of two limbs is taken as a `u128`, which compilers lower to the
//! target's widening multiplication.

use crate::math::numbers::gcd::limbs::{Matrix1, gcd_11};

/// Sets `(a, b)` to `(u11 a - u01 b, u00 b - u10 a)` over the `n` limbs at `a` and `b`, in place:
/// the pair reduced by the inverse of the matrix `m`, whose entries are below `2^63`. Both results
/// must be non-negative and fit in `n` limbs, as they do when `m` is a reduction of the pair.
///
/// Each position takes the two products of each result in an `i128`, whose magnitude stays below
/// `2^127`: a product of a limb and an entry is below `2^127 - 2^64`, and the carry that the
/// position passes on, the high limb of the difference, is within `2^63 + 1` of zero.
///
/// # Safety
///
/// `n >= 1`; the `n` limbs at `a` and at `b` must be valid for reads and writes, and apart.
#[inline]
pub(in crate::math::numbers) unsafe fn reduce_pair(a: *mut u64, b: *mut u64, n: usize, m: Matrix1) {
    let (mut ca, mut cb) = (0i128, 0i128);
    // SAFETY: as the caller guarantees; each position is read before it is written.
    unsafe {
        for i in 0..n {
            let (x, y) = (u128::from(*a.add(i)), u128::from(*b.add(i)));
            let p = (u128::from(m.u11) * x) as i128 - (u128::from(m.u01) * y) as i128 + ca;
            let q = (u128::from(m.u00) * y) as i128 - (u128::from(m.u10) * x) as i128 + cb;
            *a.add(i) = p as u64;
            *b.add(i) = q as u64;
            ca = p >> 64;
            cb = q >> 64;
        }
    }
    debug_assert!(ca == 0 && cb == 0, "a reduced pair is negative or too long");
}

/// Sets `(x, y)` to `(u00 x + u10 y, u01 x + u11 y)` over the `n` limbs at `x` and `y`, in place:
/// the row `(x, y)` times the matrix `m`, whose entries are below `2^63`. Returns the limbs that
/// carry out of the two results, each below `2^64`: two products of a limb and an entry and a carry
/// add up to less than `2^128`.
///
/// # Safety
///
/// `n >= 1`; the `n` limbs at `x` and at `y` must be valid for reads and writes, and apart.
#[inline]
pub(in crate::math::numbers) unsafe fn extend_row(
    x: *mut u64,
    y: *mut u64,
    n: usize,
    m: Matrix1,
) -> (u64, u64) {
    let (mut cx, mut cy) = (0u64, 0u64);
    // SAFETY: as the caller guarantees; each position is read before it is written.
    unsafe {
        for i in 0..n {
            let (u, v) = (u128::from(*x.add(i)), u128::from(*y.add(i)));
            let p = u128::from(m.u00) * u + u128::from(m.u10) * v + u128::from(cx);
            let q = u128::from(m.u01) * u + u128::from(m.u11) * v + u128::from(cy);
            *x.add(i) = p as u64;
            *y.add(i) = q as u64;
            cx = (p >> 64) as u64;
            cy = (q >> 64) as u64;
        }
    }
    (cx, cy)
}

/// Adds the `n` limbs at `x` times `y` to the `n` limbs at `r`, and returns the limb that carries
/// out of them.
///
/// # Safety
///
/// `n >= 1`; the `n` limbs at `x` must be valid for reads and those at `r` for reads and writes,
/// and apart from them.
#[inline]
pub(in crate::math::numbers) unsafe fn addmul_1(
    r: *mut u64,
    x: *const u64,
    n: usize,
    y: u64,
) -> u64 {
    let mut carry = 0;
    // SAFETY: as the caller guarantees. `x * y + r + carry` is at most `(2^64 - 1)^2 + 2 (2^64 -
    // 1) = 2^128 - 1`, a `u128`.
    unsafe {
        for i in 0..n {
            let p =
                u128::from(*x.add(i)) * u128::from(y) + u128::from(*r.add(i)) + u128::from(carry);
            *r.add(i) = p as u64;
            carry = (p >> 64) as u64;
        }
    }
    carry
}

/// The quotient and the remainder of `n` by `d`, for `n < 8 d` and `d < 2^61`: three conditional
/// subtractions of `4 d`, `2 d` and `d`, each a select that the compiler keeps without a branch, as
/// whether it is taken is as good as random.
#[inline(always)]
pub(in crate::math::numbers) fn small_q1(n: u64, d: u64) -> (u64, u64) {
    let mut r = n;
    let mut q = 0;
    for k in [2, 1, 0] {
        let (t, borrow) = r.overflowing_sub(d << k);
        r = core::hint::select_unpredictable(borrow, r, t);
        q = (q << 1) | u64::from(!borrow);
    }
    (q, r)
}

/// [`small_q1`] on two limbs, for `n < 8 d` and `d < 2^125`.
#[inline(always)]
pub(in crate::math::numbers) fn small_q2(n: u128, d: u128) -> (u64, u128) {
    let mut r = n;
    let mut q = 0;
    for k in [2, 1, 0] {
        let (t, borrow) = r.overflowing_sub(d << k);
        r = core::hint::select_unpredictable(borrow, r, t);
        q = (q << 1) | u64::from(!borrow);
    }
    (q, r)
}

/// The greatest common divisor of the odd `u` and `v` of up to two limbs: [`gcd_11`]'s steps on
/// two limbs while either number has a second one, then on one.
///
/// [`gcd_11`]: crate::math::numbers::gcd::limbs::gcd_11
#[inline]
pub(in crate::math::numbers) fn gcd_22(mut u: u128, mut v: u128) -> u128 {
    debug_assert!(u & v & 1 == 1);
    while (u | v) >> 64 != 0 {
        let (d, below) = u.overflowing_sub(v);
        if d == 0 {
            return u;
        }
        let zeros = d.trailing_zeros();
        if below {
            v = u;
        }
        u = if below { d.wrapping_neg() } else { d } >> zeros;
    }
    u128::from(gcd_11(u as u64, v as u64))
}