// MASH Simulation
// http://www.holmea.demon.co.uk/Frac2/Mash.htm

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <minmax.h>
#include <float.h>

#undef NDEBUG
#include <assert.h>


const double PI = 3.14159265358979323846;

#define MODULO (1<<24)
#define HANNING

#ifdef HANNING
#define SEQLEN (1<<20)      // Hanning Window
#define SEQTYP double
#else
#define SEQLEN (4*MODULO)   // 2* when F even; 4* if odd
#define SEQTYP char
#endif

// 2L = Sequence length (in ref cycles)
const double L = SEQLEN/2;

// Parameters
const int F=1;


class MASH
{
    int q[4];       // Latch outputs
    int c[4][4];    // Carry flip flops

public:
    MASH() // Constructor
    { 
        // Initial values
        memset(c, 0, sizeof c);
        q[0] = 1;
        q[1] = 0;
        q[2] = 0;
        q[3] = 0;

        // Fill flip-flops
        for (int i=0; i<4; i++) Clock();
    }

    int Clock()
    {
        int d[4]; // Adder output / latch data input

        // Carry flip-flops
        c[3][3] = c[3][2]; c[3][2] = c[3][1]; c[3][1] = c[3][0];
                           c[2][2] = c[2][1]; c[2][1] = c[2][0];
                                              c[1][1] = c[1][0];

        // Adders
        d[0] = (q[0]+   F) % MODULO;
        d[1] = (q[1]+d[0]) % MODULO;
        d[2] = (q[2]+d[1]) % MODULO;
        d[3] = (q[3]+d[2]) % MODULO;

        // Carries
        c[0][0] = (q[0]+   F) / MODULO;
        c[1][0] = (q[1]+d[0]) / MODULO;
        c[2][0] = (q[2]+d[1]) / MODULO;
        c[3][0] = (q[3]+d[2]) / MODULO;

        // Latch
        for (int i=0; i<4; i++) q[i] = d[i];

        return  c[0][0]
            +   c[1][0] -   c[1][1]
            +   c[2][0] - 2*c[2][1] +   c[2][2]
            +   c[3][0] - 3*c[3][1] + 3*c[3][2] - c[3][3]
            ;
    }

    bool operator == (MASH& m) // Debug
    {
        return !memcmp(this, &m, sizeof(MASH));
    }
};


void main(int argc, char *argv[])
{
    _control87(_PC_64, MCW_PC);

    double n, r, start, stop, span;
    int i, steps;

    if (argc!=4) return;

    start  = atof(argv[1]);
    stop   = atof(argv[2]);
    steps  = int(atof(argv[3]));
    span = log10(stop/start);
    r = pow(10, span/steps);

    for (i=0, n=start; i<steps && n<stop; n=max(floor(n*r), n+1), i++);
    steps = i;

    double *a = new double[steps];
    double *b = new double[steps];

    SEQTYP *s = new SEQTYP[SEQLEN];
    MASH m;
//  for (i=0; i<SEQLEN; i++) s[i] = i & 1024? -1:1;
    for (i=0; i<SEQLEN; i++) s[i] = m.Clock();
#ifdef HANNING
    for (i=0; i<SEQLEN; i++) s[i] *= (1 - cos(2*PI/double(SEQLEN) * double(i)))/2; // Hanning Window
#else
    assert(m==MASH());
#endif

    time_t tt = time(0);

    for (i=0, n=start; i<steps; n=max(floor(n*r), n+1), i++)
    {
        double sin_t2=0, sin_t1, re=0, recur_a = 2*pow(sin(n*PI/L/2), 2);
        double cos_t2=1, cos_t1, im=0, recur_b =       sin(n*PI/L  );

        SEQTYP *f = s;

        for (int t=SEQLEN; t--; f++)
        {
            cos_t1 = cos_t2;
            sin_t1 = sin_t2;

            re += (*f) * cos_t1;
            im -= (*f) * sin_t1;

            cos_t2 = cos_t1 - (recur_a*cos_t1 + recur_b*sin_t1);
            sin_t2 = sin_t1 - (recur_a*sin_t1 - recur_b*cos_t1);
        }

        double cos_x = cos(n*PI/L);
        double sin_x = sin(n*PI/L);

        a[i] = re * (1.0 - cos_x) - im * sin_x;
        b[i] = im * (1.0 - cos_x) + re * sin_x;

        fprintf(stderr, "\r%0.1f%%  ", 100.0*(i+1)/steps);
        __int64 rem = __int64(time(0) - tt) * __int64(steps-i-1) / (i+1);
        fprintf(stderr, "%02I64d:%02I64d:%02I64d ", rem/3600, (rem/60)%60, rem%60);
    }

    fprintf(stderr, "\nTime %d\n", time(0) - tt);

    for (i=0, n=start; i<steps; n=max(floor(n*r), n+1), i++)
    {
        double re = a[i]/PI/n;
        double im = b[i]/PI/n;
        printf("%g,%g\n", log10(double(n)), 10*log10(re*re + im*im));
//      printf("%g,%g,%g\n", log10(double(n)), re, im);
    }
}