#include <iostream>
#include <math.h>

using namespace std;



typedef struct
{
    unsigned int chain;
    int          slots; // to calculate frequency
    double       slot_frequency;
    int          length; // the effective chain length which will be <
			 // slots , for < d
    double       length_frequency;
    int          length2; // the effective chain length which will be <
			  // slots, for == d
    double       length_frequency2;
} Rec;

int
slots (unsigned int s)
{
    int n = 0;

    for (int i = 0; i < 32; i++)
    {
        if (s & (0x1 << i)) n++;
    }

    return n;
}

int
length (unsigned int s, int L, int *p)
{
    int n = 0;
    
    for (int i = L - 1 - *p; i >= 0; i--, (*p)++)
    {
	if (s & (1 << i))
	    n++;
	else if (n)
	    break;
	else
	    continue;
    }

    return n;
}

void
factorial (double f [] /* d + 1 elements */, int d)
{
    for (int i = 0; i <= d; i++)
    {
        if (i == 0)
        {
            f [i] = 1;
        }
        else
        {
            f [i] = f [i-1] * i;
        }
    }
}

char*
fmtbits (unsigned int chain, int L, char* buf)
{
    for (int i = L - 1; i >= 0; i--)
    {
        buf [L - i - 1] = (chain & (1 << i)) ? '1' : '0';
    }

    buf [L] = '\0';

    return buf;
}

Rec r [1024*1024];

int
main (int argc, char* argv [])
{
    FILE* fp1;
    FILE* fp2;
    FILE* fp3;
    int m;
    int d;
    int L;
    int s;
    double f;
    double mc [64 + 1] = {0};
    double mc2 [64 + 1] = {0};
    double mclf [64 + 1] = {0};
    double mc2lf [64 + 1] = {0};
    double mcsf [64 + 1] = {0};
    double pf [64 + 1] = {0};
    double frequency;
    char mfile [64];
    char cfile [64];
    char pfile [64];
    int m1, d1;
    double fac [1024 + 1] = {0};

    if (argc != 4)
    {
        cerr << "chains <L> <d> <f>" << endl;
        return -1;
    }

    L = atoi (argv [1]);
    d = atoi (argv [2]);
    f = atof (argv [3]);
    m = (int) L*d*f;

    if (m > 64)  // size of mcsf
    {
        cerr << "m can be upto 64" << endl;
        return -1;
    }

    if (L*d > 1024) // size of factorial
    {
        cerr << "L*d can be upto 1024" << endl;
        return -1;
    }

    if ((1 << L) > 1024*1024)
    {
        cerr << "1 << L can be upto 2^20" << endl;
        return -1;
    }

    sprintf (cfile, "L%dd%df%g", L, d, f); 
    fp1 = fopen (cfile, "w");
    cout << "generating L = " << L << " d = " << d << " f = " << f << " m = " << m << " in file " << cfile << endl;

    fprintf (fp1, "L=%d d=%d f=%g m=%d\n", L, d, f, m);

    sprintf (mfile, "../arrangements/m%dd%d", m, d);
    fp2 = fopen (mfile, "r");

    if (!fp2)
    {
        cerr << "file " << mfile << " not found" << endl;
        return -1;
    }

    fscanf (fp2, "m=%d d=%d\n", &m1, &d1);

    if (m1 != m ||
        d1 != d)
    {
        cerr << "m, d mismatch" << endl;
        return -1;
    }

    // compute factorials
    factorial (fac, 1024);

    // load the m's arrangement frequencies

    for (int i = 0; i <= m; i++)
    {
        int i1;
        double mci;
	double mc2i;

        fscanf (fp2, "%d = %lg %lg", &i1, &mci, &mc2i);

        mc [i] = mci;
	mc2 [i] = mc2i;

        if (i1 != i)
        {
            cerr << " i, i1 mismatch" << endl;
        }
    }


    // now for each possible arrangement of m,
    // find it's effective chain length, and load its frequency

    for (unsigned int i = 0; i < (1 << L); i++)
    {
        char buf [64];

        r [i].chain = i;
        r [i].slots = slots (i);
        r [i].slot_frequency = mc [r [i].slots] + mc2 [r [i].slots];

	int lmax = 0;
	int lmax2 = 0;
	
	for (int p = 0; p < L; p++)
	{
	    int l1;

	    l1 = length (i, L, &p);
	    
	    if (l1 > lmax)
	    {
		lmax = l1;
	    }

	    int l2;

	    l2 = L - p + l1;
	    
	    if (l2 > lmax2)
	    {
		lmax2 = l2;
	    }
	}

        r [i].length = lmax;
	// if slot freq is 0, len freq has to be 0
	r [i].length_frequency = r [i].slot_frequency ? mc [r [i].slots] : 0;
        r [i].length2 = lmax2;
	// if slot freq is 0, len freq has to be 0
	r [i].length_frequency2 = r [i].slot_frequency ? mc2 [r [i].slots] : 0;

        // add up frequency of chains of a given length
        mcsf [r [i].slots] += r [i].slot_frequency;
        mclf [r [i].length] += r [i].length_frequency;
	mc2lf [r [i].length2] += r [i].length_frequency2;

        fprintf (fp1, "0x%x (%s)  %d %lg  %d %lg  %d %lg\n",
                 r [i].chain,
                 fmtbits (r [i].chain, L, buf),
                 r [i].slots,
                 r [i].slot_frequency,
                 r [i].length,
                 r [i].length_frequency,
                 r [i].length2,
                 r [i].length_frequency2);
    }

    // sanity check
    frequency = 0;
    for (unsigned int i = 0; i < (1 << L); i++)
    {
        frequency += r [i].slot_frequency;
    }

    frequency *= fac [m];
    double facld = fac [L*d]/fac [L*d - m];

    if (frequency != facld)
    {
        cerr << "facld " << facld << " frequency " << frequency << endl;
    }


    sprintf (pfile, "S:L%dd%df%g", L, d, f); 
    fp3 = fopen (pfile, "w");

    
    double Ldm = 1;

    for (int i = 0; i < m; i++)
    {
        Ldm *= L*d - i;
    }


    double p = 0;
    double pfL = 0;

    for (unsigned int i = 0; i <= L; i++)
    {
        pf [i] = (mclf [i] / Ldm / (L - i + 1)) * fac [m];

        p += pf [i];
    }

    pfL += (mc2lf [L] / Ldm) * fac [m];

    p += pfL;


#define TOLERANCE 0.00001
#define FLOAT_EQ(x, y) ((x) > (y) ? (x) - (y) < TOLERANCE : (y) - (x) < TOLERANCE)

    // if very near 1, assume p is 1
    if (FLOAT_EQ (1, p)) p = 1;

    double anonymity = 1 - p;
    
    cout << "anonymity " << anonymity << endl;

    if (anonymity < 0 || anonymity > 1)
    {
	cerr << "bad anonymity " << anonymity << endl;
    }

    fprintf (fp3, "L=%d d=%d f=%g m=%d p=%lg a=%lg\n", L, d, f, m, p, anonymity);

    for (unsigned int i = 0; i <= L; i++)
    {
        fprintf (fp3, "%d %lg    %lg\n",
                 i, mclf [i], pf [i]);
    }

    fprintf (fp3, "B %lg    %lg\n",
	     mc2lf [L], pfL);


    fclose (fp1);
    fclose (fp2);
    fclose (fp3);
}