# source:git/Singular/LIB/aksaka.lib@0dd77c2

spielwiese
Last change on this file since 0dd77c2 was 0dd77c2, checked in by Hans Schoenemann <hannes@…>, 13 years ago
• Property mode set to `100644`
File size: 12.1 KB
Line
2///////////////////////////////////////////////////////////////////////////////
3version="\$Id\$";
4category="Teaching";
5info="
6LIBRARY: aksaka.lib     Procedures for primality testing after Agrawal, Saxena, Kayal
7AUTHORS: Christoph Mang
8
9OVERVIEW:
10 Algorithms for primality testing in polynomial time
11 based on the ideas of Agrawal, Saxena and  Kayal.
12
13PROCEDURES:
14
15fastExpt(a,m,n)         a^m for numbers a,m; if a^k>n n+1 is returned
16log2(n)                    logarithm to basis 2 of n
17PerfectPowerTest(n)        checks if there are a,b>1, so that a^b=n
18wurzel(r)                  square root of number r
19euler(r)                   phi-function of Euler
20coeffmod(f,n)              polynomial f modulo number n (coefficients mod n)
21powerpolyX(q,n,a,r)        (polynomial a)^q modulo (poly r,number n)
23";
24
25LIB "crypto.lib";
26LIB "ntsolve.lib";
27
28///////////////////////////////////////////////////////////////
29//                                                           //
30//   FAST (MODULAR) EXPONENTIATION                           //
31//                                                           //
32///////////////////////////////////////////////////////////////
33proc fastExpt(number a,number m,number n)
34"USAGE: fastExpt(a,m,n); a, m, n = number;
35RETURN: the m-th power of a; if a^m>n the procedure returns n+1
36NOTE:   uses fast exponentiation
37EXAMPLE:example fastExpt; shows an example
38"
39{
40  number b,c,d;
41  c=m;
42  b=a;
43  d=1;
44  while(c>=1)
45  {
46    if(b>n)
47    {
48      return(n+1);
49    }
50    if((c mod 2)==1)
51    {
52      d=d*b;
53      if(d>n)
54      {
55        return(n+1);
56      }
57    }
58    b=b^2;
59    c=intPart(c/2);
60  }
61  return(d)
62}
63example
64{ "EXAMPLE:"; echo = 2;
65   ring R = 0,x,dp;
66   fastExpt(2,10,1022);
67}
68////////////////////////////////////////////////////////////////////////////
69proc coeffmod(poly f,number n)
70"USAGE: coeffmod(f,n);
71ASSUME: poly f depends on at most var(1) of the basering
72RETURN: poly f modulo number n
73NOTE:   at first the coefficients of the monomials of the polynomial f are
74        determined, then their remainder modulo n is calculated,
75        after that the polynomial 'is put together' again
76EXAMPLE:example coeffmod; shows an example
77"
78{
79  matrix M=coeffs(f,var(1));         //Bestimmung der Polynomkoeffizienten
80  int i=1;
81  poly g;
82  int o=nrows(M);
83
84  while(i<=o)
85  {
87                   // umwandeln der Koeffizienten in numbers,
88                   // Berechnung der Reste dieser numbers modulo n,
89                   // Polynom mit neuen Koeffizienten wieder zusammensetzen
90    i++;
91  }
92  return(g);
93}
94example
95{ "EXAMPLE:"; echo = 2;
96   ring R = 0,x,dp;
97   poly f=2457*x4+52345*x3-98*x2+5;
98   number n=3;
99   coeffmod(f,n);
100}
101//////////////////////////////////////////////////////////////////////////
102proc powerpolyX(number q,number n,poly a,poly r)
103"USAGE: powerpolyX(q,n,a,r);
104RETURN: the q-th power of poly a modulo poly r and number n
105EXAMPLE:example powerpolyX; shows an example
106"
107{
108  ideal I=r;
109
110  if(q==1){return(coeffmod(reduce(a,I),n));}
111  if((q mod 5)==0){return(coeffmod(reduce(powerpolyX(q/5,n,a,r)^5,I),n));}
112  if((q mod 4)==0){return(coeffmod(reduce(powerpolyX(q/4,n,a,r)^4,I),n));}
113  if((q mod 3)==0){return(coeffmod(reduce(powerpolyX(q/3,n,a,r)^3,I),n));}
114  if((q mod 2)==0){return(coeffmod(reduce(powerpolyX(q/2,n,a,r)^2,I),n));}
115
116  return(coeffmod(reduce(a*powerpolyX(q-1,n,a,r),I),n));
117}
118example
119{ "EXAMPLE:"; echo = 2;
120   ring R=0,x,dp;
121   poly a=3*x3-x2+5;
122   poly r=x7-1;
123   number q=123;
124   number n=5;
125   powerpolyX(q,n,a,r);
126}
127///////////////////////////////////////////////////////////////
128//                                                           //
129//   GENERAL PROCEDURES                                      //
130//                                                           //
131///////////////////////////////////////////////////////////////
132proc log2(number x)
133"USAGE:  log2(x);
134RETURN:  logarithm to basis 2 of x
135NOTE:    calculates the natural logarithm of x with a power-series
136         of the ln, then the basis is changed to 2
137EXAMPLE: example log2; shows an example
138"
139{
140  number b,c,d,t,l;
141  int k;
142                                         // log2=logarithmus zur basis 2,
143                                         // log=natuerlicher logarithmus
144  b=100000000000000000000000000000000000000000000000000;
145  c=141421356237309504880168872420969807856967187537695;    // c/b=sqrt(2)
146  d=69314718055994530941723212145817656807550013436026;     // d/b=log(2)
147
148  //bringen x zunaechst zwischen 1/sqrt(2) und sqrt(2), so dass Reihe schnell
149  //konvergiert, berechnen dann Reihe bis 30. Summanden und letztendlich
150  //log2(x)=log(x)/log(2)=(log(x/2^j)+j*log(2))/log(2)  fuer grosse x
151  //log2(x)=log(x)/log(2)=(log(x*2^j)-j*log(2))/log(2)  fuer kleine x
152
153  number j=0;
154  if(x<(b/c))
155  {
156    while(x<(b/c))
157    {
158      x=x*2;
159      j=j+1;
160    }
161    t=(x-1)/(x+1);
162    k=0;
163    l=0;
164    while(k<30)       //fuer  x*2^j wird Reihe bis k-tem Summanden berechnet
165    {
166      l=l+2*(t^(2*k+1))/(2*k+1);      //l=log(x*2^j) nach k Summanden
167      k=k+1;
168    }
169    return((l*b/d)-j);         //log2(x)=log(x*2^j)/log(2)-j wird ausgegeben
170  }
171  while(x>(c/b))
172  {
173    x=x/2;
174    j=j+1;
175  }
176  t=(x-1)/(x+1);
177  k=0;
178  l=0;
179  while(k<30)         //fuer  x/2^j wird Reihe bis k-tem Summanden berechnet
180  {
181    l=l+2*(t^(2*k+1))/(2*k+1);       //l=log(x/2^j) nach k Summanden
182    k=k+1;
183  }
184  return((l*b/d)+j);     //hier wird log2(x)=log(x/2^j)/log(2)+j  ausgegeben
185}
186example
187{ "EXAMPLE:"; echo = 2;
188   ring R = 0,x,dp;
189   log2(1024);
190}
191//////////////////////////////////////////////////////////////////////////
192proc wurzel(number r)
193"USAGE: wurzel(r);
194ASSUME: characteristic of basering is 0, r>=0
195RETURN: number, square root of r
196EXAMPLE:example wurzel; shows an example
197"
198{
199  poly f=var(1)^2-r;             //Wurzel wird als Nullstelle eines Polys
200                                 //mit proc nt_solve genaehert
201  def B=basering;
202  ring R=(real,40),var(1),dp;
203  poly g=imap(B,f);
204  list l=nt_solve(g,1.1);
206  setring B;
207  return(imap(R,m));
208}
209example
210{ "EXAMPLE:"; echo = 2;
211   ring R = 0,x,dp;
212   wurzel(7629412809180100);
213}
214//////////////////////////////////////////////////////////////////////////
215proc euler(number r)
216"USAGE: euler(r);
217RETURN: number phi(r), where phi is Eulers phi-function
218NOTE:   first r is factorized with proc PollardRho, then phi(r) is
219        calculated with the help of phi(p) of every factor p;
220EXAMPLE:example euler; shows an example
221"
222{
223  list l=PollardRho(r,5000,1);           //bestimmen der Primfaktoren von r
224  int k;
225  number phi=r;
226  for(k=1;k<=size(l);k++)
227  {
228    phi=phi-phi/l[k];       //berechnen phi(r) mit Hilfe der
229  }                         //Primfaktoren und Eigenschaften der phi-Fktn
230  return(phi);
231}
232example
233{ "EXAMPLE:"; echo = 2;
234   ring R = 0,x,dp;
235   euler(99991);
236}
237///////////////////////////////////////////////////////////////
238//                                                           //
239//   PERFECT POWER TEST                                      //
240//                                                           //
241///////////////////////////////////////////////////////////////
242proc PerfectPowerTest(number n)
243"USAGE: PerfectPowerTest(n);
244RETURN: 0 if there are numbers a,b>1 with a^b=n;
245        1 if there are no numbers a,b>1 with a^b=n;
246        if printlevel>=1 and there are a,b>1 with a^b=n,
247        then a,b are printed
248EXAMPLE:example PerfectPowerTest; shows an example
249"
250{
251  number a,b,c,e,f,m,l,p;
252  b=2;
253  l=log2(n);
254  e=0;
255  f=1;
256
257  while(b<=l)
258  {
259    a=1;
260    c=n;
261
262    while(c-a>=2)
263    {
264      m=intPart((a+c)/2);
265      p=fastExpt(m,b,n);
266
267      if(p==n)
268      {
269        if(printlevel>=1){"es existieren Zahlen a,b>1 mit a^b=n";
270                      "a=";m;"b=";b;pause();}
271        return(e);
272      }
273
274      if(p<n)
275      {
276        a=m;
277      }
278      else
279      {
280        c=m;
281      }
282    }
283    b=b+1;
284  }
285  if(printlevel>=1){"es existieren keine Zahlen a,b>1  mit a^b=n";pause();}
286  return(f);
287}
288example
289{ "EXAMPLE:"; echo = 2;
290   ring R = 0,x,dp;
291   PerfectPowerTest(887503681);
292}
293///////////////////////////////////////////////////////////////
294//                                                           //
296//                                                           //
297///////////////////////////////////////////////////////////////
300ASSUME: n>1
301RETURN: 0 if n is composite;
302        1 if n is prime;
303        if printlevel>=1, you are informed what the procedure will do
304        or has calculated
305NOTE:   ASK-algorithm; uses proc powerpolyX for step 5
307"
308{
309  number c,p;
310  c=0;
311  p=1;
312
313  if(n==2) { return(p); }
314  if(printlevel>=1){"Beginn: Test ob a^b=n fuer a,b>1";pause();}
315  if(0==PerfectPowerTest(n))             // Schritt1 des ASK-Alg.
316  { return(c); }
317  if(printlevel>=1)
318  {
319    "Beginn: Berechnung von r, so dass ord(n) in Z/rZ groesser log2(n)^2 ";
320    pause();
321  }
322  number k,M,M2,b;
323  int r;                          // Zeile 526-542: Schritt 2
324  M=log2(n);                      // darin sind die Schritte
325  M2=intPart(M^2);                // 3 und 4 integriert
326  r=2;
327
328  if(gcdN(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
329                    // Schritt 3 des ASK-Alg.
330  {
331    if(printlevel>=1){"Zahl r gefunden mit r|n";"r=";r;pause();}
332    return(c);
333  }
334  if(r==n-1)              // falls diese Bedingung erfuellt ist, ist
335                         // ggT(a,n)=1 fuer 0<a<=r, schritt 4 des ASK-Alg.
336  {
337    if(printlevel>=1){"ggt(r,n)=1 fuer alle 1<r<n";pause();}
338    return(p);
339  }
340  k=1;
341  b=1;
342
343  while(k<=M2)         //Beginn des Ordnungstests fuer aktuelles r
344  {
345    b=((b*n) mod r);
346    if(b==1)   //tritt Bedingung ein so gilt fuer aktuelles r,k:
347               //n^k=1 mod r
348               //tritt Bedingung ein, wird wegen k<=M2=intPart(log2(n)^2)
349               //r erhoeht,k=0 gesetzt und Ordnung erneut untersucht
350               //tritt diese Bedingung beim groessten k=intPart(log2(n)^2)
351               //nicht ein, so ist ord_r_(n)>log2(n)^2 und Schleife
352               //wird nicht mehr ausgefuehrt
353    {
354      if(r==2147483647)
355      {
356        string e="error: r ueberschreitet integer Bereich und darf
357                  nicht als Grad eines Polynoms verwendet werden";
358        return(e);
359      }
360      r=r+1;
361      if(gcdN(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
362                         //wird aufgrund von Schritt 3 des ASK-Alg. fuer
363                         //jeden Kandidaten r getestet
364      {
365        if(printlevel>=1){"Zahl r gefunden mit r|n";"r=";r;pause();}
366        return(c);
367      }
368      if(r==n-1)         //falls diese Bedingung erfuellt ist,
369                         //ist n wegen Zeile 571 prim
370                         //wird aufgrund von Schritt 4 des ASK-Alg. fuer
371                         //jeden Kandidaten r getestet
372      {
373        if(printlevel>=1){"ggt(r,n)=1 fuer alle 1<r<n";pause();}
374        return(p);
375      }
376
377      k=0;               //fuer neuen Kandidaten r, muss k fuer den
378                         //Ordnungstest zurueckgesetzt werden
379    }
380    k=k+1;
381  }
382  if(printlevel>=1)
383  {
384    "Zahl r gefunden, so dass ord(n) in Z/rZ groesser log2(n)^2 ";
385    "r=";r;pause();
386  }
387  number l=1;                         // Zeile 603-628: Schritt 5
388  poly f,g;                           // Zeile 618 ueberprueft Gleichungen fuer
389                                      // das jeweilige a
390  g=var(1)^r-1;
391  number grenze=intPart(wurzel(euler(r))*M);
392
393  if(printlevel>=1)
394  {"Beginn: Ueberpruefung ob (x+a)^n kongruent x^n+a mod (x^r-1,n)
395        fuer 0<a<=intPart(wurzel(phi(r))*log2(n))=";grenze;pause();
396  }
397  while(l<=grenze)          //
398  {
399    if(printlevel>=1){"a=";l;}
400    f=var(1)+l;
401    if(powerpolyX(n,n,f,g)!=(var(1)^(int((n mod r)))+l))
402    {
403      if(printlevel>=1)
404      {"Fuer a=";l;"ist (x+a)^n nicht kongruent x^n+a mod (x^r-1,n)";
405        pause();
406      }
407      return(c);
408    }
409    l=l+1;
410  }
411  if(printlevel>=1)
412  {"(x+a)^n kongruent x^n+a mod (x^r-1,n) fuer 0<a<wurzel(phi(r))*log2(n)";
413    pause();
414  }
415  return(p);       // Schritt 6
416}
417example
418{ "EXAMPLE:"; echo = 2;
419   ring R = 0,x,dp;