Changeset 511f2d in git

Timestamp:

Jan 23, 2013, 5:57:20 PM (10 years ago)

Author:

Oleksandr Motsak <malex984@…>

Branches:

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

Children:

4def554e55ce8709cb931c050a80a944fe9a7633

Parents:

8b57f5e1c589185d90380ede5c9b3a0b2a44ff97fd85113e6317e9d03532a7ec8affbc338a502ea3

Message:

Merge pull request #256 from mmklee/sw_fix

Sw fix

Files:

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

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

@@ -37 +37 @@
 
 
-void divTimes(mpz_t n, mpz_t d, int* times)
-{
-  *times = 0;
-  mpz_t r; mpz_init(r);
-  mpz_t q; mpz_init(q);
-  mpz_fdiv_qr(q, r, n, d);
-  while (mpz_cmp_ui(r, 0) == 0)
-  {
-    (*times)++;
-    mpz_set(n, q);
-    mpz_fdiv_qr(q, r, n, d);
-  }
-  mpz_clear(r);
-  mpz_clear(q);
-}
-
-void divTimes_ui(mpz_t n, unsigned long d, int* times)
-{
-  *times = 0;
-  mpz_t r; mpz_init(r);
-  mpz_t q; mpz_init(q);
-  mpz_fdiv_qr_ui(q, r, n, d);
-  while (mpz_cmp_ui(r, 0) == 0)
-  {
-    (*times)++;
-    mpz_set(n, q);
-    mpz_fdiv_qr_ui(q, r, n, d);
-  }
-  mpz_clear(r);
-  mpz_clear(q);
-}
-
-static inline void divTimes_ui_ui(unsigned long *n, unsigned long d, int* times)
-{
-  *times = 0;
-  unsigned long q=(*n) / d;
-  unsigned long r=(*n) % d;
-  while (r==0)
-  {
-    (*times)++;
-    (*n)=q;
-    q=(*n)/d; r=(*n)%d;
-  }
-}
-
 void setListEntry(lists L, int index, mpz_t n)
 { /* assumes n > 0 */
@@ -114 +69 @@
 }
 
