Print("# findPQ(SIZE) finds two distinct primes p, q\n");
Print("# which are congruent to 3 mod 4, about the square root of SIZE\n");

Print("# \n");

Print("# findRootsBF(r, n) tries to find (half) of the square roots of r mod n\n");
Print("# by brute force, without even checking whether r is a square mod n\n");

Print("# \n");

Print("# findRoots := function (r, p, q) or (r, [p, q])\n");
Print("# finds the four square roots of r modulo n = p * q,\n");
Print("# where p and q are distincts primes\n");
Print("# which are congruent to 3 mod 4. Checks (about) everything...\n");

findPQ := function (SIZE)
# generates suitable p, q, about the size of square root of SIZE
local sq, delta, p, q, MAX;
sq := RootInt(SIZE);
MAX := Minimum([sq,2^60-1]);
delta := Random([2..MAX]);
p := sq-delta;
repeat
p := NextPrimeInt(p);;
until p mod 4 = 3;
q := Maximum([Int(SIZE/p),p]);
repeat
q := NextPrimeInt(q);;
until q mod 4 = 3;
return [p , q];
end;

findN := function (SIZE)
# generates suitable n, about the size of SIZE
local sq, delta, p, q, MAX;
sq := RootInt(SIZE);
MAX := Minimum([sq,2^60-1]);
delta := Random([2..MAX]);
p := sq-delta;
repeat
p := NextPrimeInt(p);;
until p mod 4 = 3;
q := Maximum([Int(SIZE/p),p]);
repeat
q := NextPrimeInt(q);;
until q mod 4 = 3;
return p * q;
end;

findSquare := function(r, n)
return r^2 mod n;
end;

isSquare := function (arg)
local r, p, q, x, res;
if Length(arg) < 2 or Length(arg) > 3
then
Print("Either isSquare(x, p) or isSquare(x, p, q)\n");
return;
fi;
r := arg[1]; p := arg[2];
if Length(arg) = 2
then
res := PowerMod(r, (p-1)/2, p);
return res = 1;
else
q := arg[3];
res := List([p, q], x -> PowerMod(r, (x-1)/2, x));
return res = [1,1];
fi;
end;

findRootsBF := function(r, n)
local i, noOfRoots, roots;
noOfRoots := 0;
roots := [];
for i in [1..Int(n/2)] do
  if i mod 1000000 = 0 then 
    Print("Tried up to ", EuclideanQuotient(i, 1000000),  
      " millions\n"); 
  fi;
  if PowerMod(i, 2, n) = r then 
    noOfRoots := noOfRoots + 1;
    Print("Root ", noOfRoots, " is ", i, "\n");
    Add(roots, i); 
  fi;
od;
return roots;
end;

findRoots := function (arg)
# arg should be r, p, q
# Assumes r is a square modulo p and q...
# Assumes p and q primes and congruent to 3 mod 4
# OK, we will check both

local r, p, q, a, x, temp;
if Length(arg) = 3 
then 
r := arg[1]; p := arg[2]; q := arg[3];
elif Length(arg) = 2 and IsList(arg[2]) and Length(arg[2]) = 2 then
r := arg[1]; p := arg[2][1]; q := arg[2][2];
else
Print("# Error: either findRoots(r, p, q) or findRoots(r,[p,q])\n");
return;
fi;


if not IsPrime(p) then Print("# Error: ", p , " is not prime\n"); return; fi;
if not IsPrime(q) then Print("# Error: ", q , " is not prime\n"); return; fi;
if p = q then Print("# Error: different primes, please\n"); return; fi;
if not (p mod 4 = 3) then 
Print("# Error: ", p , " is not congruent to 3 mod 4\n"); 
return; 
fi; 
if not (q mod 4 = 3) then 
Print("# Error: ", q , " is not congruent to 3 mod 4\n"); 
return; fi; 
if Gcd(r, p * q) > 1 then
Print("# Error: ", r , " not coprime to ", p, " * ", q, "\n"); 
return; 
fi;
if not isSquare(r, p, q)  then 
Print("# Error: ", r , " is not a square\n"); 
return; 
fi;

a := List([p,q], x -> PowerMod(r, (x+1)/4, x));

temp := List([[1,1],[1,-1],[-1,1],[-1,-1]
             ],
   x -> ChineseRem([p,q],[x[1]*a[1],x[2]*a[2]]));
Sort(temp);
return temp;
end;

squareRootModP := function(a, p)
local z, tmp, even, odd;
if not IsPrimeInt(p) then
Print(p, " is not prime\n");
return;
fi;
if Legendre (a, p ) = -1 then
Print(a, " is not a square modulo ", p, "\n");
return;
fi;
z := PrimitiveRootMod(p);

tmp := p-1; even := 1;
while (tmp mod 2) = 0 do
even := 2 * even;
tmp := EuclideanQuotient(tmp, 2);
od;
odd := (p-1)/even;
return [even,odd];
end;