佩尔方程实例讲解

hresh 755 0

佩尔方程实例讲解

前言

佩尔方程是一种不定二次方程。
下面的不定方程称为佩尔(Pell)方程:

x^2-d*y^2= 1 ........①

其中 d 为正整数,若 d 是完全平方数,则这个方程式只有平凡解(1,-1,0)。若 d 是非平方数。佩尔方程存在无穷多个解。

若佩尔方程的最小特解(最小正整数解)是(x1,y1),那么可有迭代公式:

佩尔方程

除了上述常见的佩尔方程,还有一种形式。下面的不定方程称为第II型佩尔(Pell)方程:

x^2-d*y^2= -1 ........②

如果②有正整数解,设(a,b)是②的正整数解中使 x+y√d 最小的解(称(a,b)为②的基本解),则②的全部正整数解可以表示为:

x+y√d=(a+b√d)^(2n+1) (n为任意正整数)

而且记x0+y0√d=(a+b√d)^2,则(x0,y0)为①的基本解。

关于公式②的通解,可以前往阅读数论笔记 · 佩尔方程,该文章对此进行了详细讲解。这里我们就直接拿来使用。

佩尔方程

如上图所示,当 K 为 -1 时,即和公式②一致,此时 p,q 为公式②的解,(r,s)为公式①的基本解。将公式①和公式②进行相乘,最终得到(x,y)的通解。假设 D 为 2,公式①的基本解为(3,2),则(x,y) = (3p±4q,2p±3q)。

实例分析

案例一

记得有一次全班去唱K, 其中有个活动是情歌对唱. 具体操作流程是这样的:
准备好 21 个阄(我们班 15 男 6 女), 其中只有两个是有标记的, 每人随意抓取一个, 最后取到有标记的阄的两个人去点首情歌对唱.
旁边一哥们儿幽幽地对我说, 看来搅基真是神的安排啊, 你看我们班的男女人数, 搅基的几率 C(15,2)/C(21,2) 刚好是 1/2.
给跪了, 这哥们儿对数字太敏感了, 简直是拉马努金转世啊. 不过我随之想到一个问题: (21, 15) 真的是神的唯一安排吗? 其实不是的,
神还有很多类似的安排. 比如 (4, 3), 显然 C(4,2)/C(3,2) 也等于 1/2, 当然还有 (120, 85) 等等等等.
神的安排太多太多了, 如果我们定义 (n, m) 是一个安排(其中 1 < m < n), 而如果 C(m,2)/C(n,2) = 1/2, 它就是神的安排.
现在的问题是, 给你一个不大于 10^9 的正整数 N, 有多少组神的安排 (n, m) 满足 n <= N 呢?

分析:
根据题意列出方程 2*m^2-2*m=n^2-n,转换为佩尔方程形式:(2*n-1)^2-2*(2*m-1)^2=-1
记 x = 2n-1,y=2m-1,原方程即为:x^2-2*y^2=-1。为了求解该方程的通解,需要先计算x^2-2*y^2=1的基本解,显然(3,2)为基本解,套用上述得到的通解公式,可以得到 xn,yn 的迭代公式:

x(n+1) = 3xn±4yn
y(n+1) = 2xn±3yn

#由于题目解要求为正整数且 x(n+1) 大于 xn,(1,1)为方程基本解,因此选择乘积相加,而非相减。
x(n+1) = 3xn+4yn
y(n+1) = 2xn+3yn

代码实现:

N = 54
def cal_probability2():
    global N

    x = 1
    y = 1
    n = 0

    xy_list = []
    while (x + 1) / 2 <= N:
        a = x
        b = y
        if x>y:
            xy_list.append(((x + 1) // 2, (y + 1) // 2))
            n += 1
        x = 3 * a + 4 * b
        y = 2 * a + 3 * b
    print(n)
    print(xy_list)

输出结果为:

2
[(4, 3), (21, 15)]

案例二

数1,3,6,10,....构成数列{a_n},数1,4,9,16,....构成数列{b_n}。将既是数列{a_n}中的项,又是数列{b_n}中的项按照从小到大的顺序
构成新的数列{c_n},试写出数列{c_n}的前4项。

分析:

数列{a_n}的通项公式是 a_n = n(n+1)/2,数列{b_n}的通项公式是 b_n = n^2。要找出两个数列中的公共项,可以让两个通项公式相等。此时,令 a_n = b_m,即 n(n+1)/2 = m**2,因为不确定 m 和 n 是否完全相同,因此分开来区别表示。

由 n(n+1)/2 = m2 可得 (2n+1)2-8*m^2 = 1,令 x = 2n+1,y = m,则方程变为:x^2 - 8*y^2 = 1,其中(3,1)为该方程的基本解。

代码实现:

N = 4

n = 1
a = 3   #方程最小解(x1,y1)
b = 1

nm_list = []
while n<=4:
    x = a
    y = b
    nm_list.append(((x-1)//2,y))
    a = 3*x+8*y
    b = x+y*3
    n += 1

print(nm_list)

result = []
for w in nm_list:
    b_n = w[1]**2
    result.append(b_n)
print(result)

输出结果为:

[(1, 1), (8, 6), (49, 35), (288, 204)]
[1, 36, 1225, 41616]

后续会继续更新相关知识点。

发表评论 取消回复
表情 图片 链接 代码

分享