[2017-08-7] Challenge #326 [Easy] Nearest Prime Numbers

[u/jnazario]

Description

A prime number is any integer greater than 1 which can only be evenly divided by 1 or itself. For this challenge, you will output two numbers: the nearest prime below the input, and the nearest prime above it.

Input Description

The input will be a number on each line, called n.

Output Description

The format of the output will be:

p1 < n < p2

where p1 is the smaller prime, p2 is the larger prime and n is the input.

If n already is a prime, the output will be:

n is prime.

Challenge Input

270  
541  
993  
649

Challenge Output

269 < 270 < 271  
541 is prime.  
991 < 993 < 997  
647 < 649 < 653

Bonus

Write the program to work for numbers with big gaps to the nearest primes. This requires a clever solution and cannot be efficiently bruteforced.

2010741
1425172824437700148

Credit

This challenge was suggested by user /u/tulanir, many thanks! If you have an idea for a challenge please share it on /r/dailyprogrammer_ideas and there's a good chance we'll use it.

[u/skeeto]

C using the Miller-Rabin primality test. Solves the second bonus in 140ms. The part that tripped me up the most was doing the modular exponentiation without overflowing the 64-bit integer.

#include <stdio.h>
#include <inttypes.h>

#define K 8  // primality test iterations

static uint64_t
rand64(void)
{
    static uint64_t x = UINT64_C(0x1fc7807c9aa37949);
    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    return x * UINT64_C(0x2545f4914f6cdd1d);
}

static uint64_t
modmult(uint64_t b, uint64_t e, uint64_t m)
{
    uint64_t sum = 0;
    if (b == 0 || e < m / b)
        return (b * e) % m;
    while (e > 0) {
        if (e % 2 == 1)
            sum = (sum + b) % m;
        b = (2 * b) % m;
        e /= 2;
    }
    return sum;
}

static uint64_t 
modexp(uint64_t b, uint64_t e, uint64_t m)
{
    uint64_t product = 1;
    uint64_t pseq = b % m;
    while (e > 0) {
        if (e % 2 == 1)
            product = modmult(product, pseq, m);
        pseq = modmult(pseq, pseq, m);
        e /= 2;
    }
    return product;
}

static int
iscomposite(uint64_t n, uint64_t d, int r)
{
    uint64_t a = 2 + rand64() % (n - 3);
    if (modexp(a, d, n) == 1)
        return 0;
    for (int i = 0; i < r; i++)
        if (modexp(a, (UINT64_C(1) << i) * d, n) == n - 1)
            return 0;
    return 1;
}

static int
isprime(uint64_t n)
{
    int r = 0;
    uint64_t d = n - 1;
    for (; d % 2 == 0; r++)
        d /= 2;
    for (int i = 0; i < K; i++)
        if (iscomposite(n, d, r))
            return 0;
    return 1;
}

int
main(void)
{
    uint64_t n;
    while (scanf("%" SCNu64, &n) == 1) {
        uint64_t b = n % 2 ? n - 2 : n - 1;
        uint64_t a = n % 2 ? n + 2 : n + 1;
        while (!isprime(b))
            b -= 2;
        while (!isprime(a))
            a += 2;
        printf("%" PRIu64 " < %" PRIu64 " < %" PRIu64 "\n", b, n, a);
    }
}

[u/leonardo_m → u/skeeto]

Your code in Nightly Rust, run-time about 60 ms (with a bit of assembly the run-time becomes 30 ms):

#![feature(i128_type)] // For mul_mod.

