source: git/Singular/LIB/aksaka.lib @ d2f488

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