lambdaspeech :: horner

_h1 Horner & Newton {sup (for polys)}

_p Some previous attempts can be seen in [[horner]] and [[newton]]. Here, we will follow a [[W. Kahan|http://www.cs.berkeley.edu/~wkahan/Math128/Poly.pdf]]'s paper giving a synthesis of the two methods for polynomial functions. Given the coefficients {code a{sub j}} of the polynomial {code A(x) = Σ{sub 0}{sup N} a{sub j}x{sup N-j}}, and a numerical value {code z}, we can compute both {code p := A(z)} and the derivative {code q := A'(z)} by means of {b Horner}'s recurrence:
     
{pre
q := 0 ; 
p := a{sub 0} ;
for j = 1 to N do 
    q := zq + p ;
    p := zp + a{sub j} ;
}

_p We translate this "pseudo-code" in lambdatalk:

{pre
'{def horner
 {def horner.rec
  {lambda {:l :x :q :p}
   {if {equal? :l nil}
    then {cons :p {/ :p :q}}
    else {horner.rec {cdr :l}
                     :x 
                     {+ {* :x :q} :p}
                     {+ {* :x :p} {car :l}}}}}}
 {lambda {:l :x}
  {horner.rec :l :x 0 0}}} 
-> {def horner
 {def horner.rec {lambda {:l :x :q :p}
   {if {equal? :l nil}
    then {cons :p {/ :p :q}}
    else {horner.rec {cdr :l} :x 
                   {+ {* :x :q} :p}
                   {+ {* :x :p} {car :l}}}}}}
 {lambda {:l :x} {horner.rec :l :x 0 0}}}
}
_p Note that the {code horner} function returns the pair {code [f(x),f(x)/f'(x)]}. The first term will be used to evaluate {code f(x)} for a given value {code x}. For instance:
{pre
'{car {horner {list.new 9 8 7 6 5 4 3 2 1} 10}}
-> {car {horner {list.new 9 8 7 6 5 4 3 2 1} 10}}    // as you might guess!
}

_p The second will be used in the Newton's method:

{pre
'{def newton
 {def epsilon 1.0e-15}
 {lambda {:l :x}
  {let { {:l :l} {:x :x} {:h {horner :l :x}} }
  {if {< {abs {car :h}} {epsilon}}
   then :x
   else {newton :l {- :x {cdr :h}}}}}}}
-> {def newton
 {lambda {:l :x}
  {let { {:l :l} {:x :x} {:h {horner :l :x}} }
  {if {< {abs {car :h}} 1.0e-15}
   then :x
   else {newton :l {- :x {cdr :h}}}}}}}
}

_h5 Some tests

_ul 1) Compute the root of {code x{sup 2}-2=0}  which is known to be  √2 = {code {sqrt 2}}:
{pre
'{newton {list.new 1 0 -2} 1} 
-> x = {newton {list.new 1 0 -2} 1}                              // Identical!
}

_ul 2) Compute the root of {code x{sup 2}-1-1=0}  which is known to be Φ = {code {/ {+ 1 {sqrt 5}} 2}}:
{pre
'{newton {list.new 1 -1 -1} 1.5} 
-> x = {newton {list.new 1 -1 -1} 1.5}                               // Identical!
}

_ul 3) Compute the root of {code cos(x)=0} which is known to be π/2. We will use the {b Mac Laurin's} development:

{pre {center cos(x) = Σ{sub k=0}{sup k=10} (-1){sup 2k}x{sup 2k}/{sub (2k)!}}}

_p So let's define the factorial {code !} and build the serie {code cosx} with {code 10} terms, {code k ∈ [10,1]}:

{pre

'{def ! {lambda {:a :b}
  {if {< :b 1} then :a else {! {* :a :b} {- :b 1}}}}}
-> {def !
 {lambda {:a :b}
  {if {< :b 1} then :a else {! {* :a :b} {- :b 1}}}}}

'{def cosx {list.new {map {lambda {:n}
  {/ {pow -1 :n} {! 1 {* 2 :n}}} 0} {serie 10 1 -1}} 1}}
-> {def cosx {list.new 
 {map {lambda {:n} {/ {pow -1 :n} {! 1 {* 2 :n}}} 0} {serie 10 1 -1}}
 1}}
}

_p Note: we could be 15% faster - but less pure - replacing the factorial {code '{! a b}} by {code '{* {serie a b}}}.

_p Finally we can call {code newton} on {code cosx} function with {code 1} as initial value:

{pre
'{* 2 {newton {cosx} 1}}    // which is {PI} 
-> x = {* 2 {newton {cosx} 1}}                              // Identical!

}
_p {b It's simple, quick and accurate!} {i This page is computed in 15ms on a recent laptop}.