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

fieker-DuValspielwiese
Last change on this file since a96a75 was 6391eb, checked in by Hans Schoenemann <hannes@…>, 5 years ago
version numbers
  • Property mode set to 100644
File size: 12.1 KB
Line 
1////////////////////////////////////////////////////////////////////////////////
2version="version aksaka.lib 4.1.2.0 Feb_2019 "; // $Id$
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(bigint a,bigint m,bigint 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  bigint 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=c div 2;
59  }
60  return(d)
61}
62example
63{ "EXAMPLE:"; echo = 2;
64   fastExpt(2,10,1022);
65}
66////////////////////////////////////////////////////////////////////////////
67proc coeffmod(poly f,bigint n)
68"USAGE: coeffmod(f,n);
69ASSUME: poly f depends on at most var(1) of the basering
70RETURN: poly f modulo number n
71NOTE:   at first the coefficients of the monomials of the polynomial f are
72        determined, then their remainder modulo n is calculated,
73        after that the polynomial 'is put together' again
74EXAMPLE:example coeffmod; shows an example
75"
76{
77  matrix M=coeffs(f,var(1));         //Bestimmung der Polynomkoeffizienten
78  int i=1;
79  poly g;
80  int o=nrows(M);
81
82  while(i<=o)
83  {
84    g=g+(bigint(leadcoef(M[i,1])) mod n)*var(1)^(i-1) ;
85                   // umwandeln der Koeffizienten in bigint,
86                   // Berechnung der Reste dieser bigint modulo n,
87                   // Polynom mit neuen Koeffizienten wieder zusammensetzen
88    i++;
89  }
90  return(g);
91}
92example
93{ "EXAMPLE:"; echo = 2;
94   ring R = 0,x,dp;
95   poly f=2457*x4+52345*x3-98*x2+5;
96   bigint n=3;
97   coeffmod(f,n);
98}
99//////////////////////////////////////////////////////////////////////////
100proc powerpolyX(bigint q,bigint n,poly a,poly r)
101"USAGE: powerpolyX(q,n,a,r);
102RETURN: the q-th power of poly a modulo poly r and number n
103EXAMPLE:example powerpolyX; shows an example
104"
105{
106  ideal I=r;
107
108  if(q==1){return(coeffmod(reduce(a,I),n));}
109  if((q mod 5)==0){return(coeffmod(reduce(powerpolyX(q div 5,n,a,r)^5,I),n));}
110  if((q mod 4)==0){return(coeffmod(reduce(powerpolyX(q div 4,n,a,r)^4,I),n));}
111  if((q mod 3)==0){return(coeffmod(reduce(powerpolyX(q div 3,n,a,r)^3,I),n));}
112  if((q mod 2)==0){return(coeffmod(reduce(powerpolyX(q div 2,n,a,r)^2,I),n));}
113
114  return(coeffmod(reduce(a*powerpolyX(q-1,n,a,r),I),n));
115}
116example
117{ "EXAMPLE:"; echo = 2;
118   ring R=0,x,dp;
119   poly a=3*x3-x2+5;
120   poly r=x7-1;
121   bigint q=123;
122   bigint n=5;
123   powerpolyX(q,n,a,r);
124}
125///////////////////////////////////////////////////////////////
126//                                                           //
127//   GENERAL PROCEDURES                                      //
128//                                                           //
129///////////////////////////////////////////////////////////////
130proc log2(bigint xx)
131"USAGE:  log2(x);
132RETURN:  logarithm to basis 2 of x
133NOTE:    calculates the natural logarithm of x with a power-series
134         of the ln, then the basis is changed to 2
135EXAMPLE: example log2; shows an example
136"
137{
138  ring r=0,@x,dp;
139  number b,c,d,t,l,x;
140  x=number(xx);
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 (intPart((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(intPart(((l*b/d)+j)));     //hier wird log2(x)=log(x/2^j)/log(2)+j  ausgegeben
185}
186example
187{ "EXAMPLE:"; echo = 2;
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(bigint r)
215"USAGE: euler(r);
216RETURN: bigint 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  bigint phi=r;
225  for(k=1;k<=size(l);k++)
226  {
227    phi=phi-phi div 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   euler(99991);
234}
235///////////////////////////////////////////////////////////////
236//                                                           //
237//   PERFECT POWER TEST                                      //
238//                                                           //
239///////////////////////////////////////////////////////////////
240proc PerfectPowerTest(bigint n)
241"USAGE: PerfectPowerTest(n);
242RETURN: 0 if there are numbers a,b>1 with a^b=n;
243        1 if there are no numbers a,b>1 with a^b=n;
244        if printlevel>=1 and there are a,b>1 with a^b=n,
245        then a,b are printed
246EXAMPLE:example PerfectPowerTest; shows an example
247"
248{
249  bigint a,b,c,e,f,m,l,p;
250  b=2;
251  l=log2(n);
252  e=0;
253  f=1;
254
255  while(b<=l)
256  {
257    a=1;
258    c=n;
259
260    while(c-a>=2)
261    {
262      m=(a+c) div 2;
263      p=fastExpt(m,b,n);
264
265      if(p==n)
266      {
267        if(printlevel>=1){"es existieren Zahlen a,b>1 mit a^b=n";
268                      "a=";m;"b=";b;pause();}
269        return(e);
270      }
271
272      if(p<n)
273      {
274        a=m;
275      }
276      else
277      {
278        c=m;
279      }
280    }
281    b=b+1;
282  }
283  if(printlevel>=1){"es existieren keine Zahlen a,b>1  mit a^b=n";pause();}
284  return(f);
285}
286example
287{ "EXAMPLE:"; echo = 2;
288   PerfectPowerTest(887503681);
289}
290///////////////////////////////////////////////////////////////
291//                                                           //
292//   ASK PRIMALITY TEST                                      //
293//                                                           //
294///////////////////////////////////////////////////////////////
295proc ask(bigint n)
296"USAGE: ask(n);
297ASSUME: n>1
298RETURN: 0 if n is composite;
299        1 if n is prime;
300        if printlevel>=1, you are informed what the procedure will do
301        or has calculated
302NOTE:   ASK-algorithm; uses proc powerpolyX for step 5
303EXAMPLE:example ask; shows an example
304"
305{
306  bigint c,p;
307  c=0;
308  p=1;
309
310  if(n==2) { return(p); }
311  if(printlevel>=1){"Beginn: Test ob a^b=n fuer a,b>1";pause();}
312  if(0==PerfectPowerTest(n))             // Schritt1 des ASK-Alg.
313  { return(c); }
314  if(printlevel>=1)
315  {
316    "Beginn: Berechnung von r, so dass ord(n) in Z/rZ groesser log2(n)^2 ";
317    pause();
318  }
319  bigint k,M,M2,b;
320  int r;                          // Zeile 526-542: Schritt 2
321  M=log2(n);                      // darin sind die Schritte
322  M2=M^2;                         // 3 und 4 integriert
323  r=2;
324
325  if(gcd(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
326                    // Schritt 3 des ASK-Alg.
327  {
328    if(printlevel>=1){"Zahl r gefunden mit r|n";"r=";r;pause();}
329    return(c);
330  }
331  if(r==n-1)              // falls diese Bedingung erfuellt ist, ist
332                         // ggT(a,n)=1 fuer 0<a<=r, schritt 4 des ASK-Alg.
333  {
334    if(printlevel>=1){"ggt(r,n)=1 fuer alle 1<r<n";pause();}
335    return(p);
336  }
337  k=1;
338  b=1;
339
340  while(k<=M2)         //Beginn des Ordnungstests fuer aktuelles r
341  {
342    b=((b*n) mod r);
343    if(b==1)   //tritt Bedingung ein so gilt fuer aktuelles r,k:
344               //n^k=1 mod r
345               //tritt Bedingung ein, wird wegen k<=M2=intPart(log2(n)^2)
346               //r erhoeht,k=0 gesetzt und Ordnung erneut untersucht
347               //tritt diese Bedingung beim groessten k=intPart(log2(n)^2)
348               //nicht ein, so ist ord_r_(n)>log2(n)^2 und Schleife
349               //wird nicht mehr ausgefuehrt
350    {
351      if(r==2147483647)
352      {
353        string e="error: r ueberschreitet integer Bereich und darf
354                  nicht als Grad eines Polynoms verwendet werden";
355        return(e);
356      }
357      r=r+1;
358      if(gcd(n,r)!=1)   //falls ggt ungleich eins, Teiler von n gefunden,
359                         //wird aufgrund von Schritt 3 des ASK-Alg. fuer
360                         //jeden Kandidaten r getestet
361      {
362        if(printlevel>=1){"Zahl r gefunden mit r|n";"r=";r;pause();}
363        return(c);
364      }
365      if(r==n-1)         //falls diese Bedingung erfuellt ist,
366                         //ist n wegen Zeile 571 prim
367                         //wird aufgrund von Schritt 4 des ASK-Alg. fuer
368                         //jeden Kandidaten r getestet
369      {
370        if(printlevel>=1){"ggt(r,n)=1 fuer alle 1<r<n";pause();}
371        return(p);
372      }
373
374      k=0;               //fuer neuen Kandidaten r, muss k fuer den
375                         //Ordnungstest zurueckgesetzt werden
376    }
377    k=k+1;
378  }
379  if(printlevel>=1)
380  {
381    "Zahl r gefunden, so dass ord(n) in Z/rZ groesser log2(n)^2 ";
382    "r=";r;pause();
383  }
384  ring @R=0,@x,dp;
385  bigint l=1;                         // Zeile 603-628: Schritt 5
386  poly f,g;                           // Zeile 618 ueberprueft Gleichungen fuer
387                                      // das jeweilige a
388  g=var(1)^r-1;
389  bigint grenze=intPart(wurzel(euler(r))*M);
390
391  if(printlevel>=1)
392  {"Beginn: Ueberpruefung ob (x+a)^n kongruent x^n+a mod (x^r-1,n)
393        fuer 0<a<=intPart(wurzel(phi(r))*log2(n))=";grenze;pause();
394  }
395  while(l<=grenze)          //
396  {
397    if(printlevel>=1){"a=";l;}
398    f=var(1)+l;
399    if(powerpolyX(bigint(n),bigint(n),f,g)!=(var(1)^(int((n mod r)))+l))
400    {
401      if(printlevel>=1)
402      {"Fuer a=";l;"ist (x+a)^n nicht kongruent x^n+a mod (x^r-1,n)";
403        pause();
404      }
405      return(c);
406    }
407    l=l+1;
408  }
409  if(printlevel>=1)
410  {"(x+a)^n kongruent x^n+a mod (x^r-1,n) fuer 0<a<wurzel(phi(r))*log2(n)";
411    pause();
412  }
413  return(p);       // Schritt 6
414}
415example
416{ "EXAMPLE:"; echo = 2;
417   //ask(100003);
418   ask(32003);
419}
Note: See TracBrowser for help on using the repository browser.