Changeset f3a11a in git


Ignore:
Timestamp:
Dec 14, 2006, 1:02:03 PM (17 years ago)
Author:
Gerhard Pfister <pfister@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
3703444c997bb04b388cb5ed654c3e719223bd29
Parents:
19ffafbfb9dbbdff4d14a0f428ddeb4bc349ab7f
Message:
Prozeduren schnellexp,schnellexppoly entfernt und die anderen diesbezueglich
angepasst.


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/ASK.lib

    r19ffafb rf3a11a  
    11//CM, last modified 10.12.06
    22///////////////////////////////////////////////////////////////////////////////
    3 version="$Id: ASK.lib,v 1.1 2006-12-11 12:52:42 Singular Exp $";
     3version="$Id: ASK.lib,v 1.2 2006-12-14 12:02:03 pfister Exp $";
    44category="Teaching";
    55info="
    … …  
    1212 based on the ideas of Agrawal, Saxena and  Kayal.
    1313
    14 
    15 
    16 
    1714PROCEDURES:
    1815
    19 
    20 schnellexp(a,m)            a^m for numbers a,m
    2116schnellexpt(a,m,n)         a^m for numbers a,m; if a^k>n n+1 is returned
    2217log2(n)                    logarithm to basis 2 of n
    2318PerfectPowerTest(n)        checks if there are a,b>1, so that a^b=n
    24 schnellexppoly(a,m)        a^m for number m,poly a
    2519wurzel(r)                  square root of number r
    2620euler(r)                   phi-function of euler
    … …  
    3327LIB "ntsolve.lib";
    3428
    35 
    36 
    37 
    38 
    39 ///////////////////////////////////////////////////////////////
    40 //                                                           //
    41 //   Fast Exponentiation                                     //
    42 //                                                           //
    43 ///////////////////////////////////////////////////////////////
    44 
    45 
    46 
    47 
    48 proc schnellexp(number a,number m)
    49 "USAGE: schnellexp(a,m);
    50 RETURN: the m-th power of a
    51 NOTE:   uses fast exponentiation
    52 EXAMPLE:example schnellexp; shows an example
    53 "
    54 {
    55  number b,c,d;
    56  c=m;
    57  b=a;
    58  d=1;
    59 
    60   while(c>=1)
    61   {
    62     if((c mod 2)==1)
    63      {
    64        d=d*b;
    65      }
    66 
    67    b=b^2;
    68    c=intPart(c/2);
    69 
    70   }
    71  return(d)
    72 }
    73 example
    74 { "EXAMPLE:"; echo = 2;
    75    ring R = 0,x,dp;
    76    schnellexp(24,1234);
    77 }
    78 
    79 
    80 
    81 ////////////////////////////////////////////////////////////////////////
    82 
    83 
     29///////////////////////////////////////////////////////////////
     30//                                                           //
     31//   FAST (MODULAR) EXPONENTIATION                           //
     32//                                                           //
     33///////////////////////////////////////////////////////////////
    8434proc schnellexpt(number a,number m,number n)
    8535"USAGE: schnellexpt(a,m,n);
    … …  
    12171   schnellexpt(2,10,1022);
    12272}
    123 
    124 
    125 
    126 
    127 
    128 ////////////////////////////////////////////////////////////////////////
    129 
    130 
    131 proc schnellexppoly(poly a,number m)
    132 "USAGE: schnellexppoly(a,m);
    133 RETURN: the m-th power of poly a
    134 NOTE:   uses fast exponentiation;
    135         proc schnellexp for a poly a instead of a number a;
    136 EXAMPLE:example schnellexppoly; shows an example
    137 "
    138 {
    139  number c;
    140  poly b,d;
    141 
    142  c=m;
    143  b=a;
    144  d=1;
    145 
    146   while(c>=1)
    147   {
    148     if((c mod 2)==1)
    149      {
    150        d=d*b;
    151      }
    152 
    153    b=b^2;
    154    c=intPart(c/2);
    155 
    156   }
    157  return(d)
    158 }
    159 example
    160 { "EXAMPLE:"; echo = 2;
    161    ring R = 0,x,dp;
     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;
    162126   poly a=3*x3-x2+5;
    163    number m=21;
    164    schnellexppoly(a,m);
    165 }
    166 
    167 
    168 
    169 
    170 
    171 
     127   poly r=x7-1;
     128   number q=123;
     129   number n=5;
     130   powerpolyX(q,n,a,r);
     131}
    172132///////////////////////////////////////////////////////////////
    173133//                                                           //
    … …  
    175135//                                                           //
    176136///////////////////////////////////////////////////////////////
    177 
    178 
    179 
    180 
    181137proc log2(number x)
    182138"USAGE:  log2(x);
    … …  
    187143"
    188144{
    189  number b,c,d,t,l,k;
     145   number b,c,d,t,l;
     146   int k;
    190147                                         // log2=logarithmus zur basis 2,
    191148                                         // log=natuerlicher logarithmus
    192  b=100000000000000000000000000000000000000000000000000;
    193  c=141421356237309504880168872420969807856967187537695;    // c/b=sqrt(2)
    194  d=69314718055994530941723212145817656807550013436026;     // d/b=log(2)
     149   b=100000000000000000000000000000000000000000000000000;
     150   c=141421356237309504880168872420969807856967187537695;    // c/b=sqrt(2)
     151   d=69314718055994530941723212145817656807550013436026;     // d/b=log(2)
    195152
    196153//bringen x zunaechst zwischen 1/sqrt(2) und sqrt(2), so dass Reihe schnell
    … …  
    199156//log2(x)=log(x)/log(2)=(log(x*2^j)-j*log(2))/log(2)  für kleine x
    200157
    201 number j=0;
    202 
    203 if(x<(b/c))
     158 number j=0;
     159
     160 if(x<(b/c))
    204161 {
    205162
    … …  
    215172  while(k<30)       //fuer  x*2^j wird Reihe bis k-tem Summanden berechnet
    216173   {
    217     l=l+2*schnellexp(t,2*k+1)/(2*k+1);      //l=log(x*2^j) nach k Summanden
     174    l=l+2*(t^(2*k+1))/(2*k+1);      //l=log(x*2^j) nach k Summanden
    218175
    219176    k=k+1;
    … …  
    237194 while(k<30)         //fuer  x/2^j wird Reihe bis k-tem Summanden berechnet
    238195  {
    239    l=l+2*schnellexp(t,2*k+1)/(2*k+1);       //l=log(x/2^j) nach k Summanden
     196   l=l+2*(t^(2*k+1))/(2*k+1);       //l=log(x/2^j) nach k Summanden
    240197   k=k+1;
    241198  }
    … …  
    248205   log2(1024);
    249206}
    250 
    251 
    252 
    253207//////////////////////////////////////////////////////////////////////////
    254 
    255 
    256208proc wurzel(number r)
    257209"USAGE: wurzel(r);
    … …  
    281233   wurzel(7629412809180100);
    282234}
    283 
    284 
    285 
    286235//////////////////////////////////////////////////////////////////////////
    287 
    288 
    289236proc euler(number r)
    290237"USAGE: euler(r);
    … …  
    313260   euler(99991);
    314261}
    315 
    316 
    317 
    318 
    319 
    320262///////////////////////////////////////////////////////////////
    321263//                                                           //
    … …  
    323265//                                                           //
    324266///////////////////////////////////////////////////////////////
    325 
    326 
    327267proc PerfectPowerTest(number n)
    328268"USAGE: PerfectPowerTest(n);
    … …  
    376316   PerfectPowerTest(887503681);
    377317}
    378 
    379 
    380 
    381 
    382 
    383 
    384 
    385 
    386 ///////////////////////////////////////////////////////////////
    387 //                                                           //
    388 //   FAST MODULAR EXPONENTIATION OF POLYNOMIALS              //
    389 //                                                           //
    390 ///////////////////////////////////////////////////////////////
    391 
    392 
    393 
    394 proc coeffmod(poly f,number n)
    395 "USAGE: coeffmod(f,n);
    396 ASSUME: poly f depends on at most var(1) of the basering
    397 RETURN: poly f modulo number n
    398 NOTE:   at first the coefficients of the monomials of the poly f are
    399         determined, then their remainder modulo n is calculated,
    400         after that the poly 'is put together' again
    401 EXAMPLE:example coeffmod; shows an example
    402 "
    403 {
    404   matrix M=coeffs(f,var(1));         //Bestimmung der Polynomkoeffizienten
    405   int i=1;
    406   poly g;
    407   number o=nrows(M);
    408 
    409   while(i<=o)
    410     {
    411       g=g+(leadcoef(M[i,1]) mod n)*var(1)^(i-1) ;
    412                      // umwandeln der Koeffizienten in numbers,
    413                      // Berechnung der Reste dieser numbers modulo n,
    414                      // poly mit neuen Koeffizienten wieder zusammensetzen
    415       i=i+1;
    416     }
    417   return(g);
    418 }
    419 example
    420 { "EXAMPLE:"; echo = 2;
    421    ring R = 0,x,dp;
    422    poly f=2457*x4+52345*x3-98*x2+5;
    423    number n=3;
    424    coeffmod(f,n);
    425 }
    426 
    427 
    428 
    429 
    430 
    431 
    432 
    433 
    434 //////////////////////////////////////////////////////////////////////////
    435 
    436 
    437 proc powerpolyX(number q,number n,poly a,poly r)
    438 "USAGE: powerpolyX(q,n,a,r);
    439 RETURN: the q-th power of poly a modulo poly r and number n
    440 EXAMPLE:example powerpolyX; shows an example
    441 "
    442 {
    443   ideal I=std(r);
    444 
    445    if(q==1){return(coeffmod(reduce(a,I),n));}
    446    if((q mod 5)==0){return(coeffmod(reduce(powerpolyX(q/5,n,a,r)^5,I),n));}
    447    if((q mod 4)==0){return(coeffmod(reduce(powerpolyX(q/4,n,a,r)^4,I),n));}
    448    if((q mod 3)==0){return(coeffmod(reduce(powerpolyX(q/3,n,a,r)^3,I),n));}
    449    if((q mod 2)==0){return(coeffmod(reduce(powerpolyX(q/2,n,a,r)^2,I),n));}
    450 
    451   return(coeffmod(reduce(a*powerpolyX(q-1,n,a,r),I),n));
    452 }
    453 example
    454 { "EXAMPLE:"; echo = 2;
    455    ring R=0,x,dp;
    456    poly a=3*x3-x2+5;
    457    poly r=x7-1;
    458    number q=123;
    459    number n=5;
    460    powerpolyX(q,n,a,r);
    461 }
    462 
    463 
    464 
    465 
    466318///////////////////////////////////////////////////////////////
    467319//                                                           //
    … …  
    469321//                                                           //
    470322///////////////////////////////////////////////////////////////
    471 
    472 
    473 
    474323proc ask(number n)
    475324"USAGE: ask(n);
    … …  
    596445     f=var(1)+l;
    597446
    598       if(powerpolyX(n,n,f,g)!=(schnellexppoly(var(1),(n mod r))+l))
     447       if(powerpolyX(n,n,f,g)!=(var(1)^(int((n mod r)))+l))
    599448       {
    600449         if(printlevel>=1)
    … …  
    617466{ "EXAMPLE:"; echo = 2;
    618467   ring R = 0,x,dp;
    619    ask(100003);
    620 }
    621 
    622 
     468   //ask(100003);
     469   ask(32003);
     470}
     471
     472
Note: See TracChangeset for help on using the changeset viewer.