\\% (slow) Fourier and Hartley transform and cyclic convolution
\\ Author: Joerg Arndt
\\ online at http://www.jjj.de/pari/
\\ version: 2010-May-18 (13:01)


dft(a, is=+1)=
/* Complex Fourier transform (direct computation, n^2 operations) */
{
    local(n, f, s, ph0, ph);
    n = length(a);
    f = vector(n);
    ph0 = is*2*Pi*I/n;
    for (k=0, n-1,
        ph = ph0 * k;
        f[k+1] = sum (x=0, n-1,  a[x+1] * exp(ph*x) );
    );
    return( f );
} /* ----- */

dht(v)=
/* Hartley transform (direct computation, n^2 operations) */
{
    local(n, phi, h);
    n = length(v);
    h = vector(n);
    for (k=0, n-1,
        phi = 2*Pi*k/n;
        h[k+1] = sum(j=0,n-1, v[j+1]*(cos(j*phi)+sin(j*phi)) );
    );
    return( h );
} /* ----- */


cconv(a, b)=
/* Cyclic convolution (direct computation, n^2 operations) */
/* Example: cconv([a,b],[c,d]) ==> [b*d + c*a, a*d + c*b] */
{
    local(n, f, s, k, k2);
    n = length(a);
    f = vector(n);
    for (tau=0, n-1,  \\ tau = k + k2
        s0 = 0;  k = 0;  k2 = tau;
        while (k<=tau, s0 += (a[k+1]*b[k2+1]);  k++; k2--);
        s1 = 0;  k2 = n-1; \\ k=tau+1
        while (k<n,    s1 += (a[k+1]*b[k2+1]);  k++; k2--);
        f[tau+1] = s0 + s1;
    );
    return( f );
} /* ----- */

wcconv(a, b, w)=
/* Weighted cyclic convolution (direct computation, n^2 operations) */
/* w==+1 ==> usual cyclic convolution */
/* Example: wcconv([a,b],[c,d]) ==> [b*d + c*a, a*d + c*b] */
/* Example: wcconv([a,b],[c,d],-1) ==> [-b*d + c*a, a*d + c*b] */
{
    local(n, f, s, k, k2);
    n = length(a);
    f = vector(n);
    for (tau=0, n-1,  \\ tau = k + k2
        s0 = 0;  k = 0;  k2 = tau;
        while (k<=tau, s0 += (a[k+1]*b[k2+1]);  k++; k2--);
        s1 = 0;  k2 = n-1; \\ k=tau+1
        while (k<n,    s1 += (a[k+1]*b[k2+1]);  k++; k2--);
        f[tau+1] = s0 + w * s1;
    );
    return( f );
} /* ----- */


\\ ==== end of file ====