On Nov 10, 2007 11:53 AM, Steffen <[EMAIL PROTECTED]> wrote: > Original: > Time: CPU 9.94 s, Wall: 10.09 s > Number of roots found: > 0 > Tonelli: > Time: CPU 0.00 s, Wall: 0.00 s > Number of roots found: > 0
Wow, that's an improvement! I've made this trac #1138 http://trac.sagemath.org/sage_trac/ticket/1138 > This shows only that the old implementation was inefficient for large > p. The results for the original implementation for 256 and 512 bit are > not contained in the following results, since sqrt(a) returned an > error for the 256 bit and 512 bit case. Here are the results for 100 > random values for a. Find x: x^2==a mod p > > Results for 128 bit prime: > Original: > Time: CPU 5.36 s, Wall: 5.41 s > Number of roots found: > 47 > Tonelli: > Time: CPU 0.07 s, Wall: 0.08 s > Number of roots found: > 47 > Results for 256 bit prime: > Tonelli: > Time: CPU 0.10 s, Wall: 0.10 s > Number of roots found: > 48 > Results for 512 bit prime: > Tonelli: > Time: CPU 0.30 s, Wall: 0.33 s > Number of roots found: > 50 > > > Here is the complete code I used: > > > def step3(b,p,r,x): > # Step 3: Find exponent > if GF(p)(b) == GF(p)(1): > return b,r,x,0 > m = 0 > while GF(p)(b**(2**m)) != 1: > m = m + 1 > if m == r: > return b,r,0,0 > return b,r,x,m > > def s_root(a,p): > # Step 0: Determine q: > q = 0 > e = 0 > while q % 2 != 1: > e = e+1 > q = (p-1) / 2**e > # Step 1: Find generator > n = ZZ.random_element() > while kronecker(n,p) != -1: > n = ZZ.random_element() > n = GF(p)(n) > z = GF(p)(n**q) > # Step 2: Initialize > y = z > r = e > a = GF(p)(a) > x = GF(p)(a**((q-1)/2)) > b = GF(p)(a*(x**2)) > x = GF(p)(a*x) > # Step 3: > b,r,x,m = step3(b,p,r,x) > # Step 4: Reduce exponent > while ZZ(m) != ZZ(0): > t = GF(p)(y**(2**(r-m-1))) > y = GF(p)(t**2) > r = GF(p)(m) > x = GF(p)(x*t) > b = GF(p)(b*y) > b,r,x,m = step3(b,p,r,x) > return x > > A = {} > > def go1(p): > X = {} > for i in range(0,len(A)): > X[i] = sqrt(GF(p)(A[i])) > return X > > def go2(p): > X ={} > for i in range(0,len(A)): > X[i] = s_root(A[i], p) > return X > > def analyseResult(A,X,p): > numberRoots = 0 > for i in range(0,len(A)): > x = X[i] > a = A[i] > if x in GF(p): > if (GF(p)(x) != GF(p)(0)) & (GF(p)(x**2) != GF(p)(a)): > print 'Error in implemenation' > elif GF(p)(x) != GF(p)(0): > numberRoots = numberRoots + 1 > print 'Number of roots found:' > print numberRoots > > n = 100 > > print 'Results for 128 bit prime:' > p = next_prime(2^(128)) > while p % 4 != 1: > p = next_prime(p+1) > > for i in range(0,n): > A[i] = GF(p).random_element() > > print 'Original:' > time X = go1(p) > analyseResult(A,X,p) > > print 'Tonelli:' > time X = go2(p) > analyseResult(A,X,p) > > print 'Results for 256 bit prime:' > p = next_prime(2^(256)) > while p % 4 != 1: > p = next_prime(p+1) > > for i in range(0,n): > A[i] = GF(p).random_element() > > #print 'Original:' > #time X = go1(p) > #analyseResult(A,X,p) > > print 'Tonelli:' > time X = go2(p) > analyseResult(A,X,p) > > print 'Results for 512 bit prime:' > p = next_prime(2^(512)) > while p % 4 != 1: > p = next_prime(p+1) > > for i in range(0,n): > A[i] = GF(p).random_element() > > #print 'Original:' > #time X = go1(p) > #analyseResult(A,X,p) > > print 'Tonelli:' > time X = go2(p) > analyseResult(A,X,p) > > -- Steffen > > > > > > -- William Stein Associate Professor of Mathematics University of Washington http://wstein.org --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~----------~----~----~----~------~----~------~--~---
