Changeset 4c4340 in git
 Timestamp:
 Oct 2, 2012, 12:25:11 PM (10 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
 Children:
 65b813a15472ad024e1b30ff11d4f966ec98b53a
 Parents:
 8b57f5e1c589185d90380ede5c9b3a0b2a44ff97
 gitauthor:
 Hans Schoenemann <hannes@mathematik.unikl.de>20121002 12:25:11+02:00
 gitcommitter:
 Martin Lee <martinlee84@web.de>20130123 11:27:17+01:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/misc_ip.cc
r8b57f5 r4c4340 37 37 38 38 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 84 39 void setListEntry(lists L, int index, mpz_t n) 85 40 { /* assumes n > 0 */ … … 114 69 } 115 70 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 */ 73 static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6}; 74 75 void 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[index1]++; 152 } 153 last_f=f; 154 failures = 0; 155 } 156 } 157 158 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) 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[index1]++; 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[index1]++; 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 301 void 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 } 177 334 /* n and pBound are assumed to be bigint numbers */ 178 335 lists primeFactorisation(const number n, const number pBound) 179 336 { 337 int i; 338 int index=0; 180 339 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;185 340 lists primes = (lists)omAllocBin(slists_bin); primes>Init(1000); 186 int* multiplicities = new int[1000];341 int* multiplicities = (int*)omAlloc0(1000*sizeof(int)); 187 342 int positive=1; int probTest = 0; 188 343 … … 194 349 mpz_neg(nn,nn); 195 350 } 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); 300 352 } 301 353 302 354 lists primesL = (lists)omAllocBin(slists_bin); 303 355 primesL>Init(index); 304 for (i nt i= 0; i < index; i++)356 for (i = 0; i < index; i++) 305 357 { 306 358 primesL>m[i].rtyp = primes>m[i].rtyp; … … 312 364 lists multiplicitiesL = (lists)omAllocBin(slists_bin); 313 365 multiplicitiesL>Init(index); 314 for (i nt i= 0; i < index; i++)366 for (i = 0; i < index; i++) 315 367 { 316 368 multiplicitiesL>m[i].rtyp = INT_CMD; … … 326 378 setListEntry(L, 2, nn); 327 379 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); 329 381 330 382 return L;
Note: See TracChangeset
for help on using the changeset viewer.