project euler 223 and 224

發表於 2023-02-04 更新於 2024-10-23 分類於 Project Euler Disqus：文章字數： 3.9k 所需閱讀時間 ≈ 13 分鐘

How many barely acute triangles are there with perimeter ≤ 25,000,000?

如果Latex在解密后没有渲染，刷新网页即可, 相关Issue

题目

统计整数对$(a,b,c)$个数满足

$1\le a \le b \le c$

$a+b+c \le 25,000,000$

$a^2+b^2=c^2+1$

我的思路

$(a-1)(a+1) = (c-b)(c+b)$

$a \le \frac{25,000,000}{3}$

提前线性筛的方式处理可能的$a$的一个质因子

然后暴力枚举a, 这样可以log 级别得到 a-1和a+1的分解

然后 dfs(分解) 去凑 $c+b$和$c-b$, 注意奇偶性可以剪枝，然后$c+b>c-b$也可以剪枝， $(c+b)-(c-b) = 2b \ge 2a$也可以剪枝

这样大概是 O(n * 分解的幂次之积) 的时间复杂度

i7-7700HQ 单核python3 跑了 810.686072 s, 十多分钟有点慢

定义

PT(Pythagorean triple): $a^2+b^2=c^2$

PPT(primitive Pythagorean triple): $a^2+b^2=c^2$且$gcd(a,b,c) = 1$

APT(Almost Pythagorean triple): $a^2+b^2=c^2+1$

NPT(Nearly Pythagorean triple): $a^2+b^2=c^2-1$

Doraki

大佬说

from a triplet (a,b,c), I get the next 3 triplets

1
2
3

(2c + b - 2a, 2c + 2b - a, 3c + 2b - 2a)
(2c + b + 2a, 2c + 2b + a, 3c + 2b + 2a)
(2c - 2b + a, 2c - b + 2a, 3c - 2b + 2a)

Start from (1,1,1) and (1,2,2), be careful about duplicates when a=b, and this is it.

据说有人按照这个算法用delphi实现，只跑了3.8s

为啥就刚好是这三个啊??? 看下面的PT_matrix, 本质上是通过正向和反向矩阵建立三元组的三叉关系, 让c正向单调递增, 反向唯一(父节点) 且为正, 因此能作为根的只有”自己是自己的父节点”的节点

对于a^2+b^2=c^2+k

https://mathworld.wolfram.com/PythagoreanTriple.html

如果有 (a,b,c) 满足

那么上面的容易证明 $[1]^2+[2]^2-[3]^2 = a^2+b^2-c^2 = k$

想法是

$\begin{cases} a=np+mq \\ b=mp-nq \\ c=mp+nq \\ k = (np-mq)^2 \end{cases}$

来源也是和上面 $(a-1)(a+1) = (b-c)(b+c)$ 一样的方法，就是把来源的因数给了新的符号而已

然后配合上 Brahmagupta-fibonacci identity,

${\displaystyle {\begin{aligned}\left(a^{2}+b^{2}\right)\left(c^{2}+d^{2}\right)&{}=\left(ac-bd\right)^{2}+\left(ad+bc\right)^{2}&&(1)\\&{}=\left(ac+bd\right)^{2}+\left(ad-bc\right)^{2}.&&(2)\end{aligned}}}$

cut-the-knot pt-matrix

$a^2+b^2=c^2, a,b,c>0$

参数化$m,n,q>0$

$a=q+m$

$b=q+n$

$c=q+m+n$

带入

$(q+m)^2+(q+n)^2=(q+m+n)^2$

$q=\sqrt{2mn}$

因此选$m,n$ 让$2mn$是平方数, 就有$q$和对应的$a,b,c$

一些性质

如果$m,n$互质,那么$a,b,c$ 是PPT
$q$一直是偶
$q$大于$m$或$n$,或两个
$q$小于$m+n$

$m = c-b \ge 0$

$n = c-a \ge 0$

$q = a+b-c \ge 0$

注意到它们是成对的

如果 $q=\pm \sqrt{2mn}$ 那么会得到不同的$(a,b,c)$

因此有效$(m,n)$ 能生成两个不同的$(a,b,c)$

例如$(m,n) = (25,8), q =\pm \sqrt{2mn} = \pm 20$

PPT $(a,b,c) = (q+m,q+n,q+m+n) = 45,28,23$

S(signed)PPT $(a’,b’,c’) = (q’+m,q’+n,q’+m+n) = 5,-12,13$

$0 < c’ < c$

树结构

观察PPT,SPPT程度出现, 且PPT可以且仅可以生成3个SPPT, 通过修改a,b的符号

-	A	B	C	M	N	Q	Q’	A’	B’	C’	-
PPT	3	4	5	1	2	2	-2	-1	0	1	(later)
SPPT	-3	4	5	1	8	-4	4	5	12	13	PPT
SPPT	3	-4	5	9	2	-6	6	15	8	17	PPT
SPPT	-3	-4	5	9	8	-12	12	21	20	29	PPT

同样$(5,12,13)$ 也可以生成,

-	A	B	C	M	N	Q	Q’	A’	B’	C’	-
PPT	5	12	13	1	8	4	-4	-3	4	5	(later)
SPPT	-5	12	13	1	18	-6	6	7	24	25	PPT
SPPT	5	-12	13	25	8	-20	20	45	28	53	PPT
SPPT	-5	-12	13	25	18	-30	30	55	48	73	PPT

注意到第一行就是向回的指针

$C \to C’$ 是单调递增的

按照这个想法$(1,0,1)$ 是三叉树的根(或者说$(0,1,1)$)

-	A	B	C	M	N	Q	Q’	A’	B’	C’	-
～PPT	1	0	1	1	0	0	0	1	0	1	～PPT
～SPPT	-1	0	1	1	2	-2	2	3	4	5	PPT

去掉 $M,N,Q,Q’$ 直接看$(A,B,C), (A’,B’,C’)$ 之间的转换

也就是三个正向矩阵 $\pmatrix{1&2&2 \\ -2&-1&-2 \\ 2&2&3},\pmatrix{-1&-2&-2 \\ 2&1&2 \\ 2&2&3},\pmatrix{1&2&2 \\ 2&1&2 \\ 2&2&3}$

和一个逆向矩阵$\pmatrix{-1&-2&-2 \\ -2&-1&-2 \\ 2&2&3}$

意味着PPT可以正向逆向生成SPPT,其中正向C增加,逆向C减少

The trinary tree covers the entire set of PPTs completely and uniquely.

唯一性: 因为每个点在树上有唯一的父节点, 所以如果点在树上，它到根的路径唯一

completely: 注意到$c=m+n+q,c’=m+n+q’, q’=-q, q\ge 0, m+n > q$

所以 $0 \le c’ \le c$

相等情况需要$q = 0$ 也就是上面的$(m,n,q) = (1,0,0)$的情况,是根(这种self parent)

而$q>0$时$c=m+n+q \to c’=m+n-q$ 是单调递减且为正的

对于上面的转换矩阵

这里还给了通过$(a,b,c)\to(m,n,q) \to (m,n,-q) \to (a’,b’,c’)$ 的变化

$T=\pmatrix{-1&-2&-2 \\ -2&-1&-2 \\ 2&2&3} = \pmatrix{0&-1&1 \\ -1&0&1 \\ 1&1&-1}\cdot\pmatrix{1&0&0 \\ 0&1&0 \\ 0&0&-1}\cdot\pmatrix{1&0&1 \\ 0&1&1 \\ 1&1&1}$

显然 $(\pmatrix{0&-1&1 \\ -1&0&1 \\ 1&1&-1}\cdot\pmatrix{1&0&0 \\ 0&1&0 \\ 0&0&-1}\cdot\pmatrix{1&0&1 \\ 0&1&1 \\ 1&1&1})^2 = I$

