Demasiado rápido, demasiado Fourier: FFT Code Golf

48

Implemente la Transformada rápida de Fourier en la menor cantidad de caracteres posible.

Reglas:

  • La solución más corta gana
  • Se puede suponer que la entrada es una matriz 1D cuya longitud es una potencia de dos.
  • Puede usar el algoritmo de su elección, pero la solución debe ser una Transformada rápida de Fourier, no solo una ingenua Transformada discreta de Fourier (es decir, debe tener un costo de cálculo asintótico de )O(NlogN)

Editar:

  • el código debe implementar la Transformación rápida rápida de Fourier estándar, cuya forma se puede ver en la ecuación (3) de este artículo de Wolfram ,

    ingrese la descripción de la imagen aquí

  • El uso de una función FFT de una biblioteca estándar preexistente o un paquete de estadísticas no está permitido. El reto aquí es sucinta aplicar el algoritmo FFT en sí.
jakevdp
fuente
3
Esto está poco especificado. Como mínimo, debe definir los factores de normalización, y también debe ser consciente de que cualquier ambigüedad será malinterpretada deliberadamente. Por ejemplo, ¿está satisfecho "Implementar" con la respuesta " FFT(3 caracteres): está en la biblioteca estándar"? Algunos casos de prueba también serían buenos.
Peter Taylor
¿Importa el orden de los elementos de salida, es decir, necesitamos implementar la descodificación de bits invertidos o podemos dejar la salida en orden codificado?
Paul R
Ver las ediciones de las reglas. La salida debe ser una lista / matriz con valores ordenados de acuerdo con los índices en la expresión DFT estándar, mencionada anteriormente.
jakevdp
2
¿Puedes publicar algunas entradas y salidas de ejemplo para que podamos probar nuestras implementaciones?
FUZxxl
2
El título debería haber sido "Rápido y Fourier-s" (Rápido y Furioso).
clismique

Respuestas:

12

Mathematica, 95 bytes

Otra implementación de Cooley – Tukey FFT con la ayuda de @ chyaong .

