Date: Mon, 27 Jan 1997 00:52:32 -0800 From: hbaker@netcom.com (Henry Baker) To: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol In article <5ch4g9$fbh@gannett.math.niu.edu>, rusin@vesuvius.math.niu.edu (Dave Rusin) wrote: > In article , > Henry Baker wrote: > >Q: Are there any algorithms for which the speed of evaluation of the > >Jacobi (a/p) symbol is a bottleneck/on the critical path/in the inner loop?? > > I've been doing some calculations on elliptic curves which make frequent > use of maple's Legendre function L(a,b). (I also make calls of the > form isolve(a*x^2+b*y^2+c*z^2=0) which I think use L ). > Using quadratic reciprocity, it's fairly easy to compute things like > this _if_ you can factor the arguments into prime powers. > But try for example > a:=expand(nextprime(10^21)*nextprime(10^31)): > b:=expand(nextprime(10^20)*nextprime(10^30)): > L(a,b): > This takes quite long... Of course, the > actual source of the problem is the slow speed of factorization. I don't know the implementation of maple's Legendre function, but if it is the same Legendre/Jacobi function that I know, it does _not_ require factorization, and should be of the same or slightly easier complexity than gcd(p,q). If maple is factorizing, then something is dreadfully wrong... ============================================================================== Date: Mon, 27 Jan 1997 01:01:19 -0800 From: hbaker@netcom.com (Henry Baker) To: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol Newsgroups: sci.crypt,sci.math,sci.math.symbolic (A copy of this message has also been posted to the following newsgroups: sci.crypt, sci.math,sci.math.symbolic) In article <5ch4g9$fbh@gannett.math.niu.edu>, rusin@vesuvius.math.niu.edu (Dave Rusin) wrote: > In article , > Henry Baker wrote: > >Q: Are there any algorithms for which the speed of evaluation of the > >Jacobi (a/p) symbol is a bottleneck/on the critical path/in the inner loop?? > > I've been doing some calculations on elliptic curves which make frequent > use of maple's Legendre function L(a,b). (I also make calls of the > form isolve(a*x^2+b*y^2+c*z^2=0) which I think use L ). > Using quadratic reciprocity, it's fairly easy to compute things like > this _if_ you can factor the arguments into prime powers. > But try for example > a:=expand(nextprime(10^21)*nextprime(10^31)): > b:=expand(nextprime(10^20)*nextprime(10^30)): > L(a,b): > This takes quite long... Of course, the > actual source of the problem is the slow speed of factorization. Something is dreadfully wrong here. I tried exactly the same example on PARI/GP on my Mac, and computing a took 2.7 seconds computing b took 2.0 seconds computing kro(a,b) took 16 msec. (kro is pari/gp's name for legendre/jacobi) computing kro(b,a) took 16 msec. (answer was -1 in both cases). ============================================================================== Date: Mon, 27 Jan 97 08:29:28 CST From: rusin (Dave Rusin) To: hbaker@netcom.com Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol This is very interesting. When I noticed maple having trouble with some calls to L() and isolve() I factored the arguments and noticed (1) comparable delays and (2) large factors. Since my gut reaction to computing these things myself would be to work mod p for each relevant p, I sort of assumed maple was doing the same and taking time to factor. (For the numbers given in the post, this is essentially an interminable task since I think maple used the pollard-rho as its default factorizer.) Certainly after 15 minutes this morning L(a,b) has not terminated. Perhaps I should switch to pari for this purpose. Now you've piqued my curiosity. 1. How _does_ one compute these symbols quickly? Is there a reciprocity formula which does not require prime factorization? 2. What prompted you to ask about speed if you knew about pari's speed but were unaware of maple's sluggishness? Inquiring minds etc etc dave ============================================================================== From: hbaker@netcom.com (Henry G. Baker) Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol To: rusin@math.niu.edu (Dave Rusin) Date: Mon, 27 Jan 1997 07:59:22 -0800 (PST) > This is very interesting. When I noticed maple having trouble with some > calls to L() and isolve() I factored the arguments and noticed > (1) comparable delays and (2) large factors. Since my gut reaction to > computing these things myself would be to work mod p for each relevant p, > I sort of assumed maple was doing the same and taking time to factor. > (For the numbers given in the post, this is essentially an interminable > task since I think maple used the pollard-rho as its default factorizer.) > Certainly after 15 minutes this morning L(a,b) has not terminated. > Perhaps I should switch to pari for this purpose. pari is pretty good, but only within its limited scope -- it is _not_ a symbolic algebra program, although it is very good with integers, rationals, polynomials, etc. I have found at least one bug -- gcd() doesn't work correctly over the Gaussian integers (and probably not over other algebraic integral domains, either, I suspect). > Now you've piqued my curiosity. > > 1. How _does_ one compute these symbols quickly? Is there a reciprocity > formula which does not require prime factorization? Jacobi blithely goes ahead to formally compute (a/p), even when p isn't prime. This ruins the meaning of (a/p) for non-primes p, but preserves them for primes p. The rules: Jac(a,n): n must be odd, a,n must be relatively prime ( (a,n)=1 ). Jac(1,n) = 1 Jac(a,n) = Jac(a mod n,n) Jac(-a,n) = Jac(-1,n) Jac(a,n) = (-1)^((n-1)/2) Jac(a,n) Jac(2a,n) = Jac(2,n) Jac(a,n) = (-1)^((n^2-1)/8) Jac(a,n) Jac(a,n) = (-1)^((a-1)(n-1)/4) Jac(n,a) **** a must be odd **** = (-1)^((a-1)(n-1)/4) Jac(n mod a,a) **** a must be odd **** So Jac(a,n) is roughly equivalent to a gcd calculation in complexity. Lehmer made good use of this--I think he (they?) had a number of articles about it down through the years. Don't ask me for references, though. Also, see Knuth, V.II. > 2. What prompted you to ask about speed if you knew about pari's speed > but were unaware of maple's sluggishness? I'm writing a brief article about it for ACM Sigplan Notices. > Inquiring minds etc etc > dave -- Henry Baker www/ftp directory: ftp.netcom.com:/pub/hb/hbaker/home.html ============================================================================== Date: Mon, 27 Jan 97 15:21:08 CST From: rusin (Dave Rusin) To: hbaker@netcom.com Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol >(kro is pari/gp's name for legendre/jacobi) I guess this is why maple is factoring. Pari's kro(2,15)=+1, which is as you point out easy to compute but not nearly as informative as maple's L(2,15)=-1, which is the information I need for my purposes. dave ============================================================================== From: hbaker@netcom.com (Henry G. Baker) Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol To: rusin@math.niu.edu (Dave Rusin) Date: Mon, 27 Jan 1997 14:39:42 -0800 (PST) > >(kro is pari/gp's name for legendre/jacobi) > I guess this is why maple is factoring. Pari's kro(2,15)=+1, which is as > you point out easy to compute but not nearly as informative as maple's > L(2,15)=-1, which is the information I need for my purposes. > > dave Riesel's book points out exactly this example: "However, since [(ab/n)=(a/n)(b/n)] always holds, the Jacobi symbol (a/n) occasionally assumes, in such a case, the value +1 for a quadratic non-residue a (mod n). As an example we give (2/15)=(2/3)(2/5)=(-1)(-1)=+1, although 2 is a quadratic non-residue (mod 15)." Of course, if true quadratic residues were this easy, so would factoring! Maple's problem is in not clearly distinguishing between their function L(a,b) and Legendre's symbol (a/b), which is defined _only_ for prime b. Jacobi's symbol is the generalization of (a/b) to all odd b, but it loses some of the usefulness for non-prime b, since it no longer correctly identifies non-residues. -- Henry Baker www/ftp directory: ftp.netcom.com:/pub/hb/hbaker/home.html ============================================================================== Date: Tue, 28 Jan 1997 08:17:58 -0600 From: Daniel Lichtblau To: Henry Baker , rusin@vesuvius.math.niu.edu Subject: Re: Q: Jacobi (Legendre/Kronecker) Symbol Henry Baker wrote: > > In article <5ch4g9$fbh@gannett.math.niu.edu>, rusin@vesuvius.math.niu.edu > (Dave Rusin) wrote: > > > In article , > > Henry Baker wrote: > > >Q: Are there any algorithms for which the speed of evaluation of the > > >Jacobi (a/p) symbol is a bottleneck/on the critical path/in the inner loop?? > > > > I've been doing some calculations on elliptic curves which make frequent > > use of maple's Legendre function L(a,b). (I also make calls of the > > form isolve(a*x^2+b*y^2+c*z^2=0) which I think use L ). > > Using quadratic reciprocity, it's fairly easy to compute things like > > this _if_ you can factor the arguments into prime powers. > > But try for example > > a:=expand(nextprime(10^21)*nextprime(10^31)): > > b:=expand(nextprime(10^20)*nextprime(10^30)): > > L(a,b): > > This takes quite long... Of course, the > > actual source of the problem is the slow speed of factorization. > > Something is dreadfully wrong here. I tried exactly the same example > on PARI/GP on my Mac, and > > computing a took 2.7 seconds > computing b took 2.0 seconds > computing kro(a,b) took 16 msec. (kro is pari/gp's name for legendre/jacobi) > computing kro(b,a) took 16 msec. > (answer was -1 in both cases). I tried this, and it is indeed unrelated to integer factorization. A glance at the Mathematica code indicates it is done mostly via mod and bit-shifting. Not unlike integer gcd algorithms. I expect it will be fast up to several thousands of digits; a timing for arguments near 500 digits was still a fraction of a second (granted, a Pentium Pro 200 is a faster machine than I deserve). Daniel Lichtblau Wolfram Research In[2]:= ??JacobiSymbol JacobiSymbol[n, m] gives the Jacobi symbol, an integer function of n and m. Attributes[JacobiSymbol] = {Listable, Protected} In[3]:= np[a_Integer] := Module[{j=a}, While[!PrimeQ[j], j++]; j] In[4]:= np[10^21] Out[4]= 1000000000000000000117 In[5]:= x1 = %; In[6]:= x2 = np[10^31]; In[7]:= y1 = np[10^20]; In[8]:= y2 = np[10^30]; In[9]:= e1 = x1*x2 Out[9]= 10000000000000000001170000000033000000000000000003861 In[10]:= e2 = y1*y2 Out[10]= 100000000000000000039000000005700000000000000002223 In[11]:= Timing[JacobiSymbol[e1,e2]] Out[11]= {0. Second, -1} In[12]:= Timing[JacobiSymbol[e2,e1]] Out[12]= {0. Second, -1} In[21]:= a1 = Random[Integer, 2^(16*100)]; In[22]:= a2 = Random[Integer, 2^(16*100)]; In[23]:= EvenQ[a1] Out[23]= True In[24]:= a1 += 1; In[25]:= EvenQ[a2] Out[25]= False In[26]:= Log[10., a1] Out[26]= 481.157 In[27]:= Log[10., a2] Out[27]= 480.792 In[28]:= Timing[JacobiSymbol[a1,a2]] Out[28]= {0.02 Second, 1}