#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define CMIN	(-2.0)
#define CMAX	(2.0)

#define XMIN	(-2.0)
#define XMAX	(2.0)


/*
To generate a normal distribution, use the Box-Muller Transformation:

http://mathworld.wolfram.com/Box-MullerTransformation.html

http://math.nist.gov/cgi-bin/gams-serve/GET/Source_of_ZUFALL_in_RANDOM_from_NETLIB/get-module-component/Source/12984/NETLIB.html

*/

#define GP_PI            3.14159265358979323846

double randNormal ()
        {

        static int cycle = 0;

        static double x1, x2;
        static double z1, z2;

        double temp1, temp2;
        double rtn;

        if (cycle == 0)
                {
                x1 = drand48 ();
                x2 = drand48 ();

                temp1 = sqrt (-2.0*log(x1));
                temp2 = 2.0*GP_PI*x2;

                z1 = temp1 * cos (temp2);
                z2 = temp1 * sin (temp2);

                rtn = z1;
                }
        else
                {
                rtn = z2;
                }

        cycle = (cycle + 1) % 2;

        return  rtn;
        }       /* end of randNormal */

double 
getRand (double min, double max)
	{
	return (drand48()*(max - min)) + min;
	}


dbl_compar (const void *a, const void *b)
        {
        return ( *((double *)a) < *((double *)b) ) ? -1 : 1;
        }

main (int argc, char **argv)
	{
	double coeff[1024];
	double sd;
	double x[1024], y[1024], noise[1024];
	double xp;
	int i, j, degree, numPts;
	long lseed;


	degree = atoi (argv[1]);
	numPts = atoi (argv[2]);
	sd = atof (argv[3]);

	if (argc >= 5)
		{
		lseed = atoi (argv[4]);
		srand48 (lseed);
		}

	printf ("num_coeffs: %d\n", degree+1);
	printf ("coeffs_low_to_high:\n");
	for (j = 0; j <= degree; j++)
		{
		coeff[j] = getRand (CMIN, CMAX);
		printf (" %20.12e\n", coeff[j]);
		}

	printf ("noise_standard_deviation: %20.12e\n", sd);
	printf ("num_points: %d\n", numPts);
	printf ("points:\n");

	for (i = 0; i < numPts; i++)
		{
		x[i] = getRand (XMIN, XMAX);
		}



	qsort (x, numPts, sizeof (x[0]), dbl_compar);


	for (i = 0; i < numPts; i++)
		{

		y[i] = 0; 
		xp = 1;
		for (j = 0; j <= degree; j++)
			{
			y[i] += xp * coeff[j];
			xp *= x[i];
			}

		}



	if (sd > 0.0)
		{
		for (i = 0; i < numPts; i++)
			{
			noise[i] = sd * randNormal ();
			y[i] += noise[i];
			}
		}


	for (i = 0; i < numPts; i++)
		{
		printf ("%20.12e %20.12e %20.12e\n", x[i], y[i], noise[i]);
		}



	exit (0);
	}