मुझे C में एक प्रोग्राम मिला है जिसने मिलर-राबिन प्रारंभिक परीक्षण का एक नियतात्मक संस्करण लागू किया है यहां . हालांकि, संशोधित कोड (जिसे नीचे देखा जा सकता है) 2^32 से बड़ी संख्याओं के साथ काम करते समय काम नहीं करता है, भले ही मैं अपने नंबरों को संग्रहीत करने के लिए अहस्ताक्षरित लंबे लंबे डेटा प्रकार का उपयोग करता हूं। जो 2^64 तक के पूर्णांकों को धारण करने में सक्षम होना चाहिए। यह कहाँ गलत हो जाता है?

संक्षेप में: मेरी समस्या यह है कि मेरा कोड सही ढंग से निर्धारित करता है कि कोई संख्या प्रमुख है या नहीं, लेकिन केवल अगर यह 2 ^ 32 से छोटी है, तो ऐसा नहीं होना चाहिए क्योंकि मैं संख्याओं को 2 ^ 64 तक स्टोर कर सकता हूं


unsigned long long power(unsigned long long a, unsigned long long n, unsigned long long mod)
{
    unsigned long long power = a;
    unsigned long long result = 1;

    while (n)
    {
        if (n & 1)
            result = (result * power) % mod;
        power = (power * power) % mod;
        n >>= 1;
    }
    return result;
}

// n−1 = 2^s * d with d odd by factoring powers of 2 from n−1
unsigned long long witness(unsigned long long n, unsigned long long s, unsigned long long d, unsigned long long a)
{
    unsigned long long x = power(a, d, n);
    unsigned long long y;

    while (s) {
        y = (x * x) % n;
        if (y == 1 && x != 1 && x != n-1)
            return 0;
        x = y;
        --s;
    }
    if (y != 1)
        return 0;
    return 1;
}

/*
 * if n < 1,373,653, it is enough to test a = 2 and 3;
 * if n < 9,080,191, it is enough to test a = 31 and 73;
 * if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
 * if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
 * if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
 * if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
 * if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
 */

int is_prime_mr(unsigned long long n)
{
    if (((!(n & 1)) && n != 2 ) || (n < 2) || (n % 3 == 0 && n != 3))
        return 0;
    if (n <= 3)
        return 1;

    unsigned long long d = n / 2;
    unsigned long long s = 1;
    while (!(d & 1)) {
        d /= 2;
        ++s;
    }
    unsigned long long a1 = 4759123141;
    unsigned long long a2 = 1122004669633;
    unsigned long long a3 = 2152302898747;
    unsigned long long a4 = 3474749660383;
    if (n < 1373653)
        return witness(n, s, d, 2) && witness(n, s, d, 3);
    if (n < 9080191)
        return witness(n, s, d, 31) && witness(n, s, d, 73);
    if (n < a1)
        return witness(n, s, d, 2) && witness(n, s, d, 7) && witness(n, s, d, 61);
    if (n < a2)
        return witness(n, s, d, 2) && witness(n, s, d, 13) && witness(n, s, d, 23) && witness(n, s, d, 1662803);
    if (n < a3)
        return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11);
    if (n < a4)
        return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13);
    return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13) && witness(n, s, d, 17);
}

int main()
{
  unsigned long long in1 = 4294967291;
  unsigned long long in2 = 4294967311;
  int a = is_prime_mr(in1); //4294967291 is the last prime below 2^32, should return 1 and does so correctly
  printf("%d\n",a);
  int b = is_prime_mr(in2); //4294967311 is the first prime after 2^32, should also return 1 but returns 0
  printf("%d",b);
  return 0;
}
1
Michael 13 अप्रैल 2020, 15:47

1 उत्तर

सबसे बढ़िया उत्तर

जैसा कि @Michael Dorgan ने सुझाव दिया था, समस्या दो 64-बिट पूर्णांकों के गुणन में निहित थी जिसे powermod के एक अलग कार्यान्वयन का उपयोग करके हल किया गया था।

//computes (a*b) (mod n)
unsigned long long mul_mod(unsigned long long a, unsigned long long b, unsigned long long m)
{
   unsigned long long d = 0, mp2 = m >> 1;
   if (a >= m) a %= m;
   if (b >= m) b %= m;
   for (int i = 0; i < 64; ++i)
   {
       d = (d > mp2) ? (d << 1) - m : d << 1;
       if (a & 0x8000000000000000ULL)
           d += b;
       if (d >= m) d -= m;
       a <<= 1;
   }
   return d;
}
//computes (a^b) (mod m)
unsigned long long powermod(unsigned long long a, unsigned long long b, unsigned long long m)
{
    unsigned long long r = m==1?0:1;
    while (b > 0) {
        if (b & 1)
            r = mul_mod(r, a, m);
        b = b >> 1;
        a = mul_mod(a, a, m);
    }
    return r;
}

मेरा कोड अब 2^64 तक के पूर्णांकों के लिए काम करता है।

0
Michael 17 अप्रैल 2020, 10:21