Changeset 4c4340 in git for Singular/misc_ip.cc


Ignore:
Timestamp:
Oct 2, 2012, 12:25:11 PM (12 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
65b813a15472ad024e1b30ff11d4f966ec98b53a
Parents:
8b57f5e1c589185d90380ede5c9b3a0b2a44ff97
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2012-10-02 12:25:11+02:00
git-committer:
Martin Lee <martinlee84@web.de>2013-01-23 11:27:17+01:00
Message:
fix: primeFactorization redone with trial dvivision/Pollard-Rho from gmp

Conflicts:

	Singular/misc_ip.cc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/misc_ip.cc

    r8b57f5 r4c4340  
    3737
    3838
    39 void divTimes(mpz_t n, mpz_t d, int* times)
    40 {
    41   *times = 0;
    42   mpz_t r; mpz_init(r);
    43   mpz_t q; mpz_init(q);
    44   mpz_fdiv_qr(q, r, n, d);
    45   while (mpz_cmp_ui(r, 0) == 0)
    46   {
    47     (*times)++;
    48     mpz_set(n, q);
    49     mpz_fdiv_qr(q, r, n, d);
    50   }
    51   mpz_clear(r);
    52   mpz_clear(q);
    53 }
    54 
    55 void divTimes_ui(mpz_t n, unsigned long d, int* times)
    56 {
    57   *times = 0;
    58   mpz_t r; mpz_init(r);
    59   mpz_t q; mpz_init(q);
    60   mpz_fdiv_qr_ui(q, r, n, d);
    61   while (mpz_cmp_ui(r, 0) == 0)
    62   {
    63     (*times)++;
    64     mpz_set(n, q);
    65     mpz_fdiv_qr_ui(q, r, n, d);
    66   }
    67   mpz_clear(r);
    68   mpz_clear(q);
    69 }
    70 
    71 static inline void divTimes_ui_ui(unsigned long *n, unsigned long d, int* times)
    72 {
    73   *times = 0;
    74   unsigned long q=(*n) / d;
    75   unsigned long r=(*n) % d;
    76   while (r==0)
    77   {
    78     (*times)++;
    79     (*n)=q;
    80     q=(*n)/d; r=(*n)%d;
    81   }
    82 }
    83 
    8439void setListEntry(lists L, int index, mpz_t n)
    8540{ /* assumes n > 0 */
    … …  
    11469}
    11570
    116 /* true iff p is prime */
    117 /*
    118 bool isPrime(mpz_t p)
    119 {
    120   if (mpz_cmp_ui(p, 2) == 0) return true;
    121   if (mpz_cmp_ui(p, 3) == 0) return true;
    122   if (mpz_cmp_ui(p, 5) < 0)  return false;
    123 
    124   mpz_t d; mpz_init_set_ui(d, 5); int add = 2;
    125   mpz_t sr; mpz_init(sr); mpz_sqrt(sr, p);
    126   mpz_t r; mpz_init(r);
    127   while (mpz_cmp(d, sr) <= 0)
    128   {
    129     mpz_cdiv_r(r, p, d);
    130     if (mpz_cmp_ui(r, 0) == 0)
    131     {
    132       mpz_clear(d); mpz_clear(sr); mpz_clear(r);
    133       return false;
    134     }
    135     mpz_add_ui(d, d, add);
    136     add += 2; if (add == 6) add = 2;
    137   }
    138   mpz_clear(d); mpz_clear(sr); mpz_clear(r);
    139   return true;
    140 }
    141 */
    142 
    143 /* finds the next prime q, bound >= q >= p;
    144    in case of success, puts q into p;
    145    otherwise sets q = bound + 1;
    146    e.g. p = 24; nextPrime(p, 30) produces p = 29 (success),
    147         p = 24; nextPrime(p, 29) produces p = 29 (success),
    148         p = 24; nextPrime(p, 28) produces p = 29 (no success),
    149         p = 24; nextPrime(p, 27) produces p = 28 (no success) */
    150 /*
    151 void nextPrime(mpz_t p, mpz_t bound)
    152 {
    153   int add;
    154   mpz_t r; mpz_init(r); mpz_cdiv_r_ui(r, p, 6); // r = p mod 6, 0 <= r <= 5
    155   if (mpz_cmp_ui(r, 0) == 0) { mpz_add_ui(p, p, 1); add = 4; }
    156   if (mpz_cmp_ui(r, 1) == 0) {                      add = 4; }
    157   if (mpz_cmp_ui(r, 2) == 0) { mpz_add_ui(p, p, 3); add = 2; }
    158   if (mpz_cmp_ui(r, 3) == 0) { mpz_add_ui(p, p, 2); add = 2; }
    159   if (mpz_cmp_ui(r, 4) == 0) { mpz_add_ui(p, p, 1); add = 2; }
    160   if (mpz_cmp_ui(r, 5) == 0) {                      add = 2; }
    161 
    162   while (mpz_cmp(p, bound) <= 0)
    163   {
    164     if (isPrime(p)) { mpz_clear(r); return; }
    165     mpz_add_ui(p, p, add);
    166     add += 2; if (add == 6) add = 2;
    167   }
    168   mpz_set(p, bound);
    169   mpz_add_ui(p, p, 1);
    170   mpz_clear(r);
    171   return;
    172 }
    173 */
    174 
    175 
    176 
     71
     72/* Factoring with Pollard's rho method. stolen from GMP/demos */
     73static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
     74
     75void factor_using_division (mpz_t t, unsigned int limit,lists primes, int *multiplicities,int &index)
     76{
     77  mpz_t q, r;
     78  unsigned long int f;
     79  int ai;
     80  unsigned *addv = add;
     81  unsigned int failures;
     82
     83  mpz_init (q);
     84  mpz_init (r);
     85
     86  f = mpz_scan1 (t, 0);
     87  mpz_div_2exp (t, t, f);
     88  if (f>0)
     89  {
     90    setListEntry_ui(primes, index, 2);
     91    multiplicities[index++] = f;
     92  }
     93
     94  f=0;
     95  loop
     96  {
     97    mpz_tdiv_qr_ui (q, r, t, 3);
     98    if (mpz_cmp_ui (r, 0) != 0)
     99        break;
     100    mpz_set (t, q);
     101    f++;
     102  }
     103  if (f>0)
     104  {
     105    setListEntry_ui(primes, index, 3);
     106    multiplicities[index++] = f;
     107  }
     108  f=0;
     109  loop
     110  {
     111    mpz_tdiv_qr_ui (q, r, t, 5);
     112    if (mpz_cmp_ui (r, 0) != 0)
     113        break;
     114    mpz_set (t, q);
     115    f++;
     116  }
     117  if (f>0)
     118  {
     119    setListEntry_ui(primes, index, 5);
     120    multiplicities[index++] = f;
     121  }
     122
     123  failures = 0;
     124  f = 7;
     125  ai = 0;
     126  unsigned last_f=0;
     127  while (mpz_cmp_ui (t, 1) != 0)
     128  {
     129    mpz_tdiv_qr_ui (q, r, t, f);
     130    if (mpz_cmp_ui (r, 0) != 0)
     131    {
     132      f += addv[ai];
     133      if (mpz_cmp_ui (q, f) < 0)
     134        break;
     135      ai = (ai + 1) & 7;
     136      failures++;
     137      if (failures > limit)
     138        break;
     139    }
     140    else
     141    {
     142      mpz_swap (t, q);
     143      if (f!=last_f)
     144      {
     145        setListEntry_ui(primes, index, f);
     146        multiplicities[index]++;
     147        index++;
     148      }
     149      else
     150      {
     151        multiplicities[index-1]++;
     152      }
     153      last_f=f;
     154      failures = 0;
     155    }
     156  }
     157
     158  mpz_clears (q, r, NULL);
     159}
     160
     161void factor_using_pollard_rho (mpz_t n, unsigned long a, unsigned long p,lists primes, int * multiplicities,int &index)
     162{
     163  mpz_t x, x1, y, P;
     164  mpz_t t1, t2;
     165  mpz_t last_f;
     166  unsigned long long k, l, i;
     167
     168  mpz_init (t1);
     169  mpz_init (t2);
     170  mpz_init_set_si (last_f, 0);
     171  mpz_init_set_si (y, 2);
     172  mpz_init_set_si (x, 2);
     173  mpz_init_set_si (x1, 2);
     174  mpz_init_set_ui (P, 1);
     175  k = 1;
     176  l = 1;
     177
     178  while (mpz_cmp_ui (n, 1) != 0)
     179  {
     180    loop
     181    {
     182      do
     183      {
     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
     196        mpz_sub (t1, x1, x);
     197        mpz_mul (t2, P, t1);
     198        mpz_mod (P, t2, n);
     199
     200        if (k % 32 == 1)
     201        {
     202          mpz_gcd (t1, P, n);
     203          if (mpz_cmp_ui (t1, 1) != 0)
     204            goto factor_found;
     205          mpz_set (y, x);
     206        }
     207      }
     208      while (--k != 0);
     209
     210      mpz_gcd (t1, P, n);
     211      if (mpz_cmp_ui (t1, 1) != 0)
     212        goto factor_found;
     213
     214      mpz_set (x1, x);
     215      k = l;
     216      l = 2 * l;
     217      for (i = 0; i < k; i++)
     218      {
     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        }
     230      }
     231      mpz_set (y, x);
     232    }
     233
     234  factor_found:
     235    do
     236    {
     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      }
     247      mpz_sub (t1, x1, y);
     248      mpz_gcd (t1, t1, n);
     249    }
     250    while (mpz_cmp_ui (t1, 1) == 0);
     251
     252    mpz_divexact (n, n, t1);        /* divide by t1, before t1 is overwritten */
     253
     254    if (!mpz_probab_prime_p (t1, 10))
     255    {
     256      do
     257      {
     258        mp_limb_t a_limb;
     259        mpn_random (&a_limb, (mp_size_t) 1);
     260        a = a_limb;
     261      }
     262      while (a == 0);
     263
     264      factor_using_pollard_rho (t1, a, p, primes,multiplicities,index);
     265    }
     266    else
     267    {
     268      if (mpz_cmp(t1,last_f)==0)
     269      {
     270        multiplicities[index-1]++;
     271      }
     272      else
     273      {
     274        mpz_set(last_f,t1);
     275        setListEntry(primes, index, t1);
     276        multiplicities[index++] = 1;
     277      }
     278    }
     279    mpz_mod (x, x, n);
     280    mpz_mod (x1, x1, n);
     281    mpz_mod (y, y, n);
     282    if (mpz_probab_prime_p (n, 10))
     283    {
     284      if (mpz_cmp(n,last_f)==0)
     285      {
     286        multiplicities[index-1]++;
     287      }
     288      else
     289      {
     290        mpz_set(last_f,n);
     291        setListEntry(primes, index, n);
     292        multiplicities[index++] = 1;
     293      }
     294      break;
     295    }
     296  }
     297
     298  mpz_clears (P,t2,t1,x1,x,y,last_f,NULL);
     299}
     300
     301void factor (mpz_t t,lists primes,int *multiplicities,int &index)
     302{
     303  unsigned int division_limit;
     304
     305  if (mpz_sgn (t) == 0)
     306    return;
     307
     308  /* Set the trial division limit according the size of t.  */
     309  division_limit = mpz_sizeinbase (t, 2);
     310  if (division_limit > 1000)
     311    division_limit = 1000 * 1000;
     312  else
     313    division_limit = division_limit * division_limit;
     314
     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);
     333}
    177334/* n and pBound are assumed to be bigint numbers */
    178335lists primeFactorisation(const number n, const number pBound)
    179336{
     337  int i;
     338  int index=0;
    180339  mpz_t nn; number2mpz(n, nn);
    181   mpz_t pb; number2mpz(pBound, pb);
    182   mpz_t b; number2mpz(pBound, b);
    183   mpz_t p; mpz_init(p); int tt;
    184   mpz_t sr; mpz_init(sr); int index = 0; int add;
    185340  lists primes = (lists)omAllocBin(slists_bin); primes->Init(1000);
    186   int* multiplicities = new int[1000];
     341  int* multiplicities = (int*)omAlloc0(1000*sizeof(int));
    187342  int positive=1; int probTest = 0;
    188343
    … …  
    194349      mpz_neg(nn,nn);
    195350    }
    196     divTimes_ui(nn, 2, &tt);
    197     if (tt > 0)
    198     {
    199       setListEntry_ui(primes, index, 2);
    200       multiplicities[index++] = tt;
    201     }
    202 
    203     divTimes_ui(nn, 3, &tt);
    204     if (tt > 0)
    205     {
    206       setListEntry_ui(primes, index, 3);
    207       multiplicities[index++] = tt;
    208     }
    209 
    210     unsigned long p_ui=5; add = 2;
    211     BOOLEAN b_is_0=(mpz_cmp_ui(b, 0) == 0);
    212     BOOLEAN sr_sets_pb=FALSE;
    213     mpz_sqrt(sr, nn);
    214     // there are 3 possible limits, we take the minimum:
    215     // - argument pBound (if >0)
    216     // - sr = sqrt(nn)
    217     // - 1<<31
    218     unsigned long  limit=~(0L);
    219     if (b_is_0 || (mpz_cmp(pb, sr) > 0))
    220     {
    221       mpz_set(pb, sr);
    222       sr_sets_pb=TRUE;
    223     }
    224     if (mpz_cmp_ui(pb, limit)<0)
    225     {
    226      limit=mpz_get_ui(pb);
    227     }
    228     else
    229     {
    230       mpz_set_ui(pb,limit);
    231     }
    232     while (p_ui <=limit)
    233     {
    234       divTimes_ui(nn, p_ui, &tt);
    235       if (tt > 0)
    236       {
    237         setListEntry_ui(primes, index, p_ui);
    238         multiplicities[index++] = tt;
    239         //mpz_sqrt(sr, nn);
    240         //if ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
    241         if (mpz_size1(nn)<=2)
    242         {
    243           mpz_sqrt(sr, nn);
    244           if (sr_sets_pb || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
    245           unsigned long l=mpz_get_ui(sr);
    246           if (l<limit) { limit=l; }
    247           if (mpz_size1(nn)<=1)
    248           {
    249             unsigned long nn_ui=mpz_get_ui(nn);
    250             while (p_ui <=limit)
    251             {
    252               divTimes_ui_ui(&nn_ui, p_ui, &tt);
    253               if (tt > 0)
    254               {
    255                 setListEntry_ui(primes, index, p_ui);
    256                 multiplicities[index++] = tt;
    257                 if (nn_ui==1) break;
    258                 if (nn_ui<(limit/6)) { limit=nn_ui/6;}
    259               }
    260               p_ui +=add;
    261               //add += 2; if (add == 6) add = 2;
    262               add =2+2*(add==2);
    263             }
    264             mpz_set_ui(nn,nn_ui);
    265             break;
    266           }
    267         }
    268       }
    269       p_ui +=add;
    270       //add += 2; if (add == 6) add = 2;
    271       add =2+2*(add==2);
    272     }
    273     mpz_set_ui(p, p_ui);
    274     mpz_sqrt(sr, nn);
    275     if (b_is_0 || sr_sets_pb || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
    276     while (mpz_cmp(pb, p) >= 0)
    277     {
    278       divTimes(nn, p, &tt);
    279       if (tt > 0)
    280       {
    281         setListEntry(primes, index, p);
    282         multiplicities[index++] = tt;
    283         if (mpz_cmp_ui(nn,1)==0) break;
    284         mpz_sqrt(sr, nn);
    285         if (b_is_0 || sr_sets_pb || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
    286       }
    287       mpz_add_ui(p, p, add);
    288       //add += 2; if (add == 6) add = 2;
    289       add =2+2*(add==2);
    290     }
    291     if ((mpz_cmp_ui(nn, 1) > 0) &&
    292         (b_is_0 || (mpz_cmp(nn, b) <= 0)))
    293     {
    294       setListEntry(primes, index, nn);
    295       multiplicities[index++] = 1;
    296       mpz_set_ui(nn, 1);
    297     }
    298     if ((mpz_cmp_ui(nn, 1) > 0) && (mpz_probab_prime_p(nn, 25) != 0))
    299       probTest = 1;
     351    factor(nn,primes,multiplicities,index);
    300352  }
    301353
    302354  lists primesL = (lists)omAllocBin(slists_bin);
    303355  primesL->Init(index);
    304   for (int i = 0; i < index; i++)
     356  for (i = 0; i < index; i++)
    305357  {
    306358    primesL->m[i].rtyp = primes->m[i].rtyp;
    … …  
    312364  lists multiplicitiesL = (lists)omAllocBin(slists_bin);
    313365  multiplicitiesL->Init(index);
    314   for (int i = 0; i < index; i++)
     366  for (i = 0; i < index; i++)
    315367  {
    316368    multiplicitiesL->m[i].rtyp = INT_CMD;
    … …  
    326378  setListEntry(L, 2, nn);
    327379  L->m[3].rtyp =  INT_CMD; L->m[3].data = (void*)probTest;
    328   mpz_clear(nn); mpz_clear(pb); mpz_clear(b); mpz_clear(p); mpz_clear(sr);
     380  mpz_clear(nn);
    329381
    330382  return L;
Note: See TracChangeset for help on using the changeset viewer.