44 lines
1.9 KiB
Text
44 lines
1.9 KiB
Text
{def pell
|
|
{lambda {:n}
|
|
{let { {:n :n}
|
|
{:x {BN.intPart {BN.sqrt :n}}} // x=int(sqrt(n))
|
|
} {pell.r :n :x :x 1 {* 2 :x} 1 0 0 1}
|
|
}}}
|
|
-> pell
|
|
|
|
{def pell.r
|
|
{lambda {:n :x :y :z :r :e1 :e2 :f1 :f2}
|
|
{let { {:n :n} {:x :x} {:z :z} {:r :r} // no closure ->
|
|
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} // must reassign :(
|
|
{:y {BN.- {BN.* :r :z} :y}} // y=rz-y
|
|
} {let { {:n :n} {:x :x} {:y :y} {:r :r}
|
|
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
|
|
{:z {BN.intPart
|
|
{BN./ {BN.- :n {BN.* :y :y}} :z}}} // z=(n-y*y)//z
|
|
} {let { {:n :n} {:x :x} {:y :y} {:z :z}
|
|
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
|
|
{:r {BN.intPart {BN./ {BN.+ :x :y} :z}}} // r= (x+y)//z
|
|
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
|
|
{:e1 :e2} // e1=e2
|
|
{:e2 {BN.+ {BN.* :r :e2} :e1}} // e2=r*e2+e1
|
|
{:f1 :f2} // f1=f2
|
|
{:f2 {BN.+ {BN.* :r :f2} :f1}} // f2=r*f2+f1
|
|
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
|
|
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
|
|
{:a {BN.+ :e2 {BN.* :x :f2}}} // a=e2+x*f2
|
|
{:b :f2} // b=f2
|
|
} {if {= {BN.compare {BN.- {BN.* :a :a}
|
|
{BN.* :n {BN.* :b :b}}}
|
|
1}
|
|
0} // a*a-n*b*b == 1
|
|
then {div}x{sup 2} - n*y{sup 2} = 1 for n=:n, x=:a, y=:b
|
|
else {pell.r :n :x :y :z :r :e1 :e2 :f1 :f2} // do it again
|
|
}}}}}}}}
|
|
-> pell.r
|
|
|
|
{S.map pell 61 109 181 277}
|
|
->
|
|
x^2 - n*y^2 = 1 for n=61, x=1766319049, y=226153980
|
|
x^2 - n*y^2 = 1 for n=109, x=158070671986249, y=15140424455100
|
|
x^2 - n*y^2 = 1 for n=181, x=2469645423824185801, y=183567298683461940
|
|
x^2 - n*y^2 = 1 for n=277, x=159150073798980475849, y=9562401173878027020
|