Lagrange Interpolation: Difference between revisions
(Added Wren) |
m (→{{header|jq}}) |
||
(12 intermediate revisions by 7 users not shown) | |||
Line 1: | Line 1: | ||
{{draft task}} |
{{draft task}} |
||
The task is to implement the |
The task is to implement the Lagrange Interpolation formula and use it to solve the example problem to find a polynomial |
||
P of degree<4 satisfying P(1)=1 P(2)=4 P(3)=1 P(4)=5 as described at |
P of degree<4 satisfying P(1)=1 P(2)=4 P(3)=1 P(4)=5 as described at https://brilliant.org/wiki/lagrange-interpolation |
||
or https://en.wikipedia.org/wiki/Lagrange_polynomial |
|||
;Related task [[Curve that touches three points]] |
;Related task [[Curve that touches three points]] |
||
Line 19: | Line 20: | ||
<pre> |
<pre> |
||
-21 + 215/6*x - 16*x^2 + 13/6*x^3 |
-21 + 215/6*x - 16*x^2 + 13/6*x^3 |
||
</pre> |
|||
=={{header|J}}== |
|||
From the J wiki [[j:Phrases/Polynomials|polynomial phrases]] page, we have a couple implementations: |
|||
<syntaxhighlight lang=J>lagrange=: ] +/ .* [ (] % p.~) 1 p.@;~&1\. [ |
|||
lagrange1=: %.@(^/ i.@#)@[ +/ .* ]</syntaxhighlight> |
|||
Task example: |
|||
<syntaxhighlight lang=J> 1 2 3 4 lagrange 1 4 1 5 |
|||
_21 215r6 _16 13r6 |
|||
1 2 3 4 lagrange1 1 4 1 5 |
|||
_21 35.8333 _16 2.16667</syntaxhighlight> |
|||
(The J representation of a polynomial is the list of coefficients of the polynomial such that the index of the coefficient is the power applied to the polynomial's variable.) |
|||
=={{header|jq}}== |
|||
{{works with|jq}} |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
'''With minor modifications related to the module system, the program presented here also works with jaq, the Rust implementation of jq''' |
|||
'''The function `lagrange/1` is adapted from [[#Wren|Wren]].''' |
|||
This entry uses the [[:Category:jq/polynomials.jq|"polynomials" module]] in which a polynomial is |
|||
represented by the JSON array comprised of the polynomial's |
|||
coefficients, with the entry at .[i] corresponding to the coefficient |
|||
of x^i. The canonical form of a polynomial of degree n is represented |
|||
by an array of length n+1. |
|||
<syntaxhighlight lang="jq"> |
|||
include "polynomials" {search: "./"}; |
|||
# Returns the Lagrange interpolating polynomial which passes through a list of points. |
|||
def lagrange($pts): |
|||
($pts|length) as $c |
|||
| reduce range(0;$c) as $i ({polys: []}; |
|||
.poly = [1] |
|||
| reduce range(0;$c) as $j (.; |
|||
if ($i != $j) |
|||
then .poly = (multiply(.poly; [-$pts[$j][0], 1])) |
|||
else . |
|||
end ) |
|||
| (.poly | eval($pts[$i][0])) as $d |
|||
| .polys[$i] = (.poly | scalarDivide($d)) ) |
|||
| reduce range(0;$c) as $i (.sum = [0]; |
|||
.polys[$i] = (.polys[$i] | scalarMultiply($pts[$i][1])) |
|||
| .sum = add(.sum; .polys[$i]) ) |
|||
| .sum ; |
|||
def pts: [ |
|||
[1, 1], |
|||
[2, 4], |
|||
[3, 1], |
|||
[4, 5] |
|||
]; |
|||
lagrange(pts) | pp |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
<pre> |
|||
2.1666666666666665x³-16.0x²+35.83333333333333x-21.0 |
|||
</pre> |
|||
=={{header|Julia}}== |
|||
Note the Polynomials module prints polynomials in order from degree 0 to n rather than n to zero degree. |
|||
<syntaxhighlight lang="julia">using Polynomials, SpecialPolynomials |
|||
const pts = [[1, 1], [2, 4], [3, 1], [4, 5]] |
|||
const xs = first.(pts) |
|||
const ys = last.(pts) |
|||
const p = Lagrange(xs, ys) |
|||
@show p Polynomial(p) |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
p = Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 1⋅ℓ_2(x) + 5⋅ℓ_3(x)) |
|||
Polynomial(p) = Polynomial(-21.0 + 35.83333333333333*x - 16.0*x^2 + 2.1666666666666665*x^3) |
|||
</pre> |
|||
=== Rational base polynomial answer === |
|||
<syntaxhighlight lang="julia">using Polynomials |
|||
function Lagrange(pts::Vector{Vector{Int}}) |
|||
xs = first.(pts) |
|||
cs = last.(pts) |
|||
n = length(xs) |
|||
n == 1 && return Polynomial(cs[1]) |
|||
arr = ones(Int, n) |
|||
for j in 2:n |
|||
for k in 1:j |
|||
arr[k] = (xs[k] - xs[j]) * arr[k] |
|||
end |
|||
arr[j] = prod(xs[j] - xs[i] for i in 1:(j - 1)) |
|||
end |
|||
ws = 1 .// arr |
|||
q = Polynomial(0 // 1) |
|||
x = variable(q) |
|||
for i in eachindex(ws) |
|||
m = prod(x - xs[j] for j in eachindex(xs) if j != i) |
|||
q += m * ws[i] * cs[i] |
|||
end |
|||
return q |
|||
end |
|||
const testpoints = [[1, 1], [2, 4], [3, 1], [4, 5]] |
|||
@show Lagrange(testpoints) |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
Lagrange(testpoints) = Polynomial(-21//1 + 215//6*x - 16//1*x^2 + 13//6*x^3) |
|||
</pre> |
|||
=={{header|Phix}}== |
|||
{{trans|Wren}} |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">add_poly</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">l1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">l2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">l2</span><span style="color: #0000FF;">></span><span style="color: #000000;">l1</span> <span style="color: #008080;">then</span> <span style="color: #000000;">p1</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">l2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">l1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l2</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">p1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">p2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">p1</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">mul_poly</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">p1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p1</span><span style="color: #0000FF;">),</span> |
|||
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">j</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">p1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">p2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">scalar_multiply</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">sclar_divide</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">eval_poly</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">[$]</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">)-</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">res</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">show_poly</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">l</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">poly</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;"><</span><span style="color: #000000;">l</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span><span style="color: #0000FF;">?</span><span style="color: #008000;">" "</span><span style="color: #0000FF;">:</span><span style="color: #008000;">" +"</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%g"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"x"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">></span><span style="color: #000000;">2</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"^%d"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">show</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">show_poly</span> |
|||
<span style="color: #008080;">enum</span> <span style="color: #000000;">X</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">Y</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">lagrange</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">pts</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000080;font-style:italic;">// Returns the Lagrange interpolating polynomial which passes through a list of points.</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">polys</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #004600;">null</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">c</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">poly</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">c</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">j</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #000000;">poly</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mul_poly</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">,{-</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">X</span><span style="color: #0000FF;">],</span><span style="color: #000000;">1</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">eval_poly</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">X</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #000000;">polys</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sclar_divide</span><span style="color: #0000FF;">(</span><span style="color: #000000;">poly</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">c</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #000000;">polys</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">scalar_multiply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">polys</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">Y</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">add_poly</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">polys</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">pts</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">}}</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">lip</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lagrange</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pts</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">show</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lip</span><span style="color: #0000FF;">)</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
|||
<pre> |
|||
2.16667x^3 -16x^2 +35.8333x -21 |
|||
</pre> |
|||
=={{header|Raku}}== |
|||
{{trans|Wren}} |
|||
<syntaxhighlight lang="raku" line># 20231020 Raku programming solution |
|||
class Point { has ($.x, $.y) } |
|||
sub add(@p1, @p2) { return @p1 <<+>> @p2 } # Add two polynomials. |
|||
sub multiply(@p1, @p2) { # Multiply two polynomials. |
|||
my ($m,$n) = (@p1,@p2)>>.elems; |
|||
my @prod = [ 0 xx $m + $n - 1 ]; |
|||
for ^$m X ^$n -> ($i,$j) { @prod[$i+$j] += @p1[$i] * @p2[$j] } |
|||
return @prod; |
|||
} |
|||
# Multiply a polynomial by a scalar. |
|||
sub scalarMultiply(@poly, $x) { return @poly.map: * × $x } |
|||
# Divide a polynomial by a scalar. |
|||
sub scalarDivide(@poly, $x) { return scalarMultiply(@poly, 1/$x) } |
|||
# rosettacode.org/wiki/Horner%27s_rule_for_polynomial_evaluation#Raku |
|||
sub evalPoly(@coefs, $x) { return ([o] map { $_ + $x × * }, @coefs.reverse)(0) } |
|||
sub lagrange(@pts) { |
|||
my ($c, @polys) = @pts.elems; |
|||
for ^$c -> $i { |
|||
my @poly = 1; |
|||
for ^$c -> $j { |
|||
next if $i == $j; |
|||
@poly = multiply @poly, [ -@pts[$j].x, 1 ] |
|||
} |
|||
@polys[$i] = scalarDivide @poly, evalPoly(@poly.reverse, @pts[$i].x) |
|||
} |
|||
my @sum = 0; |
|||
for ^$c -> $i { @sum = add @sum, scalarMultiply @polys[$i], @pts[$i].y } |
|||
return @sum; |
|||
} |
|||
my @pts = [<1 1>,<2 4>,<3 1>,<4 5>].map: { Point.new: x => .[0], y => .[1] }; |
|||
say [~] lagrange(@pts).kv.rotor(2).reverse.map: -> ($expo,$coef) { |
|||
"{ '+' if $coef ≥ 0 }$coef" ~ do given $expo { |
|||
when 0 { " " } |
|||
when 1 { "𝑥 " } |
|||
default { "𝑥^$_ " } |
|||
} |
|||
}</syntaxhighlight> |
|||
You may [https://ato.pxeger.com/run?1=jVS9btswEJ66-CkOiYrYsaJYbooWsSW4QIYsAYJOBQTHYCzaYSKJAinZEgJl79ipQ1EgQ_MkfYuOfYI-Qo-kZCtuhsqAQd7Px7vvPvLbD0Hu8sfHpzxbHL3_9epsHhEp4ZKzJIN7uCESupZT2GA5ZQ-qTkfm10DCsDtJXRsm6bCHUYJmuUhw58J43Pd9ZYcK9uFDGEK25pDyqEx4zEgkHQMR51HG0qh8hrMPF7X53yzALy6xmNi2kh54oDNVou87NKKxHDUxk1TwEJceBDCAogArhj5YCRyBC9MRgA5ccAFX6PmE_-jyEZrZ1q2qQwMEFutbt1Poe6ox3E3hUBUaKGOlITZ9Y_iog-S0GiCt8uFa7eWcREQ4un2zvtiSgLHIcfGMTbQ5MUlP8dyfX9EJ-oQztmIh_S98E_oi-ssFuMcqSB8juKRZRuY8pA4Xy-M1u2PH51wkVLwevpMzkUd0hhzOtmXM6IpEOckYT_Y_oqp0Jcp2ydURc04XcqeMbsCngD2ixZqpIRWq1UOoUBM63hF0RYWkve5gI7-ILAVJlqqxTCJYp5HG3DasSaUP5WwJw4x7rgZtMZPTqAUzMN4d1bZ25O0mEr-EFhmwhcr3PPSNtq4GpJE11IQGcKTqUJpRlwjlV-dUnVai1Orynk2tQdjyp_VQ02GDgWUI29N6rmoWJjKPEWrwUtONE--vXto7MmhV0zqgNNiNLDFPa11Tl0l1y8YuuL49HsIJ_r_R6xN460-Neu_Na-IkdH0KBXg-OMEADyjN0sXbNMLBkhKCh-nOcJ27lSN4xkV32Gt6N6j6wtIi5balhFLLYO8eDvoHekjKCr8_P-ETUOnNHjxAyGHJVjQBnYql1WNY36BtAGjYw1_Vtrra-uf7l6eWJ6QLgqQ1nisUb-2skBvzmNZvavO2_gU Attempt This Online!] |
|||
=={{header|RPL}}== |
|||
{{works with|HP|49}} |
|||
[[1 2 3 4][1 4 1 5]] LAGRANGE |
|||
{{out}} |
|||
<pre> |
|||
1: '(13*X^3-96*X^2+215*X-126)/6' |
|||
</pre> |
</pre> |
||
Line 26: | Line 289: | ||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
A polynomial is represented here by a list of coefficients from the lowest to the highest degree. However, the library methods which deal with polynomials expect the coefficients to be presented from highest to lowest degree so we therefore need to reverse the list before calling these methods. |
A polynomial is represented here by a list of coefficients from the lowest to the highest degree. However, the library methods which deal with polynomials expect the coefficients to be presented from highest to lowest degree so we therefore need to reverse the list before calling these methods. |
||
<syntaxhighlight lang=" |
<syntaxhighlight lang="wren">import "./dynamic" for Tuple |
||
import "./math" for Math |
import "./math" for Math |
||
import "./fmt" for Fmt |
import "./fmt" for Fmt |
||
Line 56: | Line 319: | ||
var scalarMultiply = Fn.new { |poly, x| poly.map { |coef| coef * x }.toList } |
var scalarMultiply = Fn.new { |poly, x| poly.map { |coef| coef * x }.toList } |
||
// Divide a |
// Divide a polynomial by a scalar. |
||
var scalarDivide = Fn.new { |poly, x| scalarMultiply.call(poly, 1/x) } |
var scalarDivide = Fn.new { |poly, x| scalarMultiply.call(poly, 1/x) } |
||
// Returns the Lagrange interpolating polynomial which passes through a list of points. |
|||
var lagrange = Fn.new { |pts| |
var lagrange = Fn.new { |pts| |
||
var c = pts.count |
var c = pts.count |
Latest revision as of 19:58, 28 January 2024
The task is to implement the Lagrange Interpolation formula and use it to solve the example problem to find a polynomial P of degree<4 satisfying P(1)=1 P(2)=4 P(3)=1 P(4)=5 as described at https://brilliant.org/wiki/lagrange-interpolation or https://en.wikipedia.org/wiki/Lagrange_polynomial
- Related task Curve that touches three points
F#
// Lagrange Interpolation. Nigel Galloway: September 5th., 2023
let symbol=MathNet.Symbolics.SymbolicExpression.Variable
let qi=MathNet.Symbolics.SymbolicExpression.FromInt32
let eval (g:MathNet.Symbolics.SymbolicExpression) x=let n=Map["x",MathNet.Symbolics.FloatingPoint.Real x] in MathNet.Symbolics.Evaluate.evaluate n g.Expression
let fN g=let x=symbol "x" in g|>List.fold(fun z c->(x-c)*z)(qi 1)
let fG(n,g)=let n,g=n|>List.map qi,g|>List.map qi in List.map2(fun i g->i,g,n|>List.except [i]) n g
let LIF n=fG n|>List.sumBy(fun(ci,bi,c)->bi*(fN c)/(c|>List.fold(fun z c->(ci-c)*z)(qi 1)))
printfn $"%s{LIF([1;2;3;4],[1;4;1;5]).Expand().ToString()}"
- Output:
-21 + 215/6*x - 16*x^2 + 13/6*x^3
J
From the J wiki polynomial phrases page, we have a couple implementations:
lagrange=: ] +/ .* [ (] % p.~) 1 p.@;~&1\. [
lagrange1=: %.@(^/ i.@#)@[ +/ .* ]
Task example:
1 2 3 4 lagrange 1 4 1 5
_21 215r6 _16 13r6
1 2 3 4 lagrange1 1 4 1 5
_21 35.8333 _16 2.16667
(The J representation of a polynomial is the list of coefficients of the polynomial such that the index of the coefficient is the power applied to the polynomial's variable.)
jq
Works with gojq, the Go implementation of jq
With minor modifications related to the module system, the program presented here also works with jaq, the Rust implementation of jq
The function `lagrange/1` is adapted from Wren.
This entry uses the "polynomials" module in which a polynomial is represented by the JSON array comprised of the polynomial's coefficients, with the entry at .[i] corresponding to the coefficient of x^i. The canonical form of a polynomial of degree n is represented by an array of length n+1.
include "polynomials" {search: "./"};
# Returns the Lagrange interpolating polynomial which passes through a list of points.
def lagrange($pts):
($pts|length) as $c
| reduce range(0;$c) as $i ({polys: []};
.poly = [1]
| reduce range(0;$c) as $j (.;
if ($i != $j)
then .poly = (multiply(.poly; [-$pts[$j][0], 1]))
else .
end )
| (.poly | eval($pts[$i][0])) as $d
| .polys[$i] = (.poly | scalarDivide($d)) )
| reduce range(0;$c) as $i (.sum = [0];
.polys[$i] = (.polys[$i] | scalarMultiply($pts[$i][1]))
| .sum = add(.sum; .polys[$i]) )
| .sum ;
def pts: [
[1, 1],
[2, 4],
[3, 1],
[4, 5]
];
lagrange(pts) | pp
- Output:
2.1666666666666665x³-16.0x²+35.83333333333333x-21.0
Julia
Note the Polynomials module prints polynomials in order from degree 0 to n rather than n to zero degree.
using Polynomials, SpecialPolynomials
const pts = [[1, 1], [2, 4], [3, 1], [4, 5]]
const xs = first.(pts)
const ys = last.(pts)
const p = Lagrange(xs, ys)
@show p Polynomial(p)
- Output:
p = Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 1⋅ℓ_2(x) + 5⋅ℓ_3(x)) Polynomial(p) = Polynomial(-21.0 + 35.83333333333333*x - 16.0*x^2 + 2.1666666666666665*x^3)
Rational base polynomial answer
using Polynomials
function Lagrange(pts::Vector{Vector{Int}})
xs = first.(pts)
cs = last.(pts)
n = length(xs)
n == 1 && return Polynomial(cs[1])
arr = ones(Int, n)
for j in 2:n
for k in 1:j
arr[k] = (xs[k] - xs[j]) * arr[k]
end
arr[j] = prod(xs[j] - xs[i] for i in 1:(j - 1))
end
ws = 1 .// arr
q = Polynomial(0 // 1)
x = variable(q)
for i in eachindex(ws)
m = prod(x - xs[j] for j in eachindex(xs) if j != i)
q += m * ws[i] * cs[i]
end
return q
end
const testpoints = [[1, 1], [2, 4], [3, 1], [4, 5]]
@show Lagrange(testpoints)
- Output:
Lagrange(testpoints) = Polynomial(-21//1 + 215//6*x - 16//1*x^2 + 13//6*x^3)
Phix
with javascript_semantics function add_poly(sequence p1, p2) integer l1 = length(p1), l2 = length(p2) if l2>l1 then p1 &= repeat(0,l2-l1) end if for i=1 to l2 do p1[i] += p2[i] end for return p1 end function function mul_poly(sequence p1,p2) integer m = length(p1), n = length(p2) sequence res = repeat(0,m+n-1) for i=1 to m do for j=1 to n do res[i+j-1] += p1[i]*p2[j] end for end for return res end function function scalar_multiply(sequence poly, atom x) return sq_mul(poly,x) end function function sclar_divide(sequence poly, atom x) return sq_div(poly,x) end function function eval_poly(sequence poly, atom x) integer c = length(poly) atom res = poly[$] for i=length(poly)-1 to 1 by -1 do res = res*x + poly[i] end for return res end function procedure show_poly(sequence poly) integer l = length(poly) for i=l to 1 by -1 do atom p = poly[i] if i<l then puts(1,iff(p<0?" ":" +")) end if printf(1,"%g",p) if i>1 then puts(1,"x") if i>2 then printf(1,"^%d",i-1) end if end if end for puts(1,"\n") end procedure constant show = show_poly enum X, Y function lagrange(sequence pts) // Returns the Lagrange interpolating polynomial which passes through a list of points. integer c = length(pts) sequence polys = repeat(null,c) for i=1 to c do sequence poly = {1} for j=1 to c do if i!=j then poly = mul_poly(poly,{-pts[j][X],1}) end if end for atom d = eval_poly(poly,pts[i][X]) polys[i] = sclar_divide(poly,d) end for sequence res = {0} for i=1 to c do polys[i] = scalar_multiply(polys[i],pts[i][Y]) res = add_poly(res,polys[i]) end for return res end function constant pts = {{1,1},{2,4},{3,1},{4,5}} sequence lip = lagrange(pts) show(lip)
- Output:
2.16667x^3 -16x^2 +35.8333x -21
Raku
# 20231020 Raku programming solution
class Point { has ($.x, $.y) }
sub add(@p1, @p2) { return @p1 <<+>> @p2 } # Add two polynomials.
sub multiply(@p1, @p2) { # Multiply two polynomials.
my ($m,$n) = (@p1,@p2)>>.elems;
my @prod = [ 0 xx $m + $n - 1 ];
for ^$m X ^$n -> ($i,$j) { @prod[$i+$j] += @p1[$i] * @p2[$j] }
return @prod;
}
# Multiply a polynomial by a scalar.
sub scalarMultiply(@poly, $x) { return @poly.map: * × $x }
# Divide a polynomial by a scalar.
sub scalarDivide(@poly, $x) { return scalarMultiply(@poly, 1/$x) }
# rosettacode.org/wiki/Horner%27s_rule_for_polynomial_evaluation#Raku
sub evalPoly(@coefs, $x) { return ([o] map { $_ + $x × * }, @coefs.reverse)(0) }
sub lagrange(@pts) {
my ($c, @polys) = @pts.elems;
for ^$c -> $i {
my @poly = 1;
for ^$c -> $j {
next if $i == $j;
@poly = multiply @poly, [ -@pts[$j].x, 1 ]
}
@polys[$i] = scalarDivide @poly, evalPoly(@poly.reverse, @pts[$i].x)
}
my @sum = 0;
for ^$c -> $i { @sum = add @sum, scalarMultiply @polys[$i], @pts[$i].y }
return @sum;
}
my @pts = [<1 1>,<2 4>,<3 1>,<4 5>].map: { Point.new: x => .[0], y => .[1] };
say [~] lagrange(@pts).kv.rotor(2).reverse.map: -> ($expo,$coef) {
"{ '+' if $coef ≥ 0 }$coef" ~ do given $expo {
when 0 { " " }
when 1 { "𝑥 " }
default { "𝑥^$_ " }
}
}
You may Attempt This Online!
RPL
[[1 2 3 4][1 4 1 5]] LAGRANGE
- Output:
1: '(13*X^3-96*X^2+215*X-126)/6'
Wren
A polynomial is represented here by a list of coefficients from the lowest to the highest degree. However, the library methods which deal with polynomials expect the coefficients to be presented from highest to lowest degree so we therefore need to reverse the list before calling these methods.
import "./dynamic" for Tuple
import "./math" for Math
import "./fmt" for Fmt
var Point = Tuple.create("Point", ["x", "y"])
// Add two polynomials.
var add = Fn.new { |p1, p2|
var m = p1.count
var n = p2.count
var sum = List.filled(m.max(n), 0)
for (i in 0...m) sum[i] = p1[i]
for (i in 0...n) sum[i] = sum[i] + p2[i]
return sum
}
// Multiply two polynmials.
var multiply = Fn.new { |p1, p2|
var m = p1.count
var n = p2.count
var prod = List.filled(m + n - 1, 0)
for (i in 0...m) {
for (j in 0...n) prod[i+j] = prod[i+j] + p1[i] * p2[j]
}
return prod
}
// Multiply a polynomial by a scalar.
var scalarMultiply = Fn.new { |poly, x| poly.map { |coef| coef * x }.toList }
// Divide a polynomial by a scalar.
var scalarDivide = Fn.new { |poly, x| scalarMultiply.call(poly, 1/x) }
// Returns the Lagrange interpolating polynomial which passes through a list of points.
var lagrange = Fn.new { |pts|
var c = pts.count
var polys = List.filled(c, null)
for (i in 0...c) {
var poly = [1]
for (j in 0...c) {
if (i == j) continue
poly = multiply.call(poly, [-pts[j].x, 1])
}
var d = Math.evalPoly(poly[-1..0], pts[i].x)
polys[i] = scalarDivide.call(poly, d)
}
var sum = [0]
for (i in 0...c) {
polys[i] = scalarMultiply.call(polys[i], pts[i].y)
sum = add.call(sum, polys[i])
}
return sum
}
var pts = [
Point.new(1, 1),
Point.new(2, 4),
Point.new(3, 1),
Point.new(4, 5)
]
var lip = lagrange.call(pts)
Fmt.pprint("$f", lip[-1..0], "", "x")
- Output:
2.166667x³ - 16.000000x² + 35.833333x - 21.000000