0

I tried using the Miller-Rabin algorithm, but it can't detect very large numbers.

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

long long mulmod(long long a, long long b, long long mod)
{
    long long x = 0,y = a % mod;
    while (b > 0)
    {
        if (b % 2 == 1)
        {    
            x = (x + y) % mod;
        }
        y = (y * 2) % mod;
        b /= 2;
    }
    return x % mod;
}

long long modulo(long long base, long long exponent, long long mod)
{
    long long x = 1;
    long long y = base;
    while (exponent > 0)
    {
        if (exponent % 2 == 1)
        {
            x = (x * y) % mod;  
        }
        y = (y * y) % mod;
        exponent = exponent / 2;
    }
    return x % mod;
}

int Miller(unsigned long long int p, int iteration)
{
    int i;
    long long s;
    if (p < 2)
    {
        return 0;
    }
    if (p != 2 && p % 2==0)
    {
        return 0;
    }
    s = p - 1;
    while (s % 2 == 0)
    {
        s /= 2;
    }
    for (i = 0; i < iteration; i++)
    {
        long long a = rand() % (p - 1) + 1, temp = s;
        long long mod = modulo(a, temp, p);
        while (temp != p - 1 && mod != 1 && mod != p - 1)
        {
            mod = mulmod(mod, mod, p);
            temp *= 2;
        }
        if (mod != p - 1 && temp % 2 == 0)
        {
            return 0;
        }
    }
    return 1;
}

int main()
{
    int iteration = 5, cases;
    unsigned long long int num;
    scanf("%d", &cases);
    for(int i = 0; i < cases; i++)
    {
        scanf("%llu", &num);
        if(Miller(num, iteration))
        {
            printf("YES\n");    
        } 
        else
        {
            printf("NO\n");
        }   
    }
    return 0;
}

Output examples:

10 //cases
1
NO
2
YES
3
YES
4
NO
5
YES
1299827
YES
1951
YES
379
YES
3380
NO
12102
NO

I am trying to do my homework by creating a program that tells if a number is prime or not, print out YES if prime, NO if not. However, every time I submit the code to the online judge it just says "Wrong answer", when even my last attempt to do the assignment was without any efficient algorithm it says "Time limit exceeded".

Is there any way to determine if N is a prime number or not when N is [2 <= N <= 2^63-1]?

AB3
  • 67
  • 6
  • 1
    please provide an *example* of such input – Antti Haapala -- Слава Україні Sep 25 '20 at 14:22
  • `return x % mod;` can be written as `return x;` as `x` is always updated `% mod` (maybe the compiler optimizes it ... but why trust the compiler?). – pmg Sep 25 '20 at 14:23
  • You get the largest possible prime number P less than 2^63 and for a given N you compute N^(P-1) mod P and see if it's 1. – alinsoar Sep 25 '20 at 14:25
  • 2
    On Linux [Debian](http://debian.org/) you can find several open source prime number generators (e.g. the `primesieve-bin` package). Read also wikipage on [primality test](https://en.wikipedia.org/wiki/Primality_test). Read [*Modern C*](https://modernc.gforge.inria.fr/) and the documentation of your C compiler (e.g. [GCC](http://gcc.gnu.org/)...) and debugger (e.g. [GDB](https://www.gnu.org/software/gdb/)...) – Basile Starynkevitch Sep 25 '20 at 14:31
  • 2
    Miller-Rabin's whole point is to separate *very* large primes from composites; numbers less than 64 bits aren't "very large" at all. Miller-Rabin will absolutely work, if you implement it correctly. – ShadowRanger Sep 25 '20 at 14:34
  • Take care to keep a list of Carmichael numbers as well as special cases. – alinsoar Sep 25 '20 at 14:41
  • If your maximal number is only 2^63 then it should work fast the mr test using simple methods. as commented, the things become complex when you deal with big integers. – alinsoar Sep 25 '20 at 14:41
  • Using `rand()` is a bit random (well, pseudorandom). See the Wikipedia article for the lists of "witnesses" to test with to determine primality up to certain values guaranteed. For 64-bit numbers it can be done by using all prime numbers up to 37 as witnesses. – Ian Abbott Sep 25 '20 at 15:29

1 Answers1

1

OP's code has many possibilities of overflowing 63-bit math. e.g. x * y in x = (x * y) % mod;

At a minimum, recommend to go to unsigned math. e.g.: long long --> unsigned long long or simply uintmax_t.

For a mulmod() that does not overflow: Modular exponentiation without range restriction.

I'll look more into this later. GTG.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256