Interpolación polinómica

12

Escriba un programa que realice la interpolación polinómica usando números racionales de precisión arbitraria verdadera. La entrada se ve así:

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

Puede suponer que hay exactamente un espacio en blanco antes y después del =signo, todos los números son fracciones o enteros. También puede suponer que todas las fracciones en la entrada ya son irreductibles.

No es necesario verificar errores, puede suponer que la entrada es válida y que no x se duplica en f (x).

La salida debe estar en una forma compatible con LaTeX, el código emitido de LaTeX debe producir la misma representación gráfica que la salida dada aquí.

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

La fracción debe reducirse tanto como sea posible, por ejemplo. Algo así \frac{2}{4} no está permitido. Si el número es entero, no use una fracción.

Reglas especiales:

Tu programa debería ...

  • trabajar para polinomios hasta grado 12
  • completar en menos de 1 minuto para una entrada razonable
  • no use ninguna función que haga todo el cálculo por usted
  • salida del polinomio de menor grado posible

Casos de prueba:

Los casos de prueba dados son solo para aclaración. Su programa debe producir el resultado correcto para todas las entradas correctas.

Entrada

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

Salida

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

Entrada

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

Salida

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

Entrada

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

Salida

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

Entrada

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

Salida

f (x) = 2x
     + 5

Entrada

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

Salida

f (x) = - \ frac {1} {2}
FUZxxl
fuente
¿Por qué estás hablando de números reales si todo lo que estás usando son números racionales?
Joey
Lo siento. Mi Inglés es pobre. Sí, solo usa números racionales. Los resultados tienen que ser exactos.
FUZxxl
En el primer caso de prueba, ¿son los puntos ( ...) realmente parte de la entrada?
Eelvex
@Eelvex: No. Fijo.
FUZxxl
La salida para el tercer caso de prueba es incorrecta. La respuesta correcta es -\frac{37745}{14592}x^4 - \frac{853249}{43776}x^3 + \frac{57809}{7296}x^2 + \frac{225205}{2736}x + \frac{23}{3}. Sospecho que la entrada tenía la intención de ser algo diferente :)
Timwi

Respuestas:

3

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/'

Ejecute el script sh:

./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}
Eelvex
fuente
No tiene que crear exactamente el mismo formato de código fuente. En la salida de LaTeX. Simplemente debería producir la misma representación gráfica después de ejecutar LaTeX. Siéntase libre de guardar algunos caracteres.
FUZxxl
No puedo leer J, pero por la corta longitud que tomo, ¿significa que J tiene una función incorporada para la forma escalonada de matriz?
Timwi
@Timwi: No, pero uso la "matriz inversa" incorporada. Sin embargo, J es muy conciso: incluso si tuviera que implementar "matriz invertida", tendría unos pocos caracteres de longitud.
Eelvex
3

Perl (569 caracteres)

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}

Explicación detallada:

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.

Comentarios

  • Estoy seguro de que hay un módulo para la manipulación matricial que proporciona una función para la forma escalonada. Específicamente no utilicé eso (ni siquiera busqué uno) porque supongo que es el objetivo de este concurso hacerlo usted mismo. También es más interesante de esa manera. Por supuesto, se podría decir lo mismo sobre BigInt, pero sospecho que nadie intentará este desafío ...

Ediciones

  • (630 → 585) Me di cuenta de que puedo hacer la forma escalonada en un bucle en lugar de dos. Agregue explicaciones como comentarios en el código.

  • (585 → 583) Acabo de descubrir la sintaxis del paquete que me permite usar en 'lugar de ::.

  • (583 → 573) Algunos microgolfing más

  • (573 → 569) Expresión regular más corta para analizar la entrada

Timwi
fuente
Sigo recibiendo errores de compilación: ideone.com/LoB2T
FUZxxl
@FUZxxl: Gracias por señalar eso. Solo faltaba un espacio. Corregido ahora.
Timwi
3

TI-Basic (83/84): 109 caracteres

Técnicamente 109 caracteres, TI-Basic cuenta dim (, For (, ->, rref (, [A], y se enumera como "un carácter".

La entrada se formatea en L1 y L2, en pares (x, y) [ex 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
beary605
fuente
1
Esto no usa racionales ni forma LaTeX.
lirtosiast
1

Método Lagrange, Python, 199 bytes

Un poco tarde, pero ...

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))])
Fred Frey
fuente
1
Probablemente no necesite todo ese espacio en blanco alrededor de los operadores, ¿verdad?
0
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))

Solo una versión abreviada del código de Fred Freys. Tenga en cuenta que probablemente se podría omitir pasar D a l, ya que solo puede extraerlo del alcance externo. Como probablemente pueda hacer lo mismo con yo aquí, incluso podríamos afeitarnos una lambda. Lo probaré algún día.

Teck-freak
fuente