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

Lagrange Method, Python, 199 236 202 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))])

Edit: To adhere to general rules, we have to remove references to non-builtin functions (at least without importing them). I have rewritten it to only use standard library modules and also included the imports (thanks to @TheFifthMarshal & @Tech-Freak), while deleting extraneous whitespace, bumping it up to 208 bytes. To elaborate, the input d takes the form of a 2D iterable (an nx2 matrix), where each row is a measurement, the first column is \$t\$ and the second column is \$x\$. The output is an anonymous function which can take any floating point or integer number and give the approximate result.

from functools import*
mul=lambda a,b:a*b
def l(d):
 l=lambda i,x:reduce(mul,[(x-d[j][0])/(d[i][0]-d[j][0])for j in range(len(d))if i!=j],1)
 return lambda x:sum([l(i,x)*d[i][1]for i in range(len(d))])

Usage Example:

>>> t = np.arange(0, 5, 0.5)
>>> t2 = np.arange(0, 5, 0.01)
>>> x = t**3
>>> dp = np.vstack((t, x)).T
>>> poly = l(dp)
>>> plt.plot(t2, np.vectorize(poly)(t2), label="output")
>>> plt.scatter(t, x, label="input")
>>> plt.legend()