Polynomial Interpolation

Question

Write a program that performs Polynomial Interpolation using true arbitrary precision rational numbers. The input looks like this:

f(1) = 2/3
f(2) = 4/5
f(3) = 6/7
...

You may assume that there's exactly one whitespace before and after the = sign, all the numbers are either fractions or integers. You may also assume, that all fraction in the input are already irreducible.

No error checking is needed, you may assume, that the input is valid and no x is doubled in the f(x).

The output should be in a LaTeX compatible form, the emitted LaTeX code should yield the same graphical representation as the output given here.

f(x) = 123x^2 + \frac{45}{2}x + \frac{7}{4}

The fraction must be reduced as much as possible, eg. something like \frac{2}{4} is not allowed. If the number is integer, don't use a fraction.

Special rules:

Your program should...

• work for polynomials up to degree 12
• complete in less then 1 minute for reasonable input
• not use any functions that do the whole calculation for you
• output the polynomial of smallest possible degree

Testcases:

The given testcases are just for clarification. Your program should yield correct result for all correct inputs.

Input

f(1) = 2/3
f(2) = 4/5
f(3) = 6/7

Output

f(x) = - \frac{4}{105}x^2
       + \frac{26}{105}x
       + \frac{16}{35}

Input

f(-12) = 13/2
f(5/3) = 3/5
f(13) = -6
f(1/5) = -3/4

Output

f(x) = - \frac{2186133}{239455744}x^3
       + \frac{2741731}{149659840}x^2
       + \frac{26720517}{29201920}x
       - \frac{279464297}{299319680}

Input

f(4/3) = 617/81
f(2) = 20/3
f(-8/3) = 6749/81
f(-5) = 7367/12
f(0) = 23/3

Output

f(x) = \frac{1}{2}x^4
     - 2x^3
     + \frac{7}{4}x^2
     + \frac{23}{3}

Input

f(0) = 5
f(1) = 7
f(2) = 9
f(3) = 11
f(4) = 13

Output

f(x) = 2x
     + 5

Input

f(1/2) = -1/2
f(-25) = -1/2
f(-54/12) = -1/2

Output

f(x) = -\frac{1}{2}

Answer 1

J + sh

J script:

