#From limd@cs.rpi.edu  Thu Sep 26 16:00:52 1996
#Received: from cs.rpi.edu by eos03a.eos.ncsu.edu (8.7.3/ES20May96)
#        id QAA28355; Thu, 26 Sep 1996 16:00:52 -0400 (EDT)
#From: limd@cs.rpi.edu
#Received: from fornax.cs.rpi.edu by cs.rpi.edu (5.67a/1.4-RPI-CS-Dept)
#	id AA17342; Thu, 26 Sep 1996 16:03:58 -0400 (limd from fornax.cs.rpi.edu)
#Date: Thu, 26 Sep 96 15:59:32 EDT
#Received: by fornax.cs.rpi.edu (4.1/2.3-RPI-CS-client)
#	id AA19737; Thu, 26 Sep 96 15:59:32 EDT
#Message-Id: <9609261959.AA19737@fornax.cs.rpi.edu>
#Apparently-To: kaltofen@pams.ncsu.edu
#Status: OR
#
#     It looks as though I sent you one of my older versions of the program.
#
#THe following file has the Solve procedure.
#
#						Darren
#
#From limd@cs.rpi.edu  Thu Sep 26 16:02:08 1996
#Received: from cs.rpi.edu by eos03a.eos.ncsu.edu (8.7.3/ES20May96)
#        id QAA28717; Thu, 26 Sep 1996 16:02:08 -0400 (EDT)
#From: limd@cs.rpi.edu
#Received: from juicer.cs.rpi.edu by cs.rpi.edu (5.67a/1.4-RPI-CS-Dept)
#	id AA17451; Thu, 26 Sep 1996 16:05:15 -0400 (limd from juicer.cs.rpi.edu)
#Date: Thu, 26 Sep 96 16:02:03 EDT
#Received: by juicer.cs.rpi.edu (4.1/2.3-RPI-CS-client)
#	id AA02922; Thu, 26 Sep 96 16:02:03 EDT
#Message-Id: <9609262002.AA02922@juicer.cs.rpi.edu>
#Apparently-To: kaltofen@pams.ncsu.edu
#Status: OR
#
#  Description:  This Maple program implements the multiple polynomial
#quadratic sieve method for factoring integers.
#  To run, type read `mqs2.mpl` at the Maple prompt

with(numtheory):
with(linalg):
readlib(ilog):
#Note: ilog[2] truncates

TS := proc(num)
#Description:  This procedure calculates the trucated square root of a number
local result, temp;

#  if num <4 then RETURN(2); fi;
  temp := evalf(sqrt(num),ilog[10](num)+1);
  result := trunc(temp);
  result
end:
######################################################################

Solver := proc(f,plist,pcount,rel_list,scount, N)
#Description:  This procedure solves the linear system modulo 2 for
#the Quadratic Sieve Method
#This function returns a non-trivial factor or the value 1
#Since linsolve doesn't work modulo two, we develope our answer via
#Gauss-Jordan elimination of an augmented matrix
local result, relation_index, bnd, sys, prime_index, temp_f, sys_index, newsys,
      y_squared, x_value, y_value, ntemp, zcheck;

  bnd := scount + pcount;
  sys := array(sparse, 1..scount, 1..bnd);
  for relation_index from 1 to scount do
    temp_f := subs(x = rel_list[relation_index], f);
    for prime_index from 1 to pcount do
      while temp_f mod plist[prime_index] = 0 do
        sys[relation_index, prime_index]:=1-sys[relation_index, prime_index];
        temp_f := temp_f / plist[prime_index]
      od
    od
  od;
  for sys_index from 1 to scount do
    sys[sys_index,sys_index + pcount] := 1
  od;

#print(sys);
#For some inexplicable reason, Gausselim has to be capitalized for using mod
  newsys := modp(Gausselim(sys),2);

  result := 1;
  y_squared := 1; 
  x_value := 1;
#ntemp := trunc(sqrt(N));
  ntemp := TS(N);

  relation_index := scount;
  zcheck := 0;
  for sys_index from 1 to pcount do
    zcheck := zcheck + newsys[relation_index,sys_index]
  od;
  while zcheck = 0 and result <= 1 do
    for sys_index from pcount+1 to bnd do
      if newsys[relation_index, sys_index] = 1
        then
          x_value := x_value * (rel_list[sys_index-pcount] + ntemp);
          y_squared := y_squared * (subs(x = rel_list[sys_index-pcount], f))
      fi 
    od;
    y_value := sqrt(y_squared);
    result := gcd(x_value - y_value,N) mod N;
 
    relation_index := relation_index - 1;
    zcheck := 0;
    for sys_index from 1 to pcount do
      zcheck := zcheck + newsys[relation_index,sys_index]
    od
  od;

#print(x_value);
#print(y_squared);
#print(y_value);
  
#print(result);
  result
end:
######################################################################

Sieve := proc(f,s1,s_limit, mplb, Tbound, B, N)
#Description: Finds values of x and f(x) where f(x) is smooth using
#Pomerance's improvement.  Temp is a two dimensional array, where
#column 1 is the x value, and column two is the f(x) value.
local result, temp, t_index, sieve_prime, s_index, primelog, smooth_count,
      answer, m_index;

  temp := array(sparse, 1..Tbound, 1..2);

  for t_index from 1 to Tbound do
    temp[t_index,1] := ilog[2](subs(x = t_index, f))

####temp[t_index,2] := subs(x = t_index, f)
  od;

  s_index := 1;
  while s_index <= s_limit do
#For all multiplicities
    sieve_prime := s1[s_index, 2*mplb + 1];
    primelog := ilog[2](sieve_prime);
    for m_index from 1 to mplb do
      t_index := s1[s_index,m_index * 2 - 1];
      if t_index = 0
        then t_index := sieve_prime ^ m_index;
      fi;
      while t_index <= Tbound do
        temp[t_index,1] := temp[t_index,1] - primelog;
        t_index := t_index + sieve_prime ^ m_index
      od;
      t_index := s1[s_index,m_index * 2];
      if t_index = 0
        then t_index := sieve_prime ^ m_index
      fi;
      while t_index <= Tbound do
        temp[t_index,1] := temp[t_index,1] - primelog;
        t_index := t_index + sieve_prime ^ m_index
      od
    od;
    s_index := s_index + 1
  od;

#adjust for odd t/even f(t)
#assumption n <> 1 mod 8

  t_index := 1;
  if t_index mod 2 = 1
    then t_index := t_index + 1
  fi;
  while t_index <= Tbound do
    temp[t_index,1] := temp[t_index,1] - 1;
    t_index := t_index + 2
  od; 
  
#print(temp);

  print(`Sieving Completed.`);
  result := [];
  smooth_count := 0;

  for t_index from 2 to Tbound do
if temp[t_index,1] >=7 and temp[t_index,1] < 9
then print(ifactor(subs(x = t_index,f)))
fi;
    if temp[t_index,1] < 9
      then
        smooth_count := smooth_count + 1;
#print(ifactor(subs(x = t_index,f)));
        result := [op(result),t_index]
    fi
  od;
  result
end:
######################################################################

Pre_Comp := proc(p, f)
#Description: This procedure precomputes the solutions to the congruence
#equation f(x) === 0(mod p) using the Maple function Roots.
local result;

###print(`Incoming function `.f);
  result := array(Roots(f) mod p);
###print(result);
  result
end:
######################################################################

QSS := proc(n,T, p_limit, aval, f)
#Description: This procedure attempts to find a non-trivial factor of n.
#Note: This program will not work if n has small factors
#TEST 113*127 = 14351
#TEST 541*547 = 295927
#TEST 1009*4243 = 4281187
#T is defined to be the sieve length
#p_limit is defined to be the number of primes in the factor base
local temp, p, msol, result, B, p_count, answer, mpower, mpb, solbound, a_f;

  p_count := 0;
  p := 3;
  mpb := 8;

  solbound := 2*mpb + 1;
  msol := array(sparse, 1..p_limit-1, 1..solbound);

#  subs(x = 2, f)

  while p_count < p_limit-1 do
    if L(n,p) = 1
      then
        p_count := p_count+1;
#The for loop takes care of all multiplicities upto p^mpb
        for mpower from 1 to mpb do
          temp := Pre_Comp(p ^ mpower, f);
          msol[p_count, (2 * mpower)-1] := temp[1,1];
          msol[p_count, mpower * 2] := temp[2,1]
        od;
        msol[p_count,solbound] := p
      fi;
    p := nextprime(p);
  od;

#print(msol);

  print(`Pre computing solutions to f(x) congruent to 0 mod p done.`);

#The Set B is the factor base
  B := [2];
  for temp from 1 to p_count do
    B := [op(B),msol[temp,solbound]]
  od;
  print(`The factor base is` . B);

  a_f := f * aval;
#print(a_f);

  result := Sieve(f,msol,p_count, mpb, T, B, n);
#print(`The x values with smooth f(x) values are` . result)
  result
end:
######################################################################

FindA := proc(N,M)
#Description:  This procedure determines a suitable a_value for the family
#of polynomials for MQS.  We try to find an a_value which is square that is
#around sqrt(2N)/M, where M is a bound on sieving.
local result;

  result := trunc(TS(2*N)/M);
  result := TS(result);
  while  L(N,result) <> 1 do
    result := result+1;
  od;
  result := result ^ 2;
  result
end:
######################################################################

MQS := proc(n)
#Description: This procedure attempts to find a non-trivial factor of n
#using the MQS method. 
#TEST 3571 * 3581 = 12787751
#TEST 463157*2465467 =  1141898299319  (22/20000);
#TEST 4356824539 * 1145661061 = 4991444223941575879
local a_value, b, b_value, c_value, i, numpoly, f, result, a_temp, T, pnum,
      t_result, r_index, smooth_count, answer;

  numpoly := 1;
  a_value := array(1..numpoly);
  b_value := array(1..numpoly);
  c_value := array(1..numpoly);
  f := array(1..numpoly);
  
  pnum := trunc(ilog[10](n) * ilog[10](n) / 2);
  T := pnum*500;

  a_temp := FindA(n,T);
  for i from 1 to numpoly do
    a_value[i] := a_temp;
  od;
  for i from 1 to numpoly do
    b := array(Roots(x^2 - n) mod a_value[i]);
    b_value[i] := b[i,1];
    c_value[i] := (b_value[i] ^ 2 - n)/ a_value[i];
    f[i] := x^2 * a_value[i] + x * 2 * b_value[i] + c_value[i];
  od;

  print(`Polynomials found.`);
  result := [];
  t_result := [];
  for i from 1 to numpoly do
    t_result := QSS(n, T, p_num, a_value[i], f[i]);
    for r_index from 1 to nops(t_result) do
      if not member(t_result[r_index],result)
        then result := [op(result),t_result[r_index]]
      fi
    od
  od;

  if pnum < smooth_count
    then 
      print(`We found ` . smooth_count . ` smooth values.`);
      answer := Solver(f, B, s_limit+1, result, smooth_count, N)
    else
      print(`Not enough smooth values found for our factor base.`);
      print(`The number of smooth f(x) values is ` . smooth_count);
      print(`The size of our factor base is `. pnum);
      answer := 1
  fi;
  print(`A factor found by the quadratic sieve method is `);
  answer
end: