Code Search (Find, Follow), Analysis (Metrics, Static, Dynamic, Bugs) and Change (Modernization, Migration, Generation, Optimization, Rearchitecting)

Example of C++ Smart Differencer Output

This page contains an example of the output generated by SD's C++ Smart Differencer tool
when applied to an original file and an updated version of the same file.

Output from C++ Smart Differencer on original and updated file

This page contains an example of the output generated by SD's C++ Smart Differencer tool.
Observe that the Smart Differencer produces only 20 lines of output, compared to diff's 117 lines of output; the programmer simply has to look at less code. Note the added precision of Smart Differencer output: precise deltas with line and column number (l.c) ranges, rather than simple line deltas. It notes renamed variables ('oldname' ->'newname') across a block (lines 68-104) just once rather than report many small blocks of code in which the variable name has simply changed as diff does. (With enough information, the Smart Difference will note that a name changed has occurred consistently throughout out a scope, and therefore isn't an interesting difference). Notice that source code formatting (different line breaks in line 32-54) don't fool the smart differencer, that changed comments don't confuse it, that equivalent constants (line 437) are ignored. Finally, a block of code is noted as deleted/added (moved) at line 282.
For comparison, see deltas generated from a standard diff tool.

Rename 68.42-104.10 to 62.42-98.10 with 'xtemp'->'x2'
Rename 161.1-194.1 to 155.1-188.1 with 'Vx'~>'U'
Delete 282.1-289.1 moving 282.1-289.1 to 429.1-436.1
<void DST(const ColumnVector& U, ColumnVector& V)
<{
<   // Discrete sine transform, type I
<   Tracer trace("DST");
<   REPORT
<   DST_inverse(U, V);
<   V *= (V.Nrows()-1)/2;
<}
Insert 429.1-436.1 moving 282.1-289.1 to 429.1-436.1
>void DST(const ColumnVector& U, ColumnVector& V)
>{
>   // Discrete sine transform, type I
>   Tracer trace("DST");
>   REPORT
>   DST_inverse(U, V);
>   V *= (V.Nrows()-1)/2;
>}

Output from Cygwin diff tool

Note absence of column number data, absence of identifier rename detection, failure to recognize reformatted code as unchanged (diff thinks the block of code at lines 32-54 has changed), failure to notice same numeric values is used, although it is written differently (line 437, and failure to ignore comments. You are told about the moved block of code (line 282) twice (line 443), but you aren't told it was moved or that these are related in any way.
See the actual deltas generated from C++ Smart Differencer.
See original file and updated version of the same file.

21c21
< // minimise roundoff error; SHOULD TEST THIS CODE
---
> // minimise roundoff error
32,38c32,33
<    case 0:
<        REPORT c =  cos(ratio);
<        s =  sin(ratio);
<        break;
<    case 1:
<        REPORT c = -sin(ratio); s =  cos(ratio);
<        break;
---
>    case 0: REPORT c =  cos(ratio); s =  sin(ratio); break;
>    case 1: REPORT c = -sin(ratio); s =  cos(ratio); break;
40,41c35
<    case 3: REPORT c =  sin(ratio); s =
<           -cos(ratio); break;
---
>    case 3: REPORT c =  sin(ratio); s = -cos(ratio); break;
54c48
<    Real r_arg = 1.;  Real i_arg = 0.;
---
>    Real r_arg = 1.0;  Real i_arg = 0.0;
68c62
<          Real* a1 = a++; Real* b1 = b++; Real* xtemp = x1++; Real* y2 = y1++;
---
>          Real* a1 = a++; Real* b1 = b++; Real* x2 = x1++; Real* y2 = y1++;
77c71
<                *xtemp = r_value * r_arg - i_value * i_arg + *(a2-gamma);
---
>                *x2 = r_value * r_arg - i_value * i_arg + *(a2-gamma);
80c74
<                xtemp += delta; y2 += delta;
---
>                x2 += delta; y2 += delta;
100c94
<                *xtemp = r_value; *y2 = i_value;
---
>                *x2 = r_value; *y2 = i_value;
102c96
<                xtemp += delta; y2 += delta;
---
>                x2 += delta; y2 += delta;
146,147c140,141
<    *x++ = *a + *b; *y++ = 0.;                // first complex element
<    *xn-- = *a++ - *b++; *yn-- = 0.;          // last complex element
---
>    *x++ = *a + *b; *y++ = 0.0;                // first complex element
>    *xn-- = *a++ - *b++; *yn-- = 0.0;          // last complex element
149c143
<    int j = - 1; i = n2/2;
---
>    int j = -1; i = n2/2;
161c155
< void RealFFTI(const ColumnVector& A, const ColumnVector& B, ColumnVector& Vx)
---
> void RealFFTI(const ColumnVector& A, const ColumnVector& B, ColumnVector& U)
191,192c185,186
<    Vx.ReSize(n); i = n2;
<    x = X.Store(); y = Y.Store(); Real* u = Vx.Store();
---
>    U.ReSize(n); i = n2;
>    x = X.Store(); y = Y.Store(); Real* u = U.Store();
282,290d275
< void DST(const ColumnVector& U, ColumnVector& V)
< {
<    // Discrete sine transform, type I
<    Tracer trace("DST");
<    REPORT
<    DST_inverse(U, V);
<    V *= (V.Nrows()-1)/2;
< }
< 
303c288
<    *x = *v; *y = 0.;
---
>    *x = *v; *y = 0.0;
358c343
<    *x = *(--w); *y = 0.;
---
>    *x = *(--w); *y = 0.0;
384c369
<    Real vi = *v++; *x++ = vi; *y++ = 0.;
---
>    Real vi = *v++; *x++ = vi; *y++ = 0.0;
393c378
<    vi = *v; *x = vi; *y = 0.; vi /= 2.0; sum1 += vi; sum2 += vi;
---
>    vi = *v; *x = vi; *y = 0.0; vi /= 2.0; sum1 += vi; sum2 += vi;
427c412
<    Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.;
---
>    Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.0;
430c415
<    *x = -2 * vi; *y = 0.;
---
>    *x = -2 * vi; *y = 0.0;
434c419
<    i = n2; int k = 0; *u++ = 0.; *v-- = 0.;
---
>    i = n2; int k = 0; *u++ = 0.0; *v-- = 0.0;
437c422
<       Real s = sin(15.707963267948966192e-1 * (++k) / n2);
---
>       Real s = sin(1.5707963267948966192 * (++k) / n2);
443a429,437
> void DST(const ColumnVector& U, ColumnVector& V)
> {
>    // Discrete sine transform, type I
>    Tracer trace("DST");
>    REPORT
>    DST_inverse(U, V);
>    V *= (V.Nrows()-1)/2;
> }
> 

Source Code before Change

This code is an open source Fast Fourier transform of 443 lines. A number of cleanups were made from before and after.
See updated version of the same file.

//$$ fft.cpp                         Fast fourier transform
// Copyright (C) 1991,2,3,4,8: R B Davies

#define WANT_MATH
// #define WANT_STREAM

#include "include.h"

#include "newmatap.h"

// #include "newmatio.h"

#ifdef DO_REPORT
#define REPORT { static ExeCounter ExeCount(__LINE__,19); ++ExeCount; }
#else
#define REPORT {}
#endif

static void cossin(int n, int d, Real& c, Real& s)
// calculate cos(twopi*n/d) and sin(twopi*n/d)
// minimise roundoff error; SHOULD TEST THIS CODE
{
   REPORT
   long n4 = n * 4; int sector = (int)floor( (Real)n4 / (Real)d + 0.5 );
   n4 -= sector * d;
   if (sector < 0) { REPORT sector = 3 - (3 - sector) % 4; }
   else  { REPORT sector %= 4; }
   Real ratio = 1.5707963267948966192 * (Real)n4 / (Real)d;

   switch (sector)
   {
   case 0:
       REPORT c =  cos(ratio);
       s =  sin(ratio);
       break;
   case 1:
       REPORT c = -sin(ratio); s =  cos(ratio);
       break;
   case 2: REPORT c = -cos(ratio); s = -sin(ratio); break;
   case 3: REPORT c =  sin(ratio); s =
          -cos(ratio); break;
   }
}

static void fftstep(ColumnVector& A, ColumnVector& B, ColumnVector& X,
   ColumnVector& Y, int after, int now, int before)
{
   REPORT
   Tracer trace("FFT(step)");
   // const Real twopi = 6.2831853071795864769;
   const int gamma = after * before;  const int delta = now * after;
   // const Real angle = twopi / delta;  Real temp;
   // Real r_omega = cos(angle);  Real i_omega = -sin(angle);
   Real r_arg = 1.;  Real i_arg = 0.;
   Real* x = X.Store();  Real* y = Y.Store();   // pointers to array storage
   const int m = A.Nrows() - gamma;

   for (int j = 0; j < now; j++)
   {
      Real* a = A.Store(); Real* b = B.Store(); // pointers to array storage
      Real* x1 = x; Real* y1 = y; x += after; y += after;
      for (int ia = 0; ia < after; ia++)
      {
         // generate sins & cosines explicitly rather than iteratively
         // for more accuracy; but slower
         cossin(-(j*after+ia), delta, r_arg, i_arg);

         Real* a1 = a++; Real* b1 = b++; Real* xtemp = x1++; Real* y2 = y1++;
         if (now==2)
         {
            REPORT int ib = before;
            if (ib) for (;;)
            {
               REPORT
               Real* a2 = m + a1; Real* b2 = m + b1; a1 += after; b1 += after;
               Real r_value = *a2; Real i_value = *b2;
               *xtemp = r_value * r_arg - i_value * i_arg + *(a2-gamma);
               *y2 = r_value * i_arg + i_value * r_arg + *(b2-gamma);
               if (!(--ib)) break;
               xtemp += delta; y2 += delta;
            }
         }
         else
         {
            REPORT int ib = before;
            if (ib) for (;;)
            {
               REPORT
               Real* a2 = m + a1; Real* b2 = m + b1; a1 += after; b1 += after;
               Real r_value = *a2; Real i_value = *b2;
               int in = now-1; while (in--)
               {
                  // it should be possible to make this faster
                  // hand code for now = 2,3,4,5,8
                  // use symmetry to halve number of operations
                  a2 -= gamma; b2 -= gamma;  Real temp = r_value;
                  r_value = r_value * r_arg - i_value * i_arg + *a2;
                  i_value = temp    * i_arg + i_value * r_arg + *b2;
               }
               *xtemp = r_value; *y2 = i_value;
               if (!(--ib)) break;
               xtemp += delta; y2 += delta;
            }
         }

         // temp = r_arg;
         // r_arg = r_arg * r_omega - i_arg * i_omega;
         // i_arg = temp  * i_omega + i_arg * r_omega;

      }
   }
}


void FFTI(const ColumnVector& U, const ColumnVector& V,
   ColumnVector& X, ColumnVector& Y)
{
   // Inverse transform
   Tracer trace("FFTI");
   REPORT
   FFT(U,-V,X,Y);
   const Real n = X.Nrows(); X /= n; Y /= (-n);
}

void RealFFT(const ColumnVector& U, ColumnVector& X, ColumnVector& Y)
{
   // Fourier transform of a real series
   Tracer trace("RealFFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n2), B(n2);
   Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;
   while (i--) { *a++ = *u++; *b++ = *u++; }
   FFT(A,B,A,B);
   int n21 = n2 + 1;
   X.ReSize(n21); Y.ReSize(n21);
   i = n2 - 1;
   a = A.Store(); b = B.Store();              // first els of A and B
   Real* an = a + i; Real* bn = b + i;        // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + n2; Real* yn = y + n2;      // last els of X and Y

   *x++ = *a + *b; *y++ = 0.;                // first complex element
   *xn-- = *a++ - *b++; *yn-- = 0.;          // last complex element

   int j = - 1; i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am + c * bp; Real sbpcam = s * bp - c * am;
      *x++  =  0.5 * ( ap + samcbp); *y++  =  0.5 * ( bm + sbpcam);
      *xn-- =  0.5 * ( ap - samcbp); *yn-- =  0.5 * (-bm + sbpcam);
   }
}

void RealFFTI(const ColumnVector& A, const ColumnVector& B, ColumnVector& Vx)
{
   // inverse of a Fourier transform of a real series
   Tracer trace("RealFFTI");
   REPORT
   const int n21 = A.Nrows();                     // length of arrays
   if (n21 != B.Nrows() || n21 == 0)
      Throw(ProgramException("Vector lengths unequal or zero", A, B));
   const int n2 = n21 - 1;  const int n = 2 * n2;  int i = n2 - 1;

   ColumnVector X(n2), Y(n2);
   Real* a = A.Store(); Real* b = B.Store();  // first els of A and B
   Real* an = a + n2;   Real* bn = b + n2;    // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + i;    Real* yn = y + i;     // last els of X and Y

   Real hn = 0.5 / n2;
   *x++  = hn * (*a + *an);  *y++  = - hn * (*a - *an);
   a++; an--; b++; bn--;
   int j = -1;  i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am - c * bp; Real sbpcam = s * bp + c * am;
      *x++  =  hn * ( ap + samcbp); *y++  =  - hn * ( bm + sbpcam);
      *xn-- =  hn * ( ap - samcbp); *yn-- =  - hn * (-bm + sbpcam);
   }
   FFT(X,Y,X,Y);             // have done inverting elsewhere
   Vx.ReSize(n); i = n2;
   x = X.Store(); y = Y.Store(); Real* u = Vx.Store();
   while (i--) { *u++ = *x++; *u++ = - *y++; }
}

