Crash with primefactors

[gorzel]

Crash with primefactors

Die ehemalige Version von primefactors aus general.lib war mangelhaft das aktuelle Kommando im Kernel ist nun besser, hat aber auch noch Schwachstellen:

1.) primfactors(0);

terminiert nicht: top CPU 99,9%

2.) Bei negativen Zahlen stuerzt Singular ab.

> primefactors(-2);
Singular : signal 8 (v: 3115/ 13286 ):
current line:>>primefactors(-2);<<
Segment fault/Bus error occurred at 851d4e9 because of 0 (r:1285594714)
please inform the authors
trying to restart...

3.) Wenn nur die Primzahlen <=MAXINT ermittelt werden, sollte die Rueckgabe-

liste, wie vorher aus

intvec von Faktoren, und intvec von Exponenten

und nicht aus

list von Faktoren, und list von Exponenten

bestehen.

4.) Die Reihenfolge der Elemente in der Rueckgabe liste ist nicht mehr wie

frueher. Daher wird so nicht mehr primecoefs laufen (dies wiederum

hat auch noch einen Fehler.)

Solle das derzeitige Format des Rueckgabewertes das endgueltige sein?

5.) Was besagt das 4. Slot? So wie ich verstehe, wird hiermit angezeigt,

ob die Eingabe vollstaendig zerlegt ist. Dann sollte aber in jedem Fall, die Ausgabe so sein, dass man einfach durch Multiplikation der Rueckgabewerte die urspruengliche Eingabe erhaelt.

> primefactors(37);
[1]:
   1
[2]:
   [1]:
      37
[3]:
   [1]:
      1
[4]:
   0
> primefactors(37,35);
[1]:
   37
[2]:
   empty list
[3]:
   empty list
[4]:
   1

Die leeren Listen, sollte man vielleicht durch jeweils durch 1 ersetzen.

Das Flag im 4. Slot sollte umgekehrt 1 bei erkannter vollstaendiger

0 bei unvollstanediger Zerlegung sein?

------------------------------------------------------------------------

Da ich fuer meine decomp.lib Primfaktorisierung zumindest der Exponenten benoetigte und sich die alte proc primefactors als Flaschenhals erwies, hatte ich vor ~einem Jahr primefactors fuer mich wesentliche verbessert und cgprimefactors genannt. Das Format der Rueckgabelist habe ich nicht umgedreht, so dass dies auch noch kompatibel mit primecofs ist. D.h. der dritte slot enthaelt den vorzeichenbehafteten restlichen Faktor.

Eine weitere geringfügige Beschleunigung ergibt sich daraus, dass man in der benoetigten proc primes die neuen Primzahlen hinten! am intvec anfuegt, dies habe ich in cgprimes so geandert.

Siehe die Beispiele und den Help-string meiner proc cgprimefactors.

Abgesehen, davon dass cgprimefcators nicht abstuerzt hat es noch

folgende Eigenschaften

• negative Zahlen werden faktorisiert

• konstante polynomials werden faktorisiert

• Ein Fehlermeldung bei Bruechen wird ausgegeben.

• Primzahlen<=MAXINT werden sofort als solche erkannt, d.h. das zweite Argument zur Einschraenkung der primen Kandidaten wird dann nicht beachtet.

• Zahlen |n|<=MAXINT werden stets vollstaendig zerlegt, ein zweites Argumen wird nicht beruecksichtigt.

• Das zweite Argument p, ist nicht die Beschraenkung fuer die zuletzt entdeckte Primzahl, sondern eine Schranke dafuer dass der restliche Faktor eine Primzahl <=p^{2 ist. (d.h. hiermit wird sqrt(N) fuer die noch zu untersuchende Zahl bestimmt.)}

• Es koennen so auch Primzahlen ausserhalb von MAXINT entdeckt werden. Dies kennzeichne ich mit "bigint" im 4. Slot.

• Mit dbprint wird angzeigt, welche letzte Primzahl tatsaechlich benoetigt wurde.

• Der Rueckgabetyp des restlichen Faktors im 3. Slot, ist der des Eingabetyps. D.h. die interne konvertierung auf bigint wird dort wieder rueckgaengig gemacht.

