Changeset 555b03 in git


Ignore:
Timestamp:
Aug 30, 2010, 4:29:37 PM (14 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
5e016e27b98b70a42edca584335f57145f3a8dc0
Parents:
f9bb6862e22f14387a67696733ad2d1024a69ade
Message:
primefactors via GMP, additional prob. check for unfactorized remainder

git-svn-id: file:///usr/local/Singular/svn/trunk@13153 2c84dea3-7e68-4137-9b89-c4e89433aadc
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • Singular/iparith.cc

    rf9bb68 r555b03  
    31873187static BOOLEAN jjPFAC2(leftv res, leftv u, leftv v)
    31883188{
    3189   number m = (number)u->Data();
    3190   int bound = (int)(long)v->Data();
    3191   number n;
    3192   if ((u->Typ() == BIGINT_CMD)
    3193   || ((u->Typ()==NUMBER_CMD) && rField_is_Q()))
    3194   {
    3195     n = nlCopy(m);
     3189  number n1; number n2; number temp; int i;
     3190
     3191  if ((u->Typ() == BIGINT_CMD) ||
     3192     ((u->Typ() == NUMBER_CMD) && rField_is_Q()))
     3193  {
     3194    temp = (number)u->Data();
     3195    n1 = nlCopy(temp);
     3196  }
     3197  else if (u->Typ() == INT_CMD)
     3198  {
     3199    i = (int)(long)u->Data();
     3200    n1 = nlInit(i, NULL);
    31963201  }
    31973202  else
    31983203  {
    3199     /* then, we expect some other type of number inside m;
    3200        it needs to be converted to bigint; we do this via 'int' */
    3201     int i = n_Int(m, currRing);
    3202     n = nlInit(i, NULL);
    3203   }
    3204   lists l = primeFactorisation(n, bound);
    3205   nlDelete(&n,NULL);
    3206   res->data=(char*)l;
     3204    WerrorS("wrong type: expected int, bigint, or number as 1st argument");
     3205    return TRUE;
     3206  }
     3207
     3208  if ((v->Typ() == BIGINT_CMD) ||
     3209     ((v->Typ() == NUMBER_CMD) && rField_is_Q()))
     3210  {
     3211    temp = (number)v->Data();
     3212    n2 = nlCopy(temp);
     3213  }
     3214  else if (v->Typ() == INT_CMD)
     3215  {
     3216    i = (int)(long)v->Data();
     3217    n2 = nlInit(i, NULL);
     3218  }
     3219  else
     3220  {
     3221    WerrorS("wrong type: expected int, bigint, or number as 2nd argument");
     3222    return TRUE;
     3223  }
     3224
     3225  lists l = primeFactorisation(n1, n2);
     3226  nlDelete(&n1, NULL); nlDelete(&n2, NULL);
     3227  res->data = (char*)l;
    32073228  return FALSE;
    32083229}
     
    37213742,{jjRES,       MRES_CMD,       RESOLUTION_CMD, MODUL_CMD,  INT_CMD, ALLOW_PLURAL |ALLOW_RING}
    37223743//,{nuMPResMat,  MPRES_CMD,      MODUL_CMD,      IDEAL_CMD,  INT_CMD, NO_PLURAL |ALLOW_RING}
    3723 ,{jjPFAC2,     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
    3724 ,{jjPFAC2,     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
     3744,{jjPFAC2,     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, BIGINT_CMD, ALLOW_PLURAL |ALLOW_RING}
     3745,{jjPFAC2,     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, BIGINT_CMD, ALLOW_PLURAL |ALLOW_RING}
     3746,{jjPFAC2,     PFAC_CMD,       LIST_CMD,       BIGINT_CMD, NUMBER_CMD, ALLOW_PLURAL |ALLOW_RING}
     3747,{jjPFAC2,     PFAC_CMD,       LIST_CMD,       NUMBER_CMD, NUMBER_CMD, ALLOW_PLURAL |ALLOW_RING}
    37253748#ifdef HAVE_PLURAL
    37263749,{jjPlural_num_poly, NCALGEBRA_CMD,NONE,       POLY_CMD,   POLY_CMD  , NO_PLURAL |NO_RING}
     
    45144537static BOOLEAN jjPFAC1(leftv res, leftv v)
    45154538{
     4539  /* call method jjPFAC2 with second argument = 0 (meaning that no
     4540     valid bound for the prime factors has been given) */
    45164541  sleftv tmp;
    4517   memset(&tmp,0,sizeof(tmp));
    4518   tmp.rtyp=INT_CMD;
    4519   return jjPFAC2(res,v,&tmp);
     4542  memset(&tmp, 0, sizeof(tmp));
     4543  tmp.rtyp = INT_CMD;
     4544  return jjPFAC2(res, v, &tmp);
    45204545}
    45214546static BOOLEAN jjLU_DECOMP(leftv res, leftv v)
  • Singular/misc_ip.cc

    rf9bb68 r555b03  
    2222#include <Singular/misc_ip.h>
    2323
    24 /* This works by Newton iteration, i.e.,
    25       a(1)   = n;
    26       a(i+1) = a(i)/2 + n/2/a(i), i > 0.
    27    This sequence is guaranteed to decrease monotonously and
    28    it is known to converge fast.
    29    All used numbers are bigints. */
    30 number approximateSqrt(const number n)
    31 {
    32   if (nlIsZero(n)) { number zero = nlInit(0, NULL); return zero; }
    33   number temp1; number temp2;
    34   number one = nlInit(1, NULL);
    35   number two = nlInit(2, NULL);
    36   number m = nlCopy(n);
    37   number mOld = nlSub(m, one); /* initially required to be different from m */
    38   number nHalf = nlIntDiv(n, two);
    39   bool check = true;
    40   while (!nlEqual(m, mOld) && check)
    41   {
    42     temp1 = nlIntDiv(m, two);
    43     temp2 = nlIntDiv(nHalf, m);
    44     mOld = m;
    45     m = nlAdd(temp1, temp2);
    46     nlDelete(&temp1, NULL); nlDelete(&temp2, NULL);
    47     temp1 = nlMult(m, m);
    48     check = nlGreater(temp1, n);
    49     nlDelete(&temp1, NULL);
    50   }
    51   nlDelete(&mOld, NULL); nlDelete(&two, NULL); nlDelete(&nHalf, NULL);
    52   while (!check)
    53   {
    54     temp1 = nlAdd(m, one);
    55     nlDelete(&m, NULL);
    56     m = temp1;
    57     temp1 = nlMult(m, m);
    58     check = nlGreater(temp1, n);
    59     nlDelete(&temp1, NULL);
    60   }
    61   temp1 = nlSub(m, one);
    62   nlDelete(&m, NULL);
    63   nlDelete(&one, NULL);
    64   m = temp1;
    65   return m;
    66 }
    67 
    68 /* returns the quotient resulting from division of n by the prime as many
    69    times as possible without remainder; afterwards, the parameter times
    70    will contain the highest exponent e of p such that p^e divides n
    71    e.g., divTimes(48, 4, t) = 3 with t = 2, since 48 = 4*4*3;
    72    n is expected to be a bigint; returned type is also bigint */
    73 number divTimes(const number n, const int p, int* times)
    74 {
    75   number nn = nlCopy(n);
    76   number dd = nlInit(p, NULL);
    77   number rr = nlIntMod(nn, dd);
     24void number2mpz(number n, mpz_t m)
     25{
     26  if (nlIsZero(n)) mpz_init_set_si(m, 0);
     27  else
     28  {
     29    int nAsInt = nlInt(n, NULL);
     30    if (nAsInt != 0) mpz_init_set_si(m, nAsInt);     /* n fits in an int */
     31    else             mpz_init_set(m, (mpz_ptr)n->z);
     32  }
     33}
     34
     35number mpz2number(mpz_t m)
     36{
     37  number z = (number)omAllocBin(rnumber_bin);
     38  mpz_init_set(z->z, m);
     39  mpz_init_set_ui(z->n, 1);
     40  z->s = 3;
     41  return z;
     42}
     43
     44void divTimes(mpz_t n, mpz_t d, int* times)
     45{
    7846  *times = 0;
    79   while (nlIsZero(rr))
     47  mpz_t r; mpz_init(r);
     48  mpz_t q; mpz_init(q);
     49  mpz_fdiv_qr(q, r, n, d);
     50  while (mpz_cmp_ui(r, 0) == 0)
    8051  {
    8152    (*times)++;
    82     number temp = nlIntDiv(nn, dd);
    83     nlDelete(&nn, NULL);
    84     nn = temp;
    85     nlDelete(&rr, NULL);
    86     rr = nlIntMod(nn, dd);
    87   }
    88   nlDelete(&rr, NULL); nlDelete(&dd, NULL);
    89   return nn;
     53    mpz_set(n, q);
     54    mpz_fdiv_qr(q, r, n, d);
     55  }
     56  mpz_clear(r);
     57  mpz_clear(q);
    9058}
    9159
    9260/* returns an object of type lists which contains the entries
    93    theInts[0..(length-1)] as INT_CMDs*/
     61   theInts[0..(length-1)] as INT_CMDs */
    9462lists makeListsObject(const int* theInts, int length)
    9563{
     
    10169}
    10270
    103 /* returns the i-th bit of the binary number which arises by
    104    concatenating array[length-1], ..., array[1], array[0],
    105    where array[0] contains the 32 lowest bits etc.;
    106    i is assumed to be small enough to address a valid index
    107    in the given array */
    108 bool getValue(const unsigned i, const unsigned int* ii)
    109 {
    110   if (i==2) return true;
    111   if ((i & 1)==0) return false;
    112   unsigned I= i/2;
    113   unsigned index = I / 32;
    114   unsigned offset = I % 32;
    115   unsigned int v = 1 << offset;
    116   return ((ii[index] & v) != 0);
    117 }
    118 
    119 /* sets the i-th bit of the binary number which arises by
    120    concatenating array[length-1], ..., array[1], array[0],
    121    where array[0] contains the 32 lowest bits etc.;
    122    i is assumed to be small enough to address a valid index
    123    in the given array */
    124 void setValue(const unsigned i, bool value, unsigned int* ii)
    125 {
    126   if ((i&1)==0) return; // ignore odd numbers
    127   unsigned I=i/2;
    128   unsigned index = I / 32;
    129   unsigned offset = I % 32;
    130   unsigned int v = 1 << offset;
    131   if (value) { ii[index] |= v;  }
    132   else       { ii[index] &= ~v; }
    133 }
    134 
    135 /* returns whether i is less than or equal to the bigint number n */
    136 bool isLeq(const int i, const number n)
    137 {
    138   number iN = nlInit(i - 1, NULL);
    139   bool result = nlGreater(n, iN);
    140   nlDelete(&iN, NULL);
    141   return result;
    142 }
    143 
    144 #if 0
    145 lists primeFactorisation(const number n, const int pBound)
    146 {
    147   number nn = nlCopy(n); int i;
    148   int pCounter = 0; /* for counting the number of mutually distinct
    149                        prime factors in n */
    150   /* we assume that there are at most 1000 mutually distinct prime
    151      factors in n */
    152   int* primes = new int[1000]; int* multiplicities = new int[1000];
    153 
    154   /* extra treatment for the primes 2 and 3;
    155      all other primes are equal to +1/-1 mod 6 */
    156   int e; number temp;
    157   temp = divTimes(nn, 2, &e); nlDelete(&nn, NULL); nn = temp;
    158   if (e > 0) { primes[pCounter] = 2; multiplicities[pCounter++] = e; }
    159   temp = divTimes(nn, 3, &e); nlDelete(&nn, NULL); nn = temp;
    160   if (e > 0) { primes[pCounter] = 3; multiplicities[pCounter++] = e; }
    161 
    162   /* now we simultaneously:
    163      - build the sieve of Erathostenes up to s,
    164      - divide out each prime factor of nn that we find along the way
    165        (This may result in an earlier termination.) */
    166 
    167   int s = 1<<25;       /* = 2^25 */
    168   int maxP = 2147483647; /* = 2^31 - 1, by the way a Mersenne prime */
    169   if ((pBound != 0) && (pBound < maxP))
    170   {
    171     maxP = pBound;
    172   }
    173   if (maxP< (2147483647-63) ) s=(maxP+63)/64;
    174   else                        s=2147483647/64+1;
    175   unsigned int* isPrime = new unsigned int[s];
    176   /* the lowest bit of isPrime[0] stores whether 1 is a prime,
    177      next bit is for 3, next for 5, etc. i.e.
    178      intended usage is: isPrime[0] = ...
    179      We shall make use only of bits which correspond to numbers =
    180      -1 or +1 mod 6. */
    181   //for (i = 0; i < s; i++) isPrime[i] = ~0;/*4294967295*/; /* all 32 bits set */
    182   memset(isPrime,0xff,s*sizeof(unsigned int));
    183   int p=9; while((p<maxP) && (p>0)) { setValue(p,false,isPrime); p+=6; }
    184   p = 5; bool add2 = true;
    185   /* due to possible overflows, we need to check whether p > 0, and
    186      likewise i > 0 below */
    187   while ((0 < p) && (p <= maxP) && (isLeq(p, nn)))
    188   {
    189     /* at this point, p is guaranteed to be a prime;
    190        we divide nn by the highest power of p and store p
    191        if nn is at all divisible by p */
    192     temp = divTimes(nn, p, &e);
    193     nlDelete(&nn, NULL); nn = temp;
    194     if (e > 0)
    195     { primes[pCounter] = p; multiplicities[pCounter++] = e; }
    196     /* invalidate all multiples of p, starting with 2*p */
    197     i = 2 * p;
    198     while ((0 < i) && (i < maxP)) { setValue(i, false, isPrime); i += p; }
    199     /* move on to the next prime in the sieve; we either add 2 or 4
    200        in order to visit just the numbers equal to -1/+1 mod 6 */
    201     if (add2) { p += 2; add2 = false; }
    202     else      { p += 4; add2 = true;  }
    203     while ((0 < p) && (p <= maxP) && (isLeq(p, nn)) && (!getValue(p, isPrime)))
    204     {
    205       if (add2) { p += 2; add2 = false; }
    206       else      { p += 4; add2 = true;  }
    207     }
    208   }
    209 
    210   /* build return structure and clean up */
    211   delete [] isPrime;
    212   lists primesL = makeListsObject(primes, pCounter);
    213   lists multiplicitiesL = makeListsObject(multiplicities, pCounter);
    214   delete [] primes; delete [] multiplicities;
    215   lists L=(lists)omAllocBin(slists_bin);
    216   L->Init(3);
    217   L->m[0].rtyp = BIGINT_CMD; L->m[0].data = (void *)nn;
     71void setListEntry(lists L, int index, mpz_t n)
     72{ /* assumes n > 0 */
     73  number nn = mpz2number(n);
     74  L->m[index].rtyp = BIGINT_CMD; L->m[index].data = (void*)nn;
    21875  /* try to fit nn into an int: */
    21976  int nnAsInt = nlInt(nn, NULL);
    220   if (nlIsZero(nn) || (nnAsInt != 0))
    221   {
    222     nlDelete(&nn,NULL);
    223     L->m[0].rtyp = INT_CMD;
    224     L->m[0].data = (void *)nnAsInt;
    225   }
    226   L->m[1].rtyp = LIST_CMD; L->m[1].data = (void *)primesL;
    227   L->m[2].rtyp = LIST_CMD; L->m[2].data = (void *)multiplicitiesL;
     77  if (nnAsInt != 0)
     78  {
     79    nlDelete(&nn, NULL);
     80    L->m[index].rtyp = INT_CMD; L->m[index].data = (void*)nnAsInt;
     81  }
     82}
     83
     84/* true iff p is prime */
     85/*
     86bool isPrime(mpz_t p)
     87{
     88  if (mpz_cmp_ui(p, 2) == 0) return true;
     89  if (mpz_cmp_ui(p, 3) == 0) return true;
     90  if (mpz_cmp_ui(p, 5) < 0)  return false;
     91
     92  mpz_t d; mpz_init_set_ui(d, 5); int add = 2;
     93  mpz_t sr; mpz_init(sr); mpz_sqrt(sr, p);
     94  mpz_t r; mpz_init(r);
     95  while (mpz_cmp(d, sr) <= 0)
     96  {
     97    mpz_cdiv_r(r, p, d);
     98    if (mpz_cmp_ui(r, 0) == 0)
     99    {
     100      mpz_clear(d); mpz_clear(sr); mpz_clear(r);
     101      return false;
     102    }
     103    mpz_add_ui(d, d, add);
     104    add += 2; if (add == 6) add = 2;
     105  }
     106  mpz_clear(d); mpz_clear(sr); mpz_clear(r);
     107  return true;
     108}
     109*/
     110
     111/* finds the next prime q, bound >= q >= p;
     112   in case of success, puts q into p;
     113   otherwise sets q = bound + 1;
     114   e.g. p = 24; nextPrime(p, 30) produces p = 29 (success),
     115        p = 24; nextPrime(p, 29) produces p = 29 (success),
     116        p = 24; nextPrime(p, 28) produces p = 29 (no success),
     117        p = 24; nextPrime(p, 27) produces p = 28 (no success) */
     118/*
     119void nextPrime(mpz_t p, mpz_t bound)
     120{
     121  int add;
     122  mpz_t r; mpz_init(r); mpz_cdiv_r_ui(r, p, 6); // r = p mod 6, 0 <= r <= 5
     123  if (mpz_cmp_ui(r, 0) == 0) { mpz_add_ui(p, p, 1); add = 4; }
     124  if (mpz_cmp_ui(r, 1) == 0) {                      add = 4; }
     125  if (mpz_cmp_ui(r, 2) == 0) { mpz_add_ui(p, p, 3); add = 2; }
     126  if (mpz_cmp_ui(r, 3) == 0) { mpz_add_ui(p, p, 2); add = 2; }
     127  if (mpz_cmp_ui(r, 4) == 0) { mpz_add_ui(p, p, 1); add = 2; }
     128  if (mpz_cmp_ui(r, 5) == 0) {                      add = 2; }
     129
     130  while (mpz_cmp(p, bound) <= 0)
     131  {
     132    if (isPrime(p)) { mpz_clear(r); return; }
     133    mpz_add_ui(p, p, add);
     134    add += 2; if (add == 6) add = 2;
     135  }
     136  mpz_set(p, bound);
     137  mpz_add_ui(p, p, 1);
     138  mpz_clear(r);
     139  return;
     140}
     141*/
     142
     143/* n and pBound are assumed to be bigint numbers */
     144lists primeFactorisation(const number n, const number pBound)
     145{
     146  mpz_t nn; number2mpz(n, nn);
     147  mpz_t pb; number2mpz(pBound, pb);
     148  mpz_t b; number2mpz(pBound, b);
     149  mpz_t p; mpz_init(p); int tt;
     150  mpz_t sr; mpz_init(sr); int index = 0; int add;
     151  lists primes = (lists)omAllocBin(slists_bin); primes->Init(1000);
     152  int* multiplicities = new int[1000];
     153
     154  mpz_set_ui(p, 2);
     155  if ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(p, b) <= 0))
     156  {
     157    divTimes(nn, p, &tt);
     158    if (tt > 0)
     159    {
     160      setListEntry(primes, index, p);
     161      multiplicities[index++] = tt;
     162    }
     163  }
     164
     165  mpz_set_ui(p, 3);
     166  if ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(p, b) <= 0))
     167  {
     168    divTimes(nn, p, &tt);
     169    if (tt > 0)
     170    {
     171      setListEntry(primes, index, p);
     172      multiplicities[index++] = tt;
     173    }
     174  }
     175
     176  mpz_set_ui(p, 5); add = 2;
     177  mpz_sqrt(sr, nn);
     178  if ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
     179  while (mpz_cmp(pb, p) >= 0)
     180  {
     181    divTimes(nn, p, &tt);
     182    if (tt > 0)
     183    {
     184      setListEntry(primes, index, p);
     185      multiplicities[index++] = tt;
     186      mpz_sqrt(sr, nn);
     187      if ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(pb, sr) > 0)) mpz_set(pb, sr);
     188    }
     189    mpz_add_ui(p, p, add);
     190    add += 2; if (add == 6) add = 2;
     191  }
     192  if ((mpz_cmp_ui(nn, 1) > 0) &&
     193      ((mpz_cmp_ui(b, 0) == 0) || (mpz_cmp(nn, b) <= 0)))
     194  {
     195    setListEntry(primes, index, nn);
     196    multiplicities[index++] = 1;
     197    mpz_set_ui(nn, 1);
     198  }
     199
     200  lists primesL = (lists)omAllocBin(slists_bin);
     201  primesL->Init(index);
     202  for (int i = 0; i < index; i++)
     203  {
     204    primesL->m[i].rtyp = primes->m[i].rtyp;
     205    primesL->m[i].data = primes->m[i].data;
     206  }
     207  omFreeSize((ADDRESS)primes->m, (primes->nr + 1) * sizeof(sleftv));
     208  omFreeBin((ADDRESS)primes, slists_bin);
     209
     210  lists multiplicitiesL = (lists)omAllocBin(slists_bin);
     211  multiplicitiesL->Init(index);
     212  for (int i = 0; i < index; i++)
     213  {
     214    multiplicitiesL->m[i].rtyp = INT_CMD;
     215    multiplicitiesL->m[i].data = (void*)multiplicities[i];
     216  }
     217  delete[] multiplicities;
     218
     219  lists L=(lists)omAllocBin(slists_bin);
     220  L->Init(4);
     221  setListEntry(L, 0, nn);
     222  L->m[1].rtyp = LIST_CMD; L->m[1].data = (void*)primesL;
     223  L->m[2].rtyp = LIST_CMD; L->m[2].data = (void*)multiplicitiesL;
     224  int probTest = 0;
     225  if (mpz_probab_prime_p(nn, 25) != 0) probTest = 1;
     226  L->m[3].rtyp =  INT_CMD; L->m[3].data = (void*)probTest;
     227  mpz_clear(nn); mpz_clear(pb); mpz_clear(b); mpz_clear(p); mpz_clear(sr);
     228
    228229  return L;
    229230}
    230 #else
    231 /* Factoring , from gmp-demos
    232 
    233 Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005 Free Software
    234 Foundation, Inc.
    235 
    236 This program is free software; you can redistribute it and/or modify it under
    237 the terms of the GNU General Public License as published by the Free Software
    238 Foundation; either version 3 of the License, or (at your option) any later
    239 version.
    240 
    241 This program is distributed in the hope that it will be useful, but WITHOUT ANY
    242 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
    243 PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    244 
    245 You should have received a copy of the GNU General Public License along with
    246 this program.  If not, see http://www.gnu.org/licenses/.  */
    247 
    248 
    249 #include <stdlib.h>
    250 #include <stdio.h>
    251 #include <string.h>
    252 
    253 #include <kernel/si_gmp.h>
    254 
    255 static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
    256 
    257 void factor_using_division (mpz_t t, unsigned int limit, int *L, int &L_ind, int *ex)
    258 {
    259   mpz_t q, r;
    260   unsigned long int f;
    261   int ai;
    262   unsigned *addv = add;
    263   unsigned int failures;
    264 
    265   mpz_init (q);
    266   mpz_init (r);
    267 
    268   f = mpz_scan1 (t, 0);
    269   mpz_div_2exp (t, t, f);
    270   while (f)
    271   {
    272     if ((L_ind>0) && (L[L_ind-1]==2)) ex[L_ind-1]++;
    273     else
    274     {
    275       L[L_ind]=2;
    276       L_ind++;
    277     }
    278     --f;
    279   }
    280 
    281   for (;;)
    282   {
    283     mpz_tdiv_qr_ui (q, r, t, 3);
    284     if (mpz_cmp_ui (r, 0) != 0)  break;
    285     mpz_set (t, q);
    286     if ((L_ind>0) && (L[L_ind-1]==3)) ex[L_ind-1]++;
    287     else
    288     {
    289       L[L_ind]=3;
    290       L_ind++;
    291     }
    292   }
    293 
    294   for (;;)
    295   {
    296     mpz_tdiv_qr_ui (q, r, t, 5);
    297     if (mpz_cmp_ui (r, 0) != 0) break;
    298     mpz_set (t, q);
    299     if ((L_ind>0) && (L[L_ind-1]==5)) ex[L_ind-1]++;
    300     else
    301     {
    302       L[L_ind]=5;
    303       L_ind++;
    304     }
    305   }
    306 
    307   failures = 0;
    308   f = 7;
    309   ai = 0;
    310   while (mpz_cmp_ui (t, 1) != 0)
    311   {
    312     if (f>=((unsigned long)1 <<28)) break;
    313     if (mpz_cmp_ui (t, f) < 0) break;
    314     mpz_tdiv_qr_ui (q, r, t, f);
    315     if (mpz_cmp_ui (r, 0) != 0)
    316     {
    317         f += addv[ai];
    318         ai = (ai + 1) & 7;
    319         failures++;
    320         if (failures > limit) break;
    321     }
    322     else
    323     {
    324       mpz_swap (t, q);
    325       //gmp_printf("%d: %Zd\n",f,t);
    326       // here: f in 0,,2^28-1:
    327       if ((L_ind>0) && (L[L_ind-1]==(int)f)) ex[L_ind-1]++;
    328       else
    329       {
    330         L[L_ind]=f;
    331         L_ind++;
    332       }
    333       failures = 0;
    334     }
    335   }
    336 
    337   mpz_clear (q);
    338   mpz_clear (r);
    339 }
    340 
    341 void factor_using_pollard_rho (mpz_t n, int a_int, int *L, int &L_ind, int *ex)
    342 {
    343   mpz_t x, x1, y, P;
    344   mpz_t a;
    345   mpz_t g;
    346   mpz_t t1, t2;
    347   int k, l, c, i;
    348 
    349   mpz_init (g);
    350   mpz_init (t1);
    351   mpz_init (t2);
    352 
    353   mpz_init_set_si (a, a_int);
    354   mpz_init_set_si (y, 2);
    355   mpz_init_set_si (x, 2);
    356   mpz_init_set_si (x1, 2);
    357   k = 1;
    358   l = 1;
    359   mpz_init_set_ui (P, 1);
    360   c = 0;
    361 
    362   while (mpz_cmp_ui (n, 1) != 0)
    363   {
    364 S2:
    365     mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
    366     mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
    367     c++;
    368     if (c == 20)
    369     {
    370       c = 0;
    371       mpz_gcd (g, P, n);
    372       if (mpz_cmp_ui (g, 1) != 0) goto S4;
    373       mpz_set (y, x);
    374     }
    375 S3:
    376     k--;
    377     if (k > 0) goto S2;
    378 
    379     mpz_gcd (g, P, n);
    380     if (mpz_cmp_ui (g, 1) != 0) goto S4;
    381 
    382     mpz_set (x1, x);
    383     k = l;
    384     l = 2 * l;
    385     for (i = 0; i < k; i++)
    386     {
    387       mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
    388     }
    389     mpz_set (y, x);
    390     c = 0;
    391     goto S2;
    392 S4:
    393     do
    394     {
    395       mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
    396       mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
    397     }
    398     while (mpz_cmp_ui (g, 1) == 0);
    399 
    400     mpz_div (n, n, g);        /* divide by g, before g is overwritten */
    401 
    402     if (!mpz_probab_prime_p (g, 10))
    403     {
    404       do
    405       {
    406         mp_limb_t a_limb;
    407         mpn_random (&a_limb, (mp_size_t) 1);
    408         a_int = (int) a_limb;
    409       }
    410       while (a_int == -2 || a_int == 0);
    411 
    412       factor_using_pollard_rho (g, a_int,L,L_ind,ex);
    413     }
    414     else
    415     {
    416       if ((L_ind>0) && (mpz_cmp_si(g,L[L_ind-1])==0)) ex[L_ind-1]++;
    417       else
    418       {
    419         L[L_ind]=mpz_get_si(g);
    420         L_ind++;
    421       }
    422     }
    423     mpz_mod (x, x, n);
    424     mpz_mod (x1, x1, n);
    425     mpz_mod (y, y, n);
    426     if (mpz_probab_prime_p (n, 10))
    427     {
    428       int te=mpz_get_si(n);
    429       if (mpz_cmp_si(n,te)==0) /* does it fit into an int ? */
    430       {
    431         if ((L_ind>0) && (mpz_cmp_si(n,L[L_ind-1])==0)) ex[L_ind-1]++;
    432         else
    433         {
    434           L[L_ind]=mpz_get_si(n);
    435           L_ind++;
    436         }
    437         mpz_set_si(n,1); // add n itself the list of divisors, rest is 1
    438       }
    439       break;
    440     }
    441   }
    442 
    443   mpz_clear (g);
    444   mpz_clear (P);
    445   mpz_clear (t2);
    446   mpz_clear (t1);
    447   mpz_clear (a);
    448   mpz_clear (x1);
    449   mpz_clear (x);
    450   mpz_clear (y);
    451 }
    452 
    453 void mpz_factor (mpz_t t, int *L, int & L_ind, int *ex)
    454 {
    455   unsigned int division_limit;
    456 
    457   if (mpz_sgn (t) == 0)
    458     return;
    459 
    460   /* Set the trial division limit according the size of t.  */
    461   division_limit = mpz_sizeinbase (t, 2);
    462   if (division_limit > 1000)
    463     division_limit = 1000 * 1000;
    464   else
    465     division_limit = division_limit * division_limit;
    466 
    467   factor_using_division (t, division_limit, L, L_ind, ex);
    468 
    469   if (mpz_cmp_ui (t, 1) != 0)
    470   {
    471     if (mpz_probab_prime_p (t, 10))
    472     {
    473       int tt=mpz_get_si(t);
    474       // check if t fits into int:
    475       if ((mpz_size1(t)==1)&&(mpz_cmp_si(t,tt)==0))
    476       {
    477         L[L_ind]=mpz_get_si(t);
    478         L_ind++;
    479         mpz_set_si(t,1);
    480       }
    481     }
    482     else
    483       factor_using_pollard_rho (t, 1, L,L_ind,ex);
    484   }
    485 }
    486 
    487 lists primeFactorisation(const number n, const int pBound)
    488 {
    489   mpz_t t;
    490   number nn = nlCopy(n);
    491   lists L=(lists)omAllocBin(slists_bin);
    492   L->Init(3);
    493   L->m[0].rtyp = BIGINT_CMD; L->m[0].data = (void *)nn;
    494   /* try to fit nn into an int: */
    495   int nnAsInt = nlInt(nn, NULL);
    496   if (nlIsZero(nn) || (nnAsInt != 0))
    497   {
    498     mpz_init_set_si(t,nnAsInt);
    499   }
    500   else
    501   {
    502     mpz_init_set(t,(mpz_ptr)nn->z);
    503   }
    504   int *LL=(int*)omAlloc0(1000*sizeof(int));
    505   int *ex=(int*)omAlloc0(1000*sizeof(int));
    506   int L_ind=0;
    507   mpz_factor (t,LL,L_ind,ex);
    508 
    509   nnAsInt = mpz_get_si(t);
    510   if ((mpz_size1(t)==1) && (mpz_cmp_si(t,nnAsInt)==0))
    511   {
    512     nlDelete(&nn,NULL);
    513     L->m[0].rtyp = INT_CMD;
    514     L->m[0].data = (void *)nnAsInt;
    515   }
    516   else
    517   {
    518     mpz_set(nn->z,t);
    519     L->m[0].rtyp = BIGINT_CMD;
    520     L->m[0].data = (void *)nn;
    521   }
    522   mpz_clear(t);
    523   int i;
    524   for(i=0;i<L_ind;i++) ex[i]++;
    525   L->m[1].rtyp = LIST_CMD; L->m[1].data = (void *)makeListsObject(LL,L_ind);
    526   L->m[2].rtyp = LIST_CMD; L->m[2].data = (void *)makeListsObject(ex,L_ind);
    527   return L;
    528 }
    529 #endif
    530231
    531232#include <string.h>
  • Singular/misc_ip.h

    rf9bb68 r555b03  
    2929
    3030/**
    31  * Computes an approximation of the square root of a bigint number.
     31 * Converts a non-negative bigint number into a GMP number.
    3232 *
    33  * Expects a non-negative bigint number n and returns the bigint
    34  * number m which satisfies m*m <= n < (m+1)*(m+1). This implementation
    35  * does not require a valid currRing.
     33 **/
     34void number2mpz(number n,   /**< [in]  a bigint number >= 0  */
     35                mpz_t m     /**< [out] the GMP equivalent    */
     36               );
     37               
     38/**
     39 * Converts a non-negative GMP number into a bigint number.
    3640 *
    37  * @return an approximation of the square root of n
     41 * @return the bigint number representing the given GMP number
    3842 **/
    39 number approximateSqrt(const number n   /**< [in]  a number   */
    40                       );
     43number mpz2number(mpz_t m   /**< [in]  a GMP number >= 0  */
     44                 );
     45
     46/**
     47 * Divides 'n' as many times as possible by 'd' and returns the number
     48 * of divisions (without remainder) in 'times',
     49 * e.g., n = 48, d = 4, divTimes(n, d, t) = 3 produces n = 3, t = 2,
     50 *       since 48 = 4*4*3;
     51 * assumes that d is positive
     52 **/
     53void divTimes(mpz_t n,   /**< [in]  a GMP number >= 0                      */
     54              mpz_t d,   /**< [in]  the divisor, a GMP number >= 0         */
     55              int* times /**< [out] number of divisions without remainder  */
     56             );
    4157
    4258/**
     
    4460 * than or equal to a given bound, with corresponding multiplicities.
    4561 *
    46  * The method is capable of finding only prime factors < 2^31, and it looks
    47  * only for those less than or equal the given bound, if this bound is not 0.
    48  * Therefore, after dividing out all prime factors which have been found,
    49  * some number m > 1 may survive.
    50  * The method returns a list L filled with three entries:
    51  * L[1] contains m as int or bigint, depending on the size,
    52  * L[2] a list; L[2][i] contains the i-th prime factor as int
     62 * The method finds all prime factors with multiplicities. If a non-zero
     63 * bound is given, then only the prime factors <= pBound are being found.
     64 * In this case, there may remain an unfactored portion m of n.
     65 * The method returns a list L filled with four entries:
     66 * L[1] contains the remainder m as int or bigint, depending on the size,
     67 * L[2] a list; L[2][i] contains the i-th prime factor as int or bigint
    5368 *                     (sorted in ascending order),
    5469 * L[3] a list; L[3][i] contains the multiplicity of L[2, i] in n as int
     70 * L[4] 1 iff the remainder m is probably a prime, 0 otherwise
    5571 *
    5672 * We thus have: n = L[1] * L[2][1]^L[3][1] * ... * L[2][k]^L[3][k], where
    57  * k is the number of mutually distinct prime factors < 2^31 and <= the
    58  * provided bound (if this is not zero).
     73 * k is the number of mutually distinct prime factors (<= a provided non-
     74 * zero bound).
    5975 *
    60  * @return the factorisation data in a list
     76 * @return the factorisation data in a SINGULAR-internal list
    6177 **/
    6278lists primeFactorisation(
    63        const number n,    /**< [in]  the bigint to be factorised       */
    64        const int pBound   /**< [in]  bound on the prime factors seeked */
     79       const number n,     /**< [in]  the bigint > 0 to be factorised   */
     80       const number pBound /**< [in]  bigint bound on the prime factors
     81                                      seeked                            */
    6582                        );
    6683
Note: See TracChangeset for help on using the changeset viewer.