|
Subject: aldor/axiom interoperability
From: Gregory Vanuxem
To: Peter Broadbery
Date: Thu, 13 Mar 2008 22:21:44 +0100
Dear Peter, *
I was waiting for your small Aldor patch and never see it. Did you send
it somewhere ? The main concern of this mail is if you are able to
compile ecfact.as (attached) and execute the
LenstraEllipticMethod?(Integer) function. On my system it fails, a bug
somewhere. Am I alone with this issue ? Martin ? I'm still using the
old build process and want to switch to the new one but this stops me. I
do not want to take your time but maybe you have an idea of what's going
on with this issue.
The output is :
Looking in OutputPackage() for ??349042727 with code 483270060
and an error message : "export not found"
Greg
PS : I precise that I applied one of your patchs, sorry I don't remember
which one, to be able to use recent versions of Aldor. The problem
mentioned here still remains on my machine.
aldor -- 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.
#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
aldor Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/7214143568811717751-25px001.as
using Aldor compiler and options
-O -Fasy -Fao -Flsp -lfricas -Mno-ALDOR_W_WillObsolete -DFriCAS -Y $FRICAS/algebra -I $FRICAS/algebra
Use the system command )set compiler args to change these
options.
The )library system command was not called after compilation.
fricas LenstraEllipticMethod(99)
There are no library operations named LenstraEllipticMethod
Use HyperDoc Browse or issue
)what op LenstraEllipticMethod
to learn if there is any operation containing "
LenstraEllipticMethod " in its name.
Cannot find a definition or applicable library operation named
LenstraEllipticMethod with argument type(s)
PositiveInteger
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
The problem seems to be related to trying to generate Axiom output from inside an Aldor function. E.g.
aldor #include "axiom.as"
#pile
TestOutput: with
testOutput: (Integer) -> Integer
== add
import from String
import from OutputPackage
testOutput(x:Integer):Integer ==
output("help!")
x
aldor Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/7970768026816943239-25px003.as
using Aldor compiler and options
-O -Fasy -Fao -Flsp -lfricas -Mno-ALDOR_W_WillObsolete -DFriCAS -Y $FRICAS/algebra -I $FRICAS/algebra
Use the system command )set compiler args to change these
options.
The )library system command was not called after compilation.fricas testOutput(1)
There are no library operations named testOutput
Use HyperDoc Browse or issue
)what op testOutput
to learn if there is any operation containing " testOutput " in
its name.
Cannot find a definition or applicable library operation named
testOutput with argument type(s)
PositiveInteger
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
aldor #include "axiom.as"
#pile
TestPrint: with
testPrint: (Integer) -> Integer
== add
import from String
import from OutputForm
testPrint(x:Integer):Integer ==
messagePrint("help!")$OutputForm
x
aldor Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/1420678903966076930-25px005.as
using Aldor compiler and options
-O -Fasy -Fao -Flsp -lfricas -Mno-ALDOR_W_WillObsolete -DFriCAS -Y $FRICAS/algebra -I $FRICAS/algebra
Use the system command )set compiler args to change these
options.
The )library system command was not called after compilation.fricas testPrint(1)
There are no library operations named testPrint
Use HyperDoc Browse or issue
)what op testPrint
to learn if there is any operation containing " testPrint " in
its name.
Cannot find a definition or applicable library operation named
testPrint with argument type(s)
PositiveInteger
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
Here is one way to replace Axiom's output function
with a Lisp function that can be called from Aldor.
aldor #include "axiom.as"
#pile
-- implement output for Aldor
output(x:String):Void == {
import { FORMAT: (Boolean,String,String) -> Void } from Foreign Lisp;
FORMAT(true,"~a~%",x);
}
output(x:String,y:OutputForm):Void == {
import { FORMAT: (Boolean,String,String,String) -> Void } from Foreign Lisp;
FORMAT(true,"~a ~a~%",x,unparse(convert(y)$InputForm)$InputForm);
}
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
aldor Compiling FriCAS source code from file
/var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/8598620176429485019-25px007.as
using Aldor compiler and options
-O -Fasy -Fao -Flsp -lfricas -Mno-ALDOR_W_WillObsolete -DFriCAS -Y $FRICAS/algebra -I $FRICAS/algebra
Use the system command )set compiler args to change these
options.
The )library system command was not called after compilation.The output shows that the routine works for some integers
fricas LenstraEllipticMethod(12)
There are no library operations named LenstraEllipticMethod
Use HyperDoc Browse or issue
)what op LenstraEllipticMethod
to learn if there is any operation containing "
LenstraEllipticMethod " in its name.
Cannot find a definition or applicable library operation named
LenstraEllipticMethod with argument type(s)
PositiveInteger
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
and it fails for others. (But now it gives some hints why. :-)
fricas LenstraEllipticMethod(99)
There are no library operations named LenstraEllipticMethod
Use HyperDoc Browse or issue
)what op LenstraEllipticMethod
to learn if there is any operation containing "
LenstraEllipticMethod " in its name.
Cannot find a definition or applicable library operation named
LenstraEllipticMethod with argument type(s)
PositiveInteger
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
|