{n=Length@#}~With~If[n>1,Join[+##,#-#2]&[#0@#[[;;;;2]],#0@#[[2;;;;2]]I^Array[-4#/n&,n/2,0]],#]&

Sin golf

FFT[x_] := With[{N = Length[x]},
  If[N > 1,
    With[{a = FFT[ x[[1 ;; N ;; 2]] ], 
          b = FFT[ x[[2 ;; N ;; 2]] ] * Table[E^(-2*I*Pi*k/N), {k, 0, N/2 - 1}]},
      Join[a + b, a - b]],
    x]]
millas
fuente
1
Creo #[[;;;;2]]==#[[1;;N;;2]]y [[2;;;;2]]==[[2;;N;;2]].
chyanog
1
101 caracteres:With[{L=Length@#},If[L>1,Join[+##,#-#2]&[#0@#[[;;;;2]],#0@#[[2;;;;2]]E^(-2I*Pi(Range[L/2]-1)/L)],#]]&
chyanog
Bien, puedes condensar otra función anónima dentro de ella sin entrar en conflicto con la recursiva. También aprendí que Parte completa los índices faltantes. Podemos llevarlo más lejos con Unicode.
millas
9

J, 37 bytes

_2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#

Una mejora después de unos años. Todavía usa el algoritmo Cooley-Tukey FFT.

Guardado 4 bytes usando e πi = -1, gracias a @ Leaky Nun .

Pruébalo en línea!

Uso

   f =: _2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#
   f 1 1 1 1
4 0 0 0
   f 1 2 3 4
10 _2j2 _2 _2j_2
   f 5.24626 3.90746 3.72335 5.74429 4.7983 8.34171 4.46785 0.760139
36.9894 _6.21186j0.355661 1.85336j_5.74474 7.10778j_1.13334 _0.517839 7.10778j1.13334 1.85336j5.74474 _6.21186j_0.355661

Explicación

_2&(0((+,-)]%_1^i.@#%#)&$:/@|:]\)~1<#  Input: array A
                                    #  Length
                                  1<   Greater than one?
_2&(                            )~     Execute this if true, else return A
_2                            ]\         Get non-overlapping sublists of size 2
    0                       |:           Move axis 0 to the end, equivalent to transpose
                          /@             Reduce [even-indexed, odd-indexed]
                       &$:               Call recursively on each 
                   #                     Get the length of the odd list
                i.@                      Range from 0 to that length exclusive
                    %#                   Divide each by the odd length
             _1^                         Compute (-1)^x for each x
           ]                             Get the odd list
            %                            Divide each in that by the previous
       +                                 Add the even values and modified odd values
         -                               Subtract the even values and modified odd values
        ,                                Join the two lists and return
millas
fuente
1
Vea también esto: blog.ndpar.com/2014/10/11/dft-j
FUZxxl
9

Pitón, 166 151 150 caracteres

Esto usa el algoritmo radix-2 Cooley-Tukey FFT

from math import*
def F(x):N=len(x);t=N<2or(F(x[::2]),F(x[1::2]));return N<2and x or[
a+s*b/e**(2j*pi*n/N)for s in[1,-1]for(n,a,b)in zip(range(N),*t)]

Probar el resultado

>>> import numpy as np
>>> x = np.random.random(512)
>>> np.allclose(F(x), np.fft.fft(x))
True
jakevdp
fuente
1
2 cosas: generalmente es mejor usar from x import*, y sum(([x for x in y] for y in z),[])es más largo que [x for y in z for x in y].
Boothby
1
Gracias, ¡eso ahorra 15 caracteres! 11 más y es un tweet.
jakevdp
Oh, eso es definitivamente posible. A menudo, cuando encuentras una mejora, una vieja se convierte en un obstáculo.
Boothby
5

Python 3: 140 134 113 caracteres

Versión corta: corta y dulce, cabe en un tweet (gracias a las millas ):

from math import*
def f(v):
 n=len(v)
 if n<2:return v
 a,b=f(v[::2])*2,f(v[1::2])*2;return[a[i]+b[i]/1j**(i*4/n)for i in range(n)]

(En Python 2, /se trunca la división cuando ambos lados son enteros. Por lo tanto, reemplazamos (i*4/n)por (i*4.0/n), lo que aumenta la longitud a 115 caracteres).

Versión larga: más claridad en la parte interna del clásico Cooley-Tukey FFT:

import cmath
def transform_radix2(vector):
    n = len(vector)
    if n <= 1:  # Base case
        return vector
    elif n % 2 != 0:
        raise ValueError("Length is not a power of 2")
    else:
        k = n // 2
        even = transform_radix2(vector[0 : : 2])
        odd  = transform_radix2(vector[1 : : 2])
        return [even[i % k] + odd[i % k] * cmath.exp(i * -2j * cmath.pi / n) for i in range(n)]
Nayuki
fuente
1
Acortado a 113 bytes usandoe^(-2j * pi * i / n) = (-1)^(2 * i / n) = (1j)^(4 * i / n)
millas
@miles Observación muy impresionante, ¡gracias! Después de haber implementado repetidamente DFT durante más de una década, me obsesioné con sin / cos / exp y olvidé que se pueden usar poderes simples de i. Edité mi respuesta para incorporar la nueva visión y acreditarte.
Nayuki
5

R: 142 133 99 95 bytes

¡Gracias a @Giuseppe por ayudarme a reducir 32 36 bytes!

f=function(x,n=sum(x|1),y=1:(n/2)*2)`if`(n>1,f(x[-y])+c(b<-f(x[y]),-b)*exp(-2i*(y/2-1)*pi/n),x)

Un truco adicional aquí es usar los argumentos predeterminados de la función principal para instanciar algunas variables.
El uso sigue siendo el mismo:

x = c(1,1,1,1)
f(x)
[1] 4+0i 0+0i 0+0i 0+0i

Versión de 4 años de antigüedad en 133 bytes:

f=function(x){n=length(x);if(n>1){a=Recall(x[seq(1,n,2)]);b=Recall(x[seq(2,n,2)]);t=exp(-2i*(1:(n/2)-1)*pi/n);c(a+b*t,a-b*t)}else{x}}

Con hendiduras:

f=function(x){
    n=length(x)
    if(n>1){
        a=Recall(x[seq(1,n,2)])
        b=Recall(x[seq(2,n,2)])
        t=exp(-2i*(1:(n/2)-1)*pi/n)
        c(a+b*t,a-b*t)
        }else{x}
    }

Utiliza también el algoritmo Cooley-Tukey. Los únicos trucos aquí son el uso de la función Recallque permite la recursividad y el uso de la vectorización R que acorta en gran medida el cálculo real.

Uso:

x = c(1,1,1,1)
f(x)
[1] 4+0i 0+0i 0+0i 0+0i
plannapus
fuente
1
Cuatro años después, y lo reducimos a 101 bytes . No estoy 100% seguro de por qué lo usaste, Recallya que es una función con nombre, pero bueno, ¡es fácil jugar al golf en retrospectiva! :) +1, muy agradable.
Giuseppe
Sí, Recallahora es innecesario, de hecho. Me di cuenta de eso hace unos meses pero era demasiado vago para cambiarlo :) Lo modificaré.
plannapus
¡Muy agradable! ¡Exprimí otros 4 bytes! , poniendo esto a la par con Mathematica.
Giuseppe
¡Gracias! Pensé en poner yallí, pero no me di cuenta de que también podría usarse para la exp(...)parte.
plannapus
4

Python, 134

Esto toma mucho prestado de la solución de jakevdp, por lo que configuré este en una wiki comunitaria.

from math import*
F=lambda x:x*(len(x)<2)or[a+s*b/e**(2j*pi*n/len(x))for s in(1,-1)for n,(a,b)in
enumerate(zip(F(x[::2]),F(x[1::2])))]

Cambios:

-12 caracteres: matar t.

def F(x):N=len(x);t=N<2or(F(x[::2]),F(x[1::2]));return ... in zip(range(N),*t)]
def F(x):N=len(x);return ... in zip(range(N),F(x[::2]),F(x[1::2]))]

-1 char: truco de exponente, x*y**-z == x/y**z (esto podría ayudar a otros)

...[a+s*b*e**(-2j*pi*n/N)...
...[a+s*b/e**(2j*pi*n/N)...

-2 caracteres: reemplazar andcon*

...return N<2and x or[
...return x*(N<2)or[

+1 char: lambdaize, matandoN

def F(x):N=len(x);return x*(N<2)or[a+s*b/e**(2j*pi*n/N) ... zip(range(N) ...
F=lambda x:x*(len(x)<2)or[a+s*b/e**(2j*pi*n/len(x)) ... zip(range(len(x)) ...

-2 char: usar en enumeratelugar dezip(range(len(

...for(n,a,b)in zip(range(len(x)),F(x[::2]),F(x[1::2]))]
...for n,(a,b)in enumerate(zip(F(x[::2]),F(x[1::2])))]
boothby
fuente
Sin embargo, creo que esto ya no es una transformación rápida de Fourier ... al "matar t" agregaste algunos cálculos innecesarios que lo mueven de O [N log (N)] a O [N ^ 2]
jakevdp
Parece que no puedo rechazar mi propia publicación. Tienes razón, cambié el orden de los bucles y eliminé el rendimiento. Dejaré esto por ahora, en caso de que encuentre una manera de solucionarlo.
stand
101 bytes conf=lambda x:x*(len(x)<2)or[u+v/1j**(4*i/len(x))for i,(u,v)in enumerate(zip(f(x[::2])*2,f(x[1::2])*2))]
millas
Se puede reemplazar for s in(1,-1)forcon for s in 1,-1foro incluso, si el orden no importa, for s in-1,1for.
Jonathan Frech
4

C, 259

typedef double complex cplx;
void fft(cplx buf[],cplx out[],int n,int step){
if(step < n){
fft(out, buf,n, step * 2);
fft(out+step,buf+step,n,step*2);
for(int i=0;i<n;i+=2*step){
cplx t=cexp(-I*M_PI*i/n)*out[i+step];
buf[i/2]=out[i]+t;
buf[(i+n)/2]=out[i]-t;
}}}

El problema es que tales implementaciones son inútiles, y el algoritmo directo es MUCHO más rápido.

sanaris
fuente
2
Puede eliminar algunos espacios en blanco más para obtener una menor cantidad de caracteres, por ejemplo, step < nse puede cambiar a step<ny step * 2se puede cambiar a step*2.
ProgramFOX
2
todas las variables y funciones y typedefs deben tener nombres de una letra para guardar muchos caracteres
2
Hiciste que alguien sugiriera algunas mejoras para esto. Míralos
Justin
1
Puede eliminar todas las nuevas líneas, las nuevas líneas son inútiles en C
TuxCrafting
@ TùxCräftîñg No todas las líneas nuevas son inútiles. Son necesarios para directivas como #include, #define, #if, etc.
Nayuki
3

Matlab, 128 118 107 102 101 94 93 bytes

EDIT6: ¡gracias @algmyr por otro byte!

function Y=f(Y);
n=numel(Y);
k=2:2:n;
if k;
   c=f(Y(k-1));
   d=f(Y(k)).*i.^(2*(2-k)/n);
   Y=[c+d;c-d];
end

EDIT5: Todavía se está acortando :) gracias a @sanchises

function Y=f(Y)
n=numel(Y);
k=2:2:n;
if k;
   c=f(Y(k-1));
   d=f(Y(k)).*(-1).^((2-k)/n);
   Y=[c+d;c-d];
end

EDITAR4: Sí, -1 carácter más (también podría haberlo hecho sin el k):

function Y=f(Y)
n=numel(Y);
if n>1;
   k=2:2:n;
   c=f(Y(k-1));
   d=f(Y(k)).*(-1).^((k/2-1)*2/n)';
   Y=[c+d;c-d];
end

EDIT2 / 3: ¡Gracias por @sanchises por nuevas mejoras!

function Y=f(Y)
n=numel(Y);  
if n>1;
   c=f(Y(1:2:n));
   d=f(Y(2:2:n)).*(-1).^(-(0:n/2-1)*2/n).';
   Y=[c+d;c-d]; 
end

EDITAR: podría hacer algunas mejoras, y noté que la constante de escala no es necesaria.

Esta es la versión ampliada, el recuento de caracteres es válido si elimina las nuevas líneas / espacios. (Funciona solo para vectores de columna).

function y=f(Y)
n=numel(Y);  
y=Y;
if n>1;
   c=f(Y(1:2:n));
   d=f(Y(2:2:n));
   n=n/2;
   d=d.*exp(-pi*i*(0:n-1)/n).';
   y=[c+d;c-d]; 
end
falla
fuente
Consejo: puede combinar las dos d=líneas en una sola: m=n/2;d=f(Y(2:2:n)).*exp(-pi*i*(0:m-1)/m).';. Además, considere cambiar y=f(Y)a Y=f(Y)y quitar la línea 3 (y promete que nunca se va a hacer que fuera de código de golf)
Sanchises
¡Oh gracias! ¿ function Y = f(Y)Tiene alguna desventaja aparte de la imposibilidad de leer?
flawr
Bueno, MATLAB nunca se quejará de un valor de retorno, incluso si Y nunca cambia. Sin embargo, es un poco más rápido, así que supongo que no es tan malo después de todo para algunos propósitos (es decir, una función que casi nunca cambia la variable de entrada)
Sanchises
Ahora, para afeitarse más: m=n/2podría eliminarse y mreemplazarse por n/2y n*2respectivamente. Y luego, creo firmemente, que el programa es tan corto como podría ser en MATLAB.
Sanchises
1
Y luego, creo firmemente, que el programa es tan corto como podría ser en MATLAB. - Sanchises 8 de marzo de 15 a 21:05 últimas palabras famosas ...
Sanchises
2

Gelatina, 31 30 28 26 bytes , no competitiva

LḶ÷$N-*×,N$+ḷF
s2Z߀ç/µ¹Ṗ?

Jelly fue creada después de este desafío, por lo que no es competitiva.

Esto usa el algoritmo recursivo Cooley-Tukey radix-2. Para una versión sin golf, vea mi respuesta en Mathematica.

Pruébelo en línea o Verifique múltiples casos de prueba .

Explicación

LḶ÷$N-*×,N$+ḷF  Helper link. Input: lists A and B
L               Get the length of A
   $            Operate on that length
 Ḷ                Make a range [0, 1, ..., length-1]
  ÷               Divide each by length
    N           Negate each
     -          The constant -1
      *         Compute -1^(x) for each x in that range
       ×        Multiply elementwise between that range and B, call it B'  
          $     Operate on that B'
         N        Negate each
        ,         Make a list [B', -B']
            ḷ   Get A
           +    Add vectorized, [B', -B'] + A = [A+B', A-B']
             F  Flatten that and return

s2Z߀ç/µ¹Ṗ?  Main link. Input: list X
         Ṗ   Curtail - Make a copy of X with the last value removed
          ?  If that list is truthy (empty lists are falsey)
       µ       Parse to the left as a monad
s2             Split X into sublists of length 2
  Z            Transpose them to get [even-index, odd-index]
   ߀          Call the main link recursively on each sublist
     ç/        Call the helper link as a dyad on the sublists and return
             Else
        ¹      Identity function on X and return
millas
fuente
2

C (gcc) , 188 186 184 183 bytes

#define d(a,b,c)f(a,b,c,1,0)
f(a,b,c,n,k)_Complex*a,*b;{_Complex z[c];*b=*a;if(n<c)for(f(a,z,c,n*2),f(a+n,z+n,c,n*2);k<c;k+=n*2)b[k+c>>1]=z[k]*2-(b[k/2]=z[k]+z[k+n]/cpow(1i,2.*k/c));}

Pruébalo en línea!

Ligeramente menos golf

#define d(a,b,c)f(a,b,c,1,0)
f(a,b,c,n,k)_Complex*a,*b;{
  _Complex z[c];
  *b=*a;
  if(n<c)
    for(f(a,z,c,n*2),f(a+n,z+n,c,n*2);k<c;k+=n*2)
      b[k+c>>1]=z[k]*2-(b[k/2]=z[k]+z[k+n]/cpow(1i,2.*k/c));
}
techo
fuente
1

Pari / GP, 76 caracteres

X(v)=my(t=-2*Pi*I/#v,s);vector(#v,k,s=t*(k-1);sum(n=0,#v-1,v[n+1]*exp(s*n)))

Uso

X([1,1,1,1])
%2 = [4.000000000000000000000000000, 0.E-27 + 0.E-28*I, 0.E-28 + 0.E-27*I, 0.E-27 + 0.E-28*I]
P̲̳x͓L̳
fuente
3
¿No es este el ingenuo DFT? (es decir, theta (N ^ 2))
millas del
1

Octava , 109103101100 bytes

f(f=@(f)@(x,n=rows(x)){@(d=f(f)(x(k=2:2:n)).*i.^((k*2-4)/n)')[d+(c=f(f)(x(k-1)));c-d],x}{1+(n<2)}())

Pruébalo en línea!

Ooooo, ¿me sangran los ojos de esta lambda maldita recursiva ? Grandes partes de esto fueron extraídas de la respuesta de @ flawr.

f(                                          % lambda function
  f=@(f)                                    % defined in its own argument list, 
                                            % accepts itself as parameter (for recursion)
    @(x,n=rows(x)){                         % calls another lambda,
                                            % 2nd parameter is just to define a variable
      @(d=f(f)(x(k=2:2:n)).*i.^((k*2-4)/n)')% 1/4 of FFT (argument just defines a variable)
        [d+(c=f(f)(x(k-1)));                % 2/4 of FFT
         c-d                                % 4/4 of FFT
        ],                                  % This is in a @()[] to inhibit evaluation
                                            % unless actually called
      x                                     % FFT of length 1
    }{1+(n<2)}                              % if len(x)==1, return x
                                            % else return [d+c;c-d]
  ()                                        % this is probably important too
)
techo
fuente
No entiendo lo que hiciste pero me gusta mucho.
falla
0

Axioma, 259 , 193 , 181 , 179 bytes

L(g,n,f)==>[g for i in 1..n|f]
h(a)==(n:=#a;n=1=>a;c:=h(L(a.i,n,odd? i));d:=h(L(a.i,n,even? i));n:=n/2;t:=1>0;v:=L(d.i*%i^(-2*(i-1)/n),n,t);append(L(c.i+v.i,n,t),L(c.i-v.i,n,t)))

Incluso si h (a) pudiera pasar toda la prueba y estaría bien como entrada para esta 'competencia', uno debe llamar h () o hlp () a través de fft () a continuación, para verificar los argumentos . No sé si este software puede estar bien porque solo había visto lo que otros escribieron, y busqué la forma en que podría ejecutarse en Axiom para devolver algún posible resultado correcto. Debajo del código sin golf con algunos comentarios:

-- L(g,n,f)==>[g for i in 1..n|f]
-- this macro L, build one List from other list, where in g, there is the generic element of index i
-- (as a.i, or a.i*b.i or a.i*4), n build 1..n that is the range of i, f is the condition 
-- for insert the element in the list result.

hlp(a)==
    n:=#a;n=1=>a
    -- L(a.i,n,odd? i)  it means build a list getting "even indices i of a.i as starting from index 0" [so even is odd and odd is even]
    -- L(a.i,n,even? i) it means build a list getting "odd  indices i of a.i as starting from index 0"
    c:=hlp(L(a.i,n,odd? i));d:=hlp(L(a.i,n,even? i))
    n:=n/2;t:=1>0
    v:=L(d.i*%i^(-2*(i-1)/n),n,t)
    append(L(c.i+v.i,n,t),L(c.i-v.i,n,t))

-- Return Fast Fourier transform of list a, in the case #a=2^n
fft(a)==(n:=#a;n=0 or gcd(n,2^30)~=n=>[];hlp(a))

(5) -> h([1,1,1,1])
   (5)  [4,0,0,0]
                                    Type: List Expression Complex Integer
(6) -> h([1,2,3,4])
   (6)  [10,- 2 + 2%i,- 2,- 2 - 2%i]
                                    Type: List Expression Complex Integer
(7) -> h([5.24626,3.90746,3.72335,5.74429,4.7983,8.34171,4.46785,0.760139])
   (7)
   [36.989359, - 6.2118552150 341603904 + 0.3556612739 187363298 %i,
    1.85336 - 5.744741 %i, 7.1077752150 341603904 - 1.1333387260 812636702 %i,
    - 0.517839, 7.1077752150 341603904 + 1.1333387260 812636702 %i,
    1.85336 + 5.744741 %i,
    - 6.2118552150 341603904 - 0.3556612739 187363298 %i]
                                      Type: List Expression Complex Float
(8) -> h([%i+1,2,%i-2,9])
   (8)  [10 + 2%i,3 + 7%i,- 12 + 2%i,3 - 7%i]
                                    Type: List Expression Complex Integer

en los pocos que había visto h () o fft () devolvería la solución exacta, pero si la simplificación no es buena como en:

(13) -> h([1,2,3,4,5,6,7,8])
   (13)
                    +--+                                   +--+
        (- 4 + 4%i)\|%i  - 4 + 4%i             (- 4 - 4%i)\|%i  - 4 + 4%i
   [36, --------------------------, - 4 + 4%i, --------------------------, - 4,
                    +--+                                   +--+
                   \|%i                                   \|%i
            +--+                                   +--+
    (- 4 + 4%i)\|%i  + 4 - 4%i             (- 4 - 4%i)\|%i  + 4 - 4%i
    --------------------------, - 4 - 4%i, --------------------------]
                +--+                                   +--+
               \|%i                                   \|%i
                                    Type: List Expression Complex Integer

de lo que es suficiente cambiar el tipo de un solo elemento de la lista, como se muestra a continuación en la escritura 8. (Flotante) para encontrar la solución aproximada:

(14) -> h([1,2,3,4,5,6,7,8.])
   (14)
   [36.0, - 4.0000000000 000000001 + 9.6568542494 923801953 %i, - 4.0 + 4.0 %i,
    - 4.0 + 1.6568542494 92380195 %i, - 4.0, - 4.0 - 1.6568542494 92380195 %i,
    - 4.0 - 4.0 %i, - 4.0 - 9.6568542494 923801953 %i]
                                      Type: List Expression Complex Float

Lo escribí, vi todas las otras respuestas porque en el enlace, la página era demasiado difícil, así que no sé si este código puede ser correcto. No soy un experto experto, así que todo esto puede (es probable) estar equivocado.

RosLuP
fuente
0

APL (NARS), 58 caracteres, 116 bytes

{1≥k←≢⍵:⍵⋄(∇⍵[y∼⍨⍳k])(+,-)(∇⍵[y←2×⍳t])×0J1*t÷⍨2-2×⍳t←⌊k÷2}

prueba

  f←{1≥k←≢⍵:⍵⋄(∇⍵[y∼⍨⍳k])(+,-)(∇⍵[y←2×⍳t])×0J1*t÷⍨2-2×⍳t←⌊k÷2}
  f 1 1 1 1
4J0 0J0 0J0 0J0 
  f 1 2 3 4
10J0 ¯2J2 ¯2J0 ¯2J¯2 
  f 1J1 2 ¯2J1  9
10J2 3J7 ¯12J2 3J¯7 
  f 5.24626,3.90746,3.72335,5.74429,4.7983,8.34171,4.46785,0.760139
36.989359J0 ¯6.211855215J0.3556612739 1.85336J¯5.744741 7.107775215J¯1.133338726 ¯0.517839J0 
  7.107775215J1.133338726 1.85336J5.744741 ¯6.211855215J¯0.3556612739 
RosLuP
fuente