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


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


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


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


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


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


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


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


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


f(x) = 2x
     + 5


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


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

3 Answers 3


J + sh

J script:

c=:|.(%.(x:((i.#i)^~])"0({."1 i)))(+/ .*)(|:{:"1 i)

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:

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

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


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{57809}{7296}x^2
     + \frac{225205}{2736}x
     +  \frac{23}{3}
  • \$\begingroup\$ 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. \$\endgroup\$
    – FUZxxl
    Commented Mar 21, 2011 at 21:08
  • \$\begingroup\$ 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? \$\endgroup\$
    – Timwi
    Commented Mar 21, 2011 at 21:36
  • \$\begingroup\$ @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. \$\endgroup\$
    – Eelvex
    Commented Mar 21, 2011 at 22:06

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.

# Initialises the $a+1 × $a matrix with all the values.
# $j{$x,$y} contains the numerator, $k{$x,$y} the denominator.
    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;

# Generates the matrix echelon form.
# Basically, it works like this:
    # 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.
        # 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.
        for 0..$a

# Outputs the result
    # 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.

    # 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.
        ($z='+',$w/=$s)-1?"\\frac{$u}{$w}":$u,$p>1?"x^$p":x x$p

# 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...


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

  • \$\begingroup\$ I keep getting compilation errors: ideone.com/LoB2T \$\endgroup\$
    – FUZxxl
    Commented Mar 21, 2011 at 21:24
  • \$\begingroup\$ @FUZxxl: Thanks for pointing that out. There was just a missing space. Fixed now. \$\endgroup\$
    – Timwi
    Commented Mar 21, 2011 at 21:30

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

  • 2
    \$\begingroup\$ You probably don't need all that whitespace around the operators, do you? \$\endgroup\$
    – user62131
    Commented Mar 14, 2017 at 4:31
  • \$\begingroup\$ This answer appears to be invalid for reasons I explained at length in the comments on a golfed version of it \$\endgroup\$ Commented Nov 6, 2022 at 1:12
  • \$\begingroup\$ Updating now. @TheFifthMarshal \$\endgroup\$
    – fpf3
    Commented Nov 7, 2022 at 18:13
  • \$\begingroup\$ 208 bytes \$\endgroup\$ Commented Nov 7, 2022 at 18:26
  • \$\begingroup\$ Awesome. I've updated the answer and added thanks to you & tech-freak. \$\endgroup\$
    – fpf3
    Commented Nov 7, 2022 at 18:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.