Haskell i Pell

Jak mówiłem tak zrobiłem, tj. sprawdziłem algorytm Wildbergera w praniu. Poniżej kod funkcji rozwiązującej równanie Pella w oparciu o ów algorytm w Haskellu. Poza rekurencją i "pattern matching" nie demonstruje on pazura Haskella, ale mam nadzieję, że jest równie przyjemny w czytaniu jak był w pisaniu...

solvePell :: Integer -> (Integer, Integer)

solvePell d = (b11, b21) where 
   (b11, b12, b21, b22) = searchFixedPointMtr (1, 1, 0, 1) where
      searchFixedPointMtr p | (o11, o12, o21, o22) == (1, 0, 0 , -d) = r
               | o11 > 0 && o22 < 0 = searchFixedPointMtr r
               |otherwise = searchFixedPointMtr tr where
         r  = mult p (1, 1, 0, 1)
         tr = mult p (1, 0, 1, 1)
         (o11, o12, o21, o22) = mult( mult (transpose r) (1,0,0,-d) ) r
         transpose (a11 ,a12, a21, a22) = (a11, a21, a12, a22)
         mult (a11, a12, a21, a22) (b11, b12, b21, b22)=
            (a11*b11+a12*b21, a11*b12+a12*b22, 
             a21*b11+a22*b21, a21*b12+a22*b22)

Funkcja solvePell, bierze d z rówania Pella jako argument i zwraca parę liczb będąca pewnym rozwiązaniem równiania. Trzeba pamiętać, by d nie było kwadratem - tego funkcja nie sprawdza. Na kwadratach wchodzi zwykle w nieskończoną rekurencję ... Niestety, są i smutniejsze wieści. Algorym nie zawsze znajduje rozwiązanie fundamentalne. Np. dla d = 13, rozwiązaniem fundamentalnym jest x=18 y=5, podczas gdy zaimplementowany algorytm daje rozwiązanie: x=649 y=180, które jest następnym z kolei (bo 649=18^2+13*5^2 i 180 = 2*5*18). Tym samym zrozumiałe staje się, dlaczego w pracy Wildbergera nie było stosownego dowodu ;). Warto pomyśleć, czy nie da się tej metody jakoś zmodyfikować (nie zmieniając jednak jej ducha) by algorytm zwracał zawsze rozwiazanie fundamentalne.

Do badań porównawczych użyłem sprawdzonej już w bojach, on-linowej maszyny do poszukiwania fundamentalnych rozwiązań równania Pella którą niniejszym polecam: Diophantus Quadraticus. Dla miłośników numerologii, przykładowe wynik działania funkcji (pod Hugsem):

Main> solvePell 12
(7,2)
Main> solvePell 123
(122,11)
Main> solvePell 1234
(586327869067265,16691023073856)
Main> solvePell 12345
(1196823028442576899590849641,10771703481902106796084652)
Main> solvePell 123456
(32153667637494049,91511235212695)
Main> solvePell 1234567
(20371567825887579087969922203933358793493846332810110697
4127231916998110712447355624,1833441773536251588833840127
75408990796176026949965279326746283914164614149841725)
Main> solvePell 12345678
(84798178834254206501069091776063801828858638442491294173
283622693331157941169053721682146500840778040795681426799
746373774600679194577190578795925022609963251697702932038
103575417349220734335472506700750376403558084067176117304
353961009320426319261060613871588088528284393528998330022
695959115757272835324559892633223467646161202348129574120
811322088113694783,24133985779040703487042001168546778233
886602226457226201732495439539842758227586390869258263482
721919445559004333460853375915877499010684711766338764074
438034949056494509471323897949013056855661657290867499070
288510685596491236418010367399292614822499208206228574471
914923864025707862772212062023636759525773725225355147035
172813431107542245055578655655064)