Can this be made cleaner and more efficient?
- To: mathgroup at smc.vnet.net
- Subject: [mg29807] Can this be made cleaner and more efficient?
- From: Flip at safebunch.com
- Date: Wed, 11 Jul 2001 01:26:56 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Hi All,
finally I got this algorithm working, but it is "ugly" and "slow" code.
Can an experienced Mathematica user help clean this up and make it more efficient.
By the way, this is an implementation of the Lucas Primality Test.
Notes:
1. Is it possible to move the sequence inside of the module and calculate it one
step at time? (it usually converges on the 2nd or 3rd try in finding D).
2. A check should be done on "n" and verify it is not a perfect odd sqaure
(i.e., n = 9 or n = 25, for example, will cause no solution for finding D and
this should really be checked). If n is a perfect square, terminate the test.
3. If one encounters a JacobiSymbol = 0 while searching for a d with
JacobiSymbol = -1, terminate the test as this means that n has factors (i.e., it
is composite).
Can someone here help? Thank you for any inputs ... Wilson
*** Here is the code ***
In[206]:=
Clear[a]
In[207]:=
a[1] =5;
In[208]:=
a[n_]:= a[n] = (-1)^(n+1) (Abs[a[n-1]]+2)
In[209]:=
(*find the first D in a[n] for which JS[a,n] = -1 *)
In[210]:=
findD[a_,b_]:=Module[{n=a, lst=b},
i=1;
While[JacobiSymbol[lst[[i]],n]
!= -1, i++];lst[[i]]]
In[211]:=
lpt1[a_,b_,c_]:=Module[{d=a, p=b, n=c},
Clear[u,v,k,w,r];
k = n+1;w = Reverse[IntegerDigits[k,2]];
r = Length[w]-2;q = (1-d)/4;
Print["n, k, r, d, q = ",{n,k,r,d,q}];
Print["w = ", w];t= 2n;
u = 1;v = p;
Print["u, v = ", {u,v}];
While[r > -1 ,
x = Mod[u*v, t];
y = Mod[(v^2+ d*u^2)/2, t];
{u,v}={x,y};
Print["u, v = ", {u,v}]; Print["r = ",r];
Print["Kr = ",w[[r+1]]];
If[w[[r+1]]==1,
u = Mod[(p*x+y)/2, t];
v = Mod[(p*y+ d*x)/2, t]];
(*Print["u, v = ", {u,v}]*);{u,v}={u,v};
r=r-1];
Print["u, v = ", {u,v}];
Print["For n = ",n];
If[Mod[u,n]==0, Print["LPT says Prime"],Print["LPT says Not Prime"]]]
In[212]:=
c=Table[a[n],{n,1,100}];
In[213]:=
p=1;
In[214]:=
n=5;
In[215]:=
d = findD[n,c]
Out[215]=
-7
In[216]:=
lpt1[d,p,n]