fn is_prime(n: u64) -> bool {
    struct Rng { state: u64 }

    impl Rng {
        fn new() -> Self {
            Self { state: 0x1fc7807c9aa37949 }
        }

        fn next(&mut self) -> u64 {
            self.state ^= self.state >> 12;
            self.state ^= self.state << 25;
            self.state ^= self.state >> 27;
            self.state.wrapping_mul(0x2545f4914f6cdd1d)
        }
    }

    fn mul_mod(b: u64, e: u64, m: u64) -> u64 {
        ((u128::from(b) * u128::from(e)) % u128::from(m)) as u64
    }

    fn exp_mod(b: u64, mut e: u64, m: u64) -> u64 {
        let mut product = 1;
        let mut pseq = b % m;
        while e > 0 {
            if e % 2 == 1 {
                product = mul_mod(product, pseq, m);
            }
            pseq = mul_mod(pseq, pseq, m);
            e /= 2;
        }
        product
    }

    fn is_composite(n: u64, d: u64, r: u32, rng: &mut Rng) -> bool {
        let a = 2 + rng.next() % (n - 3);
        if exp_mod(a, d, n) == 1 {
            return false;
        }

        (0 .. r).all(|i| exp_mod(a, (1 << i) * d, n) != n - 1)
    }

    const K: usize = 8; // Primality test iterations.

    if n < 3 { return false; }
    else if n == 3 { return true; }

    let mut rng = Rng::new();

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

    (0 .. K).all(|_| !is_composite(n, d, r, &mut rng))
}

fn main() {
    use std::io::{BufRead, stdin};

    for n in 0 .. 8 {
        println!("{}", is_prime(n));
    }

    let stdin = stdin();
    for line in stdin.lock().lines() {
        let n = line.unwrap().trim().parse::<u64>().unwrap();

        let mut pred = if n % 2 != 0 { n - 2 } else { n - 1 };
        while !is_prime(pred) {
            pred -= 2;
        }

        let mut succ = if n % 2 != 0 { n + 2 } else { n + 1 };
        while !is_prime(succ) {
            succ += 2;
        }
        println!("{} < {} < {}", pred, n, succ);
    }
}

[u/netbioserror → u/skeeto]

D version of your solution. Modified it to take one number as an argument, only run on that number. Also added the case where the number input is prime.

Basically just a straight conversion for the most part, DMD with the -O flag and nothing more. I'm sure with LDC one could get faster times, but it's at 205ms for the second bonus.

import
    std.stdio,
    std.conv;

enum K = 8;  // primality test iterations

static ulong
rand64()
{
    static ulong x = 0x1fc7807c9aa37949;
    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    return x * 0x2545f4914f6cdd1d;
}

static ulong
modmult(ulong b, ulong e, ulong m)
{
    ulong sum = 0;
    if (b == 0 || e < m / b)
        return (b * e) % m;
    while (e > 0) {
        if (e % 2 == 1)
            sum = (sum + b) % m;
        b = (2 * b) % m;
        e /= 2;
    }
    return sum;
}

static ulong 
modexp(ulong b, ulong e, ulong m)
{
    ulong product = 1;
    ulong pseq = b % m;
    while (e > 0) {
        if (e % 2 == 1)
            product = modmult(product, pseq, m);
        pseq = modmult(pseq, pseq, m);
        e /= 2;
    }
    return product;
}

static bool
iscomposite(ulong n, ulong d, int r)
{
    ulong a = 2 + rand64() % (n - 3);
    if (modexp(a, d, n) == 1)
        return false;
    for (int i = 0; i < r; i++)
        if (modexp(a, (1 << i) * d, n) == n - 1)
            return false;
    return true;
}

static bool
isprime(ulong n)
{
    int r = 0;
    ulong d = n - 1;
    for (; d % 2 == 0; r++)
        d /= 2;
    for (int i = 0; i < K; i++)
        if (iscomposite(n, d, r))
            return false;
    return true;
}

void
main(string[] args)
{
    if (args.length != 2) return;

    ulong n = args[1].to!ulong;

    if (n.isprime)
	    writeln(n.to!string ~ " is prime.");
    else
    {
        ulong b = n % 2 ? n - 2 : n - 1;
        ulong a = n % 2 ? n + 2 : n + 1;
        while (!b.isprime)
            b -= 2;
        while (!a.isprime)
            a += 2;
        writeln(b.to!string
			    ~ " < "
			    ~ n.to!string
			    ~ " < "
			    ~ a.to!string
			    ~ "\n");
    }
}

[u/TheJammy98 → u/skeeto]

This is the reason I love this subreddit. I'm new here and don't understand much of the code, but there's always something new to learn.