Added solutions for N=120-150

120_square_remainders.py

@@ -0,0 +1,6 @@
+#!/usr/bin/env python3
+
+def main(n):
+    return sum(a*(a-2+a%2) for a in range(3, n+1))
+
+print(main(1000))

121_disc_game_prize_fund.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python3
+
+from fractions import Fraction
+
+def search(turns, wins=0, discs=2, cache={}):
+    if discs - 2  == turns:
+        return 2* wins > turns
+    if  turns // 2 - wins > turns - discs + 1:
+        return 0
+    if (discs, wins) in cache:
+        return cache[(discs, wins)]
+    r = Fraction(1, discs) * search(turns, wins+1, discs+1)
+    r += Fraction(discs-1, discs) * search(turns, wins, discs+1)
+    cache[(discs, wins)] = r
+    return r
+
+N = 15
+print(int(search(N)**(-1)))

122_efficient_exponentiation.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+
+from itertools import combinations_with_replacement as comb
+
+bound = 200
+sieve = [i-1 for i in range(bound+1)]
+dup = []
+
+def rec(chain, sieve, bound, id=0):
+    c = 1
+    for m in chain:
+        n = m+chain[0]
+        if n > bound:
+            continue
+        if len(chain) >= max(sieve[n:n+3]): #<-- trial and error magic
+            continue
+        if len(chain) < sieve[n]:
+            sieve[n] = len(chain)
+        c += rec([n] + chain, sieve, bound, id+1)
+    return c
+
+rec([1], sieve, bound)
+print(sum(sieve[2:]))

123_prime_square_remainders.py

@@ -0,0 +1,20 @@
+#!/usr/bin/env python3
+
+def main(m, lim=10**7):
+    for n, p in enumerate(primes(lim), 1):
+        r = ((2*n)%p)*p if n%2 else 2
+        if r > m:
+            return n, p, r
+    return False
+
+def primes(n):
+    sieve = [True]*n
+    for p in range(2, len(sieve)):
+        if not sieve[p]:
+            continue
+        for i in range(2*p, len(sieve), p):
+            sieve[i] = False
+        yield p
+
+#print(main(10**9, 10**6))
+print(main(10**10, 10**6)[0])

124_ordered_radicals.py

@@ -0,0 +1,14 @@
+#!/usr/bin/env python3
+
+def radicals(n):
+    sieve = [1]*(n+1)
+    for p in range(2,n+1):
+        if sieve[p] > 1:
+            continue
+        for m in range(p, n+1, p):
+            sieve[m] *= p
+    return zip(range(n+1), sieve)
+
+rads = list(radicals(10**5))
+rads.sort(key = lambda x: x[1])
+print(rads[10**4][0])

125_palindromic_sums.py

@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+
+from math import sqrt
+
+def is_palindrome(n):
+    return str(n) == str(n)[::-1]
+
+lim = 10**8
+nums = []
+
+for n in range(2, int(sqrt(lim))+1):
+    s = (n-1)**2 + n**2
+    while s < lim:
+        if is_palindrome(s) and s not in nums:
+            nums.append(s)
+        n += 1
+        s += n**2
+
+print(sum(nums))

126_cuboid_layers.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+
+def fst_layer(a, b, c):
+    return 2*(a*b + a*c + b*c)
+
+def abc(n):
+    for a in range(1, int((n//6)**0.5)+1):
+        for b in range(a, int((a**2 + n//2)**0.5) - a):
+            for c in range(b, (n - 2*a*b) // (2*a + 2*b)):
+                yield a,b,c
+
+n = 2*10**4
+sieve = [0]*n
+
+for a,b,c in abc(n):
+    m = 0
+    lyn = fst_layer(a, b, c)
+    while lyn < n:
+        sieve[lyn] += 1
+        lyn += 4*(2*m + a + b + c)
+        m += 1
+
+print(sieve.index(1000))

127_abc-hits.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python3
+
+from math import gcd
+from collections import defaultdict
+
+def getreps(n):
+    sieve = [1] * n
+    yield 1, 1
+    for p in range(2, n):
+        if sieve[p] == 1:
+            for i in range(p, n, p):
+                sieve[i] *= p
+        yield p, sieve[p]
+
+
+N = 120000
+
+repunits = list(getreps(N))
+sortedreps = sorted(set(list(zip(*repunits))[1]))
+repmap = defaultdict(list)
+for n, rep in repunits:
+    repmap[rep].append(n)
+
+count = 0
+s = 0
+
+for c, repC in repunits:
+    root = c // repC
+    if root <= 2:
+        continue
+    for repA in sortedreps:
+        if repA >= root:
+            break
+        if gcd(repA, repC) != 1:
+            continue
+        for a in repmap[repA]:
+            if 2*a >= c:
+                break
+            b = c - a
+            repB = repunits[b-1][1]
+            if repA*repB >= root:
+                continue
+            count += 1
+            s += c
+
+print(s)

128_hexagonal_tile_diff.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python3
+
+def primes(n):
+    sieve = [True] * n
+    for p in range(2, int(n**0.5)):
+        if not sieve[p]:
+            continue
+        for i in range(p**2, n, p):
+            sieve[i] = False
+    return [False, False] + sieve[2:]
+
+def is_prime(p):
+    if p >= len(sieve):
+        raise LimitError
+    return sieve[p]
+
+def getnum(n, c):
+    m = n//6
+    if c:
+        return 3*(m-1)*m + 2
+    return 3*m*(m+1) + 1
+
+sieve = primes(10**6)
+target = 2000
+
+c = 2
+n = 6
+
+while c < target:
+    n += 6
+    if not is_prime(n-1):
+        continue
+    c += is_prime(n+1) and is_prime(2*n + 5)
+    if c == target:
+        flag = True
+        break
+    c += is_prime(n+5) and is_prime(2*n-7)
+else:
+    flag = False
+
+print(getnum(n, flag))
+

129_repunit_divisibility.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python3
+
+def minrep(n):
+    base = 10 % (9*n)
+    m = base
+    c = 1
+    while m != 1:
+        m = (m * base) % (9*n)
+        c += 1
+    return c
+
+target = 10**6
+n = target + (not target%2)
+
+while n%5 == 0 or minrep(n) < target:
+    n += 2
+
+print(n)

130_composites_with_prime_repunit_property.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python3
+
+def divisors(divs):
+    r = []
+    for p in divs:
+        r += [p*d for d in r] + [p]
+    return r
+
+def add_divisors(n, st=2):
+    for i in range(st*n, len(sieve), st*n):
+        sieve[i].append(n)
+    pw = 2
+    while n**pw < len(sieve):
+        m = n**pw
+        for i in range(st*m, len(sieve), st*m):
+            sieve[i][-1] *= n
+        pw += 1
+
+target = 25
+count = 0
+s = 0
+
+lim = 10**5
+sieve = [x for i in range(lim//2) for x in ([], True)]
+
+add_divisors(2, st=1)
+
+for p in range(3, len(sieve)):
+    if not p%2 and not sieve[p+1]:
+        for d in divisors(sieve[p]):
+            if pow(10, d, 9*p+9) == 1:
+                break
+        else:
+            continue
+        count += 1
+        s += p + 1
+        #print(count, p+1)
+        if count >= target:
+            break
+    if not p%2 or not sieve[p]:
+        continue
+    for i in range(p**2, len(sieve), 2*p):
+         sieve[i] = False
+    add_divisors(p)
+
+if count < target:
+    print("Limit reached. Count = {}/{}".format(count,target))
+
+print(s)

131_prime_cube_partnership.py

@@ -0,0 +1,25 @@
+#!/usr/bin/env python3
+
+def primes(n):
+    sieve = [False] * n
+    for p in range(2, n):
+        if sieve[p]:
+            continue
+        for i in range(p**2, n, p):
+            sieve[i] = True
+    return sieve
+
+bound = 10**6
+sieve = primes(bound)
+c = 0
+
+n = 1
+cube_diff = 7
+
+while cube_diff < bound:
+    if not sieve[cube_diff]:
+        c += 1
+    n += 1
+    cube_diff = 3*n**2 + 3*n + 1
+
+print(c)

132_large_repunit_factors.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+
+def primes(n):
+    sieve = [True] * n
+    for p in range(2, int(n**0.5)):
+        if not sieve[p]:
+            continue
+        for i in range(p**2, n, p):
+            sieve[i] = False
+    return [i for i, p in enumerate(sieve) if p][2:]
+
+target = 40
+repw = 10**9
+lim = 10**6
+nums = primes(lim)
+c = 0
+s = 0
+
+for p in nums[3:]:
+    if p%4 != 1 and p%5 != 1:
+        continue
+    if pow(10, repw, 9*p) == 1:
+        c += 1
+        s += p
+        #print c, p
+        if c == target:
+            break
+
+if c < target:
+    print("Limit reached. Count = {}/{}".format(c, target))
+
+print(s)

133_repunit_nonfactors.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+def add_divisors(n):
+    for k in range(1, (len(sieve) - 1) // n):
+        sieve[k*n] *= n * (sieve[k] // sieve[k*n])
+
+def primes(n):
+    for p in range(3, int(n**0.5), 2):
+        if sieve[p] != 1:
+            continue
+        for i in range(p**2, n, 2*p):
+            sieve[i] = 0
+
+n = 10**5
+sieve = [1] * n
+
+add_divisors(2)
+add_divisors(5)
+primes(n)
+
+s = 0
+
+for p in range(7, n, 2):
+    if sieve[p] == 1 and pow(10, sieve[p-1], 9*p) != 1:
+        s += p
+
+print(s + 2 + 3 + 5)

134_prime_pair_connection.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+def primes(n):
+    sieve = [False]*n
+    for p in range(2, int(n**0.5)):
+        if sieve[p]:
+            continue
+        for i in range(p**2, n, p):
+            sieve[i] = True
+    return [p for p, t in enumerate(sieve) if not t][2:]
+
+def solve(p1, p2, k):
+    t = extended_euclid(p2, k)
+    return ((t * p1) % k) * p2
+
+def extended_euclid(p2, k):
+    a = k
+    b = p2
+    t1, t2 = 0, 1
+    while b != 0:
+        k = a // b
+        a, b = b, a % b
+        t1, t2 = t2, t1 - k*t2
+    return t1
+
+
+n = 10**6
+
+
+nums = primes(n + 1000)
+nums = [nums[i] for i in range(1, len(nums)) if nums[i-1] <= n]
+
+k = 10
+s = 0
+
+for p1, p2 in zip(nums[1:], nums[2:]):
+    if p1 > k:
+        k *= 10
+    s += solve(p1, p2, k)
+
+print(s)

135_same_differences.py

@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+
+N = 10 ** 6
+target = 10
+
+sieve = [0] * N
+
+m = 3
+while 2*m - 1 < 3*N:
+    n = m - 1 - (not m % 2)
+    i = (m - n) * (m + n)
+    while i < 3*N and n > 0:
+        if not i % 3 and not (m - (n//2)) % 3:
+            sieve[i//3] += 1
+        n -= 2
+        i = (m - n) * (m + n)
+    m += 1
+
+print(sieve.count(target))

136_singleton_differences.py

@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+
+N = 50 * 10 ** 6
+target = 1
+
+sieve = [0] * N
+
+m = 3
+while 2*m - 1 < 3*N:
+    n = m - 1 - (not m % 2)
+    i = (m - n) * (m + n)
+    while i < 3*N and n > 0:
+        if not i % 3 and not (m - (n//2)) % 3:
+            sieve[i//3] += 1
+        n -= 2
+        i = (m - n) * (m + n)
+    m += 1
+
+print(sieve.count(target))

137_fibonacci_golden_nuggets.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+
+def sols():
+    """yields int solutions of a**2 - 5 b**2 = 1"""
+    a, b = 9, 4
+    while True:
+        yield a, b
+        a, b = 9*a + 20*b, 4*a + 9*b
+
+def Cs():
+    """yields int solutions of a in a**2 - 5 b**2 = -4"""
+    for a, b in sols():
+        yield 25*b - 11*a
+        yield 10*b - 4*a
+        yield 5*b - a
+
+gen = Cs()
+next(gen)
+count = 0
+
+N = 15
+
+while count < N:
+    c = next(gen)
+    if (c-1) % 5 != 0:
+        continue
+    count += 1
+
+n = (c-1) // 5
+print(n)
+

138_special_isosceles_triangles.py

@@ -0,0 +1,20 @@
+#!/usr/bin/env python3
+
+def sols(n):
+    """yields int solutions of a**2 - 5 b**2 = 1"""
+    a, b = 9, 4
+    for i in range(n):
+        yield a, b
+        a, b = 9*a + 20*b, 4*a + 9*b
+
+N = 12
+
+Ls = []
+
+for a, b in sols(N // 2):
+    a0, b0 = 5*b - 2*a, a - 2*b
+    m0, n0 = 2*b0 + a0, b0
+    m, n   = 2*b  + a,  b
+    Ls += [m0**2 + n0**2, m**2 + n**2]
+
+print(sum(Ls))

139_pythagorean_tiles.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+
+from math import gcd
+
+def sols():
+    """yields solutions to a**2 - 2b**2 = -1"""
+    a, b = 3, 2
+    while True:
+        yield 2*b - a, a - b
+        a, b = 3*a + 4*b, 2*a + 3*b
+
+
+def mn_from_ab(a, b, sgn):
+    M, N = a + b, b - sgn
+    div = gcd(M, N)
+    m, n = M // div, N // div
+    return m, n
+
+def gen_mn(N):
+    for a, b in sols():
+        if b > N//2:
+            break
+        yield mn_from_ab(a, b, 1)
+        yield mn_from_ab(a, b, -1)
+
+def is_primitive(m, n):
+    return (m + n) % 2 == 1 and n != 0
+
+def count_triangles(m, n, circ):
+    sides = (m**2 - n**2, 2*m*n, m**2 + n**2)
+    return (circ-1) // sum(sides)
+
+N = 100 * 10**6
+
+R = sum(count_triangles(m, n, N) for m, n in gen_mn(N) if is_primitive(m, n))
+print(R)

140_modified_fibonacci_golden_nuggets.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+
+def sols():
+    """yields int solutions of a**2 - 5 b**2 = 1"""
+    a, b = 9, 4
+    while True:
+        yield a, b
+        a, b = 9*a + 20*b, 4*a + 9*b
+
+def cs():
+    """yields int solutions of a in a**2 - 5 b**2 = 44"""
+    base_sols = [(7, 1), (8, 2), (13, 5), (17, 7), (32, 14), (43, 19)]
+    base_sols.reverse()
+    for a, b in sols():
+        for a0, b0 in base_sols:
+            yield a0 * a - 5 * b0 * b
+
+
+N = 30
+gen = cs()
+next(gen)
+
+nuggs = []
+
+while len(nuggs) < N:
+    c = next(gen)
+    if (c - 7) % 5 != 0:
+        continue
+    nuggs.append((c-7) // 5)
+
+print(sum(nuggs))

141_progressive_numbers.py

@@ -0,0 +1,29 @@
+#!/usr/bin/env python3
+
+from math import gcd
+
+def is_square(n):
+    return int(n**0.5)**2 == n
+
+def sol(p, q, n):
+    return n*q * (p**3 * n + q)
+
+N = 10**12
+a = 1
+r = 0
+
+for p in range(2, int(N**(1.0/3))):
+    for q in range(1, p):
+        if p**3 * q + q**2 >= N:
+            break
+        if gcd(p, q) > 1:
+            continue
+        n = 1
+        s = sol(p, q, n)
+        while s < N:
+            if is_square(s):
+                r += s
+            n += 1
+            s = sol(p, q, n)
+
+print(r)

142_perfect_square_collection.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+from math import sqrt
+from itertools import combinations
+
+def is_square(n):
+    return int(sqrt(n))**2 == n
+
+zbound = 10**6
+
+zs = [[] for i in range(zbound)]
+
+c = 1
+while c - 2 < len(zs):
+    d = c-2
+    while d > 0 and (c**2 - d**2) // 2 < len(zs):
+        zs[(c**2 - d**2) // 2].append((c**2 + d**2) // 2)
+        d -= 2
+    c += 1
+
+max_z = None
+min_sum = None
+for z in range(1, len(zs)):
+    if max_z != None and z > max_z:
+        #print("found")
+        break
+    nums = sorted(zs[z])[::-1]
+    if len(nums) < 2:
+        continue
+
+    for cd, ef in combinations(nums, 2):
+        if is_square(cd + ef) and is_square(cd - ef):
+            x = cd
+            y = ef
+            #print(x, y, z)
+            if min_sum == None or min_sum > x + y + z:
+                min_sum = x + y + z
+            if max_z == None or x < max_z:
+                max_z = x
+
+print(min_sum)

143_torricelli_point.py

@@ -0,0 +1,58 @@
+#!/usr/bin/env python3
+
+from collections import defaultdict
+from math import gcd
+
+def sols(N):
+    """yields solutions of x**2 + 3y**2 = z**2 where x + y <= 2N"""
+    for l in range(1, N+1, 2):
+        for k in range(1, N//l + 1, 2):
+            if gcd(l, k) != 1:
+                continue
+            x = abs(3*l**2 - k**2) // 2
+            y = l * k
+            z = (3*l**2 + k**2) // 2
+            for n in range(1, (2*N) // (x+y) + 1):
+                yield (n*x, n*y, n*z)
+
+    for l in range(1, N//2 + 1):
+        for k in range(1 + (l%2), N//(2*l) + 1, 2):
+            if gcd(l, k) != 1:
+                continue
+            x = 2 * abs(3*l**2 - k**2)
+            y = 4 * l * k
+            z = 2 * (3*l**2 + k**2)
+            if x + y > 2*N:
+                continue
+            for n in range(1, (2*N) // (x+y) + 1):
+                yield (n*x, n*y, n*z)
+
+def pq_gen(N):
+    """yields solutions of c**2 = p**2 + q**2 + pq where c is int and p + q <= N"""
+    for x, y, z in sols(N):
+        p = y
+        if x <= p or (x - p) % 2 != 0:
+            continue
+        q = (x - p) // 2
+        yield p, q
+
+
+def find_pqr(N):
+    pqs = defaultdict(lambda : set())
+    for p, q in pq_gen(N-1):
+        pqs[p].add(q)
+        pqs[q].add(p)
+    sums = set()
+    for p, qrs in dict(pqs).items():
+        for q in qrs:
+            for r in qrs:
+                if p + q + r > N or r not in pqs[q]:
+                    continue
+                sums.add(p + q + r)
+    return sum(sums)
+
+
+N = 120 * 1000
+
+print(find_pqr(N))
+

144_reflections_of_laser.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+
+from math import sqrt, pi
+import matplotlib.pyplot as plt
+import numpy as np
+
+class Ellipse:
+    def __init__(self, a, b):
+        self.a = float(a)
+        self.b = float(b)
+
+    def normal(self, vec):
+        x, y = vec
+        return norm([self.b**2 * x, self.a**2 * y])
+
+def quadratic(a, b, c):
+    D = b**2 - 4*a*c
+    if D < 0:
+        return False
+    r1 = (-b + sqrt(D)) / (2*a)
+    r2 = (-b - sqrt(D)) / (2*a)
+    return r1, r2
+
+def norm(v):
+    l = sqrt(sum(x**2 for x in v))
+    return np.array([x/l for x in v])
+
+def intersections(ellipse, point, dir):
+    m = dir[1] / dir[0]
+    yint = point[1] - point[0]*m
+
+    a = ellipse.a**2 * m**2 + ellipse.b**2
+    b = 2 * ellipse.a**2 * m * yint
+    c = ellipse.a**2 * (yint**2 - ellipse.b**2)
+
+    x1, x2 = quadratic(a, b, c)
+    y1, y2 = m*x1 + yint, m*x2 + yint
+
+    return np.array([x1, y1]), np.array([x2, y2])
+
+def reflect(ellipse, A, B):
+    point = B
+    dir = B - A
+    n = ellipse.normal(point)
+    ref = 2 * (n.dot(dir) * n - dir) + dir
+    ps = intersections(ellipse, point, ref)
+    return max(ps, key=lambda v: (B-v).dot(B-v))
+
+def plot_ellipse(ellipse):
+    t = np.linspace(0, 2*pi, 100)
+    plt.plot(ellipse.a*np.cos(t), ellipse.b*np.sin(t))
+
+ellipse = Ellipse(5, 10)
+A = np.array([0.0, 10.1])
+B = np.array([1.4, -9.6])
+x = 0.01
+
+
+plt.axis('equal')
+plot_ellipse(ellipse)
+
+count = 0
+
+while not (-x <= B[0] <= x and B[1] > 0):
+    plt.plot(*zip(A, B))
+    A, B = B, reflect(ellipse, A, B)
+    count += 1
+
+
+plt.show()
+print(count)

145_reversible_numbers.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python3
+
+# closed form expression, d = digits
+# if d is even: 20 * 30**(d/2 - 1)
+# if d = 3+4n : 5 * 20**((d+1)/4) * 25**((i-3)/4)
+# iterate all digits in range
+
+def main(d): #inclusive
+    r = 0
+    for i in range(2,d+1,2):
+        r += 20*30**(i//2 - 1)
+    for i in range(3,d+1,4):
+        r += 5 * 20**((i+1)//4) * 25**((i-3)//4)
+    return r
+
+print(main(9))

146_investigating_a_prime_pattern.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+
+from random import randint
+from operator import mul
+from functools import reduce
+
+DEPTH = 5
+SEQ = [1, 3, 7, 9, 13, 27]
+NOTP = [5, 11, 15, 17, 19, 21, 23, 25]
+
+
+def is_prime(n, depth=DEPTH):
+    r = 0
+    d = n - 1
+    while not d % 2:
+        r += 1
+        d //= 2
+    for k in range(depth):
+        a = randint(2, n - 2)
+        x = pow(a, d, n)
+        if x == 1 or x == n - 1:
+            continue
+        for i in range(r-1):
+            x = pow(x, 2, n)
+            if x == n-1:
+                break
+        else:
+            return False
+    return True
+
+def test(m, nums):
+    return all((m**2 + r) % n for n in nums for r in SEQ)
+
+def getmod(nums):
+    return [m for m in range(reduce(mul, nums)) if test(m, nums)]
+
+def testseq(n):
+    for s in SEQ:
+        if not is_prime(n**2 + s):
+            return False
+    for s in NOTP:
+        if is_prime(n**2 + s):
+            return False
+    return True
+
+
+nums = [2, 3, 5, 7, 11, 13, 17]
+
+N = 150 * 10**6
+
+n = 0
+mods = getmod(nums)
+step = reduce(mul, nums)
+
+s = 0
+
+while n < N:
+    for m in mods:
+        if n + m >= N:
+            break
+        if testseq(n+m):
+            s += n + m
+    n += step
+
+print(s)

149_searching_for_maxsum_subsequence.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+
+def search_row(row):
+    row = [a for a in row]
+    for i in range(1, len(row)):
+        prev, val = row[i-1], row[i]
+        row[i] = max(prev + val, val)
+    return max(row)
+
+def search_rows(grid):
+    return max(search_row(row) for row in grid)
+
+def search_diagonals(grid):
+    ln = len(grid)
+    m = 0
+    for d in range(-ln + 1, ln):
+        r = search_row([grid[i][i+d]
+                        for i in range(ln) if 0 <= i+d < ln])
+        m = max(m, r)
+    return m
+
+def fib_gen(k, mem):
+    if len(mem) > k and mem[k] != None:
+        return mem[k]
+    elif k < 56:
+        s = 100003 - 200003*k + 300007*k**3
+        s %= 1000000
+        s -= 500000
+    else:
+        s = (fib_gen(k-24, mem) + fib_gen(k-55, mem) + 1000000) % 1000000
+        s -= 500000
+    mem += [None] * (k - len(mem) + 1)
+    mem[k] = s
+    return s
+
+
+grid2 = [[-2, 5, 3, 2],
+        [9, -6, 5, 1],
+        [3, 2, 7, 3],
+        [-1, 8, -4, 8]]
+
+mem = []
+
+side = 2000
+
+grid = [[0]*side for i in range(side)]
+
+for ri in range(side):
+    for ci in range(side):
+        grid[ri][ci] = fib_gen(side*ri + ci + 1, mem)
+
+rows = search_rows(grid)
+cols = search_rows(zip(*grid))
+dia1 = search_diagonals(grid)
+dia2 = search_diagonals([row[::-1] for row in grid])
+
+print(max(rows, cols, dia1, dia2))