Vermutlich will/soll Euer Kernel-kommando und die Beschreibung dies im wesentlichen ebenso sein von mir implementiert.

Der Unterschied ist die Interpretation der Schranke p und der Angabe im Flag und des restlichen Faktors und des restlichen Faktors diese schreibe ich

Bei mir: Bei Euch:

----------- --------------

• Faktoren sind int, * Faktoren sind ints oder ein bigint --> intvec --> list

• restlicher Faktor im Slot 1 * restlicher Faktor auuserhalb int

hat den Typ der Eingabe stets bigint

• Bildet man die Produkte aus * Es koennen empty list vorkommen der ersten drei Slots ergibt sich stets die Eingabe.

• Eine Zerlegung ist vollstaendig wenn im 3. Slot +/-1 steht

• oder im 4. Slot "bigint" auftaucht

• das Format anedert sich nict, wenn * es aendert sich die Zerlegung schon vollstaendig ist aber p zu klein gewaehlt wurde!

Es zeigt sich, dass meine interpretertierte proc immer noch ein sehr gutes Laufzeitverhalten im Vergleich zum Kernelkommando hat. Der Engpass ist nun das Kommando primes um die naechste Primzahl zu finden. Es sollte daher eigentlich das Kernel kommando *prime* erweitert werden!

Wie gesagt benoetige solch ein Kommando selber (fuer int im Bereich der Exponenten). Daher bin ich zumindest an der Festlegung des Ausgabeformat interessiert.

( Ich hatte schon ca. 2Jahre angekuendigt, eine von general.lib

eine version mit Verbesserungen zu schicken. Aber dies stammte nicht von mir und die Aenderung nochmals in Detail zu beschreiben, dauert so lange wie die procs einfach selber zu beheben... )

Gruss,

Christian

-----------------------------------------------------------------------------