void FFT(const ColumnVector& U, const ColumnVector& V,
   ColumnVector& X, ColumnVector& Y)
{
   // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
   // but first try Sande and Gentleman
   Tracer trace("FFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   if (n != V.Nrows() || n == 0)
      Throw(ProgramException("Vector lengths unequal or zero", U, V));
   if (n == 1) { REPORT X = U; Y = V; return; }

   // see if we can use the newfft routine
   if (!FFT_Controller::OnlyOldFFT && FFT_Controller::CanFactor(n))
   {
      REPORT
      X = U; Y = V;
      if ( FFT_Controller::ar_1d_ft(n,X.Store(),Y.Store()) ) return;
   }

   ColumnVector B = V;
   ColumnVector A = U;
   X.ReSize(n); Y.ReSize(n);
   const int nextmx = 8;
#ifndef ATandT
   int prime[8] = { 2,3,5,7,11,13,17,19 };
#else
   int prime[8];
   prime[0]=2; prime[1]=3; prime[2]=5; prime[3]=7;
   prime[4]=11; prime[5]=13; prime[6]=17; prime[7]=19;
#endif
   int after = 1; int before = n; int next = 0; bool inzee = true;
   int now = 0; int b1;             // initialised to keep gnu happy

   do
   {
      for (;;)
      {
	 if (next < nextmx) { REPORT now = prime[next]; }
	 b1 = before / now;  if (b1 * now == before) { REPORT break; }
	 next++; now += 2;
      }
      before = b1;

      if (inzee) { REPORT fftstep(A, B, X, Y, after, now, before); }
      else { REPORT fftstep(X, Y, A, B, after, now, before); }

      inzee = !inzee; after *= now;
   }
   while (before != 1);

   if (inzee) { REPORT A.Release(); X = A; B.Release(); Y = B; }
}

// Trigonometric transforms
// see Charles Van Loan (1992) "Computational frameworks for the fast
// Fourier transform" published by SIAM; section 4.4.

void DCT_II(const ColumnVector& U, ColumnVector& V)
{
   // Discrete cosine transform, type II, of a real series
   Tracer trace("DCT_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = *u++; }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.CleanUp();
   V.ReSize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *v = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *(++v) = xi * c + yi * s; *(--w) = xi * s - yi * c;
   }
}

void DST(const ColumnVector& U, ColumnVector& V)
{
   // Discrete sine transform, type I
   Tracer trace("DST");
   REPORT
   DST_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

void DCT_II_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete cosine transform, type II
   Tracer trace("DCT_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *v; *y = 0.;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *(++v); Real wi = *(--w);
      *(++x) = vi * c + wi * s; *(++y) = vi * s - wi * c;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = *(--b); }
}

void DST_II(const ColumnVector& U, ColumnVector& V)
{
   // Discrete sine transform, type II, of a real series
   Tracer trace("DST_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = -(*u++); }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.CleanUp();
   V.ReSize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *(--w) = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *v++ = xi * s - yi * c; *(--w) = xi * c + yi * s;
   }
}

void DST_II_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete sine transform, type II
   Tracer trace("DST_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *(--w); *y = 0.;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *v++; Real wi = *(--w);
      *(++x) = vi * s + wi * c; *(++y) = - vi * c + wi * s;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = -(*(--b)); }
}

