-rw-r--r-- 4341 high-ctidh-20210523/costisog.py
#!/usr/bin/env python3 M = 1 S = 1 from costpoly import tree1,multiprod2,multiprod2_selfreciprocal,multieval_precompute,multieval_postcompute x2 = S+S x2DBL = M+M+M+M # but save 1M if affine xDBL = x2+x2DBL xADD = M+M+S+S+M+M def isog(lmax,push,bsgs): bs,gs = bsgs assert lmax > 0 assert lmax%2 lbits = 6 while lmax>>lbits: lbits += 2 result = 0 assert bs >= 0 assert gs >= 0 if bs == 0 or gs == 0: result += xDBL result += ((lmax-5)//2)*xADD result += ((lmax-3)//2)*2*M result += push*(2*M) result += ((lmax-3)//2)*4*M*push else: # velusqrt multiples = set() multiples.add(2) for j in range(3,2*bs+1,2): multiples.add(j) for j in range(4,lmax+1-4*bs*gs,2): multiples.add(j) multiples.add(2*bs) multiples.add(4*bs) for j in range(6*bs,2*bs*(2*gs+1),4*bs): multiples.add(j) result += len(multiples)*xADD # actually xADD or xDBL result += tree1[bs] result += gs*(3*S+4*M) # biquad_curve result += push*(3*S+M) # biquad_precompute_point result += push*gs*6*M # biquad_postcompute_point result += push*multiprod2[gs] # pushing point result += 2*multiprod2_selfreciprocal[gs] # pushing curve result += multieval_precompute[bs][gs*2+1] # reciprocal of root of product tree result += 2*(push+1)*multieval_postcompute[bs][gs*2+1] # scaled remainder tree result += 2*M*(push+1)*(bs-1) # accumulating multieval results result += 2*M*(1+2*push)*((lmax-1)//2-2*bs*gs) # stray points at the end result += push*(S+S+M+M) # final point evaluation result += 2*(S+M+M+(S+S+M)*(lbits//2-1)) # powpow8mod return result def optimize(l,push): best,bestbsgs = isog(l,push,(0,0)),(0,0) # XXX: precompute more tree1 etc.; extend bs,gs limits for bs in range(2,33,2): for gs in range(1,2*bs+1): if 2*bs*gs > (l-1)//2: break if gs >= 32: break if bs > 3*gs: continue result = isog(l,push,(bs,gs)) if result < best: best,bestbsgs = result,(bs,gs) return best,bestbsgs def test1(): lmax = 587 bs = 14 gs = 10 push = 1 lbits = 6 while lmax>>lbits: lbits += 2 multiples = set() multiples.add(2) for j in range(3,2*bs+1,2): multiples.add(j) for j in range(4,lmax+1-4*bs*gs,2): multiples.add(j) multiples.add(2*bs) multiples.add(4*bs) for j in range(6*bs,2*bs*(2*gs+1),4*bs): multiples.add(j) assert len(multiples) == 37 assert tree1[bs] == 95 assert gs*(3*S+4*M)+push*(3*S+M)+push*gs*6*M == 134 assert push*multiprod2[gs] == 139 assert 2*multiprod2_selfreciprocal[gs] == 176 assert multieval_precompute[bs][2*gs+1] == 174 assert multieval_postcompute[bs][2*gs+1] == 251 assert 2*M*(push+1)*(bs-1) == 52 assert 2*M*(1+2*push)*((lmax-1)//2-2*bs*gs) == 78 assert push*(S+S+M+M)+2*(S+M+M+(S+S+M)*(lbits//2-1)) == 34 def test2(): assert optimize(587,1) == (2108,(14,10)) assert isog(587,1,(14,10)) == 2108 assert isog(587,1,(16,9)) == 2118 assert tree1[14] == 95 assert tree1[16] == 117 assert multiprod2[10] == 139 assert multiprod2[9] == 118 assert multiprod2_selfreciprocal[10] == 88 assert multiprod2_selfreciprocal[9] == 74 assert multieval_precompute[14][21] == 174 assert multieval_precompute[16][19] == 170 assert multieval_postcompute[14][21] == 251 assert multieval_postcompute[16][19] == 285 def test3(): primes = (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,587) batchsize = (3, 4, 5, 5, 6, 6, 6, 8, 7, 9, 10, 5) batchbound = (13, 19, 20, 20, 20, 20, 20, 20, 17, 17, 17, 5) B = len(batchsize) batchstart = [sum(batchsize[:j]) for j in range(B)] batchstop = [sum(batchsize[:j+1]) for j in range(B)] bsgs = [optimize(primes[batchstart[b]],1)[1] for b in range(B)] C = [[isog(primes[batchstop[b]-1],push,bsgs[b]) for b in range(B)] for push in (0,1,2)] assert C[0] == [34, 82, 170, 250, 368, 412, 532, 653, 720, 875, 1034, 1860, ] assert C[1] == [48, 120, 252, 372, 546, 643, 830, 1031, 1138, 1385, 1644, 2898, ] assert C[2] == [62, 158, 334, 494, 724, 874, 1128, 1409, 1556, 1895, 2254, 3936, ] if __name__ == '__main__': test1() test2() test3()