source: git/Singular/LIB/aksaka.lib @ 4c20ee

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