Changeset 65b813 in git

Timestamp:

Oct 3, 2012, 1:32:09 PM (10 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')

Children:

b7e5d093cbde03bb95898a886e23c3bbcd486eaf

Parents:

4c434098fb8ab226e3e557478e5aee3ce970aa25

git-author:

Hans Schoenemann <hannes@mathematik.uni-kl.de>2012-10-03 13:32:09+02:00

git-committer:

Martin Lee <martinlee84@web.de>2013-01-23 11:27:17+01:00

Message:

fix: primefactors(): result and algorithm, bound enabled again

Conflicts:

	Singular/ChangeLog
	Singular/iparith.cc

Files:

• Singular/iparith.cc (modified) (2 diffs)
• Singular/misc_ip.cc (modified) (16 diffs)
• Singular/misc_ip.h (modified) (3 diffs)
• Singular/table.h (modified) (1 diff)
• doc/reference.doc (modified) (previous)

Singular/iparith.cc

@@ -3152 +3152 @@
 static BOOLEAN jjPFAC2(leftv res, leftv u, leftv v)
 {
-  number n1; number n2; number temp; int i;
+  number n1; int i;
 
   if ((u->Typ() == BIGINT_CMD) ||
      ((u->Typ() == NUMBER_CMD) && rField_is_Q(currRing)))
   {
-    temp = (number)u->Data();
-    n1 = n_Copy(temp,coeffs_BIGINT);
+    n1 = (number)u->CopyD();
   }
   else if (u->Typ() == INT_CMD)
@@ -3167 +3166 @@
   else
   {
-    WerrorS("wrong type: expected int, bigint, or number as 1st argument");
-    return TRUE;
-  }
-
-  if ((v->Typ() == BIGINT_CMD) ||
-     ((v->Typ() == NUMBER_CMD) && rField_is_Q(currRing)))
-  {
-    temp = (number)v->Data();
-    n2 = n_Copy(temp,coeffs_BIGINT);
-  }
-  else if (v->Typ() == INT_CMD)
-  {
-    i = (int)(long)v->Data();
-    n2 = n_Init(i, coeffs_BIGINT);
-  }
-  else
-  {
-    WerrorS("wrong type: expected int, bigint, or number as 2nd argument");
-    return TRUE;
-  }
-
-  lists l = primeFactorisation(n1, n2);
-  n_Delete(&n1, coeffs_BIGINT); n_Delete(&n2, coeffs_BIGINT);
+    return TRUE;
+  }
+
+  i = (int)(long)v->Data();
+
+  lists l = primeFactorisation(n1, i);
+  n_Delete(&n1, coeffs_BIGINT);
   res->data = (char*)l;
   return FALSE;

Singular/misc_ip.cc

@@ -73 +73 @@
 static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
 
-void factor_using_division (mpz_t t, unsigned int limit,lists primes, int *multiplicities,int &index)
+int factor_using_division (mpz_t t, unsigned int limit,lists primes, int *multiplicities,int &index, unsigned long bound)
 {
   mpz_t q, r;
@@ -80 +80 @@
   unsigned *addv = add;
   unsigned int failures;
+  int bound_not_reached=1;
 
   mpz_init (q);
@@ -124 +125 @@
   f = 7;
   ai = 0;
-  unsigned last_f=0;
+  unsigned long last_f=0;
   while (mpz_cmp_ui (t, 1) != 0)
   {
@@ -131 +132 @@
     {
       f += addv[ai];
-      if (mpz_cmp_ui (q, f) < 0)
+      if (mpz_cmp_ui (t, f) < 0)
         break;
       ai = (ai + 1) & 7;
@@ -137 +138 @@
       if (failures > limit)
         break;
+      if ((bound!=0) && (f>bound))
+      {
+        bound_not_reached=0;
+        break;
+      }
     }
     else
@@ -157 +163 @@
 
   mpz_clears (q, r, NULL);
-}
-
-void factor_using_pollard_rho (mpz_t n, unsigned long a, unsigned long p,lists primes, int * multiplicities,int &index)
+  //printf("bound=%d,f=%d,failures=%d, reached=%d\n",bound,f,failures,bound_not_reached);
+  return bound_not_reached;
+}
+
+void factor_using_pollard_rho (mpz_t n, unsigned long a, lists primes, int * multiplicities,int &index)
 {
   mpz_t x, x1, y, P;
@@ -182 +190 @@
       do
       {
-        if (p != 0)
-        {
-          mpz_powm_ui (x, x, p, n);
-          mpz_add_ui (x, x, a);
-        }
-        else
-        {
-          mpz_mul (t1, x, x);
-          mpz_mod (x, t1, n);
-          mpz_add_ui (x, x, a);
-        }
-
+        mpz_mul (t1, x, x);
+        mpz_mod (x, t1, n);
+        mpz_add_ui (x, x, a);
         mpz_sub (t1, x1, x);
         mpz_mul (t2, P, t1);
@@ -217 +216 @@
       for (i = 0; i < k; i++)
       {
-        if (p != 0)
-        {
-          mpz_powm_ui (x, x, p, n);
-          mpz_add_ui (x, x, a);
-        }
-        else
-        {
-          mpz_mul (t1, x, x);
-          mpz_mod (x, t1, n);
-          mpz_add_ui (x, x, a);
-        }
+        mpz_mul (t1, x, x);
+        mpz_mod (x, t1, n);
+        mpz_add_ui (x, x, a);
       }
       mpz_set (y, x);
@@ -235 +226 @@
     do
     {
-      if (p != 0)
-      {
-        mpz_powm_ui (y, y, p, n); mpz_add_ui (y, y, a);
-      }
-      else
-      {
-        mpz_mul (t1, y, y);
-        mpz_mod (y, t1, n);
-        mpz_add_ui (y, y, a);
-      }
+      mpz_mul (t1, y, y);
+      mpz_mod (y, t1, n);
+      mpz_add_ui (y, y, a);
       mpz_sub (t1, x1, y);
       mpz_gcd (t1, t1, n);
@@ -262 +246 @@
       while (a == 0);
 
-      factor_using_pollard_rho (t1, a, p, primes,multiplicities,index);
+      factor_using_pollard_rho (t1, a, primes,multiplicities,index);
     }
     else
@@ -292 +276 @@
         multiplicities[index++] = 1;
       }
+      mpz_set_ui(n,1);
       break;
     }
@@ -299 +284 @@
 }
 
-void factor (mpz_t t,lists primes,int *multiplicities,int &index)
+void factor (mpz_t t,lists primes,int *multiplicities,int &index,unsigned long bound)
 {
   unsigned int division_limit;
@@ -313 +298 @@
     division_limit = division_limit * division_limit;
 
-  factor_using_division (t, division_limit,primes,multiplicities,index);
-
-  if (mpz_cmp_ui (t, 1) != 0)
-  {
-    //if (flag_verbose > 0)
-    //{
-    //  printf ("[is number prime?] ");
-    //  fflush (stdout);
-    //}
-    if (mpz_probab_prime_p (t, 10))
-    {
-      setListEntry(primes, index, t);
-      multiplicities[index++] = 1;
-    }
-    else
-      factor_using_pollard_rho (t, 1L, 0,primes,multiplicities,index);
-  }
-  mpz_set_ui(t,1);
+  if (factor_using_division (t, division_limit,primes,multiplicities,index,bound))
+  {
+    if (mpz_cmp_ui (t, 1) != 0)
+    {
+      if (mpz_probab_prime_p (t, 10))
+      {
+        setListEntry(primes, index, t);
+        multiplicities[index++] = 1;
+        mpz_set_ui(t,1);
+      }
+      else
+        factor_using_pollard_rho (t, 1L, primes,multiplicities,index);
+    }
+  }
 }
 /* n and pBound are assumed to be bigint numbers */
-lists primeFactorisation(const number n, const number pBound)
+lists primeFactorisation(const number n, const int pBound)
 {
   int i;
@@ -340 +321 @@
   lists primes = (lists)omAllocBin(slists_bin); primes->Init(1000);
   int* multiplicities = (int*)omAlloc0(1000*sizeof(int));
-  int positive=1; int probTest = 0;
+  int positive=1;
 
   if (!n_IsZero(n, coeffs_BIGINT))
@@ -349 +330 @@
       mpz_neg(nn,nn);
     }
-    factor(nn,primes,multiplicities,index);
+    factor(nn,primes,multiplicities,index,pBound);
   }
 
@@ -372 +353 @@
 
   lists L=(lists)omAllocBin(slists_bin);
-  L->Init(4);
+  L->Init(3);
   if (positive==-1) mpz_neg(nn,nn);
   L->m[0].rtyp = LIST_CMD; L->m[0].data = (void*)primesL;
   L->m[1].rtyp = LIST_CMD; L->m[1].data = (void*)multiplicitiesL;
   setListEntry(L, 2, nn);
-  L->m[3].rtyp =  INT_CMD; L->m[3].data = (void*)probTest;
+
   mpz_clear(nn);
 

Singular/misc_ip.h

@@ -32 +32 @@
 
 /**
- * Divides 'n' as many times as possible by 'd' and returns the number
- * of divisions (without remainder) in 'times',
- * e.g., n = 48, d = 4, divTimes(n, d, t) = 3 produces n = 3, t = 2,
- *       since 48 = 4*4*3;
- * assumes that d is positive
- **/
-void divTimes(mpz_t n,   /**< [in]  a GMP number >= 0                      */
-              mpz_t d,   /**< [in]  the divisor, a GMP number >= 0         */
-              int* times /**< [out] number of divisions without remainder  */
-             );
-
-/**
  * Factorises a given bigint number n into its prime factors less
  * than or equal to a given bound, with corresponding multiplicities.
@@ -52 +40 @@
  * Also, when n is negative, m will contain the sign. If n is zero, m will
  * be zero.
- * The method returns a list L filled with four entries:
+ * The method returns a list L filled with three entries:
  * L[1] a list; L[1][i] contains the i-th prime factor of |n| as int or
  *                      bigint (sorted in ascending order),
  * L[2] a list; L[2][i] contains the multiplicity of L[1, i] in |n| as int
  * L[3] contains the remainder m as int or bigint, depending on the size,
- * L[4] 1 iff |m| is probably a prime, 0 otherwise
  *
- * We thus have: n = L[1][1]^L[2][1] * ... * L[1][k]^L[2][k] * L[1], where
+ * We thus have: n = L[1][1]^L[2][1] * ... * L[1][k]^L[2][k] * L[3], where
  * k is the number of mutually distinct prime factors (<= a provided non-
  * zero bound).
- * Note that for n = 0, L[2] and L[3] will be emtpy lists and L[4] will be
+ * Note that for n = 0, L[1] and L[2] will be emtpy lists and L[3] will be
  * zero.
  *
@@ -69 +56 @@
 lists primeFactorisation(
        const number n,     /**< [in]  the bigint > 0 to be factorised   */
-       const number pBound /**< [in]  bigint bound on the prime factors
+       const int pBound    /**< [in]  bound on the prime factors
                                       seeked                            */
                         );

Singular/table.h

@@ -609 +609 @@
 ,{D(jjRES),       MRES_CMD,       RESOLUTION_CMD, MODUL_CMD,  INT_CMD, ALLOW_PLURAL |ALLOW_RING}
 //,{D(nuMPResMat),  MPRES_CMD,      MODUL_CMD,      IDEAL_CMD,  INT_CMD, NO_PLURAL |ALLOW_RING}
-,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, BIGINT_CMD, ALLOW_PLURAL |ALLOW_RING}
-,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, BIGINT_CMD, ALLOW_PLURAL |ALLOW_RING}
-,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, NUMBER_CMD, ALLOW_PLURAL |ALLOW_RING}
-,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, NUMBER_CMD, ALLOW_PLURAL |ALLOW_RING}
+,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
+,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
 #ifdef HAVE_PLURAL
 ,{D(jjPlural_num_poly), NCALGEBRA_CMD,NONE,       POLY_CMD,   POLY_CMD  , NO_PLURAL |NO_RING}