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}

• Joey: Why are you talking about real numbers if all you're ever using are rational numbers?
• FUZxxl: Sorry. My english is poor. Yes, only use rational numbers. The results have to be exact.
• Eelvex: In the first testcase, are dots (...) really part of the input?
• FUZxxl: @Eelvex: Nope. Fixed.
• Timwi: The output for the third testcase is wrong. The correct answer is -\frac{37745}{14592}x^4 - \frac{853249}{43776}x^3 + \frac{57809}{7296}x^2 + \frac{225205}{2736}x + \frac{23}{3}. I suspect the input was intended to be something different :)

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}

• FUZxxl: You don't have to create the exact same sourcecode formatting. In the LaTeX output. It should just yield the same graphical representation after running through LaTeX. Feel free to save some chars.
• Timwi: I can’t read J, but from the short length I take it that means J has a built-in function for matrix echelon form?
• Eelvex: @Timwi: No, but I use the built-in "inverse matrix". J is very terse though: even if I were to implement "invert matrix" it would be a few characters long.

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}=c3ドル;$k{$a,$i++}=c4ドル||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}=c3ドル;
 $k{$a,$i++}=c4ドル||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

• FUZxxl: I keep getting compilation errors: ideone.com/LoB2T
• Timwi: @FUZxxl: Thanks for pointing that out. There was just a missing space. Fixed now.

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()

Lagrange interpolation of f(x) = x^3

• user62131: You probably don't need all that whitespace around the operators, do you?
• The Fifth Marshal: This answer appears to be invalid for reasons I explained at length in the comments on a golfed version of it
• fpf3: Updating now. @TheFifthMarshal
• The Fifth Marshal: 208 bytes
• fpf3: Awesome. I've updated the answer and added thanks to you & tech-freak.