001/* 002 * $Id$ 003 */ 004 005package edu.jas.arith; 006 007 008import java.util.ArrayList; 009import java.util.BitSet; 010import java.util.Collections; 011import java.util.List; 012import java.util.Map; 013import java.util.Random; 014import java.util.SortedMap; 015import java.util.SortedSet; 016import java.util.TreeMap; 017import java.util.TreeSet; 018 019import org.apache.logging.log4j.Logger; 020import org.apache.logging.log4j.LogManager; 021 022 023/** 024 * Integer prime factorization. Code from ALDES/SAC2 and MAS module SACPRIM. 025 * 026 * See ALDES/SAC2 or MAS code in SACPRIM. 027 * See Symja <code>org/matheclipse/core/expression/Primality.java</code> for Pollard 028 * algorithm. 029 * @author Heinz Kredel 030 */ 031 032public final class PrimeInteger { 033 034 035 private static final Logger logger = LogManager.getLogger(PrimeInteger.class); 036 037 038 /** 039 * Maximal long, which can be factored by IFACT(long). Has nothing to do 040 * with SAC2.BETA. 041 */ 042 final public static long BETA = PrimeList.getLongPrime(61, 1).longValue(); 043 044 045 /** 046 * Medium prime divisor range. 047 */ 048 //final static long IMPDS_MIN = 1000; // SAC2/Aldes 049 //final static long IMPDS_MAX = 5000; // " 050 //final static long IMPDS_MIN = 2000; 051 //final static long IMPDS_MAX = 10000; 052 final static long IMPDS_MIN = 10000; 053 054 055 final static long IMPDS_MAX = 128000; 056 057 058 /** 059 * List of small prime numbers. 060 */ 061 final public static List<Long> SMPRM = smallPrimes(2, (int) (IMPDS_MIN >> 1)); 062 063 064 /** 065 * List of units of Z mod 210. 066 */ 067 final public static List<Long> UZ210 = getUZ210(); 068 069 070 /** 071 * Digit prime generator. K and m are positive beta-integers. L is the list 072 * (p(1),...,p(r)) of all prime numbers p such that m le p lt m+2*K, with 073 * p(1) lt p(2) lt ... lt p(r). 074 * See also SACPRIM.DPGEN. 075 * @param m start integer 076 * @param K number of integers 077 * @return the list L of prime numbers p with m ≤ p < m + 2*K. 078 */ 079 public static List<Long> smallPrimes(long m, int K) { 080 int k; 081 long ms; 082 ms = m; 083 if (ms <= 1) { 084 ms = 1; 085 } 086 m = ms; 087 if (m % 2 == 0) { 088 m++; 089 K--; 090 } 091 //if (kp % 2 == 0) { 092 // k = kp/2; 093 //} else { 094 // k = (kp+1)/2; 095 //} 096 k = K; 097 098 /* init */ 099 long h = 2L * ((long)k - 1L); 100 long m2 = m + h; // mp 101 BitSet p = new BitSet(k); 102 p.set(0, k); 103 //for (int i = 0; i < k; i++) { 104 // p.set(i); 105 //} 106 107 /* compute */ 108 int r, d = 0; 109 int i, c = 0; 110 while (true) { 111 switch (c) { 112 /* mark multiples of d for d=3 and d=6n-/+1 with d**2<=m2 */ 113 case 2: 114 d += 2; 115 c = 3; 116 break; 117 case 3: 118 d += 4; 119 c = 2; 120 break; 121 case 0: 122 d = 3; 123 c = 1; 124 break; 125 case 1: 126 d = 5; 127 c = 2; 128 break; 129 default: 130 throw new RuntimeException("this should not happen"); 131 } 132 if (d > (m2 / d)) { 133 break; 134 } 135 r = (int) (m % d); 136 if (r + h >= d || r == 0) { 137 if (r == 0) { 138 i = 0; 139 } else { 140 if (r % 2 == 0) { 141 i = d - (r / 2); 142 } else { 143 i = (d - r) / 2; 144 } 145 } 146 if (m <= d) { 147 i += d; 148 } 149 while (i < k) { 150 p.set(i, false); 151 i += d; 152 } 153 } 154 } 155 /* output */ 156 int l = p.cardinality(); // l = 0 157 //for (i=0; i<k; i++) { 158 // if (p.get(i)) { 159 // l++; 160 // } 161 //} 162 if (ms <= 2) { 163 l++; 164 } 165 //if (ms <= 1) { 166 //} 167 List<Long> po = new ArrayList<Long>(l); 168 if (l == 0) { 169 return po; 170 } 171 //l = 0; 172 if (ms == 1) { 173 //po.add(2); 174 //l++; 175 p.set(0, false); 176 } 177 if (ms <= 2) { 178 po.add(2L); 179 //l++; 180 } 181 long pl = m; 182 //System.out.println("pl = " + pl + " p[0] = " + p[0]); 183 //System.out.println("k-1 = " + (k-1) + " p[k-1] = " + p[k-1]); 184 for (i = 0; i < k; i++) { 185 if (p.get(i)) { 186 po.add(pl); 187 //l++; 188 } 189 pl += 2; 190 } 191 //System.out.println("SMPRM = " + po); 192 return po; 193 } 194 195 196 /** 197 * Integer small prime divisors. n is a positive integer. F is a list of 198 * primes (q(1),q(2),...,q(h)), h non-negative, q(1) le q(2) le ... lt q(h), 199 * such that n is equal to m times the product of the q(i) and m is not 200 * divisible by any prime in SMPRM. Either m=1 or m gt 1,000,000. 201 * <br /> In JAS F is a map and m=1 or m > 4.000.000. 202 * See also SACPRIM.ISPD. 203 * @param n integer to factor. 204 * @param F a map of pairs of prime numbers and multiplicities (p,e) with p**e 205 * divides n and e maximal, F is modified. 206 * @return n/F a factor of n not divisible by any prime number in SMPRM. 207 */ 208 public static long smallPrimeDivisors(long n, SortedMap<Long, Integer> F) { 209 //SortedMap<Long, Integer> F = new TreeMap<Long, Integer>(); 210 List<Long> LP; 211 long QL = 0; 212 long PL; 213 long RL = 0; 214 boolean TL; 215 216 long ML = n; 217 LP = SMPRM; //smallPrimes(2, 500); //SMPRM; 218 TL = false; 219 int i = 0; 220 do { 221 PL = LP.get(i); 222 QL = ML / PL; 223 RL = ML % PL; 224 if (RL == 0) { 225 Integer e = F.get(PL); 226 if (e == null) { 227 e = 1; 228 } else { 229 e++; 230 } 231 F.put(PL, e); 232 ML = QL; 233 } else { 234 i++; 235 } 236 TL = (QL <= PL); 237 } while (!(TL || (i >= LP.size()))); 238 //System.out.println("TL = " + TL + ", ML = " + ML + ", PL = " + PL + ", QL = " + QL); 239 if (TL && (ML != 1L)) { 240 Integer e = F.get(ML); 241 if (e == null) { 242 e = 1; 243 } else { 244 e = e + 1; 245 } 246 F.put(ML, e); 247 ML = 1; 248 } 249 //F.put(ML, 0); // hack 250 return ML; 251 } 252 253 254 /** 255 * Integer small prime divisors. n is a positive integer. F is a list of 256 * primes (q(1),q(2),...,q(h)), h non-negative, q(1) le q(2) le ... lt q(h), 257 * such that n is equal to m times the product of the q(i) and m is not 258 * divisible by any prime in SMPRM. Either m=1 or m gt 1,000,000. 259 * <br /> In JAS F is a map and m=1 or m > 4.000.000. 260 * See also SACPRIM.ISPD. 261 * @param n integer to factor. 262 * @param F a map of pairs of prime numbers and multiplicities (p,e) with p**e 263 * divides n and e maximal, F is modified. 264 * @return n/F a factor of n not divisible by any prime number in SMPRM. 265 */ 266 public static java.math.BigInteger smallPrimeDivisors(java.math.BigInteger n, SortedMap<java.math.BigInteger, Integer> F) { 267 List<Long> LP; 268 java.math.BigInteger QL = java.math.BigInteger.ZERO; 269 java.math.BigInteger PL; 270 java.math.BigInteger RL = java.math.BigInteger.ZERO; 271 boolean TL; 272 273 java.math.BigInteger ML = n; 274 LP = SMPRM; //smallPrimes(2, 500); //SMPRM; 275 TL = false; 276 int i = 0; 277 do { 278 PL = java.math.BigInteger.valueOf( LP.get(i) ); 279 java.math.BigInteger[] xx = ML.divideAndRemainder(PL); 280 QL = xx[0]; //ML.divide(PL); 281 RL = xx[1]; //ML.remainder(PL); 282 if (RL.equals(java.math.BigInteger.ZERO)) { 283 Integer e = F.get(PL); 284 if (e == null) { 285 e = 1; 286 } else { 287 e++; 288 } 289 F.put(PL, e); 290 ML = QL; 291 } else { 292 i++; 293 } 294 TL = (QL.compareTo(PL) <= 0); 295 } while (!(TL || (i >= LP.size()))); 296 //System.out.println("TL = " + TL + ", ML = " + ML + ", PL = " + PL + ", QL = " + QL); 297 if (TL && (!ML.equals(java.math.BigInteger.ONE))) { 298 Integer e = F.get(ML); 299 if (e == null) { 300 e = 1; 301 } else { 302 e = e + 1; 303 } 304 F.put(ML, e); 305 ML = java.math.BigInteger.ONE; 306 } 307 //F.put(ML, 0); // hack 308 return ML; 309 } 310 311 312 /** 313 * Integer primality test. n is a positive integer. r is true, if n is 314 * prime, else false. 315 * @param n integer to test. 316 * @return true if n is prime, else false. 317 */ 318 public static boolean isPrime(long n) { 319 java.math.BigInteger N = java.math.BigInteger.valueOf(n); 320 return isPrime(N); 321 } 322 323 324 /** 325 * Integer primality test. n is a positive integer. r is true, if n is 326 * prime, else false. 327 * @param N integer to test. 328 * @return true if N is prime, else false. 329 */ 330 public static boolean isPrime(java.math.BigInteger N) { 331 if (N.isProbablePrime(N.bitLength())) { 332 return true; 333 } 334 SortedMap<java.math.BigInteger, Integer> F = factors(N); 335 return (F.size() == 1) && F.values().contains(1); 336 } 337 338 339 /** 340 * Test prime factorization. n is a positive integer. r is true, if n = 341 * product_i(pi**ei) and each pi is prime, else false. 342 * @param n integer to test. 343 * @param F a map of pairs of prime numbers (p,e) with p**e divides n. 344 * @return true if n = product_i(pi**ei) and each pi is prime, else false. 345 */ 346 public static boolean isPrimeFactorization(long n, SortedMap<Long, Integer> F) { 347 long f = 1L; 348 for (Map.Entry<Long, Integer> m : F.entrySet()) { 349 long p = m.getKey(); 350 if (!isPrime(p)) { 351 return false; 352 } 353 int e = m.getValue(); 354 long pe = java.math.BigInteger.valueOf(p).pow(e).longValueExact(); 355 f *= pe; 356 } 357 return n == f; 358 } 359 360 361 /** 362 * Test factorization. n is a positive integer. r is true, if n = 363 * product_i(pi**ei), else false. 364 * @param n integer to test. 365 * @param F a map of pairs of numbers (p,e) with p**e divides n. 366 * @return true if n = product_i(pi**ei), else false. 367 */ 368 public static boolean isFactorization(long n, SortedMap<Long, Integer> F) { 369 long f = 1L; 370 for (Map.Entry<Long, Integer> m : F.entrySet()) { 371 long p = m.getKey(); 372 int e = m.getValue(); 373 long pe = java.math.BigInteger.valueOf(p).pow(e).longValueExact(); 374 f *= pe; 375 } 376 return n == f; 377 } 378 379 380 /** 381 * Integer factorization. n is a positive integer. F is a list (q(1), 382 * q(2),...,q(h)) of the prime factors of n, q(1) le q(2) le ... le q(h), 383 * with n equal to the product of the q(i). <br /> In JAS F is a map. 384 * See also SACPRIM.IFACT. 385 * @param n integer to factor. 386 * @return a map of pairs of numbers (p,e) with p**e divides n. 387 */ 388 public static SortedMap<Long, Integer> factors(long n) { 389 if (n > BETA) { 390 throw new UnsupportedOperationException("factors(long) only for longs less than BETA: " + BETA); 391 } 392 long ML, PL, AL, BL, CL, MLP, RL, SL; 393 SortedMap<Long, Integer> F = new TreeMap<Long, Integer>(); 394 SortedMap<Long, Integer> FP = null; 395 // search small prime factors 396 ML = smallPrimeDivisors(n, F); // , F, ML 397 if (ML == 1L) { 398 return F; 399 } 400 //System.out.println("F = " + F); 401 // search medium prime factors 402 AL = IMPDS_MIN; 403 do { 404 MLP = ML - 1; 405 RL = (new ModLong(new ModLongRing(ML), 3)).power(MLP).getVal(); //(3**MLP) mod ML; 406 if (RL == 1L) { 407 FP = factors(MLP); 408 SL = primalityTestSelfridge(ML, MLP, FP); 409 if (SL == 1) { 410 logger.info("primalityTestSelfridge: FP = {}", FP); 411 Integer e = F.get(ML); 412 if (e == null) { 413 e = 1; 414 } else { // will not happen 415 e = e + 1; 416 } 417 F.put(ML, e); 418 return F; 419 } 420 } 421 CL = Roots.sqrtInt(new BigInteger(ML)).getVal().longValue(); //SACI.ISQRT( ML, CL, TL ); 422 //System.out.println("CL = " + CL + ", ML = " + ML + ", CL^2 = " + (CL*CL)); 423 BL = Math.max(IMPDS_MAX, CL / 3L); 424 if (AL > BL) { 425 PL = 1L; 426 } else { 427 logger.info("mediumPrimeDivisorSearch: a = {}, b = {}", AL, BL); 428 PL = mediumPrimeDivisorSearch(ML, AL, BL); //, PL, ML ); 429 //System.out.println("PL = " + PL); 430 if (PL != 1L) { 431 AL = PL; 432 Integer e = F.get(PL); 433 if (e == null) { 434 e = 1; 435 } else { 436 e = e + 1; 437 } 438 F.put(PL, e); 439 ML = ML / PL; 440 } 441 } 442 } while (PL != 1L); 443 // fixed: the ILPDS should also be in the while loop, was already wrong in SAC2/Aldes and MAS 444 // seems to be okay for integers smaller than beta 445 java.math.BigInteger N = java.math.BigInteger.valueOf(ML); 446 if (N.isProbablePrime(N.bitLength())) { 447 F.put(ML, 1); 448 return F; 449 } 450 AL = BL; 451 BL = CL; 452 logger.info("largePrimeDivisorSearch: a = {}, b = {}, m = {}", AL, BL, ML); 453 // search large prime factors 454 do { 455 //ILPDS( ML, AL, BL, PL, ML ); 456 PL = largePrimeDivisorSearch(ML, AL, BL); 457 if (PL != 1L) { 458 Integer e = F.get(PL); 459 if (e == null) { 460 e = 1; 461 } else { 462 e++; 463 } 464 F.put(PL, e); 465 ML = ML / PL; 466 AL = PL; 467 CL = Roots.sqrtInt(BigInteger.valueOf(ML)).getVal().longValue(); //SACI.ISQRT( ML, CL, TL ); 468 //System.out.println("CL = " + CL + ", ML = " + ML + ", CL^2 = " + (CL*CL)); 469 BL = Math.min(BL, CL); 470 if (AL > BL) { 471 PL = 1L; 472 } 473 } 474 } while (PL != 1L); 475 //System.out.println("PL = " + PL + ", ML = " + ML); 476 if (ML != 1L) { 477 Integer e = F.get(ML); 478 if (e == null) { 479 e = 1; 480 } else { 481 e = e + 1; 482 } 483 F.put(ML, e); 484 } 485 return F; 486 } 487 488 489 /** 490 * Integer factorization, Pollard rho algorithm. n is a positive integer. F 491 * is a list (q(1), q(2),...,q(h)) of the prime factors of n, q(1) le q(2) 492 * le ... le q(h), with n equal to the product of the q(i). <br /> In 493 * JAS F is a map. 494 * See also SACPRIM.IFACT. 495 * @param n integer to factor. 496 * @return a map F of pairs of numbers (p,e) with p**e divides n and p 497 * probable prime. 498 */ 499 public static SortedMap<Long, Integer> factorsPollard(long n) { 500 if (n > BETA) { 501 throw new UnsupportedOperationException("factors(long) only for longs less than BETA: " + BETA); 502 } 503 SortedMap<Long, Integer> F = new TreeMap<Long, Integer>(); 504 factorsPollardRho(n, F); 505 return F; 506 } 507 508 509 /** 510 * Integer medium prime divisor search. n, a and b are positive integers 511 * such that a le b le n and n has no positive divisors less than a. If n 512 * has a prime divisor in the closed interval from a to b then p is the 513 * least such prime and q=n/p. Otherwise p=1 and q=n. 514 * See also SACPRIM.IMPDS. 515 * @param n integer to factor. 516 * @param a lower bound. 517 * @param b upper bound. 518 * @return p a prime factor of n, with a ≤ p ≤ b < n. 519 */ 520 public static long mediumPrimeDivisorSearch(long n, long a, long b) { 521 List<Long> LP; 522 long R, J1Y, RL1, RL2, RL, PL; 523 524 RL = a % 210; 525 LP = UZ210; 526 long ll = LP.size(); 527 int i = 0; 528 while (RL > LP.get(i)) { 529 i++; 530 } 531 RL1 = LP.get(i); 532 PL = a + (RL1 - RL); 533 //System.out.println("PL = " + PL + ", BL = " + BL); 534 while (PL <= b) { 535 R = n % PL; //SACI.IQR( NL, PL, QL, R ); 536 if (R == 0) { 537 return PL; 538 } 539 i++; 540 if (i >= ll) { 541 LP = UZ210; 542 RL2 = (RL1 - 210L); 543 i = 0; 544 } else { 545 RL2 = RL1; 546 } 547 RL1 = LP.get(i); 548 J1Y = (RL1 - RL2); 549 PL = PL + J1Y; 550 } 551 PL = 1L; //SACI.IONE; 552 //QL = NL; 553 return PL; 554 } 555 556 557 /** 558 * Integer selfridge primality test. m is an integer greater than or equal 559 * to 3. mp=m-1. F is a list (q(1),q(2),...,q(k)), q(1) le q(2) le ... le 560 * q(k), of the prime factors of mp, with mp equal to the product of the 561 * q(i). An attempt is made to find a root of unity modulo m of order m-1. 562 * If the existence of such a root is discovered then m is prime and s=1. If 563 * it is discovered that no such root exists then m is not a prime and s=-1. 564 * Otherwise the primality of m remains uncertain and s=0. 565 * See also SACPRIM.ISPT. 566 * @param m integer to test. 567 * @param mp integer m-1. 568 * @param F a map of pairs (p,e), with primes p, multiplicity e and with 569 * p**e divides mp and e maximal. 570 * @return s = -1 (not prime), 0 (unknown) or 1 (prime). 571 */ 572 public static int primalityTestSelfridge(long m, long mp, SortedMap<Long, Integer> F) { 573 long AL, BL, QL, QL1, MLPP, PL1, PL; 574 int SL; 575 //List<Long> SMPRM = smallPrimes(2, 500); //SMPRM; 576 List<Long> PP; 577 578 List<Map.Entry<Long, Integer>> FP = new ArrayList<Map.Entry<Long, Integer>>(F.entrySet()); 579 QL1 = 1L; //SACI.IONE; 580 PL1 = 1L; 581 int i = 0; 582 while (true) { 583 do { 584 if (i == FP.size()) { 585 logger.info("SL=1: m = {}", m); 586 SL = 1; 587 return SL; 588 } 589 QL = FP.get(i).getKey(); 590 i++; 591 } while (!(QL > QL1)); 592 QL1 = QL; 593 PP = SMPRM; 594 int j = 0; 595 do { 596 if (j == PP.size()) { 597 logger.info("SL=0: m = {}", m); 598 SL = 0; 599 return SL; 600 } 601 PL = PP.get(j); 602 j++; 603 if (PL > PL1) { 604 PL1 = PL; 605 AL = (new ModLong(new ModLongRing(m), PL)).power(mp).getVal(); //(PL**MLP) mod ML; 606 if (AL != 1) { 607 logger.info("SL=-1: m = {}", m); 608 SL = (-1); 609 return SL; 610 } 611 } 612 MLPP = mp / QL; 613 BL = (new ModLong(new ModLongRing(m), PL)).power(MLPP).getVal(); //(PL**MLPP) mod ML; 614 } while (BL == 1L); 615 } 616 } 617 618 619 /** 620 * Integer large prime divisor search. n is a positive integer with no prime 621 * divisors less than 17. 1 le a le b le n. A search is made for a divisor p 622 * of the integer n, with a le p le b. If such a p is found then np=n/p, 623 * otherwise p=1 and np=n. A modular version of Fermats method is used, and 624 * the search goes from a to b. 625 * See also SACPRIM.ILPDS. 626 * @param n integer to factor. 627 * @param a lower bound. 628 * @param b upper bound. 629 * @return p a prime factor of n, with a ≤ p ≤ b < n. 630 */ 631 public static long largePrimeDivisorSearch(long n, long a, long b) { // return PL, NLP ignored 632 if (n > BETA) { 633 throw new UnsupportedOperationException( 634 "largePrimeDivisorSearch only for longs less than BETA: " + BETA); 635 } 636 List<ModLong> L = null; 637 List<ModLong> LP; 638 long RL1, RL2, J1Y, r, PL, TL; 639 long RL, J2Y, XL1, XL2, QL, XL, YL, YLP; 640 long ML = 0L; 641 long SL = 0L; 642 QL = n / b; 643 RL = n % b; 644 XL1 = b + QL; 645 SL = XL1 % 2L; 646 XL1 = XL1 / 2L; // after SL 647 if ((RL != 0) || (SL != 0)) { 648 XL1 = XL1 + 1L; 649 } 650 QL = n / a; 651 XL2 = a + QL; 652 XL2 = XL2 / 2L; 653 L = residueListFermat(n); //FRESL( NL, ML, L ); // ML not returned 654 if (L.isEmpty()) { 655 return n; 656 } 657 ML = L.get(0).ring.getModul().longValue(); // sic 658 // check is okay: sort: L = SACSET.LBIBMS( L ); revert: L = MASSTOR.INV( L ); 659 Collections.sort(L); 660 Collections.reverse(L); 661 //System.out.println("FRESL: " + L); 662 r = XL2 % ML; 663 LP = L; 664 int i = 0; 665 while (i < LP.size() && r < LP.get(i).getVal()) { 666 i++; 667 } 668 if (i == LP.size()) { 669 i = 0; //LP = L; 670 SL = ML; 671 } else { 672 SL = 0L; 673 } 674 RL1 = LP.get(i).getVal(); 675 i++; 676 SL = ((SL + r) - RL1); 677 XL = XL2 - SL; 678 TL = 0L; 679 while (XL >= XL1) { 680 J2Y = XL * XL; 681 YLP = J2Y - n; 682 //System.out.println("YLP = " + YLP + ", J2Y = " + J2Y); 683 YL = Roots.sqrtInt(BigInteger.valueOf(YLP)).getVal().longValue(); // SACI.ISQRT( YLP, YL, TL ); 684 //System.out.println("YL = sqrt(YLP) = " + YL); 685 TL = YLP - YL * YL; 686 if (TL == 0L) { 687 PL = XL - YL; 688 return PL; 689 } 690 if (i < LP.size()) { 691 RL2 = LP.get(i).getVal(); 692 i++; 693 SL = (RL1 - RL2); 694 } else { 695 i = 0; 696 RL2 = LP.get(i).getVal(); 697 i++; 698 J1Y = (ML + RL1); 699 SL = (J1Y - RL2); 700 } 701 RL1 = RL2; 702 XL = XL - SL; 703 } 704 PL = 1L; 705 // unused NLP = NL; 706 return PL; 707 } 708 709 710 /** 711 * Fermat residue list, single modulus. m is a positive beta-integer. a 712 * belongs to Z(m). L is a list of the distinct b in Z(m) such that b**2-a 713 * is a square in Z(m). 714 * See also SACPRIM.FRLSM. 715 * @param m integer to factor. 716 * @param a element of Z mod m. 717 * @return Lp a list of Fermat residues for modul m. 718 */ 719 public static List<ModLong> residueListFermatSingle(long m, long a) { 720 List<ModLong> Lp; 721 SortedSet<ModLong> L; 722 List<ModLong> S, SP; 723 int MLP; 724 ModLong SL, SLP, SLPP; 725 726 ModLongRing ring = new ModLongRing(m); 727 ModLong am = ring.fromInteger(a); 728 MLP = (int) (m / 2L); 729 S = new ArrayList<ModLong>(); 730 for (int i = 0; i <= MLP; i++) { 731 SL = ring.fromInteger(i); 732 SL = SL.multiply(SL); //SACM.MDPROD( ML, IL, IL ); 733 S.add(SL); 734 } 735 L = new TreeSet<ModLong>(); 736 SP = S; 737 for (int i = MLP; i >= 0; i -= 1) { 738 SL = SP.get(i); 739 SLP = SL.subtract(am); //SACM.MDDIF( ML, SL, AL ); 740 int j = S.indexOf(SLP); 741 if (j >= 0) { // != 0 742 SLP = ring.fromInteger(i); 743 L.add(SLP); 744 SLPP = SLP.negate(); 745 if (!SLPP.equals(SLP)) { 746 L.add(SLPP); 747 } 748 } 749 } 750 Lp = new ArrayList<ModLong>(L); 751 return Lp; 752 } 753 754 755 /** 756 * Fermat residue list. n is a positive integer with no prime divisors less 757 * than 17. m is a positive beta-integer and L is an ordered list of the 758 * elements of Z(m) such that if x**2-n is a square then x is congruent to a 759 * (modulo m) for some a in L. 760 * See also SACPRIM.FRESL. 761 * @param n integer to factor. 762 * @return Lp a list of Fermat residues for different modules. 763 */ 764 public static List<ModLong> residueListFermat(long n) { 765 List<ModLong> L, L1; 766 List<Long> H, M; 767 long AL1, AL2, AL3, AL4, BL1, HL, J1Y, J2Y, KL, KL1, ML1, ML; 768 //too large: long BETA = Long.MAX_VALUE - 1L; 769 770 // modulus 2**5. 771 BL1 = 0L; 772 AL1 = n % 32L; 773 AL2 = AL1 % 16L; 774 AL3 = AL2 % 8L; 775 AL4 = AL3 % 4L; 776 if (AL4 == 3L) { 777 ML = 4L; 778 if (AL3 == 3L) { 779 BL1 = 2L; 780 } else { 781 BL1 = 0L; 782 } 783 } else { 784 if (AL3 == 1L) { 785 ML = 8L; 786 if (AL2 == 1L) { 787 BL1 = 1L; 788 } else { 789 BL1 = 3L; 790 } 791 } else { 792 ML = 16L; 793 switch ((short) (AL1 / 8L)) { 794 case (short) 0: 795 BL1 = 3L; 796 break; 797 case (short) 1: 798 BL1 = 7L; 799 break; 800 case (short) 2: 801 BL1 = 5L; 802 break; 803 case (short) 3: 804 BL1 = 1L; 805 break; 806 default: 807 throw new RuntimeException("this should not happen"); 808 } 809 } 810 } 811 L = new ArrayList<ModLong>(); 812 ModLongRing ring = new ModLongRing(ML); 813 ModLongRing ring2; 814 if (ML == 4L) { 815 L.add(ring.fromInteger(BL1)); 816 } else { 817 J1Y = ML - BL1; 818 L.add(ring.fromInteger(BL1)); 819 L.add(ring.fromInteger(J1Y)); 820 } 821 KL = L.size(); 822 823 // modulus 3**3. 824 AL1 = n % 27L; 825 AL2 = AL1 % 3L; 826 if (AL2 == 2L) { 827 ML1 = 3L; 828 ring2 = new ModLongRing(ML1); 829 KL1 = 1L; 830 L1 = new ArrayList<ModLong>(); 831 L1.add(ring2.fromInteger(0)); 832 } else { 833 ML1 = 27L; 834 ring2 = new ModLongRing(ML1); 835 KL1 = 4L; 836 L1 = residueListFermatSingle(ML1, AL1); 837 // ring2 == L1.get(0).ring 838 } 839 //L = SACM.MDLCRA( ML, ML1, L, L1 ); 840 L = ModLongRing.chineseRemainder(ring.getONE(), ring2.getONE(), L, L1); 841 ML = (ML * ML1); 842 ring = new ModLongRing(ML); // == L.get(0).ring 843 KL = (KL * KL1); 844 //System.out.println("FRESL: L = " + L + ", ring = " + ring.toScript()); 845 846 // modulus 5**2. 847 AL1 = n % 25L; 848 AL2 = AL1 % 5L; 849 if ((AL2 == 2L) || (AL2 == 3L)) { 850 ML1 = 5L; 851 ring2 = new ModLongRing(ML1); 852 J1Y = (AL2 - 1L); 853 J2Y = (6L - AL2); 854 L1 = new ArrayList<ModLong>(); 855 L1.add(ring2.fromInteger(J1Y)); 856 L1.add(ring2.fromInteger(J2Y)); 857 KL1 = 2L; 858 } else { 859 ML1 = 25L; 860 ring2 = new ModLongRing(ML1); 861 L1 = residueListFermatSingle(ML1, AL1); 862 KL1 = 7L; 863 } 864 if (ML1 >= BETA / ML) { 865 return L; 866 } 867 //L = SACM.MDLCRA( ML, ML1, L, L1 ); 868 L = ModLongRing.chineseRemainder(ring.getONE(), ring2.getONE(), L, L1); 869 ML = (ML * ML1); 870 ring = new ModLongRing(ML); 871 KL = (KL * KL1); 872 //System.out.println("FRESL: L = " + L + ", ring = " + ring.toScript()); 873 874 // moduli 7,11,13. 875 L1 = new ArrayList<ModLong>(); 876 M = new ArrayList<Long>(3); 877 H = new ArrayList<Long>(3); 878 //M = MASSTOR.COMPi( 7, MASSTOR.COMPi( 11, 13 ) ); 879 M.add(7L); 880 M.add(11L); 881 M.add(13L); 882 //H = MASSTOR.COMPi( 64, MASSTOR.COMPi( 48, 0 ) ); 883 H.add(64L); 884 H.add(48L); 885 H.add(0L); 886 int i = 0; 887 while (true) { 888 ML1 = M.get(i); 889 if (ML1 >= BETA / ML) { 890 return L; 891 } 892 ring2 = new ModLongRing(ML1); 893 AL1 = n % ML1; 894 L1 = residueListFermatSingle(ML1, AL1); 895 KL1 = L1.size(); 896 //L = SACM.MDLCRA( ML, ML1, L, L1 ); 897 L = ModLongRing.chineseRemainder(ring.getONE(), ring2.getONE(), L, L1); 898 ML = (ML * ML1); 899 ring = new ModLongRing(ML); 900 KL = (KL * KL1); 901 //System.out.println("FRESL: L = " + L + ", ring = " + ring.toScript()); 902 HL = H.get(i); 903 i++; 904 if (KL > HL) { 905 return L; 906 } 907 } 908 // return ? 909 } 910 911 912 /** 913 * Compute units of Z sub 210. 914 * See also SACPRIM.UZ210. 915 * @return list of units of Z sub 210. 916 */ 917 public static List<Long> getUZ210() { 918 List<Long> UZ = new ArrayList<Long>(); 919 java.math.BigInteger z210 = java.math.BigInteger.valueOf(210); 920 //for (int i = 209; i >= 1; i -= 2) { 921 for (long i = 1; i <= 209; i += 2) { 922 if (z210.gcd(java.math.BigInteger.valueOf(i)).equals(java.math.BigInteger.ONE)) { 923 UZ.add(i); 924 } 925 } 926 return UZ; 927 } 928 929 930 /** 931 * Integer factorization. n is a positive integer. F is a list (q(1), 932 * q(2),...,q(h)) of the prime factors of n, q(1) le q(2) le ... le q(h), 933 * with n equal to the product of the q(i). <br /> In JAS F is a map. 934 * See also SACPRIM.IFACT, uses Pollards rho method. 935 * @param n integer to factor. 936 * @return a map of pairs of numbers (p,e) with p**e divides n. 937 */ 938 public static SortedMap<java.math.BigInteger, Integer> factors(java.math.BigInteger n) { 939 java.math.BigInteger b = java.math.BigInteger.valueOf(BETA); 940 SortedMap<java.math.BigInteger, Integer> F = new TreeMap<java.math.BigInteger, Integer>(); 941 if (n.compareTo(b) > 0) { 942 n = smallPrimeDivisors(n, F); 943 if (n.compareTo(b) > 0) { 944 logger.info("run factorsPollardRho on n = {}", n); 945 factorsPollardRho(n, F); 946 return F; 947 } 948 } 949 // n <= beta 950 long s = n.longValue(); 951 SortedMap<Long, Integer> ff = factors(s); // useless 2nd smallPrimeDiv search 952 for (Map.Entry<Long, Integer> m : ff.entrySet()) { 953 java.math.BigInteger mm = java.math.BigInteger.valueOf(m.getKey()); 954 F.put(mm, m.getValue()); 955 } 956 return F; 957 } 958 959 960 /** 961 * Integer factorization using Pollards rho algorithm. n is a positive 962 * integer. F is a list (q(1), q(2),...,q(h)) of the prime factors of n, 963 * q(1) le q(2) le ... le q(h), with n equal to the product of the q(i). 964 * <br /> In JAS F is a map. 965 * @param n integer to factor. 966 * @param F a map of pairs of numbers (p,e) with p**e divides n and p is 967 * probable prime, F is modified. 968 */ 969 public static void factorsPollardRho(java.math.BigInteger n, SortedMap<java.math.BigInteger, Integer> F) { 970 java.math.BigInteger factor; 971 java.math.BigInteger temp = n; 972 int iterationCounter = 0; 973 Integer count; 974 while (!temp.isProbablePrime(32)) { 975 factor = rho(temp); 976 if (factor.equals(temp)) { 977 if (iterationCounter++ > 4) { 978 break; 979 } 980 } else { 981 iterationCounter = 1; 982 } 983 count = F.get(factor); 984 if (count == null) { 985 F.put(factor, 1); 986 } else { 987 F.put(factor, count + 1); 988 } 989 temp = temp.divide(factor); 990 } 991 count = F.get(temp); 992 if (count == null) { 993 F.put(temp, 1); 994 } else { 995 F.put(temp, count + 1); 996 } 997 } 998 999 1000 /** 1001 * Random number generator. 1002 */ 1003 //final static SecureRandom random = new SecureRandom(); 1004 final static Random random = new Random(); 1005 1006 1007 /** 1008 * Search cycle with Pollards rho algorithm x**2 + c mod n. n is a positive 1009 * integer. <br /> 1010 * @param n integer test. 1011 * @return x-y with gcd(x-y, n) = 1. 1012 */ 1013 static java.math.BigInteger rho(java.math.BigInteger n) { 1014 java.math.BigInteger divisor; 1015 java.math.BigInteger c = new java.math.BigInteger(n.bitLength(), random); 1016 java.math.BigInteger x = new java.math.BigInteger(n.bitLength(), random); 1017 java.math.BigInteger xx = x; 1018 do { 1019 x = x.multiply(x).mod(n).add(c).mod(n); 1020 xx = xx.multiply(xx).mod(n).add(c).mod(n); 1021 xx = xx.multiply(xx).mod(n).add(c).mod(n); 1022 divisor = x.subtract(xx).gcd(n); 1023 } while (divisor.equals(java.math.BigInteger.ONE)); 1024 return divisor; 1025 } 1026 1027 1028 /** 1029 * Integer factorization using Pollards rho algorithm. n is a positive 1030 * integer. F is a list (q(1), q(2),...,q(h)) of the prime factors of n, 1031 * q(1) le q(2) le ... le q(h), with n equal to the product of the q(i). 1032 * <br /> In JAS F is a map. 1033 * @param n integer to factor. 1034 * @param F a map of pairs of numbers (p,e) with p**e divides n and p is 1035 * probable prime, F is modified. 1036 */ 1037 public static void factorsPollardRho(long n, SortedMap<Long, Integer> F) { 1038 long factor; 1039 long temp = n; 1040 int iterationCounter = 0; 1041 Integer count; 1042 while (!java.math.BigInteger.valueOf(temp).isProbablePrime(32)) { 1043 factor = rho(temp); 1044 if (factor == temp) { 1045 if (iterationCounter++ > 4) { 1046 break; 1047 } 1048 } else { 1049 iterationCounter = 1; 1050 } 1051 count = F.get(factor); 1052 if (count == null) { 1053 F.put(factor, 1); 1054 } else { 1055 F.put(factor, count + 1); 1056 } 1057 temp = temp / factor; 1058 } 1059 count = F.get(temp); 1060 if (count == null) { 1061 F.put(temp, 1); 1062 } else { 1063 F.put(temp, count + 1); 1064 } 1065 //System.out.println("random = " + random.getAlgorithm()); 1066 } 1067 1068 1069 /** 1070 * Search cycle with Pollards rho algorithm x**2 + c mod n. n is a positive 1071 * integer. c is a random constant. 1072 * @param n integer test. 1073 * @return x-y with gcd(x-y, n) == 1. 1074 */ 1075 static long rho(long n) { 1076 long divisor; 1077 int bl = java.math.BigInteger.valueOf(n).bitLength(); 1078 long c = new java.math.BigInteger(bl, random).longValue(); // .abs() 1079 long x = new java.math.BigInteger(bl, random).longValue(); // .abs() 1080 ModLongRing ring = new ModLongRing(n); 1081 ModLong cm = new ModLong(ring, c); 1082 ModLong xm = new ModLong(ring, x); 1083 ModLong xxm = xm; 1084 do { 1085 xm = xm.multiply(xm).sum(cm); 1086 xxm = xxm.multiply(xxm).sum(cm); 1087 xxm = xxm.multiply(xxm).sum(cm); 1088 divisor = gcd(xm.getVal() - xxm.getVal(), n); 1089 } while (divisor == 1L); 1090 return divisor; 1091 } 1092 1093 1094 static long gcd(long a, long b) { 1095 return BigInteger.valueOf(a).gcd(BigInteger.valueOf(b)).getVal().longValue(); 1096 } 1097 1098}