# source:git/Singular/LIB/aksaka.lib@226c28

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