proc cgprimefactors (def n, list #)
"USAGE:   cgprimefactors(n [,p]); n = int, bigint, number or a constant poly, 
                                 p = integer
COMPUTE:  if |n| <= MAXINT compute its prime factorization
           uses primes <= min(p,46436)
RETURN: list,
         _[1] : intvec, primefactors of n
         _[2] : intvec, _[2][i] = multiplicity of _[1][i]
         _[3] : remaining factor such that n=product{(_[1][i]^_[2][i])} * _l[3]
                typeof(_[3])=input type, _[3] is negative if n is negative
    optionally, if a prime factor > MAXINT is detected:
         _[4] : string, \"bigprime\" 

NOTE:     cgprimefactors accepts numbers longer than 64 digits,
          hence it can so, in some sense, more than LINUX's factor command,
          but it can not always factor theses numbers since it uses a 
          limited range of primes.

    The factorisation is complete, if the input n is of the form

     (*) n = n_1^m_1* ... * n_l^m_l * n_r 
        
    where n_i, i=1,..,l  are primes <= p and n_r <=max(p^2,MAXINT)
    is also a prime number.

    It uses the builtin command prime to detect if n is a prime <= MAXINT.
     
    Since by default p = 46349, each number <= MAXINT is factored completely.
    An argument p <= 241 will be ignored since these prime numbers are 
    always used for the primality check.
    If these prime numbers are not sufficient, then the range of the primes
    is enlarged to min(p,46349) where by default p = 46349. 
    Since this takes some time, it may be useful to choose an appropriate
    bound <= sqrt(the smallest expected prime factor) 
 
    The last used prime number will be displayed, if printlevel is set to 
    a sufficient positive value.
   
EXAMPLE: example cgprimefactors; shows some examples
"
{
   int ii,jj,q;
   int neg,detected,increased;  // these int serve as flags 

   int MAXINT = 2147483647;   // 2^31 -1 
   int p = 46349;             // then p^2 > MAXINT;  
   bigint N,M;

   intvec v =               // initial test prime numbers <=241
      2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,
     197,199,211,223,227,229,233,239,241;

   intvec w1,w2=1,1;
   list l;
  
   if (size(#) > 0)
   {
     if( typeof(#[1]) != "int")
     {
       ERROR("2nd parameter must be of type int"+newline);
     }
     if (#[1]>0) 
     {               
        p=#[1]; 
     } 
     else
     {
	 "// Warning: since the 2nd parameter is <= 0, the default value "
           + string(p)+" will be used.";
	 //p = 241;
     } 
   }

// --------------------------- trivial cases --------------------
   if( n==0 )
   {
     return(list(w1[1],w2[1],0));
   }

   if( n<0 ) { n=-n; neg = 1;}

   if (n==1) { return(list(w1[1],w2[1],(-1)^neg))};

   N = bigint(n); // this cast converts number to bigint and detects a fraction
   if (n!=N)
   {
     ERROR("The number is not an integer (int resp. bigint)");
   }
   
   if(N==prime(int(N)))  // then N<=MAXINT
   {
     return(list(int(N),1,(-1)^neg));
   }
// --------------------------------------------------------------------------- 

   while(not(detected))
   {
        for( ii=1; ii<=size(v); ii++)   // check the test primes v
        {
          jj=0;
          while(1)  //divide N by v[ii] as often as possible
          {
            q  = v[ii];
            jj++;
            M  = N div q;    
            if (M*q!=N) {jj--; break;}  
            N=M;
          }
          if( jj>0 )
          {
             w1 = w1,v[ii];        //primes
             w2 = w2,jj;           //powers
          }
     // ----  check whether we have found all factors -----------------
 
	  if(N==prime(int(N)))  // this is a primality test if N<=MAXINT
          { 
	     w1 = w1,int(N);
             w2 = w2,1;
             N = 1;
          }
	  if (N==1)            // completely factored 
          { 
            detected=1; break;
          }                       
          if (bigint(q)^2>N)  // in this way detect if N is a prime > MAXINT 
          { 
            detected=2; break;
          }

     // --- increase the testprimes if all primes in v have been checked
     // --- but N is not completely factored yet

        if ((ii == size(v)) and (v[ii] < prime(p))) // all test primes used
	{                           // but still more primes to be used
	  dbprint("//-- candidates increased upto " + string(p));
          v=cgprimes(q,p); 
          ii=0;            // to change the counter of the for-loop  
                          // is bad style but I did not want to rewrite
                          // it as a while loop   
        }
        else              // not factored yet but all primes have been used 
        { 
	  detected = -1;
        }
       } 
//	dbprint("// status is ", detected);
	dbprint("//-- last used prime number was " + string(q));
   } // end while

// ----------- build return value -------------------------------------

   if (N<MAXINT)      // convert remaining factor to input type
   {                  
       n = int(N);    
   }
   else                 
   {                    
       n = N;
   }

  if( size(w1) > 1 )              // at least one prime factor was found
  {
    l[1] = intvec(w1[2..size(w1)]);
    l[2] = intvec(w2[2..size(w2)]);
    l[3] = (-1)^neg*n;
  }
  else                   // no prime factor was found or N is a big prime
  {
    l = list(1,1,(-1)^neg*n);
  }

  if (detected==2)       // N is a big prime number
  {
    l[4] = "big prime";
  }

  return(l);
}
example
{ "EXAMPLE:"; echo = 2;
   // anpassen
   printlevel = 2;
   cgprimefactors(251*257);  // enlarges the test primes upto 46349
   cgprimefactors(251*257,241);  // too small
   cgprimefactors(251*257,251);  // sufficient

   cgprimefactors(25125725345);
   cgprimefactors(7*8*121);

   
   ring r = 0,x,dp;
   cgprimefactors(number(123456789100));

   ring rZ =0,z,dp;
   cgprimefactors(number(25125725345));
}
///////////////////////////////////////////////////////////////////////////////

proc cgprimes (int n, int m)
// changed orginal proc primes 
// add the new prime at the end instead at the beginning of the intvec
// this yields some speed up for huge intvecs (large range of primes)
"USAGE:   cgprimes(n,m);  n,m integers
RETURN:  intvec, consisting of all primes p, prime(n)<=p<=m, in increasing
         order if n<=m, resp. prime(m)<=p<=n, in decreasing order if m<n.
NOTE:    prime(n); returns the biggest prime number <= min(n,2147483647)
         if n>=2, else 2
EXAMPLE: example primes; shows an example
"
{  int change;
   if ( n>m ) { change=n; n=m ; m=change; change=1; }
   int q,p = prime(m),prime(n); intvec v = q; q = q-1;
   while ( q>=p ) { q = prime(q); v = v,q; q = q-1; }
   if ( change==0 ) { v = v[size(v)..1]; }
   return(v);
}
example
{  "EXAMPLE:"; echo = 2;
    primes(50,100);"";
    intvec v = primes(37,1); v;"";
    // The larges primes known to Singular:
    int MAXINT = 2147483647;
    primes(MAXINT-100,MAXINT);
}
///////////////////////////////////////////////////////////////////////////////

Further examples for primefactors

 // See the difference for the intepretation of the second parameter p.

> int MAXINT = 2147483647;
>   printlevel =  3;
>   cgprimefactors(bigint(2)*MAXINT,23);  // completetly factored 
//- last used prime number was 2
[1]:
   2,2147483647
[2]:
   1,1
[3]:
   1
>   primefactors(bigint(2)*MAXINT,23); 
[1]:
   2147483647  See the difference for the intepretation of the second parameter p.

[2]:
   [1]:
      2
[3]:
   [1]:
      1
[4]:
   1


>   prime(263)==263;  // this is a prime number
1
>   prime(271)==271;  // this is a prime number
1


> cgprimefactors(271*263^2,263);     // completely factored
//-- candidates increased upto 263
//- last used prime number was 263
[1]:
   263,271
[2]:
   2,1
[3]:
   1

> primefactors(271*263^2,263);   // not detected as completely factored
[1]:
   271
[2]:
   [1]:
      263
[3]:
   [1]:
      2
[4]:
   1


> cgprimefactors(0);
[1]:
   1
[2]:
   1
[3]:
   0
> cgprimefactors(1);
[1]:
   1
[2]:
   1
[3]:
   1
> cgprimefactors(-4);
//- last used prime number was 2
[1]:
   2
[2]:
   2
[3]:
   -1

///---------------------------------


> bigint M3 = 8787878783;
> bigint N = (M3*2^15*3^7*5^6*bigint(51047)^2);
>  primefactors(N,94000);    // Dies wird von euch nicht vollstaendig zerlegt
[1]:
   8787878783
[2]:
   [1]:
      2
   [2]:
      3
   [3]:
      5
   [4]:
      51047
[3]:
   [1]:
      15
   [2]:
      7
   [3]:
      6
   [4]:
      2
[4]:
   1

> cgprimefactors(N,94000);                  // Ich zerlege es vollstaendig !!
//-- candidates increased upto 94000
//- last used prime number was 93761
[1]:
   2,3,5,51047
[2]:
   15,7,6,2
[3]:
   8787878783
[4]:
   big prime

>  cgprimefactors(N^2,94000);          // Das wir von mir nicht zerlegt.
//-- candidates increased upto 94000
//- last used prime number was 93997
[1]:
   2,3,5,51047
[2]:
   30,14,12,4
[3]:
   77226813504701561089



> ring r=0,x,dp;

> cgprimefactors(2/3);     // ERROR-msg
   ? The number is not an integer (int resp. bigint)
   ? leaving ::cgprimefactors

>   primefactors(2/3);   // missing ERROR-msg
[1]:
   1
[2]:
   [1]:
      2
[3]:
   [1]:
      1
[4]:
   0

> number q =  8787878783;
> cgprimefactors(q);
//-- candidates increased upto 46349
//- last used prime number was 46349
[1]:
   1
[2]:
   1
[3]:
   8787878783
> typeof(_[3]);                // restore the input type
number
> poly f =  q;
> cgprimefactors(f,46349);
//-- candidates increased upto 46349
//- last used prime number was 46349
[1]:
   1
[2]:
   1
[3]:
   8787878783
> typeof(_[3]);                // restore the input type
poly


> primefactors(q,46349);
[1]:
   8787878783
[2]:
   empty list
[3]:
   empty list
[4]:
   1
> typeof(_[1]);          // converted to bigint
bigint
> primefactors(f,46349);
   ? primefactors(`poly`,`int`) failed
   ? error occurred in or before STDIN line 31: `primefactors(f,46349);`

 // -------------------------

>  cgprimefactors(MAXINT*2,23);   // Overflow msg, OK.
// ** int overflow(*), result may be wrong
[1]:
   2
[2]:
   1
[3]:
   -1

> primefactors(MAXINT*2,23);   // Overflow msg  with crash, bad.
// ** int overflow(*), result may be wrong
Singular : signal 8 (v: 3115/ 13286 ):
current line:>>  primefactors(MAXINT*2,23);   // Overflow msg  with crash, bad.<<
Segment fault/Bus error occurred at 851d4e9 because of 0 (r:1285606729)
please inform the authors
trying to restart...

[gorzel]

Eine Teil kann man vermultich nicht lesen, da ich ihn nicht formatierte

      Bei mir:                               Bei Euch:
     -----------                          --------------
    
 * Faktoren sind int,                 * Faktoren sind ints oder ein bigint
   --> intvec                            --> list

 * restlicher Faktor im Slot 1        * restlicher Faktor auuserhalb int
     hat den Typ der Eingabe               stets bigint

 * Bildet man die Produkte aus        * Es koennen empty list vorkommen 
   der ersten drei Slots ergibt
   sich stets die Eingabe.

 * Eine Zerlegung ist vollstaendig   
   wenn im 3. Slot +/-1 steht

 * oder im  4. Slot "bigint" auftaucht

 * das Format anedert sich nict, wenn    * es aendert sich 
   die Zerlegung schon vollstaendig ist
   aber p zu klein gewaehlt wurde!

Desweitern, hier der Nachweis, dass primecoefs nicht mehr funktioniert:

> LIB "general.lib";
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/general.lib (12904,2010-06-24)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/matrix.lib (12898,2010-06-23)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/nctools.lib (12790,2010-05-14)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/random.lib (12827,2010-05-28)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/ring.lib (12231,2009-11-02)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/primdec.lib (12962,2010-07-09)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/absfact.lib (12231,2009-11-02)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/triang.lib (12231,2009-11-02)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/elim.lib (12231,2009-11-02)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/poly.lib (12443,2010-01-19)
// ** loaded /home/gorzelc/SING311-2/Singular/3-1-1/LIB/inout.lib (12541,2010-02-09)
> example primecoeffs;
// proc primecoeffs from lib general.lib
EXAMPLE:
   primecoeffs(intvec(7*8*121,7*8));"";
   ? incompatible type in list assignment
   ? error occurred in or before general.lib::primecoeffs line 1222: `     rest = rest,l[3];`
   ? leaving general.lib::primecoeffs
   skipping text from `;` error at token `)`
   ? leaving general.lib::primecoeffs

Der Fehler in primecoefs den ich oben erwaehnt ist damit noch nicht beschrieben, dies und weiter in anderen Tickets.

[seelisch]

alles gefixt: zu den obigen Fragen

1. + 2. hat Hans geändert (falls n < 0, dann steckt das Minuszeichen im nicht zerlegten Rest bzw. dieser ist dann einfach -1)
2. Faktoren müssen nun in einer Liste stehen, da diese nun auch vom Typ bigint sein können. Ob die Exponenten weiter in einem intvec sitzen sollen, kann man diskutieren. Macht aber für den SINGULAR-Benutzer kaum einen Unterschied in der Bedienung. Aber der kernel-Code wird etwas kleiner, weil man keine Extra-Kopier-Routine für intvec's braucht.
3. alte Reihenfolge bleibt (und wurde wieder hergestellt) Das war ein Versehen. Danke für den Hinweis!
4. Der 4. Rückgabewert sagt aus, ob der Betrag des nicht zerlegten Rests gemäß probabilistischem Test wahrscheinlich eine Primzahl ist oder nicht.