Changeset 65b813 in git


Ignore:
Timestamp:
Oct 3, 2012, 1:32:09 PM (12 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4188d308699580d975efd0f6cca8dcb41c396f70')
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:
5 edited

Legend:

Unmodified
Added
Removed
  • Singular/iparith.cc

    r4c4340 r65b813  
    31523152static BOOLEAN jjPFAC2(leftv res, leftv u, leftv v)
    31533153{
    3154   number n1; number n2; number temp; int i;
     3154  number n1; int i;
    31553155
    31563156  if ((u->Typ() == BIGINT_CMD) ||
    31573157     ((u->Typ() == NUMBER_CMD) && rField_is_Q(currRing)))
    31583158  {
    3159     temp = (number)u->Data();
    3160     n1 = n_Copy(temp,coeffs_BIGINT);
     3159    n1 = (number)u->CopyD();
    31613160  }
    31623161  else if (u->Typ() == INT_CMD)
     
    31673166  else
    31683167  {
    3169     WerrorS("wrong type: expected int, bigint, or number as 1st argument");
    3170     return TRUE;
    3171   }
    3172 
    3173   if ((v->Typ() == BIGINT_CMD) ||
    3174      ((v->Typ() == NUMBER_CMD) && rField_is_Q(currRing)))
    3175   {
    3176     temp = (number)v->Data();
    3177     n2 = n_Copy(temp,coeffs_BIGINT);
    3178   }
    3179   else if (v->Typ() == INT_CMD)
    3180   {
    3181     i = (int)(long)v->Data();
    3182     n2 = n_Init(i, coeffs_BIGINT);
    3183   }
    3184   else
    3185   {
    3186     WerrorS("wrong type: expected int, bigint, or number as 2nd argument");
    3187     return TRUE;
    3188   }
    3189 
    3190   lists l = primeFactorisation(n1, n2);
    3191   n_Delete(&n1, coeffs_BIGINT); n_Delete(&n2, coeffs_BIGINT);
     3168    return TRUE;
     3169  }
     3170
     3171  i = (int)(long)v->Data();
     3172
     3173  lists l = primeFactorisation(n1, i);
     3174  n_Delete(&n1, coeffs_BIGINT);
    31923175  res->data = (char*)l;
    31933176  return FALSE;
  • Singular/misc_ip.cc

    r4c4340 r65b813  
    7373static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
    7474
    75 void factor_using_division (mpz_t t, unsigned int limit,lists primes, int *multiplicities,int &index)
     75int factor_using_division (mpz_t t, unsigned int limit,lists primes, int *multiplicities,int &index, unsigned long bound)
    7676{
    7777  mpz_t q, r;
     
    8080  unsigned *addv = add;
    8181  unsigned int failures;
     82  int bound_not_reached=1;
    8283
    8384  mpz_init (q);
     
    124125  f = 7;
    125126  ai = 0;
    126   unsigned last_f=0;
     127  unsigned long last_f=0;
    127128  while (mpz_cmp_ui (t, 1) != 0)
    128129  {
     
    131132    {
    132133      f += addv[ai];
    133       if (mpz_cmp_ui (q, f) < 0)
     134      if (mpz_cmp_ui (t, f) < 0)
    134135        break;
    135136      ai = (ai + 1) & 7;
     
    137138      if (failures > limit)
    138139        break;
     140      if ((bound!=0) && (f>bound))
     141      {
     142        bound_not_reached=0;
     143        break;
     144      }
    139145    }
    140146    else
     
    157163
    158164  mpz_clears (q, r, NULL);
    159 }
    160 
    161 void factor_using_pollard_rho (mpz_t n, unsigned long a, unsigned long p,lists primes, int * multiplicities,int &index)
     165  //printf("bound=%d,f=%d,failures=%d, reached=%d\n",bound,f,failures,bound_not_reached);
     166  return bound_not_reached;
     167}
     168
     169void factor_using_pollard_rho (mpz_t n, unsigned long a, lists primes, int * multiplicities,int &index)
    162170{
    163171  mpz_t x, x1, y, P;
     
    182190      do
    183191      {
    184         if (p != 0)
    185         {
    186           mpz_powm_ui (x, x, p, n);
    187           mpz_add_ui (x, x, a);
    188         }
    189         else
    190         {
    191           mpz_mul (t1, x, x);
    192           mpz_mod (x, t1, n);
    193           mpz_add_ui (x, x, a);
    194         }
    195 
     192        mpz_mul (t1, x, x);
     193        mpz_mod (x, t1, n);
     194        mpz_add_ui (x, x, a);
    196195        mpz_sub (t1, x1, x);
    197196        mpz_mul (t2, P, t1);
     
    217216      for (i = 0; i < k; i++)
    218217      {
    219         if (p != 0)
    220         {
    221           mpz_powm_ui (x, x, p, n);
    222           mpz_add_ui (x, x, a);
    223         }
    224         else
    225         {
    226           mpz_mul (t1, x, x);
    227           mpz_mod (x, t1, n);
    228           mpz_add_ui (x, x, a);
    229         }
     218        mpz_mul (t1, x, x);
     219        mpz_mod (x, t1, n);
     220        mpz_add_ui (x, x, a);
    230221      }
    231222      mpz_set (y, x);
     
    235226    do
    236227    {
    237       if (p != 0)
    238       {
    239         mpz_powm_ui (y, y, p, n); mpz_add_ui (y, y, a);
    240       }
    241       else
    242       {
    243         mpz_mul (t1, y, y);
    244         mpz_mod (y, t1, n);
    245         mpz_add_ui (y, y, a);
    246       }
     228      mpz_mul (t1, y, y);
     229      mpz_mod (y, t1, n);
     230      mpz_add_ui (y, y, a);
    247231      mpz_sub (t1, x1, y);
    248232      mpz_gcd (t1, t1, n);
     
    262246      while (a == 0);
    263247
    264       factor_using_pollard_rho (t1, a, p, primes,multiplicities,index);
     248      factor_using_pollard_rho (t1, a, primes,multiplicities,index);
    265249    }
    266250    else
     
    292276        multiplicities[index++] = 1;
    293277      }
     278      mpz_set_ui(n,1);
    294279      break;
    295280    }
     
    299284}
    300285
    301 void factor (mpz_t t,lists primes,int *multiplicities,int &index)
     286void factor (mpz_t t,lists primes,int *multiplicities,int &index,unsigned long bound)
    302287{
    303288  unsigned int division_limit;
     
    313298    division_limit = division_limit * division_limit;
    314299
    315   factor_using_division (t, division_limit,primes,multiplicities,index);
    316 
    317   if (mpz_cmp_ui (t, 1) != 0)
    318   {
    319     //if (flag_verbose > 0)
    320     //{
    321     //  printf ("[is number prime?] ");
    322     //  fflush (stdout);
    323     //}
    324     if (mpz_probab_prime_p (t, 10))
    325     {
    326       setListEntry(primes, index, t);
    327       multiplicities[index++] = 1;
    328     }
    329     else
    330       factor_using_pollard_rho (t, 1L, 0,primes,multiplicities,index);
    331   }
    332   mpz_set_ui(t,1);
     300  if (factor_using_division (t, division_limit,primes,multiplicities,index,bound))
     301  {
     302    if (mpz_cmp_ui (t, 1) != 0)
     303    {
     304      if (mpz_probab_prime_p (t, 10))
     305      {
     306        setListEntry(primes, index, t);
     307        multiplicities[index++] = 1;
     308        mpz_set_ui(t,1);
     309      }
     310      else
     311        factor_using_pollard_rho (t, 1L, primes,multiplicities,index);
     312    }
     313  }
    333314}
    334315/* n and pBound are assumed to be bigint numbers */
    335 lists primeFactorisation(const number n, const number pBound)
     316lists primeFactorisation(const number n, const int pBound)
    336317{
    337318  int i;
     
    340321  lists primes = (lists)omAllocBin(slists_bin); primes->Init(1000);
    341322  int* multiplicities = (int*)omAlloc0(1000*sizeof(int));
    342   int positive=1; int probTest = 0;
     323  int positive=1;
    343324
    344325  if (!n_IsZero(n, coeffs_BIGINT))
     
    349330      mpz_neg(nn,nn);
    350331    }
    351     factor(nn,primes,multiplicities,index);
     332    factor(nn,primes,multiplicities,index,pBound);
    352333  }
    353334
     
    372353
    373354  lists L=(lists)omAllocBin(slists_bin);
    374   L->Init(4);
     355  L->Init(3);
    375356  if (positive==-1) mpz_neg(nn,nn);
    376357  L->m[0].rtyp = LIST_CMD; L->m[0].data = (void*)primesL;
    377358  L->m[1].rtyp = LIST_CMD; L->m[1].data = (void*)multiplicitiesL;
    378359  setListEntry(L, 2, nn);
    379   L->m[3].rtyp =  INT_CMD; L->m[3].data = (void*)probTest;
     360
    380361  mpz_clear(nn);
    381362
  • Singular/misc_ip.h

    r4c4340 r65b813  
    3232
    3333/**
    34  * Divides 'n' as many times as possible by 'd' and returns the number
    35  * of divisions (without remainder) in 'times',
    36  * e.g., n = 48, d = 4, divTimes(n, d, t) = 3 produces n = 3, t = 2,
    37  *       since 48 = 4*4*3;
    38  * assumes that d is positive
    39  **/
    40 void divTimes(mpz_t n,   /**< [in]  a GMP number >= 0                      */
    41               mpz_t d,   /**< [in]  the divisor, a GMP number >= 0         */
    42               int* times /**< [out] number of divisions without remainder  */
    43              );
    44 
    45 /**
    4634 * Factorises a given bigint number n into its prime factors less
    4735 * than or equal to a given bound, with corresponding multiplicities.
     
    5240 * Also, when n is negative, m will contain the sign. If n is zero, m will
    5341 * be zero.
    54  * The method returns a list L filled with four entries:
     42 * The method returns a list L filled with three entries:
    5543 * L[1] a list; L[1][i] contains the i-th prime factor of |n| as int or
    5644 *                      bigint (sorted in ascending order),
    5745 * L[2] a list; L[2][i] contains the multiplicity of L[1, i] in |n| as int
    5846 * L[3] contains the remainder m as int or bigint, depending on the size,
    59  * L[4] 1 iff |m| is probably a prime, 0 otherwise
    6047 *
    61  * We thus have: n = L[1][1]^L[2][1] * ... * L[1][k]^L[2][k] * L[1], where
     48 * We thus have: n = L[1][1]^L[2][1] * ... * L[1][k]^L[2][k] * L[3], where
    6249 * k is the number of mutually distinct prime factors (<= a provided non-
    6350 * zero bound).
    64  * Note that for n = 0, L[2] and L[3] will be emtpy lists and L[4] will be
     51 * Note that for n = 0, L[1] and L[2] will be emtpy lists and L[3] will be
    6552 * zero.
    6653 *
     
    6956lists primeFactorisation(
    7057       const number n,     /**< [in]  the bigint > 0 to be factorised   */
    71        const number pBound /**< [in]  bigint bound on the prime factors
     58       const int pBound    /**< [in] bound on the prime factors
    7259                                      seeked                            */
    7360                        );
  • Singular/table.h

    r4c4340 r65b813  
    609609,{D(jjRES),       MRES_CMD,       RESOLUTION_CMD, MODUL_CMD,  INT_CMD, ALLOW_PLURAL |ALLOW_RING}
    610610//,{D(nuMPResMat),  MPRES_CMD,      MODUL_CMD,      IDEAL_CMD,  INT_CMD, NO_PLURAL |ALLOW_RING}
    611 ,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, BIGINT_CMD, ALLOW_PLURAL |ALLOW_RING}
    612 ,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, BIGINT_CMD, ALLOW_PLURAL |ALLOW_RING}
    613 ,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, NUMBER_CMD, ALLOW_PLURAL |ALLOW_RING}
    614 ,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, NUMBER_CMD, ALLOW_PLURAL |ALLOW_RING}
     611,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
     612,{D(jjPFAC2),     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
    615613#ifdef HAVE_PLURAL
    616614,{D(jjPlural_num_poly), NCALGEBRA_CMD,NONE,       POLY_CMD,   POLY_CMD  , NO_PLURAL |NO_RING}
Note: See TracChangeset for help on using the changeset viewer.