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

本文主要讨论佩尔方程(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}$

为例，

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

PQa 算法解佩尔方程

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

所谓 PQa 算法，根据论文 “Solving the generalized Pell equation x^2 − Dy^2 = N”(), 即：

取正整数 $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$ 最终呈现周期性。

$\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}$

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

$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$ .

举例：

求 $x^{2} - 7y^{2} = 1$ 的最小前 $30$ 组非平凡解: https://www.weiwen.org/v1/pell?d=7&n=30
求 $x^{2} - 991y^{2} = 1$ 的最小前 $10$ 组非平凡解（小心，非常大）: https://www.weiwen.org/v1/pell?d=991&n=10

我还写了个更直观的页面，可以点进去玩一下： https://www.weiwen.org/pell

（本文完）

Solving the generalized Pell equation x^2 − Dy^2 = N, John P. Robertson, 2004