void DCT_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete cosine transform, type I
   Tracer trace("DCT_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *v++; *x++ = vi; *y++ = 0.;
   Real sum1 = vi / 2.0; Real sum2 = sum1; vi = *v++;
   int i = n2-1;
   while (i--)
   {
      Real vi2 = *v++; sum1 += vi2 + vi; sum2 += vi2 - vi;
      *x++ = vi2; vi2 = *v++; *y++ = vi - vi2; vi = vi2;
   }
   sum1 += vi; sum2 -= vi;
   vi = *v; *x = vi; *y = 0.; vi /= 2.0; sum1 += vi; sum2 += vi;
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = sum1 / n2; *v-- = sum2 / n2;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real bz = (ai - bi) / 4 / s; Real az = (ai + bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

void DCT(const ColumnVector& U, ColumnVector& V)
{
   // Discrete cosine transform, type I
   Tracer trace("DCT");
   REPORT
   DCT_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

void DST_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete sine transform, type I
   Tracer trace("DST_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.;
   int i = n2-1;
   while (i--) { *y++ = *(++v); Real vi2 = *(++v); *x++ = vi2 - vi; vi = vi2; }
   *x = -2 * vi; *y = 0.;
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = 0.; *v-- = 0.;
   while (i--)
   {
      Real s = sin(15.707963267948966192e-1 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real az = (ai + bi) / 4 / s; Real bz = (ai - bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

Source Code after Change

This code is that same file at the next checkin with a number of changes; it is now 437 lines. It is hard to see the differences by simply looking at the text, partly just from the sheer size of it. This is why you want a diff tool of some kind.
See the actual deltas generated from C++ Smart Differencer.
For comparison, see deltas generated from a standard diff tool.

//$$ fft.cpp                         Fast fourier transform
// Copyright (C) 1991,2,3,4,8: R B Davies

#define WANT_MATH
// #define WANT_STREAM

#include "include.h"

#include "newmatap.h"

// #include "newmatio.h"

#ifdef DO_REPORT
#define REPORT { static ExeCounter ExeCount(__LINE__,19); ++ExeCount; }
#else
#define REPORT {}
#endif

static void cossin(int n, int d, Real& c, Real& s)
// calculate cos(twopi*n/d) and sin(twopi*n/d)
// minimise roundoff error
{
   REPORT
   long n4 = n * 4; int sector = (int)floor( (Real)n4 / (Real)d + 0.5 );
   n4 -= sector * d;
   if (sector < 0) { REPORT sector = 3 - (3 - sector) % 4; }
   else  { REPORT sector %= 4; }
   Real ratio = 1.5707963267948966192 * (Real)n4 / (Real)d;

   switch (sector)
   {
   case 0: REPORT c =  cos(ratio); s =  sin(ratio); break;
   case 1: REPORT c = -sin(ratio); s =  cos(ratio); break;
   case 2: REPORT c = -cos(ratio); s = -sin(ratio); break;
   case 3: REPORT c =  sin(ratio); s = -cos(ratio); break;
   }
}

static void fftstep(ColumnVector& A, ColumnVector& B, ColumnVector& X,
   ColumnVector& Y, int after, int now, int before)
{
   REPORT
   Tracer trace("FFT(step)");
   // const Real twopi = 6.2831853071795864769;
   const int gamma = after * before;  const int delta = now * after;
   // const Real angle = twopi / delta;  Real temp;
   // Real r_omega = cos(angle);  Real i_omega = -sin(angle);
   Real r_arg = 1.0;  Real i_arg = 0.0;
   Real* x = X.Store();  Real* y = Y.Store();   // pointers to array storage
   const int m = A.Nrows() - gamma;

   for (int j = 0; j < now; j++)
   {
      Real* a = A.Store(); Real* b = B.Store(); // pointers to array storage
      Real* x1 = x; Real* y1 = y; x += after; y += after;
      for (int ia = 0; ia < after; ia++)
      {
         // generate sins & cosines explicitly rather than iteratively
         // for more accuracy; but slower
         cossin(-(j*after+ia), delta, r_arg, i_arg);

         Real* a1 = a++; Real* b1 = b++; Real* x2 = x1++; Real* y2 = y1++;
         if (now==2)
         {
            REPORT int ib = before;
            if (ib) for (;;)
            {
               REPORT
               Real* a2 = m + a1; Real* b2 = m + b1; a1 += after; b1 += after;
               Real r_value = *a2; Real i_value = *b2;
               *x2 = r_value * r_arg - i_value * i_arg + *(a2-gamma);
               *y2 = r_value * i_arg + i_value * r_arg + *(b2-gamma);
               if (!(--ib)) break;
               x2 += delta; y2 += delta;
            }
         }
         else
         {
            REPORT int ib = before;
            if (ib) for (;;)
            {
               REPORT
               Real* a2 = m + a1; Real* b2 = m + b1; a1 += after; b1 += after;
               Real r_value = *a2; Real i_value = *b2;
               int in = now-1; while (in--)
               {
                  // it should be possible to make this faster
                  // hand code for now = 2,3,4,5,8
                  // use symmetry to halve number of operations
                  a2 -= gamma; b2 -= gamma;  Real temp = r_value;
                  r_value = r_value * r_arg - i_value * i_arg + *a2;
                  i_value = temp    * i_arg + i_value * r_arg + *b2;
               }
               *x2 = r_value; *y2 = i_value;
               if (!(--ib)) break;
               x2 += delta; y2 += delta;
            }
         }

         // temp = r_arg;
         // r_arg = r_arg * r_omega - i_arg * i_omega;
         // i_arg = temp  * i_omega + i_arg * r_omega;

      }
   }
}


void FFTI(const ColumnVector& U, const ColumnVector& V,
   ColumnVector& X, ColumnVector& Y)
{
   // Inverse transform
   Tracer trace("FFTI");
   REPORT
   FFT(U,-V,X,Y);
   const Real n = X.Nrows(); X /= n; Y /= (-n);
}

void RealFFT(const ColumnVector& U, ColumnVector& X, ColumnVector& Y)
{
   // Fourier transform of a real series
   Tracer trace("RealFFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n2), B(n2);
   Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;
   while (i--) { *a++ = *u++; *b++ = *u++; }
   FFT(A,B,A,B);
   int n21 = n2 + 1;
   X.ReSize(n21); Y.ReSize(n21);
   i = n2 - 1;
   a = A.Store(); b = B.Store();              // first els of A and B
   Real* an = a + i; Real* bn = b + i;        // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + n2; Real* yn = y + n2;      // last els of X and Y

   *x++ = *a + *b; *y++ = 0.0;                // first complex element
   *xn-- = *a++ - *b++; *yn-- = 0.0;          // last complex element

   int j = -1; i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am + c * bp; Real sbpcam = s * bp - c * am;
      *x++  =  0.5 * ( ap + samcbp); *y++  =  0.5 * ( bm + sbpcam);
      *xn-- =  0.5 * ( ap - samcbp); *yn-- =  0.5 * (-bm + sbpcam);
   }
}

void RealFFTI(const ColumnVector& A, const ColumnVector& B, ColumnVector& U)
{
   // inverse of a Fourier transform of a real series
   Tracer trace("RealFFTI");
   REPORT
   const int n21 = A.Nrows();                     // length of arrays
   if (n21 != B.Nrows() || n21 == 0)
      Throw(ProgramException("Vector lengths unequal or zero", A, B));
   const int n2 = n21 - 1;  const int n = 2 * n2;  int i = n2 - 1;

   ColumnVector X(n2), Y(n2);
   Real* a = A.Store(); Real* b = B.Store();  // first els of A and B
   Real* an = a + n2;   Real* bn = b + n2;    // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + i;    Real* yn = y + i;     // last els of X and Y

   Real hn = 0.5 / n2;
   *x++  = hn * (*a + *an);  *y++  = - hn * (*a - *an);
   a++; an--; b++; bn--;
   int j = -1;  i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am - c * bp; Real sbpcam = s * bp + c * am;
      *x++  =  hn * ( ap + samcbp); *y++  =  - hn * ( bm + sbpcam);
      *xn-- =  hn * ( ap - samcbp); *yn-- =  - hn * (-bm + sbpcam);
   }
   FFT(X,Y,X,Y);             // have done inverting elsewhere
   U.ReSize(n); i = n2;
   x = X.Store(); y = Y.Store(); Real* u = U.Store();
   while (i--) { *u++ = *x++; *u++ = - *y++; }
}

void FFT(const ColumnVector& U, const ColumnVector& V,
   ColumnVector& X, ColumnVector& Y)
{
   // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
   // but first try Sande and Gentleman
   Tracer trace("FFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   if (n != V.Nrows() || n == 0)
      Throw(ProgramException("Vector lengths unequal or zero", U, V));
   if (n == 1) { REPORT X = U; Y = V; return; }

   // see if we can use the newfft routine
   if (!FFT_Controller::OnlyOldFFT && FFT_Controller::CanFactor(n))
   {
      REPORT
      X = U; Y = V;
      if ( FFT_Controller::ar_1d_ft(n,X.Store(),Y.Store()) ) return;
   }

   ColumnVector B = V;
   ColumnVector A = U;
   X.ReSize(n); Y.ReSize(n);
   const int nextmx = 8;
#ifndef ATandT
   int prime[8] = { 2,3,5,7,11,13,17,19 };
#else
   int prime[8];
   prime[0]=2; prime[1]=3; prime[2]=5; prime[3]=7;
   prime[4]=11; prime[5]=13; prime[6]=17; prime[7]=19;
#endif
   int after = 1; int before = n; int next = 0; bool inzee = true;
   int now = 0; int b1;             // initialised to keep gnu happy

   do
   {
      for (;;)
      {
	 if (next < nextmx) { REPORT now = prime[next]; }
	 b1 = before / now;  if (b1 * now == before) { REPORT break; }
	 next++; now += 2;
      }
      before = b1;

      if (inzee) { REPORT fftstep(A, B, X, Y, after, now, before); }
      else { REPORT fftstep(X, Y, A, B, after, now, before); }

      inzee = !inzee; after *= now;
   }
   while (before != 1);

   if (inzee) { REPORT A.Release(); X = A; B.Release(); Y = B; }
}

// Trigonometric transforms
// see Charles Van Loan (1992) "Computational frameworks for the fast
// Fourier transform" published by SIAM; section 4.4.

void DCT_II(const ColumnVector& U, ColumnVector& V)
{
   // Discrete cosine transform, type II, of a real series
   Tracer trace("DCT_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = *u++; }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.CleanUp();
   V.ReSize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *v = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *(++v) = xi * c + yi * s; *(--w) = xi * s - yi * c;
   }
}

void DCT_II_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete cosine transform, type II
   Tracer trace("DCT_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *v; *y = 0.0;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *(++v); Real wi = *(--w);
      *(++x) = vi * c + wi * s; *(++y) = vi * s - wi * c;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = *(--b); }
}

void DST_II(const ColumnVector& U, ColumnVector& V)
{
   // Discrete sine transform, type II, of a real series
   Tracer trace("DST_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = -(*u++); }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.CleanUp();
   V.ReSize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *(--w) = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *v++ = xi * s - yi * c; *(--w) = xi * c + yi * s;
   }
}

void DST_II_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete sine transform, type II
   Tracer trace("DST_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *(--w); *y = 0.0;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *v++; Real wi = *(--w);
      *(++x) = vi * s + wi * c; *(++y) = - vi * c + wi * s;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = -(*(--b)); }
}

void DCT_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete cosine transform, type I
   Tracer trace("DCT_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *v++; *x++ = vi; *y++ = 0.0;
   Real sum1 = vi / 2.0; Real sum2 = sum1; vi = *v++;
   int i = n2-1;
   while (i--)
   {
      Real vi2 = *v++; sum1 += vi2 + vi; sum2 += vi2 - vi;
      *x++ = vi2; vi2 = *v++; *y++ = vi - vi2; vi = vi2;
   }
   sum1 += vi; sum2 -= vi;
   vi = *v; *x = vi; *y = 0.0; vi /= 2.0; sum1 += vi; sum2 += vi;
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = sum1 / n2; *v-- = sum2 / n2;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real bz = (ai - bi) / 4 / s; Real az = (ai + bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

void DCT(const ColumnVector& U, ColumnVector& V)
{
   // Discrete cosine transform, type I
   Tracer trace("DCT");
   REPORT
   DCT_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

void DST_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete sine transform, type I
   Tracer trace("DST_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.0;
   int i = n2-1;
   while (i--) { *y++ = *(++v); Real vi2 = *(++v); *x++ = vi2 - vi; vi = vi2; }
   *x = -2 * vi; *y = 0.0;
   ColumnVector A; RealFFTI(X, Y, A);
   X.CleanUp(); Y.CleanUp(); U.ReSize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = 0.0; *v-- = 0.0;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real az = (ai + bi) / 4 / s; Real bz = (ai - bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

void DST(const ColumnVector& U, ColumnVector& V)
{
   // Discrete sine transform, type I
   Tracer trace("DST");
   REPORT
   DST_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}