#include  <math.h>


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


// _hi is a temporary 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 )
//
//..............................................................................



// Reciprocals of primes 3 thru 191
// A fudge factor of .0000000000000005 is added so that when
// multiplying them any rounding error is always positive.
//
const  double  Prime_Reciprocal_49[] =
{
   0x1.5555555555558p-2,    0x1.999999999999ep-3,    0x1.2492492492495p-3,
   0x1.745d1745d1749p-4,    0x1.3b13b13b13b17p-4,    0x1.e1e1e1e1e1e22p-5,
   0x1.af286bca1af2cp-5,    0x1.642c8590b2167p-5,    0x1.1a7b9611a7b98p-5,
   0x1.0842108421086p-5,    0x1.bacf914c1bad4p-6,    0x1.8f9c18f9c18fep-6,
   0x1.7d05f417d05f7p-6,    0x1.5c9882b93105ap-6,    0x1.3521cfb2b78c4p-6,
   0x1.15b1e5f75270fp-6,    0x1.0c9714fbcda3dp-6,    0x1.e9131abf0b76bp-7,
   0x1.cd85689039b0fp-7,    0x1.c0e070381c0e4p-7,    0x1.9ec8e951033ddp-7,
   0x1.8acb90f6bf3adp-7,    0x1.702e05c0b8173p-7,    0x1.51d07eae2f818p-7,
   0x1.446f86562d9fep-7,    0x1.3e22cbce4a905p-7,    0x1.323e34a2b10c2p-7,
   0x1.2c9fb4d812ca3p-7,    0x1.21fb78121fb7bp-7,    0x1.020408102040ap-7,
   0x1.f44659e4a4275p-8,    0x1.de5d6e3f8868ep-8,    0x1.d77b654b82c38p-8,
   0x1.b7d6c3dda338fp-8,    0x1.b2036406c80ddp-8,    0x1.a16d3f97a4b06p-8,
   0x1.920fb49d0e22dp-8,    0x1.886e5f0abb04dp-8,    0x1.7ad2208e0ecc6p-8,
   0x1.6e1f76b4337cap-8,    0x1.6a13cd1537293p-8,    0x1.571ed3c506b3dp-8
};



////////////////////////////////////////////////////////////////////////////////
//
uint32  Sprp_Two_Trial_49(    // Sprp 2 merged with trial division.

          uint64   N )        // A 49 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 index to the highest
//                         factor tested.
{
   uint32   Sprp = 0;
   uint64   hi;
   uint64   square;
   uint64   mask;
   uint64   d;
   uint64   r;
   uint32   c;
   uint32   z = 0;
   uint32   i = 0;

   double   pr;
   double   nd;
   double   npr;
//
//   if (sprp < 0) {}       // If N is prime or a psuedo prime,
//   else                   // Else N is proven composite,
//   {  if (sprp & 1) {}    //    If the result is a factor of N,
//      else {}             //    Else the result is twice the recipricol index.
//   }
//
// Trial division is embedded within the exponentiation loop.  It uses
// precomputed floating point recipricols to use multiplication instead
// of division.  Usually floating arithmetic can be performed in parallel
// with integer operations so trial division imposes little overhead.
//
// Note that only in very few cases does a factored value detect a composite
// that would not have otherwise determined by the full SPRP computation. 
// The main benefit is that with trial division the algorithm terminates early
// about a third of the time when the SPRP result is zero (proven composite).
//..............................................................................


pr = Prime_Reciprocal_49[ 0 ];           // Preload 1/3.
nd = (double)N;                          // N as a double; used for factoring
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 is under 64 bits,
{  r   = (1 << d) % N;                   //    2 ^ D mod N
   npr = nd * pr;                        //    N * (1/3)
}
else                                     // ELSE modular exponentiation,
{  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.
   npr  = nd * pr;                          // N * (1/3)

   for (i = 1; i <= c; i++)                 // DO over remaining bits in P,
   {  if (npr - floor(npr) <= pr) {break;}  // UNDO IF a factor was found.

      pr = Prime_Reciprocal_49[ i ];           // Preload next reciprocol.

      TimesMod( r, r, N, r, hi );              // R = R^2 mod N  where R < N
      mask >>= 1;                              // Next bit right.

      if (d & mask)                            // IF the bit is set,
      {  r <<= 1;                              //    Just double it.
      }

      npr = nd * pr;                           // N * (1/prime)
   }

   if (r >= N) { r -= N; }                  // Reduce modulus by N.  
}

if (npr - floor(npr) <= pr)              // IF a factor was found,
{  Sprp = (uint32)lround((double)(1 / pr)); // Return the factor.
}
else                                     // ELSE
{  Sprp = i << 1;                           // Return double the index.

   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 a psuedo prime,
         {  Sprp |= 0x80000000;             //       Probably prime, set sign bit.
            break;
   }  }  }
   else if ( r )                            // ELSE IF r = 1, N-1,
   {  Sprp |= 0x80000000;                   //    Probably prime, set sign bit.
}  }

return  Sprp;
}
//
////////////////////////////  sprp_two_trial_49.c  /////////////////////////////