Solve problem 94

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)
+