i=:0".(1!:1)3
i=:((-:#i),2)$i
c=:|.(%.(x:((i.#i)^~])"0({."1 i)))(+/ .*)(|:{:"1 i)
(":(-.0=c)#(c,.i.-#c))(1!:2)2

sh script:

echo -n 'f(x) = '
tr -d 'f()=' | tr /\\n- r' '_  | ./polyint.ijs | sed -e 's/^/+/;s/_/-/;s/\([0-9]*\)r\([0-9]*\)/\\frac{\1}{\2}/;s/ \([0-9]*$\)/x^\1/;s/\^1//;s/x^0//;s/+\(.*-.*\)/\1/'

Run the sh script:

./pol-int.sh
f(1/2) = -1/2
f(-25) = -1/2
f(-54/12) = -1/2

f(x) = -\frac{1}{2}

.

./pol-int.sh
f(4/3) = 617/8
f(2) = 20/3
f(-8/3) = 6749/81
f(-5) = 7367/12
f(0) = 23/3

f(x) = -\frac{37745}{14592}x^4
       -\frac{853249}{43776}x^3
     +  \frac{57809}{7296}x^2
     + \frac{225205}{2736}x
     +  \frac{23}{3}

Answer 2

Perl (569 characters)

use Math'BigInt;sub r{($u,$v)=@_;$v?r($v,$u%$v):$u}sub c{new Math'BigInt$_[0]}$a=@c=<>;for(@c){m!(-?\d+)/?(\d*). = (-?\d+)/?(\d*)!;$j{$_,$i}=$1**c$_,$k{$_,$i|=0}=($2||1)**c$_ for 0..$a;$j{$a,$i}=c$3;$k{$a,$i++}=c$4||1}for$p(0..$a-1){for$y(0..$p-1,$p+1..$a-1){$n=$j{$p,$y}*$k{$p,$p};$o=$k{$p,$y}*$j{$p,$p};$j{$_,$y}=$j{$_,$y}*$k{$_,$p}*$o-$k{$_,$y}*$j{$_,$p}*$n,$k{$_,$y}*=$k{$_,$p}*$o for 0..$a}}print"f(x)=";for(1..$a){$s=r$t=$j{$a,$p=$a-$_}*$k{$p,$p},$w=$k{$a,$p}*$j{$p,$p};$u=abs$t,print$t>0?"$z":'-',($z='+',$w/=$s)-1?"\\frac{$u}{$w}":$u,$p>1?"x^$p":x x$p if$t/=$s}

Detailed explanation:

use Math'BigInt;

# Subroutine to calculate gcd of two numbers
sub r{($u,$v)=@_;$v?r($v,$u%$v):$u}

# Subroutine to create BigInts
sub c{new Math'BigInt$_[0]}

# Read input
# Throughout, $a is the number of equations.
$a=@c=<>;

# Initialises the $a+1 × $a matrix with all the values.
# $j{$x,$y} contains the numerator, $k{$x,$y} the denominator.
for(@c)
{
    m!(-?\d+)/?(\d*). = (-?\d+)/?(\d*)!;

    # Puzzle for the reader: why is $i|=0 in the second one,
    # not the first one? Answer at the bottom!
    $j{$_,$i}=$1**c$_,$k{$_,$i|=0}=($2||1)**c$_ for 0..$a;
    $j{$a,$i}=c$3;
    $k{$a,$i++}=c$4||1
}

# Generates the matrix echelon form.
# Basically, it works like this:
for$p(0..$a-1)
{
    # For each element along the diagonal {$p,$p}, set all the values above and
    # below it to 0 by adding a multiple of row $p to each of the other rows.
    for$y(0..$p-1,$p+1..$a-1)
    {
        # So we need to multiply the row $p by the value of {$p,$y}/{$p,$p}
        # (stored in $n/$o) and then subtract that from row $y.
        $n=$j{$p,$y}*$k{$p,$p};
        $o=$k{$p,$y}*$j{$p,$p};
            $j{$_,$y}=$j{$_,$y}*$k{$_,$p}*$o-$k{$_,$y}*$j{$_,$p}*$n,
            $k{$_,$y}*=$k{$_,$p}*$o
        for 0..$a
    }
}

# Outputs the result
print"f(x)=";
for(1..$a)
{
    # Note this sets $p = $a-$_. $p is the power of x.
    # We need to divide {$a,$p} by {$p,$p}. Store the result in $t/$w.
    # We also need to put the fraction in lowest terms, so calculate the gcd.
    $s=r$t=$j{$a,$p=$a-$_}*$k{$p,$p},$w=$k{$a,$p}*$j{$p,$p};

    # Output this term only if the numerator ($t) is non-zero.
    # Output a plus sign only if this isn’t the first term.
    # Output a fraction only if the denomator ($w) isn’t 1.
        $u=abs$t,print$t>0?"$z":'-',
        ($z='+',$w/=$s)-1?"\\frac{$u}{$w}":$u,$p>1?"x^$p":x x$p
    if$t/=$s
}

# Answer to the puzzle buried in the code above:
# It’s because the second part is passed as a second argument to c,
# hence it is evaluated before the first part.

Comments

• I am sure there is a module for matrix manipulation that provides a function for echelon form. I specifically didn’t use that (didn’t even search for one) because I assume it is the point of this contest to do it yourself. It is more interesting that way, too. Of course the same could be said about BigInt, but then I suspect nobody will attempt this challenge...

Edits

• (630 → 585) Realised I can do the echelon form in one loop instead of two. Add explanation as comments in code.

• (585 → 583) Just discovered the package syntax that lets me use ' instead of ::.

• (583 → 573) Some more microgolfing

• (573 → 569) Shorter regular expression to parse input

Answer 3

TI-Basic(83/84): 109 characters

Technically 109 characters, TI-Basic counts dim(, For(, ->, rref(, [A], and lists as "one character".

The input is formatted into L1 and L2, in (x,y) pairs [e.x. L1=(1,2,3,4), L2=(2,3,5,7)].

{1,1}->dim([A]
{dim(L1),dim(L2)+1}->dim([A]
For(A,1,dim(L1)
For(B,dim(L1)-1,0,-1
L1(A)^B->[A](A,dim(L1)-B
End
L2(A->[A](A,dim(L1)+1
End
rref([A]->[A]
{0}->L3
For(A,1,dim(L1)
[A](A,dim(L1)+1->L3(A
End
Disp L3

Answer 4

Lagrange Method, Python, 199 bytes

A little late, but...

def lagrange(dp):
l = lambda i: lambda x: (prod([(x - dp[j][0]) / (dp[i][0] - dp[j][0]) for j in range(len(dp)) if i != j]))
return lambda x: sum([l(i)(x) * dp[i][1] for i in range(len(dp))])

Answer 5

l=lambda D,i:lambda x:prod((x-Dj[0])/(D[i][0]-Dj[0])for j,Dj in enumerate(D)if i!=j)
L=lambda D:lambda x:sum(l(D,i)(x)*Di[1]for i,Di in enumerate(D))

Just a shortened version of Fred Freys code. Note that one could probably skip passing D to l, since it can just pull it from the outer scope. As you can probably do the same with i here, we could even shave off one lambda. I'll test it some day.