lambdaway :: pell

6 | list | login |load

| |

_h1 the pell's equation: {code x{sup 2}-n*y{sup 2}=1}

_p Following [[https://rosettacode.org/wiki/Pell%27s_equation|https://rosettacode.org/wiki/Pell%27s_equation]].

_p Pell's equation (also called the Pell–Fermat equation) is a Diophantine equation of the form: {code x{sup 2}-n*y{sup 2} = 1} with integer solutions for x and y, where {b n} is a given non-square positive integer. The task is to find the smallest solution in positive integers to Pell's equation for n = {b [61, 109, 181, 277]}.

_p Computing big numbers requires the {b [[lib_BN]]} library.

{pre
'{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}
}}} 
-> {def pell
 {lambda {:n}
  {let { {:n :n}
         {:x {BN.intPart {BN.sqrt :n}}}
       } {pell.r :n :x :x 1 {* 2 :x} 1 0 0 1}
}}}

'{def pell.r
 {lambda {:n :x :y :z :r :e1 :e2 :f1 :f2}
  {let { {:n :n} {:x :x} {:z :z} {:r :r}
         {:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} 
         {: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                             // a*a-n*b*b == 1
             {BN.- {BN.* :a :a}
                   {BN.* :n {BN.* :b :b}}}
             1}
            0}
   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
}}}}}}}} 
-> {def pell.r
 {lambda {:n :x :y :z :r :e1 :e2 :f1 :f2}
  {let { {:n :n} {:x :x}   {:y {BN.- {BN.* :r :z} :y}}   {:z :z} {:r :r}
         {:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}                                {:i :i}    
       } {let { {:n :n} {:x :x} {:y :y}   {:z {BN.intPart {BN./ {BN.- :n {BN.* :y :y}} :z}}}   {:r :r}
                {:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}                         {:i :i}                
              } {let { {:n :n} {:x :x} {:y :y} {:z :z}   {:r {BN.intPart {BN./ {BN.+ :x :y} :z}}}
                       {:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}                  {:i :i}                     
                     } {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
                              {:e1 :e2} {:e2 {BN.+ {BN.* :r :e2} :e1}}
                              {:f1 :f2} {:f2 {BN.+ {BN.* :r :f2} :f1}}          {:i :i}
                            } {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
                                     {:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}    {:i :i}
                                     {:a {BN.+ :e2 {BN.* :x :f2}}}
                                     {:b :f2}
                                   } {if {or {= {BN.compare {BN.- {BN.* :a :a} {BN.* :n {BN.* :b :b}}} 1} 0}
                                             {> :i 100}}
                                      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}
}}}}}}}}
 {BN.DEC 50}
'{S.map pell 61 109 181 277}
->
{S.map pell 61 109 181 277} 
}

{require lib_BN}