而每次SPPT通过取绝对值变成PPT, 无非就是$\pmatrix{-1&0&0 \\ 0&1&0 \\ 0&0&1},\pmatrix{1&0&0 \\ 0&-1&0 \\ 0&0&1},\pmatrix{-1&0&0 \\ 0&-1&0 \\ 0&0&1}$ 中的一个

$T$的特征方程为$-x^3+x^2+x-1=0$

三个解为$x=1,1,-1$

对应的$(s_1,s_2,s_1+s_2)$ 和$(s_3,s_3,s_3)$

因此根有$(0,0,0),(1,0,1),(0,1,1)$

对于223,224

$a^2+b^2-c^2=k=\pm 1$ 的形式

父向 $abc\to mnq \to mnq’ \to a’b’c’ \to |a’||b’||c’|$

子向$abc\to (\pm a)(\pm b)c \to mnq \to mnq’ \to a’b’c’$

始终都有$c \ge 0, m \ge 0, n\ge 0 (a\le c, b\le c)$

当$a$或$b$ 有负数时, 显然 $q = a+b-c < 0, c’ = m+n-q > 0$ 说明子向$c$在增大

问题是$a\ge 0 , b\ge 0$时$q \ge 0$吗

$q = a+b-c = a+b-\sqrt{a^2+b^2-k} = \frac{(a+b)^2-(a^2+b^2-k)}{a+b+\sqrt{a^2+b^2-k}} = \frac{2ab+k}{a+b+\sqrt{a^2+b^2-k}}$

所以只要$a,b > 0$ 分子就有$2ab+k \ge 2\cdot 1\cdot 1 + k = 2+k > 0$ 即$q > 0$ 即父向移动时$c’$ 在减小

而$a = 0, b = 0$时就是解$a^2-c^2=k$ 可以单独解

一个问题是$c’ = m+n+q’ \ge 0$ 是否非负

即$m+n \ge q$ 是否成立

$q^2-2mn = k$

$q^2=2mn+k \le (m+n)^2 +k$ 在$m+n = 0$ 时才取等号

对于$m \ge 1$且$n \ge 1$ 时

因为 $k = \pm 1 < 2$

$q^2=2mn+k < 1+2mn+1 \le m^2+2mn+n^2 \le (m+n)^2$

因此在非边界情况$c$ 在减小

然后是找根 $(a,b,c) \ge 0, q \le 0$ 父向边 $c’ \ge c$ 的情况

上面已经有$q = \frac{2ab+k}{a+b+\sqrt{a^2+b^2-k}}$

如果要 $q\le 0$, $a = 0, k = -1$ 时， $(0,0,1)\to (2,2,3)$

P224 的情况

另一个就是$m+n < q$ 的情况, 这种无法保证$c$始终非负

就是上面讨论的当$m=0$或$n=0$时, $q^2 = 2mn+k = k = 1$, $q = 1$

$(a,b,c) = (0+1,n+1,0+n+1) = (1,b,b)$ 为根

P223 的情况

Orrin Frink

见底部链接

APT(almost pythagorean triple)

首先尝试了一些小的 (5,5,7)(4,7,8)(8,9,12) (7,11,13) (11,13,17) (10,15,18)

发现十进制下 c的末尾是2,3,7,8

2,7,12,17… 构成等差数列(差5)

8,13,18… 构成等差数列(差5)

所以搞了个(哇，咋想着靠过来的)

$x=3t+2,y=4t+1,z=5t+2$, 容易证明满足

$x=3t+1,y=4t+3,z=5t+3$, 容易证明满足

而3,4,5 本身就是 pythagorean 基础三元组， gcd(a,b,c) = 1

从三维角度看，就是(3,4,5)方向的一条线上的点

类似的思考，对于任意的 $(p,q,r)$是 APT

那么要通过它生成$(At+p,Bt+q,Ct+r)$

需要有 $A^2+B^2=C^2$

