diff --git a/python/e094.py b/python/e094.py new file mode 100644 index 0000000..2077a23 --- /dev/null +++ b/python/e094.py @@ -0,0 +1,94 @@ + +def odd(n): + return n % 2 == 1 + +def is_square(n): + """ From the internet. I should implement my own version of this. """ + ## Trivial checks + if type(n) != int: ## integer + return False + if n < 0: ## positivity + return False + if n == 0: ## 0 pass + return True + + ## Reduction by powers of 4 with bit-logic + while n&3 == 0: + n=n>>2 + + ## Simple bit-logic test. All perfect squares, in binary, + ## end in 001, when powers of 4 are factored out. + if n&7 != 1: + return False + + if n==1: + return True ## is power of 4, or even power of 2 + + ## Simple modulo equivalency test + c = n%10 + if c in {3, 7}: + return False ## Not 1,4,5,6,9 in mod 10 + if n % 7 in {3, 5, 6}: + return False ## Not 1,2,4 mod 7 + if n % 9 in {2,3,5,6,8}: + return False + if n % 13 in {2,5,6,7,8,11}: + return False + + ## Other patterns + if c == 5: ## if it ends in a 5 + if (n//10)%10 != 2: + return False ## then it must end in 25 + if (n//100)%10 not in {0,2,6}: + return False ## and in 025, 225, or 625 + if (n//100)%10 == 6: + if (n//1000)%10 not in {0,5}: + return False ## that is, 0625 or 5625 + else: + if (n//10)%4 != 0: + return False ## (4k)*10 + (1,9) + + ## Babylonian Algorithm. Finding the integer square root. + ## Root extraction. + s = (len(str(n))-1) // 2 + x = (10**s) * 4 + A = {x, n} + while x * x != n: + x = (x + (n // x)) >> 1 + if x in A: + return False + A.add(x) + return True + + +def euler_094(): + n_max = 10 ** 9 + n, total_peri = 2, 0 + while n < n_max: + n += 2 + for b2, c in [(n, n + 1), (n, n - 1)]: + b = b2 // 2 + a_squared = c * c - b * b + + if not is_square(a_squared): + continue + + peri = 2 * c + b2 + # print("b2={} c={} peri={}".format(b2, c, peri)) + + if peri > n_max: + return total_peri + + total_peri += peri + + if b2 > 10000: + n = int(n * 3.725) - 5 + if odd(n): + n -= 1 + + +if __name__ == "__main__": + solution = euler_094() + print("e094.py: " + str(solution)) + assert(solution == 518408346) +