//////////////////////////////////////////////////////////////////////////////// // // These declarations are specific to 64 bit Visual C. // #include // Shorthand for unsigned integers. // typedef unsigned long uint32; typedef unsigned long long uint64; typedef long long int64; // _hi is a temporary uint64 for the upper 64 bits of the product. // Note that a divide overflow hardware error occurs if [a, b] / mod // is over 64 bits; even though the remainder is under 64 bits. // #define TimesMod( _a, _b, _mod, _tm, _hi ) \ _udiv128( _hi, _umul128( _a, _b, &_hi ), _mod, &_tm ) #define log2_64( _u32, _u64 ) _BitScanReverse64( &_u32, _u64 ) // //.............................................................................. //////////////////////////////////////////////////////////////////////////////// // uint32 Sprp_Two_Trial( // sprp( N, 2 ) N > 2 and is odd uint64 N ) // A 50 to 64 bit number to test // Return: The sign bit is sprp( N, 2 ). // Odd: a factor of N // Even: the low bits are 2 times the highest factor tested { uint32 Sprp = 0; uint64 square; uint64 d; uint64 mask; uint64 r = 1; uint32 z = 0; uint32 c; uint32 i; uint64 hi; uint64 prime; uint64 mod = 1; // This constant is used to generate small primes up to 137. // Each two bits encodes the delta between primes (0:2, 1:4, 2:6, 3:8). // The odd, 119, is in included as the delta from 113 to 127 is big (14). // Each round the constant is shifted right. // When more than 32 rounds the low 2 bits are 0 and odd numbers are tested. // uint64 delta = 0x9e44798629189110; // //.............................................................................. if ( (N >> 1 == 0) || (N & 1 == 0) ) { return 1; // Invalid argument } d = N >> 1; // (N - 1) / 2: Remove rightmost zero. while ( !(d & 1) ) // DO until N - 1 is right justified, { d >>= 1; z += 1; // Count the zero bits. } // r = 2 ^ D mod N // if ( d <= 63 ) // IF the power will fit in a Cell, { r = (1 << d) % N; // Compute it directly. prime = 3; mod = N % prime; } else // ELSE at least one squaring, { log2_64( c, d ); // log2( d ) c -= 5; // Number of bits in P less 6. r = (uint64)1 << (d >> (uint64)c); // Starting value for top 6 bits. r = r % N; // Keep R < N or TimesMod overflows. mask = (uint64)1 << (uint64)c; // Bit past the top 6 bits. prime = 3; for (i = 1; i <= c; i++) // DO over remaining bits in P, { mod = N % prime; TimesMod( r, r, N, r, hi ); // R = R^2 mod N where R < N mask >>= 1; // Next bit right. prime += ((delta & 3) + 1) << 1; delta >>= 2; if (mod == 0) { break; }; if (d & mask) // IF the bit is set, { if ((int64)r < 0) // IF doubling will overflow, { r = (r << 1) - N; // Double and reduce by N. } else // ELSE double is over, { r <<= 1; // Just double it. if (r >= N) // IF reducable, { r -= N; // Reduce modulus by N. } } } } if (r >= N) { r -= N; } // Reduce modulus by N. } Sprp = prime; // Return any prime factor. if ( mod ) // IF not factored, { Sprp = prime << 1; if ( (r > 1) && (r != (N - 1)) ) // IF R is not 0, 1, nor N - 1, { // See if 2^(2^R) mod N = N - 1, where 0 <= R < Z. // for ( ; z > 0; z-- ) // DO over each zero bit, { TimesMod( r, r, N, r, hi ); // R = R^2 mod N if ( r == (N - 1) ) // IF prime or strong psuedo prime, { Sprp |= 0x80000000; // Not known to be composite. break; // UNDO } } } else if ( r ) // ELSE IF r = 1, N-1, { Sprp |= 0x80000000; // Probably prime, set sign bit. } } return Sprp; } // //////////////////////////// sprp_two_trial.c ////////////////////////////////