Source code

001/*
002 * $Id: PrimeInteger.java 5940 2018-10-19 08:53:13Z kredel $
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 = 2 * (k - 1);
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        if (N.isProbablePrime(N.bitLength())) {
321            return true;
322        }
323        SortedMap<Long, Integer> F = factors(n);
324        return (F.size() == 1) && F.values().contains(1);
325    }
326
327
328    /**
329     * Test prime factorization. n is a positive integer. r is true, if n =
330     * product_i(pi**ei) and each pi is prime, else false.
331     * @param n integer to test.
332     * @param F a map of pairs of prime numbers (p,e) with p**e divides n.
333     * @return true if n = product_i(pi**ei) and each pi is prime, else false.
334     */
335    public static boolean isPrimeFactorization(long n, SortedMap<Long, Integer> F) {
336        long f = 1L;
337        for (Map.Entry<Long, Integer> m : F.entrySet()) {
338            long p = m.getKey();
339            if (!isPrime(p)) {
340                return false;
341            }
342            int e = m.getValue();
343            long pe = java.math.BigInteger.valueOf(p).pow(e).longValue();
344            f *= pe;
345        }
346        return n == f;
347    }
348
349
350    /**
351     * Test factorization. n is a positive integer. r is true, if n =
352     * product_i(pi**ei), else false.
353     * @param n integer to test.
354     * @param F a map of pairs of numbers (p,e) with p**e divides n.
355     * @return true if n = product_i(pi**ei), else false.
356     */
357    public static boolean isFactorization(long n, SortedMap<Long, Integer> F) {
358        long f = 1L;
359        for (Map.Entry<Long, Integer> m : F.entrySet()) {
360            long p = m.getKey();
361            int e = m.getValue();
362            long pe = java.math.BigInteger.valueOf(p).pow(e).longValue();
363            f *= pe;
364        }
365        return n == f;
366    }
367
368
369    /**
370     * Integer factorization. n is a positive integer. F is a list (q(1),
371     * q(2),...,q(h)) of the prime factors of n, q(1) le q(2) le ... le q(h),
372     * with n equal to the product of the q(i). <br /> In JAS F is a map.
373     * See also SACPRIM.IFACT.
374     * @param n integer to factor.
375     * @return a map of pairs of numbers (p,e) with p**e divides n.
376     */
377    public static SortedMap<Long, Integer> factors(long n) {
378        if (n > BETA) {
379            throw new UnsupportedOperationException("factors(long) only for longs less than BETA: " + BETA);
380        }
381        long ML, PL, AL, BL, CL, MLP, RL, SL;
382        SortedMap<Long, Integer> F = new TreeMap<Long, Integer>();
383        SortedMap<Long, Integer> FP = null;
384        // search small prime factors
385        ML = smallPrimeDivisors(n, F); // , F, ML
386        if (ML == 1L) {
387            return F;
388        }
389        //System.out.println("F = " + F);
390        // search medium prime factors
391        AL = IMPDS_MIN;
392        do {
393            MLP = ML - 1;
394            RL = (new ModLong(new ModLongRing(ML), 3)).power(MLP).getVal(); //(3**MLP) mod ML; 
395            if (RL == 1L) {
396                FP = factors(MLP);
397                SL = primalityTestSelfridge(ML, MLP, FP);
398                if (SL == 1) {
399                    logger.info("primalityTestSelfridge: FP = " + FP);
400                    Integer e = F.get(ML);
401                    if (e == null) {
402                        e = 1;
403                    } else { // will not happen
404                        e = e + 1;
405                    }
406                    F.put(ML, e); 
407                    return F;
408                }
409            }
410            CL = Roots.sqrtInt(new BigInteger(ML)).getVal().longValue(); //SACI.ISQRT( ML, CL, TL );
411            //System.out.println("CL = " + CL + ", ML = " + ML + ", CL^2 = " + (CL*CL));
412            BL = Math.max(IMPDS_MAX, CL / 3L);
413            if (AL > BL) {
414                PL = 1L;
415            } else {
416                logger.info("mediumPrimeDivisorSearch: a = " + AL + ", b = " + BL);
417                PL = mediumPrimeDivisorSearch(ML, AL, BL); //, PL, ML );
418                //System.out.println("PL = " + PL);
419                if (PL != 1L) {
420                    AL = PL;
421                    Integer e = F.get(PL);
422                    if (e == null) {
423                        e = 1;
424                    } else {
425                        e = e + 1;
426                    }
427                    F.put(PL, e); 
428                    ML = ML / PL;
429                }
430            }
431        } while (PL != 1L);
432        // fixed: the ILPDS should also be in the while loop, was already wrong in SAC2/Aldes and MAS
433        // seems to be okay for integers smaller than beta
434        java.math.BigInteger N = java.math.BigInteger.valueOf(ML);
435        if (N.isProbablePrime(N.bitLength())) {
436            F.put(ML, 1);
437            return F;
438        }
439        AL = BL;
440        BL = CL;
441        logger.info("largePrimeDivisorSearch: a = " + AL + ", b = " + BL + ", m = " + ML);
442        // search large prime factors
443        do {
444            //ILPDS( ML, AL, BL, PL, ML );
445            PL = largePrimeDivisorSearch(ML, AL, BL);
446            if (PL != 1L) {
447                Integer e = F.get(PL);
448                if (e == null) {
449                    e = 1;
450                } else {
451                    e++;
452                }
453                F.put(PL, e);
454                ML = ML / PL;
455                AL = PL;
456                CL = Roots.sqrtInt(BigInteger.valueOf(ML)).getVal().longValue(); //SACI.ISQRT( ML, CL, TL );
457                //System.out.println("CL = " + CL + ", ML = " + ML + ", CL^2 = " + (CL*CL));
458                BL = Math.min(BL, CL);
459                if (AL > BL) {
460                    PL = 1L;
461                }
462            }
463        } while (PL != 1L);
464        //System.out.println("PL = " + PL + ", ML = " + ML);
465        if (ML != 1L) {
466            Integer e = F.get(ML);
467            if (e == null) {
468                e = 1;
469            } else {
470                e = e + 1;
471            }
472            F.put(ML, e);
473        }
474        return F;
475    }
476
477
478    /**
479     * Integer factorization, Pollard rho algorithm. n is a positive integer. F
480     * is a list (q(1), q(2),...,q(h)) of the prime factors of n, q(1) le q(2)
481     * le ... le q(h), with n equal to the product of the q(i). <br /> In
482     * JAS F is a map.
483     * See also SACPRIM.IFACT.
484     * @param n integer to factor.
485     * @return a map F of pairs of numbers (p,e) with p**e divides n and p
486     *            probable prime.
487     */
488    public static SortedMap<Long, Integer> factorsPollard(long n) {
489        if (n > BETA) {
490            throw new UnsupportedOperationException("factors(long) only for longs less than BETA: " + BETA);
491        }
492        SortedMap<Long, Integer> F = new TreeMap<Long, Integer>();
493        factorsPollardRho(n, F); 
494        return F;
495    }
496
497
498    /**
499     * Integer medium prime divisor search. n, a and b are positive integers
500     * such that a le b le n and n has no positive divisors less than a. If n
501     * has a prime divisor in the closed interval from a to b then p is the
502     * least such prime and q=n/p. Otherwise p=1 and q=n.
503     * See also SACPRIM.IMPDS.
504     * @param n integer to factor.
505     * @param a lower bound.
506     * @param b upper bound.
507     * @return p a prime factor of n, with a &le; p &le; b &lt; n.
508     */
509    public static long mediumPrimeDivisorSearch(long n, long a, long b) {
510        List<Long> LP;
511        long R, J1Y, RL1, RL2, RL, PL;
512
513        RL = a % 210;
514        LP = UZ210;
515        long ll = LP.size();
516        int i = 0;
517        while (RL > LP.get(i)) {
518            i++; 
519        }
520        RL1 = LP.get(i); 
521        PL = a + (RL1 - RL); 
522        //System.out.println("PL = " + PL + ", BL = " + BL);
523        while (PL <= b) {
524            R = n % PL; //SACI.IQR( NL, PL, QL, R );
525            if (R == 0) {
526                return PL;
527            }
528            i++; 
529            if (i >= ll) { 
530                LP = UZ210;
531                RL2 = (RL1 - 210L);
532                i = 0;
533            } else {
534                RL2 = RL1;
535            }
536            RL1 = LP.get(i); 
537            J1Y = (RL1 - RL2);
538            PL = PL + J1Y; 
539        }
540        PL = 1L; //SACI.IONE;
541        //QL = NL;
542        return PL;
543    }
544
545
546    /**
547     * Integer selfridge primality test. m is an integer greater than or equal
548     * to 3. mp=m-1. F is a list (q(1),q(2),...,q(k)), q(1) le q(2) le ... le
549     * q(k), of the prime factors of mp, with mp equal to the product of the
550     * q(i). An attempt is made to find a root of unity modulo m of order m-1.
551     * If the existence of such a root is discovered then m is prime and s=1. If
552     * it is discovered that no such root exists then m is not a prime and s=-1.
553     * Otherwise the primality of m remains uncertain and s=0.
554     * See also SACPRIM.ISPT.
555     * @param m integer to test.
556     * @param mp integer m-1.
557     * @param F a map of pairs (p,e), with primes p, multiplicity e and with
558     *            p**e divides mp and e maximal.
559     * @return s = -1 (not prime), 0 (unknown) or 1 (prime).
560     */
561    public static int primalityTestSelfridge(long m, long mp, SortedMap<Long, Integer> F) {
562        long AL, BL, QL, QL1, MLPP, PL1, PL;
563        int SL;
564        //List<Long> SMPRM = smallPrimes(2, 500); //SMPRM;
565        List<Long> PP;
566
567        List<Map.Entry<Long, Integer>> FP = new ArrayList<Map.Entry<Long, Integer>>(F.entrySet());
568        QL1 = 1L; //SACI.IONE;
569        PL1 = 1L;
570        int i = 0;
571        while (true) {
572            do {
573                if (i == FP.size()) { 
574                    logger.info("SL=1: m = " + m);
575                    SL = 1;
576                    return SL;
577                }
578                QL = FP.get(i).getKey();
579                i++; 
580            } while (!(QL > QL1));
581            QL1 = QL;
582            PP = SMPRM;
583            int j = 0;
584            do {
585                if (j == PP.size()) {
586                    logger.info("SL=0: m = " + m);
587                    SL = 0;
588                    return SL;
589                }
590                PL = PP.get(j); 
591                j++; 
592                if (PL > PL1) {
593                    PL1 = PL;
594                    AL = (new ModLong(new ModLongRing(m), PL)).power(mp).getVal(); //(PL**MLP) mod ML; 
595                    if (AL != 1) {
596                        logger.info("SL=-1: m = " + m);
597                        SL = (-1);
598                        return SL;
599                    }
600                }
601                MLPP = mp / QL; 
602                BL = (new ModLong(new ModLongRing(m), PL)).power(MLPP).getVal(); //(PL**MLPP) mod ML; 
603            } while (BL == 1L);
604        }
605    }
606
607
608    /**
609     * Integer large prime divisor search. n is a positive integer with no prime
610     * divisors less than 17. 1 le a le b le n. A search is made for a divisor p
611     * of the integer n, with a le p le b. If such a p is found then np=n/p,
612     * otherwise p=1 and np=n. A modular version of Fermats method is used, and
613     * the search goes from a to b.
614     * See also SACPRIM.ILPDS.
615     * @param n integer to factor.
616     * @param a lower bound.
617     * @param b upper bound.
618     * @return p a prime factor of n, with a &le; p &le; b &lt; n.
619     */
620    public static long largePrimeDivisorSearch(long n, long a, long b) { // return PL, NLP ignored
621        if (n > BETA) {
622            throw new UnsupportedOperationException(
623                            "largePrimeDivisorSearch only for longs less than BETA: " + BETA);
624        }
625        List<ModLong> L = null;
626        List<ModLong> LP;
627        long RL1, RL2, J1Y, r, PL, TL;
628        long RL, J2Y, XL1, XL2, QL, XL, YL, YLP;
629        long ML = 0L;
630        long SL = 0L;
631        QL = n / b;
632        RL = n % b;
633        XL1 = b + QL;
634        SL = XL1 % 2L;
635        XL1 = XL1 / 2L; // after SL
636        if ((RL != 0) || (SL != 0)) {
637            XL1 = XL1 + 1L;
638        }
639        QL = n / a;
640        XL2 = a + QL;
641        XL2 = XL2 / 2L;
642        L = residueListFermat(n); //FRESL( NL, ML, L ); // ML not returned
643        if (L.isEmpty()) {
644            return n;
645        }
646        ML = L.get(0).ring.getModul().longValue(); // sic
647        // check is okay: sort: L = SACSET.LBIBMS( L ); revert: L = MASSTOR.INV( L );
648        Collections.sort(L);
649        Collections.reverse(L);
650        //System.out.println("FRESL: " + L);
651        r = XL2 % ML;
652        LP = L;
653        int i = 0;
654        while (i < LP.size() && r < LP.get(i).getVal()) {
655            i++; 
656        }
657        if (i == LP.size()) {
658            i = 0; //LP = L;
659            SL = ML;
660        } else {
661            SL = 0L;
662        }
663        RL1 = LP.get(i).getVal(); 
664        i++; 
665        SL = ((SL + r) - RL1);
666        XL = XL2 - SL;
667        TL = 0L;
668        while (XL >= XL1) {
669            J2Y = XL * XL;
670            YLP = J2Y - n;
671            //System.out.println("YLP = " + YLP + ", J2Y = " + J2Y);
672            YL = Roots.sqrtInt(BigInteger.valueOf(YLP)).getVal().longValue(); // SACI.ISQRT( YLP, YL, TL );
673            //System.out.println("YL = sqrt(YLP) = " + YL);
674            TL = YLP - YL * YL;
675            if (TL == 0L) {
676                PL = XL - YL;
677                return PL;
678            }
679            if (i < LP.size()) {
680                RL2 = LP.get(i).getVal(); 
681                i++; 
682                SL = (RL1 - RL2);
683            } else {
684                i = 0;
685                RL2 = LP.get(i).getVal(); 
686                i++; 
687                J1Y = (ML + RL1);
688                SL = (J1Y - RL2);
689            }
690            RL1 = RL2;
691            XL = XL - SL;
692        }
693        PL = 1L; 
694        // unused NLP = NL;
695        return PL;
696    }
697
698
699    /**
700     * Fermat residue list, single modulus. m is a positive beta-integer. a
701     * belongs to Z(m). L is a list of the distinct b in Z(m) such that b**2-a
702     * is a square in Z(m).
703     * See also SACPRIM.FRLSM.
704     * @param m integer to factor.
705     * @param a element of Z mod m.
706     * @return Lp a list of Fermat residues for modul m.
707     */
708    public static List<ModLong> residueListFermatSingle(long m, long a) {
709        List<ModLong> Lp;
710        SortedSet<ModLong> L;
711        List<ModLong> S, SP;
712        int MLP;
713        ModLong SL, SLP, SLPP;
714
715        ModLongRing ring = new ModLongRing(m);
716        ModLong am = ring.fromInteger(a);
717        MLP = (int) (m / 2L);
718        S = new ArrayList<ModLong>();
719        for (int i = 0; i <= MLP; i++) {
720            SL = ring.fromInteger(i);
721            SL = SL.multiply(SL); //SACM.MDPROD( ML, IL, IL );
722            S.add(SL); 
723        }
724        L = new TreeSet<ModLong>();
725        SP = S;
726        for (int i = MLP; i >= 0; i -= 1) {
727            SL = SP.get(i); 
728            SLP = SL.subtract(am); //SACM.MDDIF( ML, SL, AL );
729            int j = S.indexOf(SLP); 
730            if (j >= 0) { // != 0
731                SLP = ring.fromInteger(i);
732                L.add(SLP); 
733                SLPP = SLP.negate(); 
734                if (!SLPP.equals(SLP)) {
735                    L.add(SLPP); 
736                }
737            }
738        }
739        Lp = new ArrayList<ModLong>(L);
740        return Lp;
741    }
742
743
744    /**
745     * Fermat residue list. n is a positive integer with no prime divisors less
746     * than 17. m is a positive beta-integer and L is an ordered list of the
747     * elements of Z(m) such that if x**2-n is a square then x is congruent to a
748     * (modulo m) for some a in L.
749     * See also SACPRIM.FRESL.
750     * @param n integer to factor.
751     * @return Lp a list of Fermat residues for different modules.
752     */
753    public static List<ModLong> residueListFermat(long n) {
754        List<ModLong> L, L1;
755        List<Long> H, M;
756        long AL1, AL2, AL3, AL4, BL1, HL, J1Y, J2Y, KL, KL1, ML1, ML;
757        //too large: long BETA = Long.MAX_VALUE - 1L; 
758
759        // modulus 2**5.
760        BL1 = 0L;
761        AL1 = n % 32L; 
762        AL2 = AL1 % 16L; 
763        AL3 = AL2 % 8L; 
764        AL4 = AL3 % 4L; 
765        if (AL4 == 3L) {
766            ML = 4L;
767            if (AL3 == 3L) {
768                BL1 = 2L;
769            } else {
770                BL1 = 0L;
771            }
772        } else {
773            if (AL3 == 1L) {
774                ML = 8L;
775                if (AL2 == 1L) {
776                    BL1 = 1L;
777                } else {
778                    BL1 = 3L;
779                }
780            } else {
781                ML = 16L;
782                switch ((short) (AL1 / 8L)) {
783                case (short) 0: 
784                    BL1 = 3L;
785                    break;
786                case (short) 1: 
787                    BL1 = 7L;
788                    break;
789                case (short) 2: 
790                    BL1 = 5L;
791                    break;
792                case (short) 3: 
793                    BL1 = 1L;
794                    break;
795                default: 
796                    throw new RuntimeException("this should not happen");
797                }
798            }
799        }
800        L = new ArrayList<ModLong>();
801        ModLongRing ring = new ModLongRing(ML);
802        ModLongRing ring2;
803        if (ML == 4L) {
804            L.add(ring.fromInteger(BL1));
805        } else {
806            J1Y = ML - BL1;
807            L.add(ring.fromInteger(BL1));
808            L.add(ring.fromInteger(J1Y));
809        }
810        KL = L.size();
811
812        // modulus 3**3.
813        AL1 = n % 27L; 
814        AL2 = AL1 % 3L; 
815        if (AL2 == 2L) {
816            ML1 = 3L;
817            ring2 = new ModLongRing(ML1);
818            KL1 = 1L;
819            L1 = new ArrayList<ModLong>();
820            L1.add(ring2.fromInteger(0));
821        } else {
822            ML1 = 27L;
823            ring2 = new ModLongRing(ML1);
824            KL1 = 4L;
825            L1 = residueListFermatSingle(ML1, AL1);
826            // ring2 == L1.get(0).ring
827        }
828        //L = SACM.MDLCRA( ML, ML1, L, L1 );
829        L = ModLongRing.chineseRemainder(ring.getONE(), ring2.getONE(), L, L1);
830        ML = (ML * ML1);
831        ring = new ModLongRing(ML); // == L.get(0).ring
832        KL = (KL * KL1);
833        //System.out.println("FRESL: L = " + L + ", ring = " + ring.toScript());
834
835        // modulus 5**2.
836        AL1 = n % 25L; 
837        AL2 = AL1 % 5L; 
838        if ((AL2 == 2L) || (AL2 == 3L)) {
839            ML1 = 5L;
840            ring2 = new ModLongRing(ML1);
841            J1Y = (AL2 - 1L);
842            J2Y = (6L - AL2);
843            L1 = new ArrayList<ModLong>();
844            L1.add(ring2.fromInteger(J1Y));
845            L1.add(ring2.fromInteger(J2Y));
846            KL1 = 2L;
847        } else {
848            ML1 = 25L;
849            ring2 = new ModLongRing(ML1);
850            L1 = residueListFermatSingle(ML1, AL1);
851            KL1 = 7L;
852        }
853        if (ML1 >= BETA / ML) {
854            return L;
855        }
856        //L = SACM.MDLCRA( ML, ML1, L, L1 );
857        L = ModLongRing.chineseRemainder(ring.getONE(), ring2.getONE(), L, L1);
858        ML = (ML * ML1);
859        ring = new ModLongRing(ML);
860        KL = (KL * KL1);
861        //System.out.println("FRESL: L = " + L + ", ring = " + ring.toScript());
862
863        // moduli 7,11,13.
864        L1 = new ArrayList<ModLong>();
865        M = new ArrayList<Long>(3);
866        H = new ArrayList<Long>(3);
867        //M = MASSTOR.COMPi( 7, MASSTOR.COMPi( 11, 13 ) );
868        M.add(7L);
869        M.add(11L);
870        M.add(13L);
871        //H = MASSTOR.COMPi( 64, MASSTOR.COMPi( 48, 0 ) );
872        H.add(64L);
873        H.add(48L);
874        H.add(0L);
875        int i = 0;
876        while (true) {
877            ML1 = M.get(i); 
878            if (ML1 >= BETA / ML) {
879                return L;
880            }
881            ring2 = new ModLongRing(ML1);
882            AL1 = n % ML1; 
883            L1 = residueListFermatSingle(ML1, AL1);
884            KL1 = L1.size(); 
885            //L = SACM.MDLCRA( ML, ML1, L, L1 );
886            L = ModLongRing.chineseRemainder(ring.getONE(), ring2.getONE(), L, L1);
887            ML = (ML * ML1);
888            ring = new ModLongRing(ML);
889            KL = (KL * KL1);
890            //System.out.println("FRESL: L = " + L + ", ring = " + ring.toScript());
891            HL = H.get(i); 
892            i++;
893            if (KL > HL) {
894                return L;
895            }
896        }
897        // return ?
898    }
899
900
901    /**
902     * Compute units of Z sub 210.
903     * See also SACPRIM.UZ210.
904     * @return list of units of Z sub 210.
905     */
906    public static List<Long> getUZ210() {
907        List<Long> UZ = new ArrayList<Long>();
908        java.math.BigInteger z210 = java.math.BigInteger.valueOf(210);
909        //for (int i = 209; i >= 1; i -= 2) {
910        for (long i = 1; i <= 209; i += 2) {
911            if (z210.gcd(java.math.BigInteger.valueOf(i)).equals(java.math.BigInteger.ONE)) {
912                UZ.add(i);
913            }
914        }
915        return UZ;
916    }
917
918
919    /**
920     * Integer factorization. n is a positive integer. F is a list (q(1),
921     * q(2),...,q(h)) of the prime factors of n, q(1) le q(2) le ... le q(h),
922     * with n equal to the product of the q(i). <br /> In JAS F is a map.
923     * See also SACPRIM.IFACT, uses Pollards rho method.
924     * @param n integer to factor.
925     * @return a map of pairs of numbers (p,e) with p**e divides n.
926     */
927    public static SortedMap<java.math.BigInteger, Integer> factors(java.math.BigInteger n) {
928        java.math.BigInteger b = java.math.BigInteger.valueOf(BETA);
929        SortedMap<java.math.BigInteger, Integer> F = new TreeMap<java.math.BigInteger, Integer>();
930        if (n.compareTo(b) > 0) {
931            n = smallPrimeDivisors(n, F);
932            if (n.compareTo(b) > 0) {
933                logger.info("run factorsPollardRho on n = " + n);
934                factorsPollardRho(n, F); 
935                return F;
936            }
937        }
938        long s = n.longValue();
939        SortedMap<Long, Integer> ff = factors(s); // useless 2nd smallPrimeDiv search
940        for (Map.Entry<Long, Integer> m : ff.entrySet()) {
941            java.math.BigInteger mm = java.math.BigInteger.valueOf(m.getKey());
942            F.put(mm, m.getValue());            
943        }
944        return F;
945    }
946
947
948    /**
949     * Integer factorization using Pollards rho algorithm. n is a positive
950     * integer. F is a list (q(1), q(2),...,q(h)) of the prime factors of n,
951     * q(1) le q(2) le ... le q(h), with n equal to the product of the q(i).
952     * <br /> In JAS F is a map.
953     * @param n integer to factor.
954     * @param F a map of pairs of numbers (p,e) with p**e divides n and p is
955     *            probable prime, F is modified.
956     */
957    public static void factorsPollardRho(java.math.BigInteger n, SortedMap<java.math.BigInteger, Integer> F) {
958        java.math.BigInteger factor;
959        java.math.BigInteger temp = n;
960        int iterationCounter = 0;
961        Integer count;
962        while (!temp.isProbablePrime(32)) {
963            factor = rho(temp);
964            if (factor.equals(temp)) {
965                if (iterationCounter++ > 4) {
966                    break;
967                }
968            } else {
969                iterationCounter = 1;
970            }
971            count = F.get(factor);
972            if (count == null) {
973                F.put(factor, 1);
974            } else {
975                F.put(factor, count + 1);
976            }
977            temp = temp.divide(factor);
978        }
979        count = F.get(temp);
980        if (count == null) {
981            F.put(temp, 1);
982        } else {
983            F.put(temp, count + 1);
984        }
985    }
986
987
988    /**
989     * Random number generator.
990     */
991    //final static SecureRandom random = new SecureRandom();
992    final static Random random = new Random();
993
994
995    /**
996     * Search cycle with Pollards rho algorithm x**2 + c mod n. n is a positive
997     * integer. <br />
998     * @param n integer test.
999     * @return x-y with gcd(x-y, n) = 1.
1000     */
1001    static java.math.BigInteger rho(java.math.BigInteger n) {
1002        java.math.BigInteger divisor;
1003        java.math.BigInteger c = new java.math.BigInteger(n.bitLength(), random);
1004        java.math.BigInteger x = new java.math.BigInteger(n.bitLength(), random);
1005        java.math.BigInteger xx = x;
1006        do {
1007            x = x.multiply(x).mod(n).add(c).mod(n);
1008            xx = xx.multiply(xx).mod(n).add(c).mod(n);
1009            xx = xx.multiply(xx).mod(n).add(c).mod(n);
1010            divisor = x.subtract(xx).gcd(n);
1011        } while (divisor.equals(java.math.BigInteger.ONE));
1012        return divisor;
1013    }
1014
1015
1016    /**
1017     * Integer factorization using Pollards rho algorithm. n is a positive
1018     * integer. F is a list (q(1), q(2),...,q(h)) of the prime factors of n,
1019     * q(1) le q(2) le ... le q(h), with n equal to the product of the q(i).
1020     * <br /> In JAS F is a map.
1021     * @param n integer to factor.
1022     * @param F a map of pairs of numbers (p,e) with p**e divides n and p is
1023     *            probable prime, F is modified.
1024     */
1025    public static void factorsPollardRho(long n, SortedMap<Long, Integer> F) {
1026        long factor;
1027        long temp = n;
1028        int iterationCounter = 0;
1029        Integer count;
1030        while (!java.math.BigInteger.valueOf(temp).isProbablePrime(32)) {
1031            factor = rho(temp);
1032            if (factor == temp) {
1033                if (iterationCounter++ > 4) {
1034                    break;
1035                }
1036            } else {
1037                iterationCounter = 1;
1038            }
1039            count = F.get(factor);
1040            if (count == null) {
1041                F.put(factor, 1);
1042            } else {
1043                F.put(factor, count + 1);
1044            }
1045            temp = temp / factor;
1046        }
1047        count = F.get(temp);
1048        if (count == null) {
1049            F.put(temp, 1);
1050        } else {
1051            F.put(temp, count + 1);
1052        }
1053        //System.out.println("random = " + random.getAlgorithm());
1054    }
1055
1056
1057    /**
1058     * Search cycle with Pollards rho algorithm x**2 + c mod n. n is a positive
1059     * integer. c is a random constant.
1060     * @param n integer test.
1061     * @return x-y with gcd(x-y, n) == 1.
1062     */
1063    static long rho(long n) {
1064        long divisor;
1065        int bl = java.math.BigInteger.valueOf(n).bitLength();
1066        long c = new java.math.BigInteger(bl, random).longValue(); // .abs()
1067        long x = new java.math.BigInteger(bl, random).longValue(); // .abs()
1068        ModLongRing ring = new ModLongRing(n);
1069        ModLong cm = new ModLong(ring, c);
1070        ModLong xm = new ModLong(ring, x);
1071        ModLong xxm = xm;
1072        do {
1073            xm = xm.multiply(xm).sum(cm);
1074            xxm = xxm.multiply(xxm).sum(cm);
1075            xxm = xxm.multiply(xxm).sum(cm);
1076            divisor = gcd(xm.getVal() - xxm.getVal(), n);
1077        } while (divisor == 1L);
1078        return divisor;
1079    }
1080
1081
1082    static long gcd(long a, long b) {
1083        return BigInteger.valueOf(a).gcd(BigInteger.valueOf(b)).getVal().longValue();
1084    }
1085
1086}