#OK to post homework
#Joseph Koutsoutis, 01-28-2024, Assignment 2

with(numtheory):
with(ArrayTools):
Help:=proc():
print(` MyIfactorSimple(n), MyIfactorRho(N) `):
end:

#1

#Here are some small examples based on pages 30-60
#of Frank Garvan's awesome Maple booklet

#a := [seq(f(i), i=1..5)]:
#f := x -> x^2:
#a[3..4];
#f := 'f': a := 'a': 

#a := {1,2,3}:
#b := {2,4,6}:
#member(3,a minus b);
#a := 'a': b := 'b':

#p := x^2+1:
#f := unapply(p,x):
#map(f, [seq(i,i=1..5)]);
#f := 'f': p := 'p':

#Limit((x+1)/(x-1), x=infinity): value(%);

#maximize(sin(x), x=0..1);

#Int(sqrt(1+x^6),x=0..1):
#value(%); evalf(%);


#2

#Here is my example for the first problem:

#n := 12345678987654321:
#for a from rand(1..n)() while gcd(n,a) <> 1 do od:
#a&^phi(n) mod n; #this outputs 1

#Here is my example for the second problem:

#p := 123456789098765432111: #isprime(p) confirms this is prime
#q := 123456789098765432111303: #isprime(q) confirms this is prime
#n := p*q:
#e := 40999169172903038692291405278497292027581:
#d := 13485563134850386244097435291778029864163861: #e*d mod phi(n) outputs 1
#m := 1234567890987654321: 
#m - (m&^e mod n)&^d mod n: #this outputs 0


#3

#MyIfactorSimple(n) factors an integer n >= 1 using a simple while loop that
#searches for small primes factoring n before inspecting larger primes,
#and outputs an array consisting of the prime factors of n.
#This never beat the built-in ifactor in my testing.
#(13 digit integers took 0.0 seconds for ifactor to compute vs >10 seconds for MyIfactorSimple)
MyIfactorSimple := proc(n) local n1,p,q,factors:
 factors := Array([]):
 p := 2:
 n1 := n:
 while (n1 <> 1) do:
  q := n1 / p:
  while type(q,integer) do:
   Append(factors,p):
   n1 := q:
   q := n1 / p:
  od:
  p := nextprime(p);
 od:
 
 RETURN(factors):
end:

#MyIfactorRho is an algorithm for factoring a prime N >= 2, with a chance of failure.
#This implementation is based on Donald Knuth's description of John Pollard's rho method
#found in TAOCP vol 2.
#During my testing, this function was faster than MyIfactorSimple for integers with
#10+ digits, but it still could not beat the built-in ifactor.
#(around 28-29 digits took the built-in ifactor roughly .1 seconds to factor vs roughly 5 seconds
#for MyIfactorRho)
MyIfactorRho := proc(N) local x1,x2,k,l,n,g,factors:
 factors := Array([]):
 x1 := 5: x2 := 2: k := 1: l := 1: n := N:

 while not isprime(n) do:
  g := gcd(x2-x1,n):
  while g = 1 do:
   k -= 1:
   if k = 0 then 
    x2 := x1:
    l *= 2:
    k := l:
   fi:
   x1 := x1*x1 + 1 mod n:
   g := gcd(x2-x1,n):
  od:
  if not isprime(g) then RETURN(FAIL): fi:
  Append(factors,g):
  n /= g:
  x1 := x1 mod n:
  x2 := x2 mod n:
 od:

 Append(factors,n):
 RETURN(sort(factors)):
end: