> interface( verboseproc = 2 );
> eval(`ifactor/lenstra`);

                                  1


  proc(n)
local
i, s, prime, f, A, X, Z, rgen, sp, a, r, curves, B1, kg, kgg;
option `Copyright (c) 1990 by the University of Waterloo. A\
ll rights reserved.`;
    if nargs < 1 or 3 < nargs then
        error "wrong number of arguments."
    end if;
    if nargs < 3 then B1 := 1000000 else B1 := args[3] end if
    ;
    if nargs < 2 then curves := 30
    else curves := args[2]
    end if;
    if modp(n, 2) = 0 then return 2 end if;
    if modp(n, 3) = 0 then return 3 end if;
    s := evalf(1/2*sqrt(5) - 1/2, 30);
    A := array(1 .. curves);
    X := array(1 .. curves);
    Z := array(1 .. curves, [1 $ curves]);
    rgen := rand(1 .. n - 1);
    for i to curves do
        a := 0;
        while modp(a*(a^2 - 1)*(9*a^2 - 1), n) = 0 do
            r := rgen();
            kg := r^2 + 6;
            kgg := igcd(kg, n);
            if kgg <> 1 then return kgg end if;
            a := modp(6*r/kg, n)
        end do;
        A[i] := modp(1/16*(-3*a^4 - 6*a^2 + 1)/a^3 + 1/2, n);
        X[i] := modp(3/4*a, n)
    end do;
    prime := 2;
    while prime <= B1 do
        sp := round(s*prime);
        `ifactor/lenstra/mulpp`(1, n, A, sp, prime, B1, X, Z)
            ;
        f := Z[1];
        for i from 2 to curves do
            `ifactor/lenstra/mulpp`(i, n, A, sp, prime, B1, X,
                Z);
            f := modp(f*Z[i], n)
        end do;
        f := igcd(f, n);
        if f <> 1 then return f end if;
        prime := nextprime(prime)
    end do;
    FAIL
end proc

> eval(`ifactor/lenstra/mulpp`);

  proc(i, n, A, mm, nn, B1, X, Z)
local pow, ax, az;
option `Copyright (c) 1990 by the University of Waterloo. A\
ll rights reserved.`;
    ax := X[i];
    az := Z[i];
    pow := nn;
    while pow <= B1 do
        `ifactor/lenstra/ellmul`(n, A[i], mm, nn, ax, az,
            'ax', 'az');
        pow := pow*nn
    end do;
    X[i] := ax;
    Z[i] := az
end proc

> eval(`ifactor/lenstra/ellmul`);

  proc(n, A, mm, nn, px, pz, aax, aaz)
local ax, az, bx, bz, cx, cz, tmpx, tmpz, d, e, t1, t2;
option `Copyright (c) 1990 by the University of Waterloo. A\
ll rights reserved.`;
    cx := px;
    cz := pz;
    e := mm;
    d := nn - mm;
    if e < d then
        `ifactor/lenstra/elldoub`(n, A, px, pz, 'bx', 'bz');
        ax := px;
        az := pz;
        d := d - e
    else
        `ifactor/lenstra/elldoub`(n, A, px, pz, 'ax', 'az');
        bx := px;
        bz := pz;
        e := e - d
    end if;
    while e <> 0 do
        if e < d then
            tmpx := bx;
            tmpz := bz;
            t1 := modp((ax - az)*(bx + bz), n);
            t2 := modp((ax + az)*(bx - bz), n);
            bx := modp(cz*modp((t1 + t2)^2, n), n);
            bz := modp(cx*modp((t1 - t2)^2, n), n);
            d := d - e
        else
            tmpx := ax;
            tmpz := az;
            t1 := modp((ax - az)*(bx + bz), n);
            t2 := modp((ax + az)*(bx - bz), n);
            ax := modp(cz*modp((t1 + t2)^2, n), n);
            az := modp(cx*modp((t1 - t2)^2, n), n);
            e := e - d
        end if;
        cx := tmpx;
        cz := tmpz
    end do;
    aax := ax;
    aaz := az
end proc

> eval(`ifactor/lenstra/elldoub`);

  proc(n, A, ax, az, cx, cz)
local t1, t2;
option `Copyright (c) 1990 by the University of Waterloo. A\
ll rights reserved.`;
    t1 := modp((ax + az)^2, n);
    t2 := modp((ax - az)^2, n);
    cx := modp(t1*t2, n);
    cz := modp((t1 - t2modp(At1 - t2) + t2, n), n)
end proc

>