dsif/

lib.rs

#[macro_use]
extern crate log;

use traceidr::TraceId;
use std::fmt;

fn u64_sqrt(n: u64) -> u64 {
    if n == 0 || n == 1 {
        return n;
    }

    let mut i = 0;
    let mut result = 0;

    while result <= n {
        i += 1;
        result = i * i;
    }

    return i - 1;
}

#[derive(Copy, Clone, PartialEq, PartialOrd, Ord, Eq, Debug)]
struct U128 {
    hi: u64,
    lo: u64,
}

fn modulo(mut a: U128, m: u64) -> u64 {
    if a.hi >= m {
        a.hi -= (a.hi / m) * m;
    }
    let mut x = a.hi;
    let mut y = a.lo;
    for _ in 0..64 {
        let t = (x as i64 >> 63) as u64;
        x = (x << 1) | (y >> 63);
        y <<= 1;
        if (x | t) >= m {
            x = x.wrapping_sub(m);
            y += 1;
        }
    }
    x
}
fn mul128(u: u64, v: u64) -> U128 {
    let u1 = u >> 32;
    let u0 = u & (!0 >> 32);
    let v1 = v >> 32;
    let v0 = v & (!0 >> 32);

    let t = u0 * v0;
    let w0 = t & (!0 >> 32);
    let k = t >> 32;

    let t = u1 * v0 + k;
    let w1 = t & (!0 >> 32);
    let w2 = t >> 32;

    let t = u0 * v1 + w1;
    let k = t >> 32;
    U128 {
        lo: (t << 32) + w0,
        hi: u1 * v1 + w2 + k,
    }
}
fn mod_mul_(a: u64, b: u64, m: u64) -> u64 {
    modulo(mul128(a, b), m)
}

fn mod_mul(a: u64, b: u64, m: u64) -> u64 {
    match a.checked_mul(b) {
        Some(r) => {
            if r >= m {
                r % m
            } else {
                r
            }
        }
        None => mod_mul_(a, b, m),
    }
}

fn mod_sqr(a: u64, m: u64) -> u64 {
    if a < (1 << 32) {
        let r = a * a;
        if r >= m {
            r % m
        } else {
            r
        }
    } else {
        mod_mul_(a, a, m)
    }
}

fn mod_exp(mut x: u64, mut d: u64, n: u64) -> u64 {
    let mut ret: u64 = 1;
    while d != 0 {
        if d % 2 == 1 {
            ret = mod_mul(ret, x, n)
        }
        d /= 2;
        x = mod_sqr(x, n);
    }
    ret
}

fn miller_rabin(n: u64) -> bool {
    const HINT: &'static [u64] = &[2];

    // we have a strict upper bound, so we can just use the witness
    // table of Pomerance, Selfridge & Wagstaff and Jeaschke to be as
    // efficient as possible, without having to fall back to
    // randomness.
    const WITNESSES: &'static [(u64, &'static [u64])] = &[
        (2_046, HINT),
        (1_373_652, &[2, 3]),
        (9_080_190, &[31, 73]),
        (25_326_000, &[2, 3, 5]),
        (4_759_123_140, &[2, 7, 61]),
        (1_112_004_669_632, &[2, 13, 23, 1662803]),
        (2_152_302_898_746, &[2, 3, 5, 7, 11]),
        (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]),
        (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]),
        (0xFFFF_FFFF_FFFF_FFFF, &[2, 3, 5, 7, 11, 13, 17, 19, 23]),
    ];

    let trace = TraceId::new();
    debug!(
        "Running Miller-Rabin test for the number {} with TraceId {}. This is an expensive task.",
        n, trace
    );

    if n % 2 == 0 {
        return n == 2;
    }
    if n == 1 {
        return false;
    }

    let mut d = n - 1;
    let mut s = 0;
    while d % 2 == 0 {
        d /= 2;
        s += 1
    }

    let witnesses = WITNESSES
        .iter()
        .find(|&&(hi, _)| hi >= n)
        .map(|&(_, wtnss)| wtnss)
        .unwrap();
    'next_witness: for &a in witnesses.iter() {
        let mut power = mod_exp(a, d, n);
        assert!(power < n);
        if power == 1 || power == n - 1 {
            continue 'next_witness;
        }

        for _r in 0..s {
            power = mod_sqr(power, n);
            assert!(power < n);
            if power == 1 {
                return false;
            }
            if power == n - 1 {
                continue 'next_witness;
            }
        }
        debug!(
            "Miller-Rabin test completed for number {} with TraceId {} with a result of false.",
            n, trace
        );
        return false;
    }

    debug!(
        "Miller-Rabin test completed for number {} with TraceId {} with a result of true.",
        n, trace
    );
    true
}

fn is_prime(num: u64) -> bool {
    let trace = TraceId::new();
    debug!(
        "Running primality tests for {} with TraceId {}.",
        num, trace
    );
    let is_prime = miller_rabin(num);
    debug!(
        "Primality test for number {} with TraceId {} produced result {}.",
        num, trace, is_prime
    );
    is_prime
}

pub struct Factors {
    inner: Vec<u64>,
}

impl Factors {
    fn new() -> Self {
        Self {
            inner: Vec::new(),
        }
    }

    fn add(&mut self, f: u64) {
        self.inner.push(f);
    }
}

impl fmt::Display for Factors {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut flag = false;
        for (i, factor) in self.inner.iter().enumerate() {
            if flag {
                write!(f, "* {}", factor)?;
            } else {
                flag = true;
                write!(f, "{}", factor)?;
            }
            if i != self.inner.len() - 1 {
                write!(f, " ")?;
            }
        }
        Ok(())
    }
}

pub fn factorize_to_primes(mut num: u64) -> Factors {
    let trace = TraceId::new();
    debug!("Factorizing number {} with TraceId {}.", num, trace);
    let mut factors = Factors::new();

    if is_prime(num) {
        debug!(
            "Number passed primality test. Returning as is. TraceId {}.",
            trace
        );
        factors.add(num);
        return factors;
    } else {
        debug!(
            "Number did not pass primality test. Continuing. TraceId {}.",
            trace
        );
    }

    while num % 2 == 0 {
        factors.add(2);
        num /= 2;
    }

    let mut i = 3;
    let e = u64_sqrt(num);
    while i <= e {
        while num % i == 0 {
            factors.add(i);
            num /= i;
        }

        i += 2;
    }

    if num > 2 {
        factors.add(num);
    }

    debug!(
        "Finished factorization of number {} with TraceId {}.",
        num, trace
    );

    factors
}