$Ap+Bq=Cr$

因为是方向向量，所以任意整数倍也是解, 所以其实就是解个二元二次方程组, 所以下面的结果可能并不互质, gcd一下就好

容易得到$(A,B,C) = (pr+q,qr-p,r^2+1)$, $(A’,B’,C’) = (pr+q,qr-p,r^2+1)$

这里我们的结论是任意合法的$(p,q,r)$ 如果延长的直线上都是合法的，且方向满足pythagorean三元组，那么方向向量有且只有两种

对于一个基础 pythagoream 三元组, $(a,b,c) = (2uv,u^2-v^2,u^2+v^2)$

令 $(x,y,z)=(at+p,bt+q,ct+r)$, $(x’,y’,z’)=(at+p’,bt+q’,ct+r’)$

其中 $p+p’=a,q+q’=b,r+r’=c$

$(x,y,z),(x’,y’,z’)$ 只要一个满足 APT,另一个也满足（根据正负对称性显然）

注意到$(p,q,r)$是合法的, 根据上面结论，$r^2 + 1$ 是 $c$的倍数(??? c的质因子性质??? TODO)

注意到$r$只用考虑在$[0,c)$ 之间取值, 因为如果不在这个范围，相当于把多的部分给t贡献就等价了

任取一个满足的r,剩下就是解关于p,q的方程ap+bq=cr, (经典的扩展欧几里德)

这里结论就是，给定a,b,c反过去找r,再找p,q 也可以得到多个解(??? 为啥说是2个 ??? TODO, c又不保证是质数, 我没有数论基础)

综上，任意(p,q,r) 对应两个 (a,b,c), 所以枚举a,b,c反找(p,q,r)可以找到所有(p,q,r)

用类似的方法，对于给定整数s, 可以解决 $x^2+y^2=z^2+s^2$

对于等腰 $a==b$, 可以用 $x^2-2y^2=-1$ 来解(pell 方程)

NPT(Nearly Pythagorean Triples) (project euler 224)

$x^2+y^2=z^2-1$

一些例子$(2,2,3) (4,8,9) (12,12,17) (8,32,33) (10,50,51) (22,46,51) (34,38,51)$

上面办法并不能用了

THEOREM: 如果$(p,q,r)$ 是APT, 那么$(2pr,2qr,2r^2+1)$ 是NPT, 反过来如果$(p,q,r)$是NPT,那么$(2p^2+1,2pq,2pr)$是APT

易证

这样也得到无限多NPT, 但并不是所有NPT能这样得到

其实也是在解 $x^2+y^2+w^2 = z^2$ 其中$w=1$

Catalan 1885年证明了解是

$x=2(pr+qs),y=2(qr-ps),w=p^2+q^2-r^2-s^2,z=p^2+q^2+r^2+s^2$

那么当 $1=w=p^2+q^2-r^2-s^2$, 令$r^2+s^2 = n$ 则$p^2+q^2=n+1$

$x=2(pr+qs),y=2(qr-ps),z=2n+1$

因此需要找 $n,n+1$且都能表示成平方和，就能得到$x,y,z$

A Note on Generateing APT

tldr; 非一般情况，没啥营养

举了个从$(3,4,5)$开始的实际例子

但对于大的PPT, 难以找$(p,q,r)$

?不太懂为什么原文要搞6个，很多是附属的方程，两个就够了

不妨设选定的PPT $(a,b,c) = (m^2-n^2,2mn,m^2+n^2)$

令 $m=i+1,n=i$

$(a,b,c) = (2i-1,2i^2-2i,2i^2-2i+1)$

要解方程组

$p^2+q^2=r^2=1$

$((2i-1)+p)^2+((2i^2-2i)+q)^2=((2i^2-2i+1)+r)^2=1$, 这里虽然没有t, 但是本质都是要 ap+bq=cr

即

$p^2+q^2=r^2=1$

$(2i-1)p+(2i^2-2i)q=(2i^2-2i+1)r$

显然有一个 q=1,p=r=2i-1 的解

对应的$(p’,q’,r’) = (2i-2,2i^2-4i+1,2i^2-4i+2)$

结论就是特殊的PPT的一个APT可以快速得到

$x^2+y^2+z^2=u^2$

Ayoub

所有除以u

$x’=\frac{x}{u},y’=\frac{y}{u},z’=\frac{z}{u}$

即$x’^2+y’^2+z’^2=1$, 解有理数方程

然后就和 pythagorean的经典方法一样，这就是个半径为1的球面，上的有理数点一样的想法

如果 $x’,y’,z’$ 其中有0, 那么就是勾股定理了 $(\frac{m^2-n^2}{m^2+n^2},\frac{2mn}{m^2+n^2},0)$

考虑经过这个点的一条线 $(x’,y’,z’) = (\frac{m^2-n^2}{m^2+n^2}+rt,\frac{2mn}{m^2+n^2}+pt,0+qt)$,其中方向为$(r,p,q)$为整数

那么和球面交就要满足$(\frac{m^2-n^2}{m^2+n^2}+rt)^2+(\frac{2mn}{m^2+n^2}+pt)^2+(qt)^2=1$

那么它是关于 t的二次方程

t = 0 是一个解

$t=\frac{-2r(m^2-n^2)-4mnp}{(m^2+n^2)(p^2+q^2+r^2)}$

所以直接无脑把分母丢给$u$, 则有

$x = (m^2-n^2)(p^2+q^2-r^2) - 4mnpr$

$y = 2mn(r^2-p^2+q^2)-2rp(m^2-n^2)$

$z = -2qr(m^2-n^2) - 4mnpq$

$u = (m^2+n^2)(p^2+q^2+r^2)$

那么文中提到的 V. A. Lebesgue, U. Bainelli, C. Gill, J. A. Euler, Japanese Matsunango, P. Cossali 的6个方程，不过是这里的特例

但是这个方程不像勾股定理的方程，它无法直接给出所有解，这里直接的意思是比如u=27

而 22^2+7^2+14^2 = 27^2 不能找到m,n,p,q,r 让u=27且得到上面的式子

def xyzu(m,n,p,q,r):
    x = (m*m-n*n)*(p*p+q*q-r*r) - 4*m*n*p*r
    y = 2*m*n*(r*r-p*p+q*q) - 2 *r *p*(m*m-n*n)
    z = -2*q*r*(m*m-n*n) - 4*m*n*p*q
    u = (m*m+n*n)*(p*p+q*q+r*r)

    print(f"x={x}\ty={y}\tz={z}\tu={u}")

for m in range(0,27):
    for n in range(m,27+1):
        m2n2=m*m+n*n
        if n == 0 or 27 % m2n2 != 0:
            continue
        for p in range(0,27+1):
            for q in range(p,27+1):
                if m2n2*(p*p+q*q) > 27:
                    break
                for r in range(q,27+1):
                    if m2n2*(p*p+q*q+r*r) > 27:
                        break
                    if m2n2*(p*p+q*q+r*r) == 27:
                        xyzu(m, n, p, q, r)
                        xyzu(m, n, p, q,-r)
                        xyzu(m, n, p,-q, r)
                        xyzu(m, n, p,-q,-r)
                        xyzu(m, n,-p, q, r)
                        xyzu(m, n,-p, q,-r)
                        xyzu(m, n,-p,-q, r)
                        xyzu(m, n,-p,-q,-r)
                        xyzu(m,-n, p, q, r)
                        xyzu(m,-n, p, q,-r)
                        xyzu(m,-n, p,-q, r)
                        xyzu(m,-n, p,-q,-r)
                        xyzu(m,-n,-p, q, r)
                        xyzu(m,-n,-p, q,-r)
                        xyzu(m,-n,-p,-q, r)
                        xyzu(m,-n,-p,-q,-r)