\\% (slow) Fourier and Hartley transform and cyclic convolution
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-December-11 (15:36)


dft(a, is=+1)=
\\ Complex Fourier transform (direct computation, n^2 operations)
\\ without normalization.
{
    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 );
} /* ----- */


lconv(a, b)=
\\ Linear convolution (direct computation, n^2 operations)
\\ Example: cconv([a,b,c],[d,e]) ==> [   ]
{
    local(na, nb, n2, f, s, k, k2);
    na = length(a);
    nb = length(b);
    \\ include final element zero (for easy comparison with DFT computation):
    n2 = na + nb;
    f = vector(n2);
    for (tau=0, n2-1,  \\ tau = k + k2
        k2 = min(nb-1, tau);  k = tau - k2;
        s = 0;
        while ( k<na && k2>=0, s += (a[k+1]*b[k2+1]);  k++;  k2-- );
        f[tau+1] = s;
    );
    return( f );
} /* ----- */

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
        k = 0;  k2 = tau;
        s0 = 0;
        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 ====