凤凰木的笔记

浅说佩尔方程(2)：PQA解法

更新于 2024-05-24 03:50

本文主要讨论佩尔方程(Pell equations)

x^{2} - Dy^{2} = 1 \tag{e1}

的解法实现。

无理算术平方根的连分式展开

我们知道，任何无理算术平方根\sqrt{N}都可以展开成循环的连分式。以 \sqrt{7} 为例，展开如下:

\begin{align} \sqrt{7} & = 2 + \sqrt{7} - 2\\\ & = 2 + \cfrac{1}{\cfrac{1}{\sqrt{7} - 2}}\\\ & = 2 + \cfrac{1}{\cfrac{\sqrt{7} + 2}{3}}\\\ & = 2 + \cfrac{1}{1+\cfrac{\sqrt{7} - 1}{3}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{\cfrac{\sqrt{7}+1}{2}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{\sqrt{7}-1}{2}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{\cfrac{\sqrt{7}+1}{3}}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{\cfrac{\sqrt{7}-2}{3}}}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{\sqrt{7}-2}}}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{1+\cfrac{\sqrt{7}-1}{3}}}}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{1+\cfrac{\sqrt{7}-1}{3}}}}}}}}}}\\\ & = 2 + \cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1 + \cfrac{1}{4+\cfrac{1}{1+\cfrac{\sqrt{7}-1}{3}}}}}}}}}}}}}}\\\ & = [2;\dot{1}, 1, 1, \dot{4}] \end{align}

这种连分式展开具有很多有意思的特性，其中最有趣的可能是生成佩尔方程(Pell equations)

x^{2} - Dy^{2} = 1 \tag{15}

的正整数解(其中 D 为正整数，且非完全平方数)。

从无理算术平方根的连分式展开生成佩尔方程的解

可以证明，无理算术平方根 \sqrt{D} 的连分式展开, 取前 n 项约分到最简分数，某些项的分子和分母恰好是佩尔方程 (e1) 的一组正整数解，并且佩尔方程 (e1) 的所有非平凡正整数解都能从这样的展开约分中得到。

以

x^{2} - 7y^{2} = 1 \tag{16}

为例，

  1. 取 [2, 1, 1, 1] 约分，得 \cfrac{8}{3}, 于是得到式 (16) 的最小非平凡正整数解 (8, 3) .
  2. 取 [2, 1, 1, 1, 4, 1, 1, 1] 约分，得 \cfrac{127}{48}, 于是得到式 (16) 的下一组非平凡正整数解 (127, 48) .
  3. 取 [2, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1] 约分，得 \cfrac{2024}{765}, 于是得到式 (16) 的下一组非平凡正整数解 (2024, 765) .
  4. 这样继续下去，我们就能得到方程 x^{2} - 7y^{2} = 1 的所有正整数解。

PQa 算法解佩尔方程

基于上述过程的 PQa 算法(PQa algorithm)，可以求出标准佩尔方程 x^{2} - Dy^{2} = 1 的所有非平凡正整数解。

所谓 PQa 算法，根据论文 “Solving the generalized Pell equation x2 − Dy2 = N”(1), 即：

取正整数 P_0, Q_0, D, 使得 Q_0 \ne 0, D > 0, D 不为完全平方数， 且 P^2_0 \equiv D(mod \quad Q_0).

对任意的 i \geqslant 0, 令 \begin{align} & a_i = \lfloor \cfrac{P_i + \sqrt{D}}{Q_i} \rfloor, \\ & A_i = a_iA_{i-1} + A_{i-2}, \\ & B_i = a_iB_{i-1} + B_{i-2} , \\ & G_i = a_iG_{i-1} + G_{i-2} , \\ \end{align}

对任意的 i \geqslant 1, 令

\begin{align} & P_i = a_{i-1}Q_{i-1} − P_{i-1}, \\ & Q_i = \cfrac{D - {P_i}^2}{Q_{i-1}} . \\ \end{align}

特别地，

\begin{align} & A_{-2} = 0, A_{-1} = 1, \\ & B_{-2} = 1, B_{-1} = 0, \\ & G_{-2} = -P_0 , G_{-1} = Q_0 . \\ \end{align}

这样构造出来的数列 {A_i}, {B_i}, {G_i}, {P_i}, {Q_i}, {a_i} 具有很多有趣的性质，论文1列出了至少 27 种。其中特别需要注意的：

A. {A_i}, {B_i}, {G_i}, {P_i}, {Q_i}, {a_i} 全部是整数列，并且 P_i, Q_i, a_i 最终呈现周期性。

B.

\begin{align} \cfrac{P_0 + \sqrt{D}}{Q_0} & = a_0+\cfrac{1}{a_1+\cfrac{1}{a_2+\cfrac{1}{a_3+\cfrac{1}{a_4+\cfrac{1}{a_5+\cfrac{1}{a_6+\cfrac{1}{a_7+\ldots}}}}}}} \\ & = [a_0;a_1,a_2,a_3,a_4,a_5,a_6,a_7\ldots] \tag{17} \end{align}

C.

\cfrac{P_0 + \sqrt{D}}{Q_0} = \lim_{i \to \infty} \cfrac{A_i}{B_i} \tag{18}

D.

G^2_{i-1} - DB^2_{i-1} = (-1)^iQ_0Q_i, for{\quad}i > 0 \tag{19}

式 {(19)} 尤为重要。只要我们取 Q_0 = 1, i = 2k, 式 {(19)} 即变为

G^2_{i-1} - DB^2_{i-1} = Q_i

这表明只要 Q_i = 1, 那么 (G_{i-1}, B_{i-1}) 即为标准型佩尔方程 (e1) 的一组解。事实上，(e1) 的所有解都可以由此生成，这里略过证明。

下面用 Haskell 实现这个算法。

首先，定义一个数据结构来存储数列 {A_i}, {B_i}, {G_i}, {P_i}, {Q_i}, {a_i} 各项:

import Prelude hiding (pi)

data PQa = PQa { ai :: Integer
               , bi :: Integer
               , gi :: Integer
               , aj :: Integer
               , pi :: Integer
               , qi :: Integer
               } deriving Show

下面这个 pqa 函数迭代出 [({PQa}_{i}, {PQa}_{i-1})] 的无穷数列：

pqas :: Integer -> Integer -> Integer -> [(PQa, PQa)]
pqas d p0 q0 = iterate pqa' (pqa0, pqa_1) where
    pqa0 = PQa { ai = c
               , bi = 1
               , gi = c * q0 - p0
               , aj = c
               , pi = p0
               , qi = q0
               }
    pqa_1 = PQa { ai = 1
                , bi = 0
                , gi = q0
                , aj = c
                , pi = 1
                , qi = 1
                }
    c = floor $ (fromInteger p0 + sqrt (fromInteger d)) / (fromInteger q0)
    pqa' :: (PQa, PQa) -> (PQa, PQa)
    pqa' (acc, acc') = (acc'', acc) where
        pi' = (aj acc) * (qi acc) - pi acc
        qi' = quot (d - pi'^2) (qi acc)
        aj' = floor $ (fromInteger pi' + sqrt (fromInteger d)) / (fromInteger qi')
        ai' = aj' * (ai acc) + ai acc'
        bi' = aj' * (bi acc) + bi acc'
        gi' = aj' * (gi acc) + gi acc'
        acc'' = PQa { ai = ai'
                    , bi = bi'
                    , gi = gi'
                    , aj = aj'
                    , pi = pi'
                    , qi = qi'
                    }

取 P_0 = 0, Q_0 = 1, 得到

pqas01 :: Integer -> [(PQa, PQa)]
pqas01 d = pqas d 0 1

滤掉 Q_i \ne 1 的项，就得到标准佩尔方程 (e1) 的所有非平凡正整数解，这个结果以 Haskell 的惰性无穷列表来表示：

pqaSolutions01 :: Integer -> [(Integer, Integer)]
pqaSolutions01 d = map (\(x, y) -> (gi y, bi y)) $ filter (\(x, y) -> qi x == 1 && bi y /= 0) $ pqas01 d

然后求最小正整数解只要取第一项就可以了:

minimumPositiveSolution :: Integer -> (Integer, Integer)                   
minimumPositiveSolution = head . pqaSolutions01

比如我们来求佩尔方程 x^{2} - 7y^{2} = 1 的最小前 30 组非平凡解，只需要在 ghci 里:

> mapM_ print $ take 40 $ pqaSolutions01 7
(8,3)
(127,48)
(2024,765)
(32257,12192)
(514088,194307)
(8193151,3096720)
(130576328,49353213)
(2081028097,786554688)
(33165873224,12535521795)
(528572943487,199781794032)
(8424001222568,3183973182717)
(134255446617601,50743789129440)
(2139663144659048,808716652888323)
(34100354867927167,12888722657083728)
(543466014742175624,205410845860451325)
(8661355881006882817,3273684811110137472)
(138038228081367949448,52173546131901748227)
(2199950293420880308351,831503053299317834160)
(35061166466652716984168,13251875306657183598333)
(558778713173022591438337,211198501853215619739168)
(8905398244301708746029224,3365924154344792732228355)
(141927593195654317345029247,53643587967663468095914512)
(2261936092886167368774438728,854931483328270696802403837)
(36049049892983023583045990401,13625260145284667680742546880)
(574522862194842209959961407688,217149230841226412195078346243)
(9156316745224492335776336532607,3460762433314337927440510993008)
(145926545061397035162461423114024,55155049702188180426853097541885)
(2325668404237128070263606433291777,879020032801696548902209049677152)
(37064767922732652089055241509554408,14009165475124956602008491697292547)
(590710618359485305354620257719578751,223267627569197609083233658107003600)
(9414305125829032233584868882003705608,3558272875632036788729730038014765053)
(150038171394905030432003281854339710977,56709098382543391010592446950129237248)
(2391196437192651454678467640787431670024,903787301245062219380749421164053030915)
(38109104823687518244423478970744567009407,14403887721538452119081398291674719257392)
(607354480741807640456097195891125640480488,229558416243370171685921623245631455087357)
(9679562587045234729053131655287265680678401,3658530772172384294855664573638428562140320)
(154265646911981948024394009288705125250373928,58306933938514778546004711554969225539157763)
(2458570788004665933661251016963994738325304447,929252412244064072441219720305869180064383888)
(39182866961162672990555622262135210687954497224,14809731661966510380513510813338937655490984445)
(624467300590598101915228705177199376268946651137,236026454179220102015774953293117133307791367232)

可以看到，解的数值以一种相当恐怖的速度增长(事实上是以指数增长的)。

简单看下对不同的 D, 佩尔方程 e1 最小非平凡正整数解的分布情况。

minimumPositiveSolutions :: Integer -> [(Integer, (Integer, Integer))]
minimumPositiveSolutions n = map (\x -> (x, minimumPositiveSolution x)) xs where
    xs = filter (\x -> let y = (round $ sqrt $ fromInteger x) in y^2 /= x) [1..n]

令 n = 100, 得到 2 \leqslant D \leqslant 99 范围内的最小非平凡正整数解，如下表：

D x y
2 3 2
3 2 1
5 9 4
6 5 2
7 8 3
8 3 1
10 19 6
11 10 3
12 7 2
13 649 180
14 15 4
15 4 1
17 33 8
18 17 4
19 170 39
20 9 2
21 55 12
22 197 42
23 24 5
24 5 1
26 51 10
27 26 5
28 127 24
29 9801 1820
30 11 2
31 1520 273
32 17 3
33 23 4
34 35 6
35 6 1
37 73 12
38 37 6
39 25 4
40 19 3
41 2049 320
42 13 2
43 3482 531
44 199 30
45 161 24
46 24335 3588
47 48 7
48 7 1
50 99 14
51 50 7
52 649 90
53 66249 9100
54 485 66
55 89 12
56 15 2
57 151 20
58 19603 2574
59 530 69
60 31 4
61 1766319049 226153980
62 63 8
63 8 1
65 129 16
66 65 8
67 48842 5967
68 33 4
69 7775 936
70 251 30
71 3480 413
72 17 2
73 2281249 267000
74 3699 430
75 26 3
76 57799 6630
77 351 40
78 53 6
79 80 9
80 9 1
82 163 18
83 82 9
84 55 6
85 285769 30996
86 10405 1122
87 28 3
88 197 21
89 500001 53000
90 19 2
91 1574 165
92 1151 120
93 12151 1260
94 2143295 221064
95 39 4
96 49 5
97 62809633 6377352
98 99 10
99 10 1

可见，最大的一组是 D = 61 时，(x, y) = (1766319049, 226153980).

> take 10 $ sortBy (\(_, (x1, _)) (_, (x2, _)) -> compare x2 x1 ) $ minimumPositiveSolutions 1000

得到 D \leqslant 1000 时，最大的前 10 组:

D x y
661 16421658242965910275055840472270471049 638728478116949861246791167518480580
541 3707453360023867028800645599667005001 159395869721270110077187138775196900
421 3879474045914926879468217167061449 189073995951839020880499780706260
769 535781868388881310859702308423201 19320788325040337217824455505160
937 480644425002415999597113107233 15701968936415353889062192632
613 464018873584078278910994299849 18741545784831997880308784340
991 379516400906811930638014896080 12055735790331359447442538767
601 38902815462492318420311478049 1586878942101888360258625080
673 4765506835465395993032041249 183696788896587421699032600
919 4481603010937119451551263720 147834442396536759781499589 
> maximumBy (\(_, (x1, _)) (_, (x2, _)) -> compare x1 x2 ) $ minimumPositiveSolutions 10000

得到 D \leqslant 10^4 时，最大的一组是当 D = 9949 时，

(x, y) = (23551019614858223475933893515741198183163217312913587552899320396564478041197360918469501097146448985821854465768234479384482435117587576296319428592757548743265811454938493105633433315887574461850060798834186249,236113054062810988826514929828649213339688520849720684015415366388626019230322623673232286474879711003505448178417385617250641629212134427833135509077013929303770208680820795381507114806491325360400076633910900)

同样的办法，可以得到 D \leqslant 10^5，最大的一组是当 D=92821 时，

(x, y) = (9138623307350640837938507246130056541284450249628186652711211526290141593088447074959889691586109808446586788487329855457947449359122812196406183471094085243300140113776722320673519989258197664615506581635543258793642137960637895993371068043922291726707457424732865511811603416709407607285667077345150397643032877987726841522351750655900253781916899660152174806370963781778224969807395941841736594668461346490187409918540644500596164105739747237839541686461147582701486262272823400803243678345216257264359841059899047944315969645240323636858831783554703663860186919961543958325467259751629633944639211198469786177829884186709831569216425531620090438319686256712276877552760998820897069126292769743115626179169131031894821449, 29995606985908278391394125599135320972040960690735594014932999302562021607917845625308781046774668965028391500344676078764117785968987025596642346687378667667268554217805689486419938846247966636235622466824853671443727025761340681500617596724921086840802447598463611192953164168684707145263201173090354218890224078617132621981270144387271894050379162839910605924620561997843353743092508084766465038567232036892276652930541497664634108271302723027950641070447197899109815123567584348859003464872979489356061130456764804964145771699318883015346111831572404557647659469760522059443590269944115798752448599110602793086535348170280145151753392282563680433826089979371590683960862942643266548255841329663100421901070962552821740)

当 D \leqslant 10^8 时，

index D x 位数 y 位数
1 9747061 8982 8979
2 9670621 8682 8678
3 9554869 8625 8621
4 9726781 8527 8524
5 9697549 8515 8512
6 9965341 8493 8489
7 9549901 8493 8490
8 9923629 8480 8477
9 9508669 8437 8433
10 8658589 8404 8401
11 9233701 8402 8399
12 9895141 8343 8340
13 9957061 8341 8337
14 9882181 8303 8300
15 9905989 8267 8264
16 8687149 8235 8232
17 9058669 8229 8226
18 8983549 8222 8219
19 9431269 8209 8205
20 9499669 8159 8155
21 8866621 8155 8152
22 9834829 8150 8147
23 9853141 8129 8125
24 8675221 8123 8119
25 8691589 8122 8118
26 9344581 8121 8117
27 9880861 8114 8110
28 9941149 8109 8105
29 9573061 8105 8101
30 9773989 8102 8098
31 9369229 8085 8082
32 8615149 8076 8072
33 9342589 8061 8057
34 9850801 8043 8039
35 9257329 8018 8014
36 9269269 8014 8010
37 9676549 8008 8005
38 9117109 8005 8002
39 9843661 8005 8001
40 9837829 7999 7996
41 9795109 7993 7989
42 8825821 7986 7983
43 9945181 7976 7972
44 9313789 7953 7950
45 8356261 7952 7949
46 9097141 7942 7939
47 9761761 7940 7936
48 9919849 7940 7937
49 8061421 7935 7931
50 8977021 7931 7928
51 7934749 7928 7925
52 9238909 7911 7908
53 9824389 7892 7888
54 9929929 7886 7883
55 9550069 7885 7882
56 9959041 7884 7880
57 9660421 7879 7876
58 8533669 7877 7873
59 9808681 7868 7865
60 9547561 7864 7861
61 9576121 7862 7859
62 8484589 7858 7855
63 9791149 7845 7841
64 9340501 7842 7838
65 8718109 7839 7836
66 9335869 7836 7833
67 9213829 7830 7826
68 9690349 7820 7817
69 9697801 7817 7813
70 9973261 7812 7808
71 9620389 7803 7800
72 8998861 7802 7798
73 9626269 7799 7795
74 9799969 7798 7794
75 9283429 7796 7793
76 9508801 7794 7791
77 9979621 7782 7779
78 9987961 7779 7775
79 8316421 7777 7774
80 9589141 7774 7771
81 9032269 7772 7768
82 9627601 7771 7768
83 8562709 7763 7760
84 9808621 7762 7758
85 9901189 7761 7758
86 9882601 7754 7750
87 8499661 7752 7749
88 9045061 7741 7737
89 7485949 7739 7736
90 9238189 7733 7730
91 9857149 7730 7726
92 9414409 7726 7723
93 9991669 7724 7720
94 9457561 7720 7717
95 8954401 7720 7717
96 7726261 7714 7710
97 9719581 7708 7705
98 9781249 7707 7704
99 8575141 7705 7702
100 8950309 7700 7696

我顺便把上面这个 Haskell 实现写成了 Web API。这个 API 接受 d, 和 n 两个参数。其中 d 即方程里的参数 D, 而 n 表示数出方程的前 n 组最小整整数解。url 是 /v1/pell.

受服务器物理限制，有范围: 0 \lt d \lt 10^9, 0 \lt n \lt 100.

举例：

（本文完）