import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.TreeMap;

public class PisanoPeriod {

    public static void main(String[] args) {
        System.out.printf("Print pisano(p^2) for every prime p lower than 15%n");
        for ( long i = 2 ; i < 15 ; i++ ) {
            if ( isPrime(i) ) {
                long n = i*i;
                System.out.printf("pisano(%d) = %d%n", n, pisano(n));
            }
        }

        System.out.printf("%nPrint pisano(p) for every prime p lower than 180%n");
        for ( long n = 2 ; n < 180 ; n++ ) {
            if ( isPrime(n) ) {
                System.out.printf("pisano(%d) = %d%n", n, pisano(n));
            }
        }

        System.out.printf("%nPrint pisano(n) for every integer from 1 to 180%n");
        for ( long n = 1 ; n <= 180 ; n++ ) {
            System.out.printf("%3d  ", pisano(n));
            if ( n % 10 == 0 ) {
                System.out.printf("%n");
            }
        }


    }

    private static final boolean isPrime(long test) {
        if ( test == 2 ) {
            return true;
        }
        if ( test % 2 == 0 ) {
            return false;
        }
        for ( long i = 3 ; i <= Math.sqrt(test) ; i += 2 ) {
            if ( test % i == 0 ) {
                return false;
            }
        }
        return true;
    }


    private static Map<Long,Long> PERIOD_MEMO = new HashMap<>();
    static {
        PERIOD_MEMO.put(2L, 3L);
        PERIOD_MEMO.put(3L, 8L);
        PERIOD_MEMO.put(5L, 20L);
    }

    //  See http://webspace.ship.edu/msrenault/fibonacci/fib.htm
    private static long pisano(long n) {
        if ( PERIOD_MEMO.containsKey(n) ) {
            return PERIOD_MEMO.get(n);
        }
        if ( n == 1 ) {
            return 1;
        }
        Map<Long,Long> factors = getFactors(n);

        //  Special cases
        //  pisano(2^k) = 3*n/2
        if ( factors.size() == 1 & factors.get(2L) != null && factors.get(2L) > 0 ) {
            long result = 3 * n / 2;
            PERIOD_MEMO.put(n, result);
            return result;
        }
        //  pisano(5^k) = 4*n
        if ( factors.size() == 1 & factors.get(5L) != null && factors.get(5L) > 0 ) {
            long result = 4*n;
            PERIOD_MEMO.put(n, result);
            return result;
        }
        //  pisano(2*5^k) = 6*n
        if ( factors.size() == 2 & factors.get(2L) != null && factors.get(2L) == 1 && factors.get(5L) != null && factors.get(5L) > 0 ) {
            long result = 6*n;
            PERIOD_MEMO.put(n, result);
            return result;
        }

        List<Long> primes = new ArrayList<>(factors.keySet());
        long prime = primes.get(0);
        if ( factors.size() == 1 && factors.get(prime) == 1 ) {
            List<Long> divisors = new ArrayList<>();
            if ( n % 10 == 1 || n % 10 == 9 ) {
                for ( long divisor : getDivisors(prime-1) ) {
                    if ( divisor % 2 == 0 ) {
                        divisors.add(divisor);
                    }
                }
            }
            else {
                List<Long> pPlus1Divisors = getDivisors(prime+1);
                for ( long divisor : getDivisors(2*prime+2) ) {
                    if ( !  pPlus1Divisors.contains(divisor) ) {
                        divisors.add(divisor);
                    }
                }
            }
            Collections.sort(divisors);
            for ( long divisor : divisors ) {
                if ( fibModIdentity(divisor, prime) ) {
                    PERIOD_MEMO.put(prime, divisor);
                    return divisor;
                }
            }
            throw new RuntimeException("ERROR 144: Divisor not found.");
        }
        long period = (long) Math.pow(prime, factors.get(prime)-1) * pisano(prime);
        for ( int i = 1 ; i < primes.size() ; i++ ) {
            prime = primes.get(i);
            period = lcm(period, (long) Math.pow(prime, factors.get(prime)-1) * pisano(prime));
        }
        PERIOD_MEMO.put(n, period);
        return period;
    }

    //  Use Matrix multiplication to compute Fibonacci numbers.
    private static boolean fibModIdentity(long n, long mod) {
        long aRes = 0;
        long bRes = 1;
        long cRes = 1;
        long aBase = 0;
        long bBase = 1;
        long cBase = 1;
        while ( n > 0 ) {
            if ( n % 2 == 1 ) {
                long temp1 = 0, temp2 = 0, temp3 = 0;
                if ( aRes > SQRT || aBase > SQRT || bRes > SQRT || bBase > SQRT || cBase > SQRT || cRes > SQRT ) {
                    temp1 = (multiply(aRes, aBase, mod) + multiply(bRes, bBase, mod)) % mod;
                    temp2 = (multiply(aBase, bRes, mod) + multiply(bBase, cRes, mod)) % mod;
                    temp3 = (multiply(bBase, bRes, mod) + multiply(cBase, cRes, mod)) % mod;
                }
                else {
                    temp1 = ((aRes*aBase % mod) + (bRes*bBase % mod)) % mod;
                    temp2 = ((aBase*bRes % mod) + (bBase*cRes % mod)) % mod;
                    temp3 = ((bBase*bRes % mod) + (cBase*cRes % mod)) % mod;
                }
                aRes = temp1;
                bRes = temp2;
                cRes = temp3;
            }
            n >>= 1L;
            long temp1 = 0, temp2 = 0, temp3 = 0;
            if ( aBase > SQRT || bBase > SQRT || cBase > SQRT ) {
                temp1 = (multiply(aBase, aBase, mod) + multiply(bBase, bBase, mod)) % mod;
                temp2 = (multiply(aBase, bBase, mod) + multiply(bBase, cBase, mod)) % mod;
                temp3 = (multiply(bBase, bBase, mod) + multiply(cBase, cBase, mod)) % mod;
            }
            else {
                temp1 = ((aBase*aBase % mod) + (bBase*bBase % mod)) % mod;
                temp2 = ((aBase*bBase % mod) + (bBase*cBase % mod)) % mod;
                temp3 = ((bBase*bBase % mod) + (cBase*cBase % mod)) % mod;
            }
            aBase = temp1;
            bBase = temp2;
            cBase = temp3;
        }
        return aRes % mod == 0 && bRes % mod == 1 && cRes % mod == 1;
    }

    private static final long SQRT = (long) Math.sqrt(Long.MAX_VALUE);

    //  Result is a*b % mod, without overflow.
    public static final long multiply(long a, long b, long modulus) {
        //System.out.println("    multiply : a = " + a + ", b = " + b + ", mod = " + modulus);
        long x = 0;
        long y = a % modulus;
        long t;
        while ( b > 0 ) {
            if ( b % 2 == 1 ) {
                t = x + y;
                x = (t > modulus ? t-modulus : t);
            }
            t = y << 1;
            y = (t > modulus ? t-modulus : t);
            b >>= 1;
        }
        //System.out.println("    multiply : answer = " + (x % modulus));
        return x % modulus;
    }

    private static final List<Long> getDivisors(long number) {
        List<Long> divisors = new ArrayList<>();
        long sqrt = (long) Math.sqrt(number);
        for ( long i = 1 ; i <= sqrt ; i++ ) {
            if ( number % i == 0 ) {
                divisors.add(i);
                long div = number / i;
                if ( div != i ) {
                    divisors.add(div);
                }
            }
        }
        return divisors;
    }

    public static long lcm(long a, long b) {
        return a*b/gcd(a,b);
    }

    public static long gcd(long a, long b) {
        if ( b == 0 ) {
            return a;
        }
        return gcd(b, a%b);
    }

    private static final Map<Long,Map<Long,Long>> allFactors = new TreeMap<Long,Map<Long,Long>>();
    static {
        Map<Long,Long> factors = new TreeMap<Long,Long>();
        factors.put(2L, 1L);
        allFactors.put(2L, factors);
    }

    public static Long MAX_ALL_FACTORS = 100000L;

    public static final Map<Long,Long> getFactors(Long number) {
        if ( allFactors.containsKey(number) ) {
            return allFactors.get(number);
        }
        Map<Long,Long> factors = new TreeMap<Long,Long>();
        if ( number % 2 == 0 ) {
            Map<Long,Long> factorsdDivTwo = getFactors(number/2);
            factors.putAll(factorsdDivTwo);
            factors.merge(2L, 1L, (v1, v2) -> v1 + v2);
            if ( number < MAX_ALL_FACTORS ) {
                allFactors.put(number, factors);
            }
            return factors;
        }
        boolean prime = true;
        long sqrt = (long) Math.sqrt(number);
        for ( long i = 3 ; i <= sqrt ; i += 2 ) {
            if ( number % i == 0 ) {
                prime = false;
                factors.putAll(getFactors(number/i));
                factors.merge(i, 1L, (v1, v2) -> v1 + v2);
                if ( number < MAX_ALL_FACTORS ) {
                    allFactors.put(number, factors);
                }
                return factors;
            }
        }
        if ( prime ) {
            factors.put(number, 1L);
            if ( number < MAX_ALL_FACTORS ) {
                allFactors.put(number, factors);
            }
        }
        return factors;
    }

}