-/* true iff p is prime */
-/*
-bool isPrime(mpz_t p)
-{
-  if (mpz_cmp_ui(p, 2) == 0) return true;
-  if (mpz_cmp_ui(p, 3) == 0) return true;
-  if (mpz_cmp_ui(p, 5) < 0)  return false;
-
-  mpz_t d; mpz_init_set_ui(d, 5); int add = 2;
-  mpz_t sr; mpz_init(sr); mpz_sqrt(sr, p);
-  mpz_t r; mpz_init(r);
-  while (mpz_cmp(d, sr) <= 0)
-  {
-    mpz_cdiv_r(r, p, d);
-    if (mpz_cmp_ui(r, 0) == 0)
-    {
-      mpz_clear(d); mpz_clear(sr); mpz_clear(r);
-      return false;
-    }
-    mpz_add_ui(d, d, add);
-    add += 2; if (add == 6) add = 2;
-  }
-  mpz_clear(d); mpz_clear(sr); mpz_clear(r);
-  return true;
-}
-*/
-
-/* finds the next prime q, bound >= q >= p;
-   in case of success, puts q into p;
-   otherwise sets q = bound + 1;
-   e.g. p = 24; nextPrime(p, 30) produces p = 29 (success),
-        p = 24; nextPrime(p, 29) produces p = 29 (success),
-        p = 24; nextPrime(p, 28) produces p = 29 (no success),
-        p = 24; nextPrime(p, 27) produces p = 28 (no success) */
-/*
-void nextPrime(mpz_t p, mpz_t bound)
-{
-  int add;
-  mpz_t r; mpz_init(r); mpz_cdiv_r_ui(r, p, 6); // r = p mod 6, 0 <= r <= 5
-  if (mpz_cmp_ui(r, 0) == 0) { mpz_add_ui(p, p, 1); add = 4; }
-  if (mpz_cmp_ui(r, 1) == 0) {                      add = 4; }
-  if (mpz_cmp_ui(r, 2) == 0) { mpz_add_ui(p, p, 3); add = 2; }
-  if (mpz_cmp_ui(r, 3) == 0) { mpz_add_ui(p, p, 2); add = 2; }
-  if (mpz_cmp_ui(r, 4) == 0) { mpz_add_ui(p, p, 1); add = 2; }
-  if (mpz_cmp_ui(r, 5) == 0) {                      add = 2; }
-
-  while (mpz_cmp(p, bound) <= 0)
-  {
-    if (isPrime(p)) { mpz_clear(r); return; }
-    mpz_add_ui(p, p, add);
-    add += 2; if (add == 6) add = 2;
-  }
-  mpz_set(p, bound);
-  mpz_add_ui(p, p, 1);
-  mpz_clear(r);
-  return;
-}
-*/
-
-
-
+/* Factoring with Pollard's rho method. stolen from GMP/demos */
+static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
+
+static int factor_using_division (mpz_t t, unsigned int limit,lists primes, int *multiplicities,int &index, unsigned long bound)
+{
+  mpz_t q, r;
+  unsigned long int f;
+  int ai;
+  unsigned *addv = add;
+  unsigned int failures;
+  int bound_not_reached=1;
+
+  mpz_init (q);
+  mpz_init (r);
+
+  f = mpz_scan1 (t, 0);
+  mpz_div_2exp (t, t, f);
+  if (f>0)
+  {
+    setListEntry_ui(primes, index, 2);
+    multiplicities[index++] = f;
+  }
+
+  f=0;
+  loop
+  {
+    mpz_tdiv_qr_ui (q, r, t, 3);
+    if (mpz_cmp_ui (r, 0) != 0)
+        break;
+    mpz_set (t, q);
+    f++;
+  }
+  if (f>0)
+  {
+    setListEntry_ui(primes, index, 3);
+    multiplicities[index++] = f;
+  }
+  f=0;
+  loop
+  {
+    mpz_tdiv_qr_ui (q, r, t, 5);
+    if (mpz_cmp_ui (r, 0) != 0)
+        break;
+    mpz_set (t, q);
+    f++;
+  }
+  if (f>0)
+  {
+    setListEntry_ui(primes, index, 5);
+    multiplicities[index++] = f;
+  }
+
+  failures = 0;
+  f = 7;
+  ai = 0;
+  unsigned long last_f=0;
+  while (mpz_cmp_ui (t, 1) != 0)
+  {
+    mpz_tdiv_qr_ui (q, r, t, f);
+    if (mpz_cmp_ui (r, 0) != 0)
+    {
+      f += addv[ai];
+      if (mpz_cmp_ui (t, f) < 0)
+        break;
+      ai = (ai + 1) & 7;
+      failures++;
+      if (failures > limit)
+        break;
+      if ((bound!=0) && (f>bound))
+      {
+        bound_not_reached=0;
+        break;
+      }
+    }
+    else
+    {
+      mpz_swap (t, q);
+      if (f!=last_f)
+      {
+        setListEntry_ui(primes, index, f);
+        multiplicities[index]++;
+        index++;
+      }
+      else
+      {
+        multiplicities[index-1]++;
+      }
+      last_f=f;
+      failures = 0;
+    }
+  }
+
+  mpz_clear (q);
+  mpz_clear (r);
+  //printf("bound=%d,f=%d,failures=%d, reached=%d\n",bound,f,failures,bound_not_reached);
+  return bound_not_reached;
+}
+
+static void factor_using_pollard_rho (mpz_t n, unsigned long a, lists primes, int * multiplicities,int &index)
+{
+  mpz_t x, x1, y, P;
+  mpz_t t1, t2;
+  mpz_t last_f;
+  unsigned long long k, l, i;
+
+  mpz_init (t1);
+  mpz_init (t2);
+  mpz_init_set_si (last_f, 0);
+  mpz_init_set_si (y, 2);
+  mpz_init_set_si (x, 2);
+  mpz_init_set_si (x1, 2);
+  mpz_init_set_ui (P, 1);
+  k = 1;
+  l = 1;
+
+  while (mpz_cmp_ui (n, 1) != 0)
+  {
+    loop
+    {
+      do
+      {
+        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);
+        mpz_mod (P, t2, n);
+
+        if (k % 32 == 1)
+        {
+          mpz_gcd (t1, P, n);
+          if (mpz_cmp_ui (t1, 1) != 0)
+            goto factor_found;
+          mpz_set (y, x);
+        }
+      }
+      while (--k != 0);
+
+      mpz_gcd (t1, P, n);
+      if (mpz_cmp_ui (t1, 1) != 0)
+        goto factor_found;
+
+      mpz_set (x1, x);
+      k = l;
+      l = 2 * l;
+      for (i = 0; i < k; i++)
+      {
+        mpz_mul (t1, x, x);
+        mpz_mod (x, t1, n);
+        mpz_add_ui (x, x, a);
+      }
+      mpz_set (y, x);
+    }
+
+  factor_found:
+    do
+    {
+      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);
+    }
+    while (mpz_cmp_ui (t1, 1) == 0);
+
+    mpz_divexact (n, n, t1);        /* divide by t1, before t1 is overwritten */
+
+    if (!mpz_probab_prime_p (t1, 10))
+    {
+      do
+      {
+        mp_limb_t a_limb;
+        mpn_random (&a_limb, (mp_size_t) 1);
+        a = a_limb;
+      }
+      while (a == 0);
+
+      factor_using_pollard_rho (t1, a, primes,multiplicities,index);
+    }
+    else
+    {
+      if (mpz_cmp(t1,last_f)==0)
+      {
+        multiplicities[index-1]++;
+      }
+      else
+      {
+        mpz_set(last_f,t1);
+        setListEntry(primes, index, t1);
+        multiplicities[index++] = 1;
+      }
+    }
+    mpz_mod (x, x, n);
+    mpz_mod (x1, x1, n);
+    mpz_mod (y, y, n);
+    if (mpz_probab_prime_p (n, 10))
+    {
+      if (mpz_cmp(n,last_f)==0)
+      {
+        multiplicities[index-1]++;
+      }
+      else
+      {
+        mpz_set(last_f,n);
+        setListEntry(primes, index, n);
+        multiplicities[index++] = 1;
+      }
+      mpz_set_ui(n,1);
+      break;
+    }
+  }
+
+  mpz_clear (P);
+  mpz_clear (t2);
+  mpz_clear (t1);
+  mpz_clear (x1);
+  mpz_clear (x);
+  mpz_clear (y);
+  mpz_clear (last_f);
+}
+
+static void factor_gmp (mpz_t t,lists primes,int *multiplicities,int &index,unsigned long bound)
+{
+  unsigned int division_limit;
+
+  if (mpz_sgn (t) == 0)
+    return;
+
+  /* Set the trial division limit according the size of t.  */
+  division_limit = mpz_sizeinbase (t, 2);
+  if (division_limit > 1000)
+    division_limit = 1000 * 1000;
+  else
+    division_limit = division_limit * division_limit;
+
+  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;
+  int index=0;
   mpz_t nn; number2mpz(n, nn);
-  mpz_t pb; number2mpz(pBound, pb);
-  mpz_t b; number2mpz(pBound, b);
-  mpz_t p; mpz_init(p); int tt;
-  mpz_t sr; mpz_init(sr); int index = 0; int add;
   lists primes = (lists)omAllocBin(slists_bin); primes->Init(1000);
-  int* multiplicities = new int[1000];
-  int positive=1; int probTest = 0;
+  int* multiplicities = (int*)omAlloc0(1000*sizeof(int));
+  int positive=1;
 
   if (!n_IsZero(n, coeffs_BIGINT))
@@ -194 +336 @@
       mpz_neg(nn,nn);
     }
-    divTimes_ui(nn, 2, &tt);
-    if (tt > 0)
-    {
-      setListEntry_ui(primes, index, 2);
-      multiplicities[index++] = tt;
-    }
-
-    divTimes_ui(nn, 3, &tt);
-    if (tt > 0)
-    {
-      setListEntry_ui(primes, index, 3);
-      multiplicities[index++] = tt;
-    }
-
-    unsigned long p_ui=5; add = 2;
-    BOOLEAN b_is_0=(mpz_cmp_ui(b, 0) == 0);
-    BOOLEAN sr_sets_pb=FALSE;
-    mpz_sqrt(sr, nn);
-    // there are 3 possible limits, we take the minimum:
-    // - argument pBound (if >0)
-    // - sr = sqrt(nn)
-    // - 1<<31
-    unsigned long  limit=~(0L);
-    if (b_is_0 || (mpz_cmp(pb, sr) > 0))
-    {
-      mpz_set(pb, sr);
-      sr_sets_pb=TRUE;
-    }
-    if (mpz_cmp_ui(pb, limit)<0)
-    {
-     limit=mpz_get_ui(pb);
-    }
-    else
-    {
-      mpz_set_ui(pb,limit);
-    }
-    while (p_ui <=limit)
-    {
-      divTimes_ui(nn, p_ui, &tt);
-      if (tt > 0)
-      {
-        setListEntry_ui(primes, index, p_ui);
-        multiplicities[index++] = tt;
-        //mpz_sqrt(sr, nn);
-        //if ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
-        if (mpz_size1(nn)<=2)
-        {
-          mpz_sqrt(sr, nn);
-          if (sr_sets_pb || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
-          unsigned long l=mpz_get_ui(sr);
-          if (l<limit) { limit=l; }
-          if (mpz_size1(nn)<=1)
-          {
-            unsigned long nn_ui=mpz_get_ui(nn);
-            while (p_ui <=limit)
-            {
-              divTimes_ui_ui(&nn_ui, p_ui, &tt);
-              if (tt > 0)
-              {
-                setListEntry_ui(primes, index, p_ui);
-                multiplicities[index++] = tt;
-                if (nn_ui==1) break;
-                if (nn_ui<(limit/6)) { limit=nn_ui/6;}
-              }
-              p_ui +=add;
-              //add += 2; if (add == 6) add = 2;
-              add =2+2*(add==2);
-            }
-            mpz_set_ui(nn,nn_ui);
-            break;
-          }
-        }
-      }
-      p_ui +=add;
-      //add += 2; if (add == 6) add = 2;
-      add =2+2*(add==2);
-    }
-    mpz_set_ui(p, p_ui);
-    mpz_sqrt(sr, nn);
-    if (b_is_0 || sr_sets_pb || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
-    while (mpz_cmp(pb, p) >= 0)
-    {
-      divTimes(nn, p, &tt);
-      if (tt > 0)
-      {
-        setListEntry(primes, index, p);
-        multiplicities[index++] = tt;
-        if (mpz_cmp_ui(nn,1)==0) break;
-        mpz_sqrt(sr, nn);
-        if (b_is_0 || sr_sets_pb || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
-      }
-      mpz_add_ui(p, p, add);
-      //add += 2; if (add == 6) add = 2;
-      add =2+2*(add==2);
-    }
-    if ((mpz_cmp_ui(nn, 1) > 0) &&
-        (b_is_0 || (mpz_cmp(nn, b) <= 0)))
-    {
-      setListEntry(primes, index, nn);
-      multiplicities[index++] = 1;
-      mpz_set_ui(nn, 1);
-    }
-    if ((mpz_cmp_ui(nn, 1) > 0) && (mpz_probab_prime_p(nn, 25) != 0))
-      probTest = 1;
+    factor_gmp(nn,primes,multiplicities,index,pBound);
   }
 
   lists primesL = (lists)omAllocBin(slists_bin);
   primesL->Init(index);
-  for (int i = 0; i < index; i++)
+  for (i = 0; i < index; i++)
   {
     primesL->m[i].rtyp = primes->m[i].rtyp;
     primesL->m[i].data = primes->m[i].data;
-  }
-  omFreeSize((ADDRESS)primes->m, (primes->nr + 1) * sizeof(sleftv));
-  omFreeBin((ADDRESS)primes, slists_bin);
+    primes->m[i].rtyp=0;
+    primes->m[i].data=NULL;
+  }
+  primes->Clean(NULL);
 
   lists multiplicitiesL = (lists)omAllocBin(slists_bin);
   multiplicitiesL->Init(index);
-  for (int i = 0; i < index; i++)
+  for (i = 0; i < index; i++)
   {
     multiplicitiesL->m[i].rtyp = INT_CMD;
     multiplicitiesL->m[i].data = (void*)multiplicities[i];
   }
-  delete[] multiplicities;
+  omFree(multiplicities);
 
   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); mpz_clear(pb); mpz_clear(b); mpz_clear(p); mpz_clear(sr);
+
+  mpz_clear(nn);
 
   return L;

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}

factory/cf_gcd_smallp.cc

@@ -4606 +4606 @@
     maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
   }
-  else if (p < 50 && algExtension && !CFFactory::gettype() == GaloisFieldDomain)
+  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
   {
     int d= degree (getMipo (a));