begin groundwork for sage testing

Makefile

@@ -55,13 +55,17 @@ ifeq ($(CC),clang)
 WARNFLAGS += -Wgcc-compat
 endif
 
+SAGE ?= sage
+SAGES= $(shell ls test/*.sage)
+BUILDPYS= $(SAGES:test/%.sage=build/%.py)
+
 ARCHFLAGS += $(XARCHFLAGS)
 CFLAGS  = $(LANGFLAGS) $(WARNFLAGS) $(INCFLAGS) $(OFLAGS) $(ARCHFLAGS) $(GENFLAGS) $(XCFLAGS)
 CXXFLAGS = $(LANGXXFLAGS) $(WARNFLAGS) $(INCFLAGS) $(OFLAGS) $(ARCHFLAGS) $(GENFLAGS) $(XCXXFLAGS) 
 LDFLAGS = $(ARCHFLAGS) $(XLDFLAGS)
 ASFLAGS = $(ARCHFLAGS) $(XASFLAGS)
 
-.PHONY: clean all test bench todo doc lib bat
+.PHONY: clean all test bench todo doc lib bat sage sagetest
 .PRECIOUS: build/%.s
 
 HEADERS= Makefile $(shell find src include test -name "*.h") $(shell find . -name "*.hxx") build/timestamp
@@ -149,6 +153,21 @@ build/%.s: src/$(FIELD)/$(ARCH)/%.c $(HEADERS)
 
 build/%.s: src/$(FIELD)/%.c $(HEADERS)
 	$(CC) $(CFLAGS) -S -c -o $@ $<
+	
+sage: $(BUILDPYS)
+
+sagetest: sage lib
+	LD_LIBRARY_PATH=build sage build/test_decaf.sage
+
+$(BUILDPYS): $(SAGES) build/timestamp
+	cp -f $(SAGES) build/
+	$(SAGE) --preparse $(SAGES:test/%.sage=build/%.sage)
+	# some sage versions compile to .sage.py
+	for f in $(SAGES:test/%.sage=build/%); do \
+		 if [ -e $$f.sage.py ]; then \
+		 	 mv $$f.sage.py $$f.py; \
+		 fi; \
+	  done
 
 doc/timestamp:
 	mkdir -p doc

test/curve.sage

@@ -0,0 +1,292 @@
+from idealized import Idealized
+from collections import namedtuple
+
+debugging = True
+def debug_print(foo):
+    if debugging: print foo
+
+checkGroupLaws = True
+checkTorsion = True
+checkIsogenies = True
+
+def memoize(f):
+    # list cache because my __hash__ hack doesn't seem to work
+    cache = []
+    def ff(*args, **kwargs):
+        key = (tuple(args),tuple(sorted(kwargs.iteritems())))
+        for key_,value in cache:
+            if key == key_: return value
+        out = f(*args,**kwargs)
+        cache.append((key,out))
+        return out   
+        
+    try:
+        ff.__name__ = f.__name__
+    except AttributeError: pass
+    return ff
+
+def EcBase(curvename,varnames,ad=()):
+    if isinstance(ad,str) or isinstance(ad[0],str):
+        ad = Idealized.vars(ad)
+    
+    class Inner(namedtuple(curvename,(v for v in varnames))):
+        params = ad
+        torsion_points = {}
+        def __new__(cls,*xy):
+            def apply_invariants(xy,x):
+                for inv in cls.invariants(*(ad+xy)):
+                    x = x.assuming(inv)
+                return x
+            
+            xy = tuple(xy)
+            if len(xy) == 0:
+                xy = Idealized.uvars(varnames)
+                xy = [apply_invariants(xy,x) for x in xy]
+            else:
+                for i,inv in enumerate(cls.invariants(*(ad + xy))):
+                    if inv != 0:
+                        raise Exception("Invariant inv[%d] not satisfied for %s: got \n%s" %
+                             (i,curvename,str(inv)))
+
+            return super(Inner,cls).__new__(cls,*xy)
+                    
+        varnames = "xy"
+        
+        @classmethod
+        def invariants(self,*args): return []
+        
+        @classmethod
+        @memoize
+        def check_group(cls):
+            if checkGroupLaws:
+                debug_print("Checking group law for %s..." % cls.__name__)
+                a,b,c,z = cls(),cls(),cls(),cls.basepoint
+                if a+z != a:
+                    raise Exception("Base point is not identity!")
+                if a-a != z:
+                    raise Exception("Subtraction doesn't work!")
+                if a+b != b+a:
+                    raise Exception("Addition is not commutative!")
+                #if a+(b+c) != (a+b)+c:
+                #    raise Exception("Addition is not associative!")
+            
+            for t,n in cls.torsion():
+                if checkTorsion:
+                    debug_print("  Checking %d-torsion..." % n)
+                    cls.check_torsion(t,n)
+                #if n not in cls.torsion_points:
+                #    cls.torsion_points[n] = set()
+                #cls.torsion_points[n].add(cls(*t(cls.basepoint)))
+                
+        @classmethod
+        def check_torsion(cls,f,n):
+            P = Q = cls()
+            good = False
+            for i in xrange(1,n+1):
+                Q = cls(*f(Q))
+                if Q == P:
+                    if i==n:
+                        good = True
+                        break
+                    raise Exception("Claimed %d-torsion, but is actually %d-torsion" % (n,i))
+            if not good: raise Exception("Claimed %d-torsion, but isn't" % n)
+            if n*P+n*cls(*f(P)) == cls.basepoint:
+                raise Exception("Torsion operation inverts element")
+        
+        @classmethod
+        def torsion(cls):
+            return []
+        
+        def __sub__(self,other):
+            return self + (-other)
+            
+        def __mul__(self,other):
+            if other==0: return self.basepoint
+            if other < 0: return -(self*-other)
+            if other==1: return self
+            if is_even(other): return (self+self)*(other//2)
+            return (self+self)*(other//2) + self
+        
+        def __rmul__(self,other):
+            return self*other
+            
+    Inner.__name__ = curvename + "_base"
+    return Inner
+
+class Isogeny(object):
+
+    isograph = DiGraph(weighted=True)
+    isomap = {}
+    
+    @classmethod
+    def generate(cls, fro, to):
+        path = cls.isograph.shortest_path(fro,to,by_weight=True)
+        if len(path):
+            iso = cls.isomap[(path[0], path[1])]
+            for i in xrange(1,len(path)-1):
+                iso = cls.isomap[(path[i],path[i+1])].compose(iso)
+            return iso
+        else:
+            return None
+    
+    def __init__(self,c1,c2,deg,fw,rv,check=True,dual=None,add=True):
+        self.c1 = c1
+        self.c2 = c2
+        self.fw = fw
+        self.rv = rv
+        self.deg = deg
+        
+        if add:
+            Isogeny.isomap[(c1,c2)] = self
+            Isogeny.isograph.add_edge(c1,c2,log(deg)/log(2) + 0.1)
+        
+        if dual is not None:
+            self.dual = dual
+        else:
+            self.dual = Isogeny(c2,c1,deg,rv,fw,False,self,add)
+        if not check: return
+        
+        
+        if not checkIsogenies: return
+        
+        debug_print("Checking isogeny %s <-%d-> %s..." % (c1.__name__,deg,c2.__name__))
+        if c2(*fw(*c1.basepoint)) != c2.basepoint:
+            raise Exception("Isogeny doesn't preserve basepoints")
+        if c1(*fw(*c2.basepoint)) != c1.basepoint:
+            raise Exception("Isogeny dual doesn't preserve basepoints")
+            
+        foo = c1()
+        bar = c2()
+        
+        c2(*fw(*foo))
+        c1(*rv(*bar))
+        
+        if c1(*rv(*c2(*fw(*foo)))) != deg*foo:
+            raise Exception("Isogeny degree is wrong")
+        if c2(*fw(*c1(*rv(*bar)))) != deg*bar:
+            raise Exception("Isogeny degree is wrong")
+        if -c2(*fw(*foo)) != c2(*fw(*(-foo))):
+            raise Exception("Isogeny uses wrong negmap")
+        if -c1(*rv(*bar)) != c1(*rv(*(-bar))):
+            raise Exception("Isogeny uses wrong negmap")
+            
+        
+    def __call__(self,ipt,**kwargs):
+        return self.c2(*self.fw(*ipt,**kwargs))
+        
+    def __repr__(self): return str(self)
+    def __str__(self):
+        out = "Isogeny %s%s <-%d-> %s%s..." %\
+            (self.c1.__name__,str(self.c1.params),self.deg,
+                self.c2.__name__,self.c2.params)
+        out += "\n  fw: %s" % str(self(self.c1()))
+        out += "\n  rv: %s" % str(self.dual(self.c2()))
+        return out
+        
+    def compose(self,other):
+        def fw(*args): return self.fw(*other.fw(*args))
+        def rv(*args): return other.rv(*self.rv(*args))
+        return Isogeny(other.c1,self.c2,self.deg*other.deg,fw,rv,False,None,False)
+
+def ec_family(defs,vars):
+    def inner1(CLS):
+        @memoize
+        def inner2(*args,**kwargs):
+            if len(args)==0 and len(kwargs)==0:
+                args = tuple(defs)
+                chk = True
+            else:
+                chk = False
+            
+            class ret(CLS,EcBase(CLS.__name__,vars,args)):
+                def __new__(cls,*args,**kwargs):
+                    return super(ret,cls).__new__(cls,*args,**kwargs)
+                
+            ret.__name__ = CLS.__name__
+            ret.basepoint = ret(*ret.get_basepoint())
+                
+            if chk: ret.check_group()
+            return ret
+            
+        inner2.__name__ = CLS.__name__ + "_family"
+        inner2()
+        return inner2
+        
+    return inner1
+
+#checkGroupLaws = checkTorsion = False
+
+@ec_family("ad","xy")
+class Edwards:
+    @classmethod
+    def invariants(cls,a,d,x,y):
+        return [y^2 + a*x^2 - 1 - d*x^2*y^2]
+        
+    def __neg__(self):
+        return self.__class__(-self.x,self.y)
+    
+    def __add__(self,other):
+        (x,y) = self
+        (X,Y) = other
+        a,d = self.params
+        dd = d*x*X*y*Y
+        return self.__class__((x*Y+X*y)/(1+dd),(y*Y-a*x*X)/(1-dd))
+
+    @classmethod
+    def get_basepoint(cls): return (0,1)
+
+    @classmethod
+    @memoize
+    def torsion(cls):
+        a,d = cls.params
+        sa = a.sqrt()
+        sd = d.sqrt()
+        sad = (a*d).sqrt()
+        def tor2_1((x,y)): return (-x,-y)
+        def tor4_1((x,y)): return (y/sa,-x*sa)
+        def tor4_2((x,y)): return (1/(sd*y),-1/(sd*x))
+        def tor2_2((x,y)): return (-1/(sad*x),-a/(sad*y))
+        
+        return [(tor2_1,2),(tor2_2,2),(tor4_1,4),(tor4_2,4)]
+
+@ec_family("eA","st")
+class JacobiQuartic:
+    @classmethod
+    def invariants(cls,e,A,s,t):
+        return [-t^2 + e*s^4 + 2*A*s^2 + 1]
+        
+    def __neg__(self):
+        return self.__class__(-self.s,self.t)
+    
+    def __add__(self,other):
+        (x,y) = self
+        (X,Y) = other
+        e,A = self.params
+        dd = e*(x*X)^2
+        YY = (1+dd)*(y*Y+2*A*x*X) + 2*e*x*X*(x^2+X^2)
+        return self.__class__((x*Y+X*y)/(1-dd),YY/(1-dd)^2)
+
+    @classmethod
+    def get_basepoint(cls): return (0,1)
+
+    @classmethod
+    @memoize
+    def torsion(cls):
+        e,A = cls.params
+        se = e.sqrt()
+        def tor2_1((s,t)): return (-s,-t)
+        def tor2_2((s,t)): return (1/(se*s),-t/(se*s^2))
+        return [(tor2_1,2),(tor2_2,2)]
+
+a,d = Idealized.vars("ad")
+def phi_iso(a,d):
+    return Isogeny(Edwards(a,d),JacobiQuartic(a^2,a-2*d),
+        2,
+        lambda x,y: (x/y, (2-y^2-a*x^2)/y^2),
+        lambda s,t: (2*s/(1+a*s^2), (1-a*s^2)/t)
+    )
+
+print phi_iso(a,d)
+print phi_iso(-a,d-a)
+
+print Isogeny.generate(Edwards(a,d),Edwards(-a,d-a))

test/idealized.sage

@@ -0,0 +1,273 @@
+class Unique(object):
+    def __init__(self,name):
+        self.name = name
+    
+    def __str__(self):
+        return self.name
+    
+    def __repr__(self):
+        return "Unique(\"%s\")" % self.name
+
+class Idealized(object):
+    UNION = ["UNION"]
+    
+    def __init__(self, R, idealMap = 0, vars = {}):
+        self.varnames = vars
+        if not isinstance(idealMap,dict):
+            idealMap = {()*R:idealMap}
+        self.idealMap = idealMap
+        self.R = R
+        self._sqrt = None
+        self._isqrt = None
+    
+    @staticmethod
+    def uvar(x):
+        return Idealized.var(Unique(x))
+    
+    @staticmethod
+    def var(x):
+        name = str(x)
+        R = PolynomialRing(QQ,[name])
+        rx = R.gens()[0]
+        return Idealized(R,rx,{x:(name,rx)})
+    
+    @staticmethod
+    def vars(xs):
+        return tuple((Idealized.var(x) for x in xs))
+    
+    @staticmethod
+    def uvars(xs):
+        return tuple((Idealized.uvar(x) for x in xs))
+    
+    def __str__(self):
+        def rep(I,x):
+            x = str(x)
+            gs = I.gens()
+            gs = [g for g in gs if g != 0]
+            if len(gs) == 0: return x
+            else:
+                g = ", ".join(["(%s)" % str(gen) for gen in gs])
+                return g + ": " + x
+        return "\n".join([rep(I,self.idealMap[I]) for I in self.idealMap])
+    
+    def __repr__(self):
+        # HACK!
+        if len(self.idealMap) == 0:
+            return "undef"
+        if len(self.idealMap) > 1:
+            return str(self)
+        for _,v in self.idealMap.iteritems():
+            return str(v)
+        
+    def prune(self):
+        self.idealMap = {I:v for I,v in self.idealMap.iteritems() if not (I*self.R).is_one()}
+        return self
+        
+    def __add__(self,other):
+        def f(x,y): return x+y
+        return self.op(other,f)
+        
+    def __radd__(self,other):
+        def f(x,y): return y+x
+        return self.op(other,f)
+        
+    def __rsub__(self,other):
+        def f(x,y): return y-x
+        return self.op(other,f)
+        
+    def __neg__(self):
+        def f(x,y): return y-x
+        return self.op(0,f)
+        
+    def __sub__(self,other):
+        def f(x,y): return x-y
+        return self.op(other,f)
+        
+    def is_square(self):
+        for _,v in self.idealMap.iteritems():
+            if not is_square(v): return False
+        return True
+        
+    def sqrt(self):
+        if self._sqrt is None:
+            s = Idealized.uvar("s")
+            self._sqrt = s.assuming(s^2 - self)
+        return self._sqrt
+        
+    def isqrt(self):
+        if self._isqrt is None:
+            s = Idealized.uvar("s")
+            z = Idealized(0).assuming(Self)
+            self._isqrt = s.assuming(s^2*self-1).union(z)
+        return self._isqrt
+        
+    def __mul__(self,other):
+        def f(x,y): return x*y
+        return self.op(other,f)
+        
+    def __rmul__(self,other):
+        def f(x,y): return y*x
+        return self.op(other,f)
+    
+    def __pow__(self,n):
+        if n < 0: return 1/self^(-n)
+        if n == 0: return 1
+        if n == 1: return self
+        if is_even(n): return (self*self)^(n//2)
+        if is_odd(n): return (self*self)^(n//2) * self
+        
+    def __div__(self,other):
+        def f(x,y): return x/y
+        return self.op(other,f)
+        
+    def __rdiv__(self,other):
+        def f(x,y): return y/x
+        return self.op(other,f)
+        
+    def union(self,other):
+        return self.op(other,Idealized.UNION)
+        
+    def __eq__(self,other):
+        return (self - other).is_zero()
+        
+    def __ne__(self,other):
+        return not (self==other)
+        
+    def __hash__(self):
+        return 0
+
+    def assume_zero(self):
+        out = {}
+        for I,J in self.idealMap.iteritems():
+            IJ = I+J.numerator()
+            if IJ.is_one(): continue
+            out[IJ] = self.R(0)
+        
+        if len(out) == 0:
+            raise Exception("Inconsistent assumption")
+        
+        return Idealized(self.R,out,self.varnames)
+    
+    def assuming(self,other):
+        return self + other.assume_zero()
+    
+    def is_zero(self):
+        for I,v in self.idealMap.iteritems():
+            if v.denominator() in I: return False
+            if v.numerator() not in I: return False
+        return True
+    
+    def op(self,other,f):
+        if not isinstance(other,Idealized):
+            other = Idealized(self.R,other,self.varnames)
+        
+        bad = False
+        for v in self.varnames:
+            if v not in other.varnames or self.varnames[v] != other.varnames[v]:
+                bad = True
+                break
+        for v in other.varnames:
+            if v not in self.varnames or self.varnames[v] != other.varnames[v]:
+                bad = True
+                break
+                
+        if bad:
+            def incrVar(v):
+                if v[-1] not in "0123456789": return v + "1"
+                elif v[-1] == 9: return incrVar(v[:-1]) + "0"
+                else: return v[:-1] + str(int(v[-1])+1)
+        
+            vars = {}
+            names = set()
+            for v,(name,_) in self.varnames.iteritems():
+                assert(name not in names)
+                names.add(name)
+                vars[v] = name
+            subMe = {n:n for n in names}
+            subThem = {}
+            for v,(name,_) in other.varnames.iteritems():
+                if v in self.varnames:
+                    subThem[name] = self.varnames[v][0]
+                else:
+                    oname = name
+                    while name in names:
+                        name = incrVar(name)
+                    names.add(name)
+                    subThem[oname] = name
+                    vars[v] = name
+            
+            R = PolynomialRing(QQ,sorted(list(names)),order="degrevlex")
+            gd = R.gens_dict()
+            subMe = {m:gd[n] for m,n in subMe.iteritems()}
+            subThem = {m:gd[n] for m,n in subThem.iteritems()}
+        
+            vars = {v:(n,gd[n]) for v,n in vars.iteritems()}
+        
+            def subIdeal(I,sub):
+                return [g(**sub) for g in I.gens()]*R
+            idealMe = {subIdeal(I,subMe):v(**subMe) for I,v in self.idealMap.iteritems()}
+            idealThem = {subIdeal(I,subThem):v(**subThem) for I,v in other.idealMap.iteritems()}
+        else:
+            R = self.R
+            idealMe = self.idealMap
+            idealThem = other.idealMap
+            vars = self.varnames
+        
+        def consist(I,x,y):
+            if (x-y).numerator() not in I:
+                raise Exception("Inconsistent: %s != %s in ideal %s" %
+                    (str(x),str(y),str(I)))
+            
+        out = {}
+        if f is Idealized.UNION:
+            for I,v in idealMe.iteritems():
+                if I in idealThem:
+                    consist(I,v,idealThem[I])
+                out[I] = v
+            for I,v in idealThem.iteritems():
+                if I in idealMe:
+                    consist(I,v,idealMe[I])
+                out[I] = v
+        
+        else:
+            for I,v in idealMe.iteritems():
+                if I in idealThem:
+                    x = f(v,idealThem[I])
+                    if I in out:
+                        consist(I,x,out[I])
+                    else: out[I] = x
+                else:
+                    for J,w in idealThem.iteritems():
+                        IJ = I+J
+                        if not IJ.is_one():
+                            x = f(v,w)
+                            if IJ in out:
+                                consist(IJ,x,out[IJ])
+                            else:
+                                out[IJ] = x
+        
+        def gb(I):
+            II = [0]*R
+            for g in I.gens():
+                if g not in II: II = II+[g]*R
+            return II
+
+        def red(I,v):
+            if I.is_zero(): return v
+            return I.reduce(R(v.numerator())) / I.reduce(R(v.denominator()))
+            
+        out = {gb(I):v for I,v in out.iteritems()}
+        out = {I:red(I,v) for I,v in out.iteritems()}
+        
+        return Idealized(R,out,vars)
+    
+    def reduce(self):
+        def red(I,v):
+            if I.is_zero(): return v
+            return I.reduce(R(v.numerator())) / I.reduce(R(v.denominator()))
+        out = {I:red(I,v) for I,v in self.idealMap.iteritems()}
+        return Idealized(self.R,out,self.vars)
+
+Idealized.INF = Idealized.uvar("inf")
+Idealized.ZOZ = Idealized.uvar("zoz")
+    

test/test_decaf.sage

@@ -0,0 +1,4 @@
+from ctypes import *
+
+decaf = CDLL("libdecaf.so")
+