From 9fcd15dba8ee78a189a7de4b3f8a06bb1b4847c8 Mon Sep 17 00:00:00 2001 From: Tibor Bizjak Date: Sun, 26 Mar 2023 17:21:44 +0200 Subject: [PATCH] Added solutions for N=120-150 --- 120_square_remainders.py | 6 ++ 121_disc_game_prize_fund.py | 18 +++++ 122_efficient_exponentiation.py | 23 ++++++ 123_prime_square_remainders.py | 20 ++++++ 124_ordered_radicals.py | 14 ++++ 125_palindromic_sums.py | 19 +++++ 126_cuboid_layers.py | 23 ++++++ 127_abc-hits.py | 46 ++++++++++++ 128_hexagonal_tile_diff.py | 42 +++++++++++ 129_repunit_divisibility.py | 18 +++++ 130_composites_with_prime_repunit_property.py | 49 +++++++++++++ 131_prime_cube_partnership.py | 25 +++++++ 132_large_repunit_factors.py | 32 +++++++++ 133_repunit_nonfactors.py | 27 +++++++ 134_prime_pair_connection.py | 41 +++++++++++ 135_same_differences.py | 19 +++++ 136_singleton_differences.py | 19 +++++ 137_fibonacci_golden_nuggets.py | 31 ++++++++ 138_special_isosceles_triangles.py | 20 ++++++ 139_pythagorean_tiles.py | 36 ++++++++++ 140_modified_fibonacci_golden_nuggets.py | 31 ++++++++ 141_progressive_numbers.py | 29 ++++++++ 142_perfect_square_collection.py | 41 +++++++++++ 143_torricelli_point.py | 58 +++++++++++++++ 144_reflections_of_laser.py | 71 +++++++++++++++++++ 145_reversible_numbers.py | 16 +++++ 146_investigating_a_prime_pattern.py | 65 +++++++++++++++++ 149_searching_for_maxsum_subsequence.py | 57 +++++++++++++++ 28 files changed, 896 insertions(+) create mode 100644 120_square_remainders.py create mode 100644 121_disc_game_prize_fund.py create mode 100644 122_efficient_exponentiation.py create mode 100644 123_prime_square_remainders.py create mode 100644 124_ordered_radicals.py create mode 100644 125_palindromic_sums.py create mode 100644 126_cuboid_layers.py create mode 100644 127_abc-hits.py create mode 100644 128_hexagonal_tile_diff.py create mode 100644 129_repunit_divisibility.py create mode 100644 130_composites_with_prime_repunit_property.py create mode 100644 131_prime_cube_partnership.py create mode 100644 132_large_repunit_factors.py create mode 100644 133_repunit_nonfactors.py create mode 100644 134_prime_pair_connection.py create mode 100644 135_same_differences.py create mode 100644 136_singleton_differences.py create mode 100644 137_fibonacci_golden_nuggets.py create mode 100644 138_special_isosceles_triangles.py create mode 100644 139_pythagorean_tiles.py create mode 100644 140_modified_fibonacci_golden_nuggets.py create mode 100644 141_progressive_numbers.py create mode 100644 142_perfect_square_collection.py create mode 100644 143_torricelli_point.py create mode 100644 144_reflections_of_laser.py create mode 100644 145_reversible_numbers.py create mode 100644 146_investigating_a_prime_pattern.py create mode 100644 149_searching_for_maxsum_subsequence.py diff --git a/120_square_remainders.py b/120_square_remainders.py new file mode 100644 index 0000000..42e234d --- /dev/null +++ b/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)) diff --git a/121_disc_game_prize_fund.py b/121_disc_game_prize_fund.py new file mode 100644 index 0000000..f2a2157 --- /dev/null +++ b/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))) diff --git a/122_efficient_exponentiation.py b/122_efficient_exponentiation.py new file mode 100644 index 0000000..01f41af --- /dev/null +++ b/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:])) diff --git a/123_prime_square_remainders.py b/123_prime_square_remainders.py new file mode 100644 index 0000000..9e8f91e --- /dev/null +++ b/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]) diff --git a/124_ordered_radicals.py b/124_ordered_radicals.py new file mode 100644 index 0000000..4c8130c --- /dev/null +++ b/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]) diff --git a/125_palindromic_sums.py b/125_palindromic_sums.py new file mode 100644 index 0000000..397fb59 --- /dev/null +++ b/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)) diff --git a/126_cuboid_layers.py b/126_cuboid_layers.py new file mode 100644 index 0000000..acb367a --- /dev/null +++ b/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)) diff --git a/127_abc-hits.py b/127_abc-hits.py new file mode 100644 index 0000000..8617cd3 --- /dev/null +++ b/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) diff --git a/128_hexagonal_tile_diff.py b/128_hexagonal_tile_diff.py new file mode 100644 index 0000000..6c55415 --- /dev/null +++ b/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)) + diff --git a/129_repunit_divisibility.py b/129_repunit_divisibility.py new file mode 100644 index 0000000..2d6db73 --- /dev/null +++ b/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) diff --git a/130_composites_with_prime_repunit_property.py b/130_composites_with_prime_repunit_property.py new file mode 100644 index 0000000..b7e22f2 --- /dev/null +++ b/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) diff --git a/131_prime_cube_partnership.py b/131_prime_cube_partnership.py new file mode 100644 index 0000000..c7b98c1 --- /dev/null +++ b/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) diff --git a/132_large_repunit_factors.py b/132_large_repunit_factors.py new file mode 100644 index 0000000..56d893e --- /dev/null +++ b/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) diff --git a/133_repunit_nonfactors.py b/133_repunit_nonfactors.py new file mode 100644 index 0000000..801f485 --- /dev/null +++ b/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) diff --git a/134_prime_pair_connection.py b/134_prime_pair_connection.py new file mode 100644 index 0000000..35213ac --- /dev/null +++ b/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) diff --git a/135_same_differences.py b/135_same_differences.py new file mode 100644 index 0000000..69652c4 --- /dev/null +++ b/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)) diff --git a/136_singleton_differences.py b/136_singleton_differences.py new file mode 100644 index 0000000..b3e06c7 --- /dev/null +++ b/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)) diff --git a/137_fibonacci_golden_nuggets.py b/137_fibonacci_golden_nuggets.py new file mode 100644 index 0000000..fc28266 --- /dev/null +++ b/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) + diff --git a/138_special_isosceles_triangles.py b/138_special_isosceles_triangles.py new file mode 100644 index 0000000..00f015c --- /dev/null +++ b/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)) diff --git a/139_pythagorean_tiles.py b/139_pythagorean_tiles.py new file mode 100644 index 0000000..8b4da1c --- /dev/null +++ b/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) diff --git a/140_modified_fibonacci_golden_nuggets.py b/140_modified_fibonacci_golden_nuggets.py new file mode 100644 index 0000000..91e0fb5 --- /dev/null +++ b/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)) diff --git a/141_progressive_numbers.py b/141_progressive_numbers.py new file mode 100644 index 0000000..3de362d --- /dev/null +++ b/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) diff --git a/142_perfect_square_collection.py b/142_perfect_square_collection.py new file mode 100644 index 0000000..ea43fb8 --- /dev/null +++ b/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) diff --git a/143_torricelli_point.py b/143_torricelli_point.py new file mode 100644 index 0000000..fe0704f --- /dev/null +++ b/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)) + diff --git a/144_reflections_of_laser.py b/144_reflections_of_laser.py new file mode 100644 index 0000000..d733724 --- /dev/null +++ b/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) diff --git a/145_reversible_numbers.py b/145_reversible_numbers.py new file mode 100644 index 0000000..eccfaed --- /dev/null +++ b/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)) diff --git a/146_investigating_a_prime_pattern.py b/146_investigating_a_prime_pattern.py new file mode 100644 index 0000000..297106d --- /dev/null +++ b/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) diff --git a/149_searching_for_maxsum_subsequence.py b/149_searching_for_maxsum_subsequence.py new file mode 100644 index 0000000..124d9b1 --- /dev/null +++ b/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))