Source code

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 &le; p &lt; 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 &gt; 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 &gt; 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 &le; p &le; b &lt; 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 &le; p &le; b &lt; 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}