////////////////////////////////////////////////////////////////////////////////
//
// These declarations are specific to 64 bit Visual C.
//
#include  <intrin.h>

// 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  ////////////////////////////////