Changeset 26508d in git


Ignore:
Timestamp:
Jul 20, 2007, 1:26:46 PM (16 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
7c38bcb5f2878e308ca07340e00bfe7b4fdd406e
Parents:
3eadab96fecbb57aa75151e22b15ce9249b408b6
Message:
*hannes: format


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/atkins.lib

    r3eadab r26508d  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: atkins.lib,v 1.3 2007-07-20 10:02:38 Singular Exp $";
     2version="$Id: atkins.lib,v 1.4 2007-07-20 11:26:46 Singular Exp $";
    33category="Teaching";
    44info="
     
    9595    bubblesort(L);
    9696}
    97 
    98 
    9997
    10098proc disc(number N, int k)
     
    193191}
    194192
    195 
    196 
    197193proc CornacchiaModified(number D, number p)
    198194"USAGE: CornacchiaModified(D,p);
     
    204200"
    205201{
    206    if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p))
    207    {// (0)[Test if assumptions well-defined]
    208       return(0);
    209       // ERROR("Parameters wrong selected!");
    210    }
    211    else
    212    {
    213       if(p==2)                          // (1)[Case p=2]
    214       {
    215          if((D+8)==intRoot(D+8)^2)
    216          {
    217             return(intRoot(D+8),1);
    218          }
    219          else
    220          {
    221             return(-1);
    222             // ERROR("The Diophantine equation has no solution!");
    223          }
     202  if((D>=0)||((D mod 4)==2)||((D mod 4)==3)||(absValue(D)>=4*p))
     203  {// (0)[Test if assumptions well-defined]
     204    return(0);
     205    // ERROR("Parameters wrong selected!");
     206  }
     207  else
     208  {
     209    if(p==2)                          // (1)[Case p=2]
     210    {
     211      if((D+8)==intRoot(D+8)^2)
     212      {
     213        return(intRoot(D+8),1);
    224214      }
    225215      else
    226216      {
    227          number k,x(0),a,b,l,r,c;
    228          k=Jacobi(D,p);                // (2)[Test if residue]
    229          if(k==-1)
    230          {
    231             return(-1);
    232             // ERROR("The Diophantine equation has no solution!");
    233          }
    234          else
    235          {
    236             x(0)=squareRoot(D,p);    // (3)[Compute square root]
    237             if((x(0) mod 2)!=(D mod 2))
    238             {
    239                x(0)=p-x(0);
    240             }
    241             a=2*p;
    242             b=x(0);
    243             l=intRoot(4*p);
    244             while(b>l)     // (4)[Euclidean algorithm]
    245             {
    246                r=a mod b;
    247                a=b;
    248                b=r;
    249             }
    250             c=(4*p-b^2)/absValue(D);// (5)[Test solution]
    251             if((((4*p-b^2) mod absValue(D))!=0)||(c!=intRoot(c)^2))
    252             {
    253                return(-1);
    254                // ERROR("The Diophantine equation has no solution!");
    255             }
    256             else
    257             {
    258                list L=b,intRoot(c);
    259                return(L);
    260             }
    261          }
    262       }
    263    }
     217        return(-1);
     218        // ERROR("The Diophantine equation has no solution!");
     219      }
     220    }
     221    else
     222    {
     223      number k,x(0),a,b,l,r,c;
     224      k=Jacobi(D,p);                // (2)[Test if residue]
     225      if(k==-1)
     226      {
     227        return(-1);
     228        // ERROR("The Diophantine equation has no solution!");
     229      }
     230      else
     231      {
     232        x(0)=squareRoot(D,p);    // (3)[Compute square root]
     233        if((x(0) mod 2)!=(D mod 2))
     234        {
     235          x(0)=p-x(0);
     236        }
     237        a=2*p;
     238        b=x(0);
     239        l=intRoot(4*p);
     240        while(b>l)     // (4)[Euclidean algorithm]
     241        {
     242          r=a mod b;
     243          a=b;
     244          b=r;
     245        }
     246        c=(4*p-b^2)/absValue(D);// (5)[Test solution]
     247        if((((4*p-b^2) mod absValue(D))!=0)||(c!=intRoot(c)^2))
     248        {
     249          return(-1);
     250          // ERROR("The Diophantine equation has no solution!");
     251        }
     252        else
     253        {
     254          list L=b,intRoot(c);
     255          return(L);
     256        }
     257      }
     258    }
     259  }
    264260}
    265261example
     
    268264    CornacchiaModified(-107,1319);
    269265}
    270 
    271 
    272266
    273267proc pFactor1(number n,int B, list P)
     
    281275"
    282276{
    283       int i;
    284       number k=1;
    285       number w;
    286       while(i<size(P))
    287          {
    288            i++;
    289            w=P[i];
    290            if(w>B) {break;}
    291            while(w*P[i]<=B)
    292               {
    293                 w=w*P[i];
    294               }
    295            k=k*w;
    296          }
    297       number a=random(2,2147483629);
    298       number d=gcdN(powerN(a,k,n)-1,n);
    299       if((d>1)&&(d<n)){return(d);}
    300       return(n);
     277  int i;
     278  number k=1;
     279  number w;
     280  while(i<size(P))
     281  {
     282    i++;
     283    w=P[i];
     284    if(w>B) {break;}
     285    while(w*P[i]<=B)
     286    {
     287      w=w*P[i];
     288    }
     289    k=k*w;
     290  }
     291  number a=random(2,2147483629);
     292  number d=gcdN(powerN(a,k,n)-1,n);
     293  if((d>1)&&(d<n)){return(d);}
     294  return(n);
    301295}
    302296example
     
    310304}
    311305
    312 
    313 
    314306proc maximum(list L)
    315307"USAGE: maximum(list L);
     
    318310"
    319311{
    320       number max=L[1];
    321 
    322       int i;
    323       for(i=2;i<=size(L);i++)
    324          {
    325            if(L[i]>max)
    326               {
    327                 max=L[i];
    328               }
    329          }
    330 
    331       return(max);
     312  number max=L[1];
     313  int i;
     314  for(i=2;i<=size(L);i++)
     315  {
     316    if(L[i]>max)
     317    {
     318      max=L[i];
     319    }
     320  }
     321  return(max);
    332322}
    333323example
     
    337327    maximum(L);
    338328}
    339 
    340 
    341329
    342330static proc cmod(number x, number y)
     
    349337"
    350338{
    351       int rest=int(x-y*int(x/y));
    352       if(rest<0)
    353          {
    354            rest=rest+int(y);
    355          }
    356 
    357       return(rest);
     339  int rest=int(x-y*int(x/y));
     340  if(rest<0)
     341  {
     342    rest=rest+int(y);
     343  }
     344  return(rest);
    358345}
    359346example
     
    365352}
    366353
    367 
    368 
    369 static proc sqr(number w, int k)
     354proc sqr(number w, int k)
    370355"USAGE: sqr(w,k);
    371356RETURN: the square root of w
     
    404389  number q=1;
    405390  number e=1;
    406 
    407391  int n;
    408392  for(n=1;n<=k;n++)
     
    432416"
    433417{
    434       number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j;
    435 
    436       number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989;
    437 
    438       tr=repart(t);
    439       ti=impart(t);
    440       if(tr==-1/2){qr1=-1;}
    441       if(tr==0){qr1=1;}
    442       if((tr!=-1/2)&&(tr!=0))
    443          {
    444            tr=tr-round(tr);
    445            qr1=expo(2*i*pi*tr,10*k);
    446          }
    447 
    448       qi1=expo(-pi*ti,10*k);
    449       q1=qr1*qi1^2;
    450       q2=q1^2;
    451 
    452       int n=1;
    453       while(n<=k)
    454          {
    455            m1=m1+(-1)^n*(q1^(n*(3*n-1)/2)+q1^(n*(3*n+1)/2));
    456            m2=m2+(-1)^n*(q2^(n*(3*n-1)/2)+q2^(n*(3*n+1)/2));
    457            n=n+1;
    458          }
    459 
    460       f=q1*((1+m2)/(1+m1))^24;
    461 
    462       j=(256*f+1)^3/f;
    463       return(j);
    464 }
    465 
     418  number q1,q2,qr1,qi1,tr,ti,m1,m2,f,j;
     419
     420  number pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989;
     421
     422  tr=repart(t);
     423  ti=impart(t);
     424  if(tr==-1/2){qr1=-1;}
     425  if(tr==0){qr1=1;}
     426  if((tr!=-1/2)&&(tr!=0))
     427  {
     428    tr=tr-round(tr);
     429    qr1=expo(2*i*pi*tr,10*k);
     430  }
     431
     432  qi1=expo(-pi*ti,10*k);
     433  q1=qr1*qi1^2;
     434  q2=q1^2;
     435
     436  int n=1;
     437  while(n<=k)
     438  {
     439    m1=m1+(-1)^n*(q1^(n*(3*n-1)/2)+q1^(n*(3*n+1)/2));
     440    m2=m2+(-1)^n*(q2^(n*(3*n-1)/2)+q2^(n*(3*n+1)/2));
     441    n++;
     442  }
     443
     444  f=q1*((1+m2)/(1+m1))^24;
     445
     446  j=(256*f+1)^3/f;
     447  return(j);
     448}
    466449example
    467450{ "EXAMPLE:"; echo = 2;
     
    470453    jOft(t,50);
    471454}
    472 
    473 
    474455
    475456proc round(number r)
     
    538519"
    539520{
    540       if(D>=0)                                                         // (0)[Test if assumptions well-defined]
    541          {
    542            ERROR("Parameter wrong selected!");
    543          }
    544 
    545       else
    546          {
    547            def S=basering;
    548            ring R=0,x,dp;
    549 
    550            string s1,s2,s3;
    551            number a1,b1,t1,g1;
    552            number D=imap(S,D);
    553            number B=intRoot(absValue(D)/3);
    554 
    555            ring C=(complex,10*k,i),x,dp;
    556            number D=imap(S,D);
    557 
    558            poly P=1;                                                   // (1)[Initialize]
    559            number b=cmod(D,2);
    560            number B=imap(R,B);
    561 
    562            number t,a,g,tau,j;
    563            list L;
    564 
    565            int step=2;
    566            while(1)
    567               {
    568                 if(step==2)                                            // (2)[Initialize a]
    569                    {
    570                      t=(b^2-D)/4;
    571                      L=b,1;
    572                      a=maximum(L);
    573                      step=3;
    574                    }
    575 
    576                 if(step==3)                                            // (3)[Test]
    577                    {
    578                      if((cmod(t,a)!=0))
    579                         {
    580                           step=4;
    581                         }
    582 
    583                      else
    584                         {
    585                           s1=string(a);
    586                           s2=string(b);
    587                           s3=string(t);
    588 
    589                           setring R;
    590                           execute("a1="+s1+";");
    591                           execute("b1="+s2+";");
    592                           execute("t1="+s3+";");
    593                           g1=gcd(gcd(a1,b1),t1/a1);
    594                           setring C;
    595                           g=imap(R,g1);
    596 
    597                           if(g!=1)
    598                              {
    599                                step=4;
    600                              }
    601 
    602                           else
    603                              {
    604                                tau=(-b+i*sqr(absValue(D),5*k))/(2*a);
    605                                j=jOft(tau,k);
    606                                if((a==b)||(a^2==t)||(b==0))
    607                                   {
    608                                     P=P*(var(1)-repart(j));
    609                                     step=4;
    610                                   }
    611 
    612                                else
    613                                   {
    614                                     P=P*(var(1)^2-2*repart(j)*var(1)+repart(j)^2+impart(j)^2);
    615                                     step=4;
    616                                   }
    617                              }
    618                         }
    619                    }
    620 
    621                 if(step==4)                                            // (4)[Loop on a]
    622                    {
    623                      a=a+1;
    624                      if(a^2<=t)
    625                         {
    626                           step=3;
    627                           continue;
    628                         }
    629 
    630                      else
    631                         {
    632                           step=5;
    633                         }
    634                    }
    635 
    636                 if(step==5)                                            // (5)[Loop on b]
    637                    {
    638                      b=b+2;
    639                      if(b<=B)
    640                         {
    641                           step=2;
    642                         }
    643 
    644                      else
    645                         {
    646                           break;
    647                         }
    648                    }
    649               }
    650 
    651            matrix M=coeffs(P,var(1));
    652 
    653            list liste;
    654            int n;
    655            for(n=1;n<=nrows(M);n++)
    656               {
    657                 liste[n]=round(repart(number(M[n,1])));
    658               }
    659 
    660            poly Q;
    661            int m;
    662            for(m=1;m<=size(liste);m++)
    663               {
    664                 Q=Q+liste[m]*var(1)^(m-1);
    665               }
    666 
    667            string s=string(Q);
    668            setring S;
    669            execute("poly Q="+s+";");
    670            return(Q);
    671          }
     521  if(D>=0)                // (0)[Test if assumptions well-defined]
     522  {
     523    ERROR("Parameter wrong selected!");
     524  }
     525  else
     526  {
     527    def S=basering;
     528    ring R=0,x,dp;
     529
     530    string s1,s2,s3;
     531    number a1,b1,t1,g1;
     532    number D=imap(S,D);
     533    number B=intRoot(absValue(D)/3);
     534
     535    ring C=(complex,10*k,i),x,dp;
     536    number D=imap(S,D);
     537
     538    poly P=1;          // (1)[Initialize]
     539    number b=cmod(D,2);
     540    number B=imap(R,B);
     541
     542    number t,a,g,tau,j;
     543    list L;
     544
     545    int step=2;
     546    while(1)
     547    {
     548      if(step==2)              // (2)[Initialize a]
     549      {
     550        t=(b^2-D)/4;
     551        L=b,1;
     552        a=maximum(L);
     553        step=3;
     554      }
     555
     556      if(step==3)               // (3)[Test]
     557      {
     558        if((cmod(t,a)!=0))
     559        {
     560          step=4;
     561        }
     562        else
     563        {
     564          s1=string(a);
     565          s2=string(b);
     566          s3=string(t);
     567
     568          setring R;
     569          execute("a1="+s1+";");
     570          execute("b1="+s2+";");
     571          execute("t1="+s3+";");
     572          g1=gcd(gcd(a1,b1),t1/a1);
     573          setring C;
     574          g=imap(R,g1);
     575
     576          if(g!=1)
     577          {
     578            step=4;
     579          }
     580          else
     581          {
     582            tau=(-b+i*sqr(absValue(D),5*k))/(2*a);
     583            j=jOft(tau,k);
     584            if((a==b)||(a^2==t)||(b==0))
     585            {
     586              P=P*(var(1)-repart(j));
     587              step=4;
     588            }
     589            else
     590            {
     591              P=P*(var(1)^2-2*repart(j)*var(1)+repart(j)^2+impart(j)^2);
     592              step=4;
     593            }
     594          }
     595        }
     596      }
     597
     598      if(step==4)                  // (4)[Loop on a]
     599      {
     600        a=a+1;
     601        if(a^2<=t)
     602        {
     603          step=3;
     604          continue;
     605        }
     606        else
     607        {
     608          step=5;
     609        }
     610      }
     611
     612      if(step==5)                 // (5)[Loop on b]
     613      {
     614        b=b+2;
     615        if(b<=B)
     616        {
     617          step=2;
     618        }
     619        else
     620        {
     621          break;
     622        }
     623      }
     624    }
     625
     626    matrix M=coeffs(P,var(1));
     627
     628    list liste;
     629    int n;
     630    for(n=1;n<=nrows(M);n++)
     631    {
     632      liste[n]=round(repart(number(M[n,1])));
     633    }
     634
     635    poly Q;
     636    int m;
     637    for(m=1;m<=size(liste);m++)
     638    {
     639      Q=Q+liste[m]*var(1)^(m-1);
     640    }
     641
     642    string s=string(Q);
     643    setring S;
     644    execute("poly Q="+s+";");
     645    return(Q);
     646  }
    672647}
    673648example
     
    677652    HilbertClassPoly(D,50);
    678653}
    679 
    680 
    681654
    682655proc rootsModp(int p, poly P)
     
    689662"
    690663{
    691       if(p<3)                                                              // (0)[Test if assumptions well-defined]
    692          {
    693            ERROR("Parameter wrong selected, since p<3!");
    694          }
    695 
    696       else
    697          {
    698            def S=basering;
    699            ring R=p,var(1),dp;
    700 
    701            poly P=imap(S,P);
    702            number d;
    703            int a;
    704            list L;
    705 
    706            poly A=gcd(var(1)^p-var(1),P);                                  // (1)[Isolate roots in Z/pZ]
    707            if(subst(A,var(1),0)==0)
    708               {
    709                 L[1]=0;
    710                 A=A/var(1);
    711               }
    712 
    713            if(deg(A)==0)                                                   // (2)[Small degree?]
    714               {
    715                 return(L);
    716               }
    717 
    718            if(deg(A)==1)
    719               {
    720                 matrix M=coeffs(A,var(1));
    721                 L[size(L)+1]=-leadcoef(M[1,1])/leadcoef(M[2,1]);
    722                 setring S;
    723                 list L=imap(R,L);
    724                 return(L);
    725               }
    726 
    727            if(deg(A)==2)
    728               {
    729                 matrix M=coeffs(A,var(1));
    730                 d=leadcoef(M[2,1])^2-4*leadcoef(M[1,1])*leadcoef(M[3,1]);
    731 
    732                 ring T=0,var(1),dp;
    733                 number d=imap(R,d);
    734                 number e=squareRoot(d,p);
    735                 setring R;
    736                 number e=imap(T,e);
    737 
    738                 L[size(L)+1]=(-leadcoef(M[2,1])+e)/(2*leadcoef(M[3,1]));
    739                 L[size(L)+1]=(-leadcoef(M[2,1])-e)/(2*leadcoef(M[3,1]));
    740                 setring S;
    741                 list L=imap(R,L);
    742                 return(L);
    743               }
    744 
    745            poly B=1;                                                       // (3)[Random splitting]
    746            poly C;
    747            while((deg(B)==0)||(deg(B)==deg(A)))
    748               {
    749                 a=random(0,p-1);
    750                 B=gcd((var(1)+a)^((p-1)/2)-1,A);
    751                 C=A/B;
    752               }
    753 
    754            setring S;                                                      // (4)[Recurse]
    755            poly B=imap(R,B);
    756            poly C=imap(R,C);
    757            list l=L+rootsModp(p,B)+rootsModp(p,C);
    758            return(l);
    759          }
     664  if(p<3)                         // (0)[Test if assumptions well-defined]
     665  {
     666    ERROR("Parameter wrong selected, since p<3!");
     667  }
     668  else
     669  {
     670    def S=basering;
     671    ring R=p,var(1),dp;
     672
     673    poly P=imap(S,P);
     674    number d;
     675    int a;
     676    list L;
     677
     678    poly A=gcd(var(1)^p-var(1),P);         // (1)[Isolate roots in Z/pZ]
     679    if(subst(A,var(1),0)==0)
     680    {
     681      L[1]=0;
     682      A=A/var(1);
     683    }
     684
     685    if(deg(A)==0)                          // (2)[Small degree?]
     686    {
     687      return(L);
     688    }
     689
     690    if(deg(A)==1)
     691    {
     692      matrix M=coeffs(A,var(1));
     693      L[size(L)+1]=-leadcoef(M[1,1])/leadcoef(M[2,1]);
     694      setring S;
     695      list L=imap(R,L);
     696      return(L);
     697    }
     698
     699    if(deg(A)==2)
     700    {
     701      matrix M=coeffs(A,var(1));
     702      d=leadcoef(M[2,1])^2-4*leadcoef(M[1,1])*leadcoef(M[3,1]);
     703
     704      ring T=0,var(1),dp;
     705      number d=imap(R,d);
     706      number e=squareRoot(d,p);
     707      setring R;
     708      number e=imap(T,e);
     709
     710      L[size(L)+1]=(-leadcoef(M[2,1])+e)/(2*leadcoef(M[3,1]));
     711      L[size(L)+1]=(-leadcoef(M[2,1])-e)/(2*leadcoef(M[3,1]));
     712      setring S;
     713      list L=imap(R,L);
     714      return(L);
     715    }
     716
     717    poly B=1;                       // (3)[Random splitting]
     718    poly C;
     719    while((deg(B)==0)||(deg(B)==deg(A)))
     720    {
     721      a=random(0,p-1);
     722      B=gcd((var(1)+a)^((p-1)/2)-1,A);
     723      C=A/B;
     724    }
     725
     726    setring S;                      // (4)[Recurse]
     727    poly B=imap(R,B);
     728    poly C=imap(R,C);
     729    list l=L+rootsModp(p,B)+rootsModp(p,C);
     730    return(l);
     731  }
    760732}
    761733example
     
    767739    rootsModp(1223,g);
    768740}
    769 
    770 
    771741
    772742proc wUnit(number D)
     
    12251195  }
    12261196}
    1227 
    12281197example
    12291198{ "EXAMPLE:"; echo = 2;
Note: See TracChangeset for help on using the changeset viewer.