::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
:
function SPRP:  Strong Probable Prime test.

 entry N     word,         &Unsignes odd number to test
       B     word          :Base to test against

  exit Sprp  Bit           :0 - Proven composite; 1 - Likely prime; maybe not

 local Z     byte,         &Number of zero bits on the right side of D
       D     word,         &N - 1; right justified
       R     word,         &Residual
       Nc    cell          :N in 64 bits

precondition
   N // 1  fault  N;        FAULT:  Cannot test 0 or 1 for primality.
   N /\ 1  fault  N|d#;     FAULT:  Cannot test even numbers for primality.
.
:...............................................................................


: Find D and Z that satisfy N - 1 = D * 2^Z, where D is odd and Z >= 0.
:
D = (N - 1) // 1;                      Always remove the rightmost zero.

DO until  D /\ 1:                      DO until N - 1 is right justified,
   D //= 1
   Z  += 1;                               Count the zero bits.
-

R = Power.Mod( B, D, N );              B ^ D mod N

IF R = 1  or  R = N - 1:               IF residual is 1 or N - 1,
   Sprp = 1;                              N is a SPRP base B (possibly prime).

ELSE IF R:                             ELSE IF checking if provably composite,
   : Otherwise, see if B^(2^R) mod N = N - 1, where 0 <= R < Z.
   :
   Nc = positive( N );                    Zero extend N to 64 bits.

   DO Z times:                            DO over each zero bit,
      R = cell{ positive( R |* R ) % Nc };   R^2 mod N

      IF R = N - 1:                          IF prime or a strong psuedo prime,
         Sprp = 1;                              Not known to be composite
   UNDO;                                  UNDO
   -  .
ELSE IF le( N, B ):                    ELSE IF R = 0  and  N is small,
   Sprp = 1;                              Not known to be composite
.

return