--* From bmt@spadserv.watson.ibm.com  Fri Dec 18 21:22:18 1992
--* Received: from spadserv.watson.ibm.com by radical.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA11080; Fri, 18 Dec 1992 21:22:18 -0500
--* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA24959; Fri, 18 Dec 1992 21:27:50 -0500
--* Date: Fri, 18 Dec 1992 21:27:50 -0500
--* From: bmt@spadserv.watson.ibm.com
--* X-External-Networks: yes
--* Message-Id: <9212190227.AA24959@spadserv.watson.ibm.com>
--* To: axc-bug@radical.watson.ibm.com
--* Subject: compiler errors seem incorrect and unintelligible

--@ Fixed  by: SMW Thu Oct 07 09:54:01 1993
--@ Tested by: <name of new or existing file in test directory>
--@ Summary:   <One line description of real problem and the fix>


--Copyright The Numerical Algorithms Group Limited 1991.

#include "aslib.as"
#include "prime.as"
#library listlib "list.aso"
import listlib

macro 
  Boolean == Bit
  B == Bit
  N == Integer
  min(aa,bb) == if aa < bb then aa else bb
  one? aa == aa = 1

ListMultiDictionary(S:with Object): with
   dictionary: () -> %
   insert!: (S,%) -> %
   insert!: (S,%,N) -> %
   empty?: % -> B
   inspect: % -> S
   count: (S, %) -> N
   remove!: (S, %) -> %
 == add

facFlags ==> Enumeration(nil, sqfr, irred, prime)

FFE(S:with Object):with
   Object
   bracket: (facFlags, S, Integer) -> %
   apply: (%, Enumeration(xpnt)) -> Integer
   set!: (%, Enumeration(xpnt), Integer) -> Integer
 == add

   Rep:=Record(flg:facFlags, fctr:S, xpnt:Integer)
   

Factored(S: with Object) : with
   makeFR: (S, List FFE(S)) -> %
   factorList: % -> List FFE(S)
   unit: % -> S
 == add


--% IntegerFactorizationPackage
-- recoded MBM Nov/87

+++ This Package contains basic methods for integer factorization.
+++ The factor operation employs trial division up to 10,000.  It
+++ then tests to see if n is a perfect power before using Pollards
+++ rho method.  Because Pollards method may fail, the result
+++ of factor may contain composite factors.  We should also employ
+++ Lenstra's elliptic curve method.

I      ==> Integer
B      ==> Bit
FF     ==> Factored I
LMI    ==> ListMultiDictionary I

IntegerFactorizationPackage: with

    factor : I -> FF
      ++ factor(n) returns the full factorization of integer n
    squareFree   : I -> FF
      ++ squareFree(n) returns the square free factorization of integer n
    BasicMethod : I -> FF
      ++ BasicMethod(n) returns the factorization
      ++ of integer n by trial division
    PollardSmallFactor: I -> Partial(I)
       ++ PollardSmallFactor(n) returns a factor
       ++ of n or "failed" if no one is found

  == add

    import Segment I
    import IntegerRoots
    import IntegerPrimesPackage
    import Partial I
    import Enumeration(nil,sqfr,irred,prime)
    import FFE(I)
    import List(FFE(I))
    import FF
    import Enumeration(xpnt)

    local concat!: (List FFE I, FFE I) -> List FFE I

    BasicSieve: (I, I) -> FF


    squareFree(n:I):FF ==
       if n<0 then (m := -n; u:I := -1)
              else (m := n; u := 1)
       (m > 1) and not(failed? (v := perfectSqrt m)) =>
          for rec in (l := factorList(sv := squareFree(v::I))) repeat
            rec.xpnt := 2 * rec.xpnt
          makeFR(u * unit sv, l)
       lim := 1 + approxNthRoot(m,3)
       lim > 100000 => makeFR(u, factorList factor m)
       x := BasicSieve(m, lim)
       y :=
         one? unit x => factorList x
         concat_!(factorList x, [sqfr,unit(x),1]$FFE(I))
       makeFR(u, y)

    -- Pfun(y: I,n: I): I == (y**2 + 5) rem n
    PollardSmallFactor(n:I):Partial(I) ==
       -- Use the Brent variation
       x0:I := random()$I
       m := 100
       y:I := x0 rem n
       r:I := 1
       q:I := 1
       G:I := 1
       while G <= 1 repeat
          x:I := y
          for i in 1..r repeat
             y := (y*y+5) rem n
             q := (q*abs(x-y)) rem n
             k:I := 0
          while not( (k>=r) or (G>1)) repeat
             ys:I := y
             for i in 1..min(m,r-k) repeat
                y := (y*y+5) rem n
                q := q*abs(x-y) rem n
             G := gcd(q,n)
             k := k+m
          r := 2*r
       if G=n then
          while G<=1 repeat
             ys := (ys*ys+5) rem n
             G := gcd(abs(x-ys),n)
       G=n => failed
       G::Partial(I)

    BasicSieve(r:I, lim:I):FF ==
       l:List(I) :=
          [generate(~1;~2;~2;~4;~2;~4;~2;~4;~6;~2;~6)]
       concat_!(l, rest rest rest l)
       d := 2
       n := r
       ls := empty()$List(FFE(I))
       for s in l repeat
          d > lim => return makeFR(n, ls)
          if n<d*d then
             if n>1 then ls := concat_!(ls, [prime,n,1]$FFE(I))
             return makeFR(1, ls)
          m:I:=0
          for mm in 0.. while zero?(n rem d) repeat (n := n quo d; m:=mm)
          if m>0 then ls := concat_!(ls, [prime,d, m]$FFE(I))
          d := d+s
       never 

    BasicMethod(n:I):FF ==
       local u:I
       if n<0 then (m := -n; u := -1)
              else (m := n; u := 1)
       x := BasicSieve(m, 1 + approxSqrt m)
       makeFR(u, factorList x)

    factor(m:I):FF ==
       local u:I
       zero? m => makeFR(m,nil)
       if negative? m then (n := -m; u := -1)
                      else (n := m; u := 1)
       bfac := BasicSieve(n, 10000)
       flb := factorList bfac
       one?(n := unit bfac) => makeFR(u, flb)
       a:LMI := dictionary() -- numbers yet to be factored
       b:LMI := dictionary() -- prime factors found
       f:LMI := dictionary() -- number which could not be factored
       insert_!(n, a)
       while not empty? a repeat
          n := inspect a; c:I := count(n, a); remove_!(n, a)
          prime?(n)$IntegerPrimesPackage => insert_!(n, b, c)
          -- test for a perfect power
          (s :Record(base:I, exponent:I) := perfectNthRoot n).exponent > 1 =>
            insert_!(s.base, a, c * s.exponent)
          not(failed? (dd := PollardSmallFactor n)) =>
             d := retract dd
             for m in 0.. while zero?(n rem d) repeat n := n quo d
             insert_!(d, a, m * c)
             if n > 1 then insert_!(n, a, c)
          -- an elliptic curve factorization attempt should be made here
          insert_!(n, f, c)
       -- insert prime factors found
       while not empty? b repeat
          n := inspect b; c := count(n, b); remove_!(n, b)
          flb := concat_!(flb, [prime,n,c]$FFE(I))
       -- insert non-prime factors found
       while not empty? f repeat
          n := inspect f; c := count(n, f); remove_!(n, f)
          flb := concat_!(flb, [nil,n,c]$FFE(I))
       makeFR(u, flb)