# MAS3008 NUMBER THEORY
> 
# SOME NUMBER-THEORETIC ALGORITHMS
> 
# 1. Binary Powering Algorithm
#  powerod(a,k,n) computes a^k mod n
> 
> 
> powermod:=proc(a,k,n) local b,c,i;
> c:=a;b:=1;i:=k;
> while i>0 do
>   if (i mod 2=1) then b:=b*c mod n; i:=i-1;fi;
>   c:=c*c mod n; i:=i/2;
> od;
> RETURN(b);
> end:# 
> 
# An example:
> powermod(15,90,91);

# 2. Miller-Rabin Primality Test.
# spsp(n,a) returns "Pass" if n passes the Miller-Rabin test to base a,
# that is, either n is prime or n is a strong pseudoprime to base a. It
# returns "Fail" otherwise, in which case n is composite. Where possible
# a non-trivial factor of n is also returned.
# 
# The algorithm as given below uses the powermod subroutine of (1)
# above. It would run slightly faster if instead we used MAPLE's
# built-in implementation of the Binary Powering Algorithm, i.e.  
# replace b:=powmod(a,m,n); in line 7 below by b:=a&^m mod n;
# 
> spsp:=proc(n,a) local g,m,s,b,b1,i;
> if a>=n or a<2 then RETURN(`Invalid input`) fi;
> if n mod 2=0 then RETURN(`Fail: factor`,2) fi;
> g:=gcd(a,n); if g>1 then RETURN(`Fail: factor`,g) fi;
> m:=n-1; s:=0;
> while m mod 2=0 do m:=m/2; s:=s+1 od;
> b:=powermod(a,m,n);
> if b=1 then RETURN(Pass) fi;
> if b=n-1 then RETURN(Pass) fi;
> for i from 1 to s-1 do
> b1:=b*b mod n;
> if b1=1 then RETURN(`Fail: factor`,gcd(b-1,n)) fi;
> if b1=n-1 then RETURN(Pass) fi;
> b:=b1;
> od;
> RETURN(Fail);
> end:
> # 
# Some examples:
# 
> spsp(561,1);spsp(561,561);
> spsp(56,2);
> spsp(567,7);
> spsp(561,2);
> 
> spsp(2047,2);# 
> spsp(2047,3);# 
> n:=1373653;
> spsp(n,2);spsp(n,3);spsp(n,5);
> n:=1001152801;
> spsp(n,2);
# 
# 3. Pollard's p-1 method.
# pollard(n,a,kmax) runs Pollard's p-1 algorithm to look for a proper
# factor of n, starting with base a, up to kmax steps.
# 
# Fail2 indicates its worth increasing kmax. Fail1 indicates that 
# increasing kmax won't help - changing a might help, but is unlikely to
# in general.
# 
# Again, we use powermod from above to carry out the Binary Powering
# Algorithm
# 
> pollard:=proc(n,a0,kmax) local a,g,k;
> a:=a0;
> for k from 2 to kmax do
>   a:=powermod(a,k,n);
>   g:=gcd(a-1,n);
>   if g=n then RETURN(Fail1); fi;
>   if g>1 then RETURN(g,n,k) fi;
> od;
> RETURN(Fail2);
> end:
> 
> 
# Some examples (as described in lectures):
# 
> pollard(403,2,50);
> pollard(1891,2,50);
> pollard(1891,3,50);
> pollard(1891,5,50);
> n:=5157437;# 
> pollard(n,2,50);
# 
> 
> n:=2^32+1;# 
> pollard(n,2,100);# 
> pollard(n,3,100);
> 
> 
# 4. Pollard's `rho' method
# 
# rho(n,c,x0,imax) looks for a proper factor of n, carrying out up to
# imax iterations of the rho algorithm, using the function x->x^2+c
# starting with the given x0.
# 
# The output is either Fail, or a factor, followed by the number of step
# taken to find it.
# 
> rho:=proc(n,c,x0,imax) local x,y,i,g;
> x:=x0;y:=x0;
> for i from 1 to imax do
>   x:=x^2+c mod n;
>   y:=y^2+c mod n; y:=y^2+c mod n;
>   g:=gcd(x-y,n);
>   if g=n then RETURN(Fail); fi;
>   if g>1 then RETURN(g,i); fi; 
> od;
> RETURN(Fail);
> end:
> 
> 
# We'll try it out on the same examples as above:
> 
> rho(403,1,1,100);
> rho(1891,1,1,100);
> rho(5157437,1,1,100);
> rho(2^32+1,1,1,100);
> rho(2^64+1,1,1,100);
> rho(2^64+1,1,1,1000);