141 half

python/141.1.py

@@ -0,0 +1,99 @@
+from functools import reduce
+from tools import number_theory
+from tools import algebra
+import sys
+
+def next_lst(lst, std):
+    for i in range(len(lst)):
+        if lst[i] < std[i]:
+            lst[i] += 1
+            lst[:i] = [0] * i
+            return True
+    else:
+        return False
+
+def implode(p_lst, e_lst):
+    return reduce(lambda x, y: x * (y[0] ** y[1]), zip(p_lst, e_lst), 1)
+
+def factor(k):
+    p = 2
+    p_lst = []
+    e_lst = []
+    while p ** 2 <= k:
+        if not k % p:
+            e = 0
+            while not k % p:
+                e += 1
+                k //= p
+            p_lst.append(p)
+            e_lst.append(e)
+        p += 1
+    if k > 1:
+        p_lst.append(k)
+        e_lst.append(1)
+    return p_lst, e_lst
+
+def pow_factor(k):
+    p_lst, e_lst = factor(k)
+    e_lst = list(map(lambda x: x * 2, e_lst))
+    g_lst = [0] * len(e_lst)
+    result = [1]
+    while next_lst(g_lst, e_lst):
+        result.append(implode(p_lst, g_lst))
+    return list(sorted(result))
+
+def simplify(x, y):
+    d = number_theory.gcd(x, y)
+    return x // d, y // d
+
+def try_root(rest, a, b):
+    if rest < a:
+        return 0
+    q, r = divmod(rest - a, b)
+    if r:
+        return 0
+    u = algebra.is_root(q, 3)
+    if u != a:
+        return u
+    return 0
+
+def calc_k(k):
+    n = k ** 2
+    pf = pow_factor(k)
+    lf = len(pf)
+    for i in range(lf):
+        t = pf[i]
+        for j in range(i, lf - i):
+            v = pf[j]
+            rest = n // t // v
+            if rest > 0:
+                if try_root(rest, t, v):
+                    return n
+                if try_root(rest, v, t):
+                    return n
+
+def search(limit):
+    for k in range(1, int(limit ** 0.5) + 1):
+        n = calc_k(k)
+        if n:
+            yield n
+
+def debug(limit):
+    total = 0
+    k = 1
+    n = 1
+    while n < limit:
+        for r in range(1, n):
+            U, V = simplify(n - r, r ** 2)
+            u = algebra.is_root(U, 3)
+            v = algebra.is_root(V, 3)
+            if u * v:
+                total += n
+                print(k, n, r, r*U//V, r*U*U//V//V)
+        k += 1
+        n = k ** 2
+    print(total)
+
+#print(calc_k(int(sys.argv[1])))
+#print(sum(search(int(sys.argv[1]))))
+debug(100000)

python/141.md

@@ -0,0 +1,26 @@
+$$
+n=qd+r
+$$
+
+显然 $r<d$ ，因此有以下三种可能
+$$
+q<r<d \Rightarrow r^2=qd \mbox{ (a)}\\
+r<q<d \Rightarrow q^2=rd \mbox{ (b)}\\
+r<d<q \Rightarrow d^2=rq \mbox{ (c)}
+$$
+ 对于 (2.a)， $n=r^2+r$ ，$n$ 不可能是平方数。
+
+无论 (2.b) 或 (2.c)，设等比数列公比为 $e$ ，则 (1) 可写为
+$$
+n=r^2e^3+r=k^2
+$$
+
+设 $e=\frac{u}{v}$ ，其中 $u, v$ 互素且 $u>v$ 。又由于 $e^2r$ 为整数，因此必有 $v^2|r$ ，不妨设 $r=tv^2$ ，由方程 (3) 可化为
+$$
+k^2=tv^2+t^2u^3v
+$$
+显然 $tv|k$ ，对 (4) 进一步转换，得
+$$
+\frac{k^2}{tv}=v+tu^3
+$$
+对 $\frac{k^2}{tv}$ 遍历 $k^2$ 的所有因子。