August 13, 2017 13:46
diff --git a/gistfile1.txt b/gistfile1.txt
 // Compile with: rustc -C opt-level=3 -C target-cpu=native -C panic=abort

 #![feature(inclusive_range_syntax, placement_in_syntax, i128_type, const_fn, iterator_step_by)]

 fn main() {
    use std::mem::size_of;
    use std::ops::RemAssign;

    fn gcd<T>(mut u: T, mut v: T) -> T
    where T: PartialEq + RemAssign + Copy + Clone + std::iter::Sum {
        let zero: T = std::iter::empty().sum();

        if v == zero { return u; }
        loop {
            u %= v;
            if u == zero { return v; }
            v %= u;
            if v == zero { return u; }
        }
    }


    struct BitVec { bits: Vec<usize>, len: usize }

    impl BitVec {
        const BPW: usize = size_of::<usize>() * 8;

        fn new(n_bits: usize, init_value: bool) -> Self {
            let val = if init_value { std::usize::MAX } else { 0 };
            let vec_len = n_bits / Self::BPW + if n_bits % Self::BPW == 0 { 0 } else { 1 };
            Self { bits: vec![val; vec_len],
                   len: n_bits }
        }

        fn is_set(&self, i: usize) -> bool {
            debug_assert!(i < self.len);
            let offset = i / Self::BPW;
            let mask = 1 << (i % Self::BPW);
            unsafe { (*self.bits.get_unchecked(offset) & mask) != 0 }
        }

        fn set(&mut self, i: usize) {
            debug_assert!(i < self.len);
            let offset = i / Self::BPW;
            let mask = 1 << (i % Self::BPW);
            unsafe { *self.bits.get_unchecked_mut(offset) |= mask; }
        }
    }


    struct Xorshift128 { state: [u32; 4] }

    impl Xorshift128 {
        fn new(mut x0: u32) -> Self {
            let mut state = [0; 4];

            for i in 0 .. state.len() {
                x0 = 1_812_433_253u32.wrapping_mul(x0 ^ (x0 >> 30)) +
                     i as u32 + 1;
                state[i] = x0;
            }

            // All state must not be 0.
            for i in 0 .. state.len() {
                if state[i] == 0 {
                    state[i] = i as u32 + 1;
                }
            }

            let mut x = Self { state };
            x.next();
            x
        }

        fn next(&mut self) {
            const A: u32 = 11;
            const B: u32 = 8;
            const C: u32 = 19;
            let temp = self.state[0] ^ (self.state[0] << A);
            self.state[0] = self.state[1];
            self.state[1] = self.state[2];
            self.state[2] = self.state[3];
            self.state[3] = self.state[3] ^ (self.state[3] >> C)
                            ^ temp ^ (temp >> B);
        }

        fn uniform_u64(&mut self, sup: u64) -> u64 {
            // This is acceptable only if sup is small.
            let part1 = u64::from(self.state[self.state.len() - 1]);
            self.next();
            let part2 = u64::from(self.state[self.state.len() - 1]);
            self.next();
            (part1 | (part2 << 32)) % sup
        }
    }


    fn e216() -> u32 {
        fn mul_mod(a: u64, b: u64, p: u64) -> u64 {
            (a * b) % p
        }

        fn pow_mod(mut a: u64, mut k: u64, p: u64) -> u64 {
            let mut res = 1;

            while k > 0 {
                if (k & 1) == 1 {
                    res = (res * a) % p;
                }

                a = (a * a) % p;
                k >>= 1;
            }

            res
        }

        fn jacobi(mut a: u64, mut m: u64) -> i32 {
            let mut t = 1;
            a %= m;

            while a != 0 {
                while (a & 1) == 0 {
                    a >>= 1;

                    if (m & 7) == 3 || (m & 7) == 5 {
                        t = -t;
                    }
                }

                std::mem::swap(&mut a, &mut m);

                if (a & 3) == 3 && (m & 3) == 3 {
                    t = -t;
                }

                a %= m;
            }

            if m == 1 { t } else { 0 }
        }

        fn sqrt_mod(mut a: u64, p: u64, rng: &mut Xorshift128) -> u64 {
            if p == 2 { return a & 1; }
            a %= p;
            if jacobi(a, p) != 1 { return 0; }

            let p8 = (p & 7) as u32;

            if p8 == 3 || p8 == 5 || p8 == 7 {
                if (p8 & 3) == 3 {
                    return pow_mod(a, (p + 1) / 4, p);
                }

                let x = pow_mod(a, (p + 3) / 8, p);
                let c = mul_mod(x, x, p);
                return if c == a { x } else { mul_mod(x, pow_mod(2, (p - 1) / 4, p), p) };
            }

            let mut alpha: u32 = 0;
            let mut s: u64 = p - 1;

            while (s & 1) == 0 {
                s >>= 1;
                alpha += 1;
            }

            let r = pow_mod(a, (s + 1) / 2, p);
            let r2a = mul_mod(r, pow_mod(a, (s + 1) / 2 - 1, p), p);

            let n = loop {
                let n = 2 + rng.uniform_u64(p - 2);
                if jacobi(n, p) == -1 { break n; }
            };

            let b = pow_mod(n, s, p);
            let mut j1 = 0;

            for i in 0 .. alpha - 1 {
                let mut c = pow_mod(b, 2 * j1, p);
                c = mul_mod(r2a, c, p);
                c = pow_mod(c, 1 << (alpha - i - 2), p);

                if c == p - 1 {
                    j1 += 1 << i;
                }
            }

            mul_mod(r, pow_mod(b, j1, p), p)
        }

        const N: u32 = 50_000_000;
        const NL: u64 = N as u64;

        let mut rng = Xorshift128::new(10);
        let max: u32 = ((2 * NL * NL + 1) as f64).sqrt() as u32;
        let lim: u32 = f64::from(max).sqrt() as u32;

        let mut sieve = BitVec::new(max as usize + 1, false);

        for i in 2 ... lim {
            if sieve.is_set(i as usize) { continue; }

            for n in (2 * i ... max).step_by(i as usize) {
                sieve.set(n as usize);
            }
        }

        let mut sieve2 = BitVec::new(N as usize + 1, false);

        for p in 3 ... max {
            if sieve.is_set(p as usize) { continue; }

            let mod8 = p % 8;
            if mod8 == 3 || mod8 == 5 { continue; }

            let n1 = sqrt_mod((u64::from(p) + 1) / 2, u64::from(p), &mut rng);

            if n1 == 0 { continue; }

            let n2 = u64::from(p) - n1;

            for x in (n1 ... N as u64).step_by(p as usize) {
                if 2 * x * x - 1 == u64::from(p) { continue; }
                sieve2.set(x as usize);
            }

            for x in (n2 ... N as u64).step_by(p as usize) {
                if 2 * x * x - 1 == u64::from(p) { continue; }
                sieve2.set(x as usize);
            }

        }

        let mut result = 0;
        for n in 2 ... N {
            if !sieve2.is_set(n as usize) {
                result += 1;
            }
        }

        result
    }
    assert_eq!(e216(), 5_437_849);


    fn e251() -> u32 {
        const LIM: i32 = 110_000_000; // Input.
        const LIMK: i32 = LIM * 2 / 3;
        let mut count = 0;

        for k in (5 ... LIMK).step_by(8) {
            let limy = f64::from(8 * LIM / 3 / k).sqrt() as i32;
            for y in (1 .. limy).step_by(2) {
                let a = (1 + 3 * y * y * k) / 8;
                let top = (3 + y * y * k) / 8;
                let lowx = (1.0 + top as f32 * y as f32 / LIM as f32) as i32;
                for x in lowx .. LIM {
                    let c = k * x * x;
                    if a + c > LIM { break; }
                    if top % x != 0 || gcd(x, y) > 1 { continue; }
                    let b = (top / x) * y;
                    if a + b + c > LIM { continue; }
                    count += 1;
                }
            }
        }
        count
    }
    assert_eq!(e251(), 18_946_051);
 } // Total run-time: 4.40 seconds.
	// Compile with: rustc -C opt-level=3 -C target-cpu=native -C panic=abort

	#![feature(inclusive_range_syntax, placement_in_syntax, i128_type, const_fn, iterator_step_by)]

	fn main() {
	use std::mem::size_of;
	use std::ops::RemAssign;

	fn gcd<T>(mut u: T, mut v: T) -> T
	where T: PartialEq + RemAssign + Copy + Clone + std::iter::Sum {
	let zero: T = std::iter::empty().sum();

	if v == zero { return u; }
	loop {
	u %= v;
	if u == zero { return v; }
	v %= u;
	if v == zero { return u; }
	}
	}


	struct BitVec { bits: Vec<usize>, len: usize }

	impl BitVec {
	const BPW: usize = size_of::<usize>() * 8;

	fn new(n_bits: usize, init_value: bool) -> Self {
	let val = if init_value { std::usize::MAX } else { 0 };
	let vec_len = n_bits / Self::BPW + if n_bits % Self::BPW == 0 { 0 } else { 1 };
	Self { bits: vec![val; vec_len],
	len: n_bits }
	}

	fn is_set(&self, i: usize) -> bool {
	debug_assert!(i < self.len);
	let offset = i / Self::BPW;
	let mask = 1 << (i % Self::BPW);
	unsafe { (*self.bits.get_unchecked(offset) & mask) != 0 }
	}

	fn set(&mut self, i: usize) {
	debug_assert!(i < self.len);
	let offset = i / Self::BPW;
	let mask = 1 << (i % Self::BPW);
	unsafe { *self.bits.get_unchecked_mut(offset) \|= mask; }
	}
	}


	struct Xorshift128 { state: [u32; 4] }

	impl Xorshift128 {
	fn new(mut x0: u32) -> Self {
	let mut state = [0; 4];

	for i in 0 .. state.len() {
	x0 = 1_812_433_253u32.wrapping_mul(x0 ^ (x0 >> 30)) +
	i as u32 + 1;
	state[i] = x0;
	}

	// All state must not be 0.
	for i in 0 .. state.len() {
	if state[i] == 0 {
	state[i] = i as u32 + 1;
	}
	}

	let mut x = Self { state };
	x.next();
	x
	}

	fn next(&mut self) {
	const A: u32 = 11;
	const B: u32 = 8;
	const C: u32 = 19;
	let temp = self.state[0] ^ (self.state[0] << A);
	self.state[0] = self.state[1];
	self.state[1] = self.state[2];
	self.state[2] = self.state[3];
	self.state[3] = self.state[3] ^ (self.state[3] >> C)
	^ temp ^ (temp >> B);
	}

	fn uniform_u64(&mut self, sup: u64) -> u64 {
	// This is acceptable only if sup is small.
	let part1 = u64::from(self.state[self.state.len() - 1]);
	self.next();
	let part2 = u64::from(self.state[self.state.len() - 1]);
	self.next();
	(part1 \| (part2 << 32)) % sup
	}
	}


	fn e216() -> u32 {
	fn mul_mod(a: u64, b: u64, p: u64) -> u64 {
	(a * b) % p
	}

	fn pow_mod(mut a: u64, mut k: u64, p: u64) -> u64 {
	let mut res = 1;

	while k > 0 {
	if (k & 1) == 1 {
	res = (res * a) % p;
	}

	a = (a * a) % p;
	k >>= 1;
	}

	res
	}

	fn jacobi(mut a: u64, mut m: u64) -> i32 {
	let mut t = 1;
	a %= m;

	while a != 0 {
	while (a & 1) == 0 {
	a >>= 1;

	if (m & 7) == 3 \|\| (m & 7) == 5 {
	t = -t;
	}
	}

	std::mem::swap(&mut a, &mut m);

	if (a & 3) == 3 && (m & 3) == 3 {
	t = -t;
	}

	a %= m;
	}

	if m == 1 { t } else { 0 }
	}

	fn sqrt_mod(mut a: u64, p: u64, rng: &mut Xorshift128) -> u64 {
	if p == 2 { return a & 1; }
	a %= p;
	if jacobi(a, p) != 1 { return 0; }

	let p8 = (p & 7) as u32;

	if p8 == 3 \|\| p8 == 5 \|\| p8 == 7 {
	if (p8 & 3) == 3 {
	return pow_mod(a, (p + 1) / 4, p);
	}

	let x = pow_mod(a, (p + 3) / 8, p);
	let c = mul_mod(x, x, p);
	return if c == a { x } else { mul_mod(x, pow_mod(2, (p - 1) / 4, p), p) };
	}

	let mut alpha: u32 = 0;
	let mut s: u64 = p - 1;

	while (s & 1) == 0 {
	s >>= 1;
	alpha += 1;
	}

	let r = pow_mod(a, (s + 1) / 2, p);
	let r2a = mul_mod(r, pow_mod(a, (s + 1) / 2 - 1, p), p);

	let n = loop {
	let n = 2 + rng.uniform_u64(p - 2);
	if jacobi(n, p) == -1 { break n; }
	};

	let b = pow_mod(n, s, p);
	let mut j1 = 0;

	for i in 0 .. alpha - 1 {
	let mut c = pow_mod(b, 2 * j1, p);
	c = mul_mod(r2a, c, p);
	c = pow_mod(c, 1 << (alpha - i - 2), p);

	if c == p - 1 {
	j1 += 1 << i;
	}
	}

	mul_mod(r, pow_mod(b, j1, p), p)
	}

	const N: u32 = 50_000_000;
	const NL: u64 = N as u64;

	let mut rng = Xorshift128::new(10);
	let max: u32 = ((2 * NL * NL + 1) as f64).sqrt() as u32;
	let lim: u32 = f64::from(max).sqrt() as u32;

	let mut sieve = BitVec::new(max as usize + 1, false);

	for i in 2 ... lim {
	if sieve.is_set(i as usize) { continue; }

	for n in (2 * i ... max).step_by(i as usize) {
	sieve.set(n as usize);
	}
	}

	let mut sieve2 = BitVec::new(N as usize + 1, false);

	for p in 3 ... max {
	if sieve.is_set(p as usize) { continue; }

	let mod8 = p % 8;
	if mod8 == 3 \|\| mod8 == 5 { continue; }

	let n1 = sqrt_mod((u64::from(p) + 1) / 2, u64::from(p), &mut rng);

	if n1 == 0 { continue; }

	let n2 = u64::from(p) - n1;

	for x in (n1 ... N as u64).step_by(p as usize) {
	if 2 * x * x - 1 == u64::from(p) { continue; }
	sieve2.set(x as usize);
	}

	for x in (n2 ... N as u64).step_by(p as usize) {
	if 2 * x * x - 1 == u64::from(p) { continue; }
	sieve2.set(x as usize);
	}

	}

	let mut result = 0;
	for n in 2 ... N {
	if !sieve2.is_set(n as usize) {
	result += 1;
	}
	}

	result
	}
	assert_eq!(e216(), 5_437_849);


	fn e251() -> u32 {
	const LIM: i32 = 110_000_000; // Input.
	const LIMK: i32 = LIM * 2 / 3;
	let mut count = 0;

	for k in (5 ... LIMK).step_by(8) {
	let limy = f64::from(8 * LIM / 3 / k).sqrt() as i32;
	for y in (1 .. limy).step_by(2) {
	let a = (1 + 3 * y * y * k) / 8;
	let top = (3 + y * y * k) / 8;
	let lowx = (1.0 + top as f32 * y as f32 / LIM as f32) as i32;
	for x in lowx .. LIM {
	let c = k * x * x;
	if a + c > LIM { break; }
	if top % x != 0 \|\| gcd(x, y) > 1 { continue; }
	let b = (top / x) * y;
	if a + b + c > LIM { continue; }
	count += 1;
	}
	}
	}
	count
	}
	assert_eq!(e251(), 18_946_051);
	} // Total run-time: 4.40 seconds.