# From limd@cs.rpi.edu  Fri Sep 20 15:06:23 1996
# Received: from cs.rpi.edu by eos03a.eos.ncsu.edu (8.7.3/ES20May96)
#         id PAA23767; Fri, 20 Sep 1996 15:06:23 -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 AA02654; Fri, 20 Sep 1996 15:09:26 -0400 (limd from juicer.cs.rpi.edu)
# Date: Fri, 20 Sep 96 15:04:51 EDT
# Received: by juicer.cs.rpi.edu (4.1/2.3-RPI-CS-client)
# 	id AA23604; Fri, 20 Sep 96 15:04:51 EDT
# Message-Id: <9609201904.AA23604@juicer.cs.rpi.edu>
# Apparently-To: kaltofen@pams.ncsu.edu
# Status: OR
# 
#       Hello, it is nice to hear from you again.  Sorry about the delay,
# but I had moved my maple files to disk and I have just retreived them today.
# The code for MPQS will be in the following email message.
# 							Darren
# 
# From limd@cs.rpi.edu  Fri Sep 20 15:07:37 1996
# Received: from cs.rpi.edu by eos03a.eos.ncsu.edu (8.7.3/ES20May96)
#         id PAA24050; Fri, 20 Sep 1996 15:07:37 -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 AA02804; Fri, 20 Sep 1996 15:10:39 -0400 (limd from juicer.cs.rpi.edu)
# Date: Fri, 20 Sep 96 15:07:31 EDT
# Received: by juicer.cs.rpi.edu (4.1/2.3-RPI-CS-client)
# 	id AA23819; Fri, 20 Sep 96 15:07:31 EDT
# Message-Id: <9609201907.AA23819@juicer.cs.rpi.edu>
# Apparently-To: kaltofen@pams.ncsu.edu
# Status: OR
# 
#  Description:  This Maple program implements the multiple quadratic sieve
#  method for factoring integers.
#  To run, type read `mqs.mpl` at the Maple prompt

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

Sieve := proc(f,s1,s_limit,s2, Tbound)
#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;

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

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

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

  s_index := 1;
  while s_index <= s_limit do
#Multiplicity 1
    sieve_prime := s1[s_index,3];
    t_index := s1[s_index,1];
    if t_index = 0
      then t_index := sieve_prime
    fi;
    while t_index <= Tbound do
      temp[t_index,2] := temp[t_index,2] - ilog[2](sieve_prime);
      t_index := t_index + sieve_prime
    od;
    t_index := s1[s_index,2];
    if t_index = 0
      then t_index := sieve_prime
    fi;
    while t_index <= Tbound do
      temp[t_index,2] := temp[t_index,2] - ilog[2](sieve_prime);
      t_index := t_index + sieve_prime
    od;
#Multiplicity 2
    sieve_prime := s2[s_index,3];
    t_index := s2[s_index,1];
    if t_index = 0
      then t_index := sieve_prime^2
    fi;
    while t_index <= Tbound do
      temp[t_index,2] := temp[t_index,2] - ilog[2](sieve_prime);
      t_index := t_index + sieve_prime^2
    od;
    t_index := s2[s_index,2];
    if t_index = 0
      then t_index := sieve_prime^2
    fi;
    while t_index <= Tbound do
      temp[t_index,2] := temp[t_index,2] - ilog[2](sieve_prime);
      t_index := t_index + sieve_prime^2
    od;
    s_index := s_index + 1;
  od;

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

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

  result := {};
  for t_index from 2 to Tbound do
    if temp[t_index,2] < 2
      then result := result union {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;

  result := array(Roots(f) mod p);
#print(p);
#print(result);
  result
end:

QS := proc(n, T, p_limit, f)
#Description: This procedure attempts to find a non-trivial factor of n
#using the QS method.  It takes as parameters n, a composite number; T, a
#bound for sieving; and p_limit, a bound on the factor base
#Note: This program will not work if n has small factors
#TEST 113*127 = 14351
#TEST 541*547 = 295927
local temp, p, p_index, msol1, msol1_index, msol2, msol2_index,
      result, B;

  p_index := 2;
  p := ithprime(p_index);

  msol1_index := 0;
  msol1 := array(sparse, 1..p_limit, 1..3);
  msol2_index := 0;
  msol2 := array(sparse, 1..p_limit, 1..3);

###Removed since a family of polynomials are found 
###f := (x + trunc(sqrt(n)))^2 - n;
#  subs(x = 2, f)

  while p < p_limit do
    if L(n,p) = 1
      then
        temp := Pre_Comp(p, f);
        msol1_index := msol1_index + 1;
        msol1[msol1_index,1] := temp[1,1];
        msol1[msol1_index,2] := temp[2,1];
        msol1[msol1_index,3] := p;
#Multiplicity of 2
        temp := Pre_Comp(p*p, f);
        msol2_index := msol2_index + 1;
        msol2[msol2_index,1] := temp[1,1];
        msol2[msol2_index,2] := temp[2,1];
        msol2[msol2_index,3] := p;
        #This term is not p*p in order to facilitate sieving later.
    fi;
    p_index := p_index + 1;
    p := ithprime(p_index);
  od;

#print(msol1);
#print(msol2);

#The Set B is the factor base
  B := [2];
  for temp from 1 to msol1_index do
    B := [op(B), msol1[temp,3]];
  od;
  print(`The factor base is` . B);
  result := Sieve(f,msol1,msol1_index,msol2, T);
  result
end:

FindA := proc(N,M,past_result)
#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(sqrt(2*N)/M);
  result := trunc(sqrt(result));
  while  L(N,result^2) <> 1 or result <= past_result do
    result := result + 1;
  od;
  result := result ^ 2;
  result
end:

MQS := proc(n, T, p_limit)
#Description: This procedure attempts to find a non-trivial factor of n
#using the MQS method.  It takes as parameters n, a composite number; T, a
#bound for sieving; and p_limit, a bound on the factor base.
#TEST: MQS(3571 * 3581, 200, 400);
local a_value, b, b_value, c_value, i, numpoly, f, result;

#number of polynomials to be used
  numpoly := 2;
  a_value := array(0..numpoly);
  b_value := array(1..numpoly);
  c_value := array(1..numpoly);
  f := array(0..numpoly);
  
  a_value[0] := 1;
  f[0] := expand((x + trunc(sqrt(n)))^2 - n);

  a_value[1] := 0;
  for i from 1 to numpoly do
    a_value[i] := FindA(n,T,a_value[i]);
    b := array(Roots(x^2 - n) mod a_value[i]);
    b_value[i] := b[2,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];
    if i <> numpoly
      then a_value[i+1] := a_value[i];
    fi;
  od;

  result := {};
  for i from 0 to numpoly do
    print(f[i]);
    # result := result union QS(n, T, p_limit, a_value[i]*f[i]);
    result := QS(n, T, p_limit, a_value[i]*f[i]);
    print(`The x values with smooth f(x) values are` , result)
  od;
end;