Changeset 226c28 in git


Ignore:
Timestamp:
Feb 21, 2009, 11:02:00 AM (14 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
59c445958ea9eeff400a858e69fea9959018518e
Parents:
d2f48833961a856dcfd3669df0f6f883f23efe6a
Message:
hannes: format, opt.


git-svn-id: file:///usr/local/Singular/svn/trunk@11420 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/aksaka.lib

    rd2f488 r226c28  
    11//CM, last modified 10.12.06
    22///////////////////////////////////////////////////////////////////////////////
    3 version="$Id: aksaka.lib,v 1.2 2009-02-20 19:02:43 Singular Exp $";
     3version="$Id: aksaka.lib,v 1.3 2009-02-21 10:02:00 Singular Exp $";
    44category="Teaching";
    55info="
    … …  
    3939"
    4040{
    41  number b,c,d;
    42  c=m;
    43  b=a;
    44  d=1;
    45 
     41  number b,c,d;
     42  c=m;
     43  b=a;
     44  d=1;
    4645  while(c>=1)
    4746  {
    48    if(b>n)
    49      {
     47    if(b>n)
     48    {
    5049      return(n+1);
    51      }
    52 
     50    }
    5351    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)
     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)
    6763}
    6864example
    … …  
    8581  int i=1;
    8682  poly g;
    87   number o=nrows(M);
     83  int o=nrows(M);
    8884
    8985  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     }
     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  }
    9793  return(g);
    9894}
    … …  
    111107"
    112108{
    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));}
     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));}
    120116
    121117  return(coeffmod(reduce(a*powerpolyX(q-1,n,a,r),I),n));
    … …  
    143139"
    144140{
    145    number b,c,d,t,l;
    146    int k;
     141  number b,c,d,t,l;
     142  int k;
    147143                                         // log2=logarithmus zur basis 2,
    148144                                         // 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;
     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;
    166175    j=j+1;
    167    }
    168 
     176  }
    169177  t=(x-1)/(x+1);
    170178  k=0;
    171179  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 
     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
    176183    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
     184  }
     185  return((l*b/d)+j);     //hier wird log2(x)=log(x/2^j)/log(2)+j  ausgegeben
    201186}
    202187example
    … …  
    213198"
    214199{
    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);
     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);
     206  number m=leadcoef(l[1][1]);
     207  setring B;
     208  return(imap(R,m));
    229209}
    230210example
    … …  
    243223{
    244224  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 
     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);
    256232}
    257233example
    … …  
    281257
    282258  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    }
     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  }
    310286  if(printlevel>=1){"es existieren keine Zahlen a,b>1  mit a^b=n";pause();}
    311287  return(f);
    … …  
    336312  p=1;
    337313
    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,
     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,
    363330                    // Schritt 3 des ASK-Alg.
    364    {
     331  {
    365332    if(printlevel>=1){"Zahl r gefunden mit r|n";"r=";r;pause();}
    366333    return(c);
    367    }
    368 
    369  if(r==n-1)              // falls diese Bedingung erfüllt ist, ist
     334  }
     335  if(r==n-1)              // falls diese Bedingung erfüllt ist, ist
    370336                         // ggT(a,n)=1 für 0<a<=r, schritt 4 des ASK-Alg.
    371337  {
    372    if(printlevel>=1){"ggt(r,n)=1 fuer alle 1<r<n";pause();}
    373    return(p);
    374   }
    375 
     338    if(printlevel>=1){"ggt(r,n)=1 fuer alle 1<r<n";pause();}
     339    return(p);
     340  }
    376341  k=1;
    377342  b=1;
    378343
    379344  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)
     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)
    460413  {"(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
     414    pause();
     415  }
     416  return(p);       // Schritt 6
    464417}
    465418example
    … …  
    469422   ask(32003);
    470423}
    471 
Note: See TracChangeset for help on using the changeset viewer.