# Polynomial Interpolation

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:

• 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}

• Why are you talking about real numbers if all you're ever using are rational numbers? – Joey Mar 21 '11 at 16:17
• Sorry. My english is poor. Yes, only use rational numbers. The results have to be exact. – FUZxxl Mar 21 '11 at 16:22
• In the first testcase, are dots (...) really part of the input? – Eelvex Mar 21 '11 at 18:01
• @Eelvex: Nope. Fixed. – FUZxxl Mar 21 '11 at 18:03
• 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 :) – Timwi Mar 21 '11 at 20:11

## 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}  • 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. – FUZxxl Mar 21 '11 at 21:08 • 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? – Timwi Mar 21 '11 at 21:36 • @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. – Eelvex Mar 21 '11 at 22:06 # Perl (569 characters) use Math'BigInt;sub r{($u,$v)=@_;$v?r($v,$u%$v):$u}sub c{new Math'BigInt$_}$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$_} # 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.


• 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

• I keep getting compilation errors: ideone.com/LoB2T – FUZxxl Mar 21 '11 at 21:24
• @FUZxxl: Thanks for pointing that out. There was just a missing space. Fixed now. – Timwi Mar 21 '11 at 21:30

# 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

• This doesn't use rationals or LaTeX form. – lirtosiast Oct 6 '15 at 4:31

# Lagrange Method, Python, 199 bytes

A little late, but...

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

• You probably don't need all that whitespace around the operators, do you? – user62131 Mar 14 '17 at 4:31
l=lambda D,i:lambda x:prod((x-Dj)/(D[i]-Dj)for j,Dj in enumerate(D)if i!=j)
L=lambda D:lambda x:sum(l(D,i)(x)*Difor 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.