Initial population.

1 files changed, 282 insertions, 0 deletions

diff --git a/src/input/ecfact.as.pamphlet b/src/input/ecfact.as.pamphlet
new file mode 100644
index 00000000..bf330973
--- /dev/null
+++ b/src/input/ecfact.as.pamphlet
@@ -0,0 +1,282 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input ecfact.as}
+\author{The Axiom Team}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{License}
+<<license>>=
+-- Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+-- All rights reserved.
+--
+-- Redistribution and use in source and binary forms, with or without
+-- modification, are permitted provided that the following conditions are
+-- met:
+--
+--     - Redistributions of source code must retain the above copyright
+--       notice, this list of conditions and the following disclaimer.
+--
+--     - Redistributions in binary form must reproduce the above copyright
+--       notice, this list of conditions and the following disclaimer in
+--       the documentation and/or other materials provided with the
+--       distribution.
+--
+--     - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--       names of its contributors may be used to endorse or promote products
+--       derived from this software without specific prior written permission.
+--
+-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+-- IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+-- TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+-- PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+-- OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+-- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+-- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+@
+<<*>>=
+<<license>>
+#include "axiom.as"
+#pile
+
+--% Elliptic curve method for integer factorization
+--  This file implements Lenstra's algorithm for integer factorization.
+--  A divisor of N is found by computing a large multiple of a rational
+--  point on a randomly generated elliptic curve in P2 Z/NZ.
+--  The Hessian model is used for the curve (1) to simplify the selection
+--  of the initial point on the random curve and (2) to minimize the
+--  cost of adding points.
+--  Ref:  IBM RC 11262, DV Chudnovsky & GV Chudnovsky
+--  SMW Sept 86.
+
+--% EllipticCurveRationalPoints
+--)abbrev domain ECPTS EllipticCurveRationalPoints
+
+EllipticCurveRationalPoints(x0:Integer, y0:Integer, z0:Integer, n:Integer): ECcat == ECdef where
+    Point   ==> Record(x: Integer, y: Integer, z: Integer)
+
+    ECcat ==> AbelianGroup with
+        double: %  -> %
+        p0:     %
+        HessianCoordinates: % -> Point
+
+    ECdef ==> add
+        Rep == Point
+        import from Rep
+        import from List Integer
+
+        Ex == OutputForm
+
+        default u, v:  %
+
+        apply(u:%,x:'x'):Integer == rep(u).x
+        apply(u:%,y:'y'):Integer == rep(u).y
+        apply(u:%,z:'z'):Integer == rep(u).z
+        import from 'x'
+        import from 'y'
+        import from 'z'
+
+        coerce(u:%): Ex        == [u.x, u.y, u.z]$List(Integer) :: Ex
+        p0:%                   == per [x0 rem n, y0 rem n, z0 rem n]
+        HessianCoordinates(u:%):Point  == rep u
+
+        0:%   ==
+            per [1, (-1) rem n, 0]
+        -(u:%):%  ==
+            per [u.y, u.x, u.z]
+        (u:%) = (v:%):Boolean ==
+            XuZv := u.x * v.z
+            XvZu := v.x * u.z
+            YuZv := u.y * v.z
+            YvZu := v.y * u.z
+            (XuZv-XvZu) rem n = 0 and (YuZv-YvZu) rem n = 0
+        (u:%) + (v:%): % ==
+            XuZv := u.x * v.z
+            XvZu := v.x * u.z
+            YuZv := u.y * v.z
+            YvZu := v.y * u.z
+            (XuZv-XvZu) rem n = 0 and (YuZv-YvZu) rem n = 0 => double u
+            XuYv := u.x * v.y
+            XvYu := v.x * u.y
+            Xw := XuZv*XuYv - XvZu*XvYu
+            Yw := YuZv*XvYu - YvZu*XuYv
+            Zw := XvZu*YvZu - XuZv*YuZv
+            per [Yw rem n, Xw rem n, Zw rem n]
+        double(u:%): % ==
+            import from PositiveInteger
+            X3 := u.x**(3@PositiveInteger)
+            Y3 := u.y**(3@PositiveInteger)
+            Z3 := u.z**(3@PositiveInteger)
+            Xw := u.x*(Y3 - Z3)
+            Yw := u.y*(Z3 - X3)
+            Zw := u.z*(X3 - Y3)
+            per [Yw rem n, Xw rem n, Zw rem n]
+        (n:Integer)*(u:%): % ==
+            n < 0 => (-n)*(-u)
+            v := 0
+            import from UniversalSegment Integer
+            for i in 0..length n - 1 repeat
+                if bit?(n,i) then v := u + v
+                u := double u
+            v
+
+
+--% EllipticCurveFactorization
+--)abbrev package ECFACT EllipticCurveFactorization
+
+EllipticCurveFactorization: with
+        LenstraEllipticMethod: (Integer)                   -> Integer
+        LenstraEllipticMethod: (Integer, Float)         -> Integer
+        LenstraEllipticMethod: (Integer, Integer, Integer) -> Integer
+        LenstraEllipticMethod: (Integer, Integer)          -> Integer
+
+        lcmLimit: Integer -> Integer
+        lcmLimit: Float-> Integer
+
+        solveBound: Float -> Float
+        bfloor:     Float -> Integer
+        primesTo:   Integer -> List Integer
+        lcmTo:      Integer -> Integer
+    == add
+        import from List Integer
+        Ex == OutputForm
+        import from Ex
+        import from String
+        import from Float
+
+        NNI==> NonNegativeInteger
+        import from OutputPackage
+        import from Integer, NonNegativeInteger
+        import from UniversalSegment Integer
+
+        blather:Boolean := true
+
+        --% Finding the multiplier
+        flabs (f: Float): Float == abs f
+        flsqrt(f: Float): Float == sqrt f
+        nthroot(f:Float,n:Integer):Float == exp(log f/n::Float)
+
+        bfloor(f: Float): Integer == wholePart floor f
+
+        lcmLimit(n: Integer):Integer ==
+            lcmLimit nthroot(n::Float, 3)
+        lcmLimit(divisorBound: Float):Integer ==
+            y := solveBound divisorBound
+            lcmLim := bfloor exp(log divisorBound/y)
+            if blather then
+                output("The divisor bound is", divisorBound::Ex)
+                output("The lcm Limit is", lcmLim::Ex)
+            lcmLim
+
+        -- Solve the bound equation using a Newton iteration.
+        --
+        -- f = y**2 - log(B)/log(y+1)
+        --
+        -- f/f' = fdf =
+        --    2                 2
+        --   y (y + 1)log(y + 1)  - (y + 1)log(y + 1) logB
+        --   ---------------------------------------------
+        --                                 2
+        --              2y(y + 1)log(y + 1)  + logB
+        --
+        fdf(y: Float, logB: Float): Float ==
+            logy  := log(y + 1)
+            ylogy := (y + 1)*logy
+            ylogy2:= y*logy*ylogy
+            (y*ylogy2 - logB*ylogy)/((2@Integer)*ylogy2 + logB)
+        solveBound(divisorBound:Float):Float ==
+            -- solve               y**2 = log(B)/log(y + 1)
+            -- although it may be  y**2 = log(B)/(log(y)+1)
+            relerr := (10::Float)**(-5)
+            logB := log divisorBound
+            y0   := flsqrt log10 divisorBound
+            y1   := y0 - fdf(y0, logB)
+            while flabs((y1 - y0)/y0) > relerr repeat
+                y0 := y1
+                y1 := y0 - fdf(y0, logB)
+            y1
+
+        -- maxpin(p, n, logn) is max d s.t. p**d <= n
+        maxpin(p:Integer,n:Integer,logn:Float): NonNegativeInteger ==
+            d: Integer := bfloor(logn/log(p::Float))
+            if d < 0 then d := 0
+            d::NonNegativeInteger
+
+        multiple?(i: Integer, plist: List Integer): Boolean ==
+            for p in plist repeat if i rem p = 0 then return true
+            false
+
+        primesTo(n:Integer):List Integer ==
+            n < 2 => []
+            n = 2 => [2]
+            plist := [3, 2]
+            i:Integer := 5
+            while i <= n repeat
+                if not multiple?(i, plist) then plist := cons(i, plist)
+                i := i + 2
+                if not multiple?(i, plist) then plist := cons(i, plist)
+                i := i + 4
+            plist
+        lcmTo(n:Integer):Integer ==
+            plist := primesTo n
+            m: Integer := 1
+            logn := log(n::Float)
+            for p in plist repeat m := m * p**maxpin(p,n,logn)
+            if blather then
+                output("The lcm of 1..", n::Ex)
+                output("            is", m::Ex)
+            m
+        LenstraEllipticMethod(n: Integer):Integer ==
+            LenstraEllipticMethod(n, flsqrt(n::Float))
+        LenstraEllipticMethod(n: Integer, divisorBound: Float):Integer ==
+            lcmLim0 := lcmLimit divisorBound
+            multer0 := lcmTo lcmLim0
+            LenstraEllipticMethod(n, lcmLim0, multer0)
+        InnerLenstraEllipticMethod(n:Integer, multer:Integer, 
+                             X0:Integer, Y0:Integer, Z0:Integer):Integer ==
+            import from EllipticCurveRationalPoints(X0,Y0,Z0,n)
+            import from Record(x: Integer, y: Integer, z: Integer)
+            p  := p0
+            pn := multer * p
+            Zn := HessianCoordinates.pn.z
+            gcd(n, Zn)
+
+        LenstraEllipticMethod(n: Integer, multer: Integer):Integer ==
+            X0:Integer := random()
+            Y0:Integer := random()
+            Z0:Integer := random()
+            InnerLenstraEllipticMethod(n, multer, X0, Y0, Z0)
+
+        LenstraEllipticMethod(n:Integer, lcmLim0:Integer, multer0:Integer):Integer ==
+            nfact: Integer := 1
+            for i:Integer in 1.. while nfact = 1 repeat
+                output("Trying elliptic curve number", i::Ex)
+                X0:Integer := random()
+                Y0:Integer := random()
+                Z0:Integer := random()
+                nfact := InnerLenstraEllipticMethod(n, multer0, X0, Y0, Z0)
+                if nfact = n then
+                    lcmLim := lcmLim0
+                    while nfact = n repeat
+                        output("Too many iterations... backing off")
+                        lcmLim := bfloor(lcmLim * 0.6)
+                        multer := lcmTo lcmLim
+                        nfact := InnerLenstraEllipticMethod(n, multer0, X0, Y0, Z0)
+            nfact
+
+
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}