100   '           =======================
  110   '           |  Proving Primality  |
  120   '           =======================
  130   '
  140   ' program: aks_c4.ub
  150   ' Aldo D., February 1, 2003
  160   '
  170   ' This program is based on AKS Conjecture 4 and on
  180   ' Michael Scott's MIRACL program. In this program
  190   ' parameter r is prime number or a power of prime.
  200   '
  210   ' Integer n having more than 4288 bits may cause
  220   ' an overflow.
  230   '
  240   emaword -268
  250   '
  260   input "Input n:";N
  270   '
  280   I=1:R=2:K=N@R
  290   while (K^2@R=1)
  300      inc R
  310      K=N@R
  320   wend
  330   '
  340   if K=0 then print "N is divisible by";R:end
  350   print "r =";R;":";"k =";K
  360   '
  370   S=0:T=1
  380   dim ema(S;R+1,0),ema(T;R+1,0),A(2)
  390   '
  400   A(0)=-1:A(1)=1:'                            A = x-1
  410   ema(T;0,0)=-1:ema(T;1,0)=1
  420   '
  430   '   Mod exp: (x-1)^n (mod x^r-1, n)
  440   '   -------------------------------
  450   I=len(N)-2
  460   while (I>=0)
  470     gosub *S_x_S:'                    (S*S)@(x^r-1)@n
  480     if bit(I,N)=1 then gosub *S_x_A:' (S*A)@(x^r-1)@n
  490     dec I
  500   wend
  510   '
  520   '   Test result
  530   '   -----------
  540   block ema(S;*,0)=0
  550   ema(S;K,0)=1
  560   ema(S;0,0)=N-1
  570   '
  580   I=0
  590   while and{I<R,ema(S;I,0)=ema(T;I,0)}
  600     inc I
  610   wend
  620   '
  630   if I=R then print "N is prime"
  640         :else print "N is composite"
  650   end
  660   ' -------------------------------------------------
  670   '
  680   '  (S(x)*S(x))@(x^r-1)@n
  690   '  ---------------------
  700   *S_x_S
  710   local I,J,K
  720   '
  730   block ema(S;*,0)=block ema(T;*,0)
  740   block ema(T;*,0)=0
  750   '
  760   for I=0 to R-1
  770     for J=0 to R-1
  780       K=(I+J)@R
  790       ema(T;K,0)=(ema(T;K,0)+ema(S;I,0)*ema(S;J,0))@N
  800     next J
  810   next I
  820   return
  830   '
  840   '  (S(x)*A(x))@(x^r-1)@n
  850   '  ---------------------
  860   *S_x_A
  870   local I,J,K
  880   '
  890   block ema(S;*,0)=block ema(T;*,0)
  900   block ema(T;*,0)=0
  910   '
  920   for I=0 to R-1
  930     for J=0 to 1
  940       K=(I+J)@R
  950       ema(T;K,0)=(ema(T;K,0)+ema(S;I,0)*A(J))@N
  960     next J
  970   next I
  980   return
  990   ' -------------------------------------------------