Pisano period: Difference between revisions
Content added Content deleted
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
|||
Line 516: | Line 516: | ||
map pisanoPeriod [1..180] == map pisano [1..180] = True |
map pisanoPeriod [1..180] == map pisano [1..180] = True |
||
map pisanoPeriod [1..180] == map pisanoConjecture [1..180] = True</pre> |
map pisanoPeriod [1..180] == map pisanoConjecture [1..180] = True</pre> |
||
=={{header|Java}}== |
|||
Use efficient algorithm to calculate period. |
|||
<lang Java> |
|||
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; |
|||
} |
|||
} |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
Print pisano(p^2) for every prime p lower than 15 |
|||
pisano(4) = 6 |
|||
pisano(9) = 24 |
|||
pisano(25) = 100 |
|||
pisano(49) = 112 |
|||
pisano(121) = 110 |
|||
pisano(169) = 364 |
|||
Print pisano(p) for every prime p lower than 180 |
|||
pisano(2) = 3 |
|||
pisano(3) = 8 |
|||
pisano(5) = 20 |
|||
pisano(7) = 16 |
|||
pisano(11) = 10 |
|||
pisano(13) = 28 |
|||
pisano(17) = 36 |
|||
pisano(19) = 18 |
|||
pisano(23) = 48 |
|||
pisano(29) = 14 |
|||
pisano(31) = 30 |
|||
pisano(37) = 76 |
|||
pisano(41) = 40 |
|||
pisano(43) = 88 |
|||
pisano(47) = 32 |
|||
pisano(53) = 108 |
|||
pisano(59) = 58 |
|||
pisano(61) = 60 |
|||
pisano(67) = 136 |
|||
pisano(71) = 70 |
|||
pisano(73) = 148 |
|||
pisano(79) = 78 |
|||
pisano(83) = 168 |
|||
pisano(89) = 44 |
|||
pisano(97) = 196 |
|||
pisano(101) = 50 |
|||
pisano(103) = 208 |
|||
pisano(107) = 72 |
|||
pisano(109) = 108 |
|||
pisano(113) = 76 |
|||
pisano(127) = 256 |
|||
pisano(131) = 130 |
|||
pisano(137) = 276 |
|||
pisano(139) = 46 |
|||
pisano(149) = 148 |
|||
pisano(151) = 50 |
|||
pisano(157) = 316 |
|||
pisano(163) = 328 |
|||
pisano(167) = 336 |
|||
pisano(173) = 348 |
|||
pisano(179) = 178 |
|||
Print pisano(n) for every integer from 1 to 180 |
|||
1 3 8 6 20 24 16 12 24 60 |
|||
10 24 28 48 40 24 36 24 18 60 |
|||
16 30 48 24 100 84 72 48 14 120 |
|||
30 48 40 36 80 24 76 18 56 60 |
|||
40 48 88 30 120 48 32 24 112 300 |
|||
72 84 108 72 20 48 72 42 58 120 |
|||
60 30 48 96 140 120 136 36 48 240 |
|||
70 24 148 228 200 18 80 168 78 120 |
|||
216 120 168 48 180 264 56 60 44 120 |
|||
112 48 120 96 180 48 196 336 120 300 |
|||
50 72 208 84 80 108 72 72 108 60 |
|||
152 48 76 72 240 42 168 174 144 120 |
|||
110 60 40 30 500 48 256 192 88 420 |
|||
130 120 144 408 360 36 276 48 46 240 |
|||
32 210 140 24 140 444 112 228 148 600 |
|||
50 36 72 240 60 168 316 78 216 240 |
|||
48 216 328 120 40 168 336 48 364 180 |
|||
72 264 348 168 400 120 232 132 178 120 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |