/* Copyright Julian D. A. Wiseman 2005 and 2010. This version 19:00 Mon 16 August 2010. */
/* This program may be distributed under the terms of the GNU General Public License,    */
/* currently available at www.gnu.org/licenses/gpl.txt. It is requested, but not         */
/* required, that authors of modifcations or improvements inform the original author     */
/* of these, contactable via www.jdawiseman.com/papers/easymath/chow_robbins.html.       */

/* Flip a coin repeatedly, stopping any time after the first toss. */
/* Score the proportion of throws that are heads: h/(h+t)          */
/* With an optimal strategy, what is the expected score?           */
/* Answer: 0.79295350640...                                        */
/* Also, is h=5 t=3 a stopping state. That depends on the value of */
/* the game at h=5 t=4, estimated here. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include <time.h>

#define arraysize (67108864L)  /* At 16 bytes per long double this uses 1G memory */

typedef long double realtype;
#define epsilonrealtype (LDBL_EPSILON)
#define longtypeprintingcode "%0.55Lf"

/* There is an occasional print to the screen, to confirm process is alive, this often */
#define logfrequency 1073741824 /* 2^30 */

#define hProbe 5
#define tProbe 4


int main()
{
	FILE *fp1, *fpProbe;
	int i, start_n;
	register signed long j,t,h;
	signed long max_t, count=0, nextmaxh, maxh, prevmaxh, prevmaxh_minus_maxh, j_plus_prevmaxh_minus_maxh;
	realtype *ev, d;

	float multiplier_n;
	realtype prevev00=0.75, prevprevev00=0.5;  /* Any values would do, though 0.75 and 0.5 are natural choices */
	realtype prevevht=0.02, prevprevevht=0.01, limitht, prevlimitht=0.03;  /* Any values would do, there not being natural choices */
	char filename1[127], filenameProbe[127];
	clock_t TimeCPU;


	printf("sizeof(float)=%lu   sizeof(double)=%lu   sizeof(long double)=%lu   used size=%lu\n", sizeof(float), sizeof(double), sizeof(long double), sizeof(realtype) );

	/* The next flag can, by setting hard-wired defaults, avoid having any user input */
#if false
	start_n = 2;  multiplier_n = 2;
#else
	do {
		printf("\nStarting value of n, >=2? ");
		scanf ("%d",&start_n);
	} while( start_n < 2 );
	do {
		printf("\nEach loop n multiplies by a number >= %f. What number? ", (float)(1.0+(1.0/((float)start_n))) );
		scanf ("%f",&multiplier_n);
	} while( ((int)(start_n * multiplier_n)) <= start_n  );
#endif
	printf("start_n = %i        multiplier_n = %f\n", start_n, multiplier_n);
	

	sprintf(filename1,     "coin_%i_%f_scores.txt", start_n, multiplier_n);
	sprintf(filenameProbe, "coin_%i_%f_h%i_t%i_scores.txt", start_n, multiplier_n, hProbe, tProbe);

	printf("Opening '%s'.\n",filename1);
	if(NULL==(fp1=fopen(filename1, "w")))
		{printf("fopen problem with '%s'.",filename1); fprintf(stderr,"fopen problem with '%s'.",filename1); exit(EXIT_FAILURE);}

	printf("Opening '%s'.\n",filenameProbe);
	if(NULL==(fpProbe=fopen(filenameProbe, "w")))
		{printf("fopen problem with '%s'.",filenameProbe); fprintf(stderr,"fopen problem with '%s'.",filenameProbe); exit(EXIT_FAILURE);}
	
	fprintf(fp1,  /* scores */
				"i"
		"\t"	"max_t"
		"\t"	"EV"
		"\t"	"DiffEV"
		"\t"	"RatioDiffEV"
		"\t"	"EstimatedLimit"
		"\n"
	);

	fprintf(fpProbe,  /* scores */
			"i"
			"\t"	"max_t"
			"\t"	"EV at h=%i t=%i"
			"\t"	"DiffEV"
			"\t"	"RatioDiffEV"
			"\t"	"EstimatedLimit"
			"\t"	"DiffEstimatedLimit"
			"\t"	"DiffEstimatedLimit*max_t"
			"\t"	"clock()"
			"\n",
		hProbe, tProbe);

	if( NULL == (ev = malloc( (arraysize) * sizeof(realtype) ) ) )
	{
		printf("malloc problem with ev"); fprintf(stderr,"malloc problem with ev");
		fclose(fp1); fclose(fpProbe); exit(EXIT_FAILURE);
	}

	max_t = start_n ;
	i = 1;  /* i used only to label output rows */
	while( i < INT_MAX )
	{
		nextmaxh = maxh = max_t;

		/* initialise rightmost column, that being t==max_t */
		t = max_t;
		for( j=0 ; j<arraysize && j<=max_t ; j++ )
		{
			ev[j] = ((realtype)0.5);
		}  /* for(j...), that being the right-column initialisation */ 

		/* Now compute all columns except the rightmost, being t<max_t */
		for( t-- ; t>=0 ; t-- )
		{
			prevmaxh = maxh; 
			maxh = nextmaxh; 
			prevmaxh_minus_maxh = prevmaxh - maxh; 

			/* last row */
			h=maxh;  /* which is not less than t. Implicit j=0; h=maxh-j; */
			ev[0] = ((realtype)h) / (h+t); 

			/* rows between the last and first */
			count += h;  /* We add the full amount to go, and then subtract those undone, calculating count without having it in the innermost loop */
			for( j=1 ;  ; j++ )  /* testing whether to stop */
			{
				if( j>= arraysize )
				{
					printf("Exiting(A) with i=%i max_t=%li t=%li h=%li j=%li>=%li\n",i,max_t,t,h,j,arraysize);
					fclose(fp1); fclose(fpProbe);
					free(ev);
					exit(EXIT_SUCCESS);
				} 

				h = maxh - j;

				if( h <= 0 )
					break;  /* don't do top row to avoid possible divide by zero */
				
				j_plus_prevmaxh_minus_maxh = j + prevmaxh_minus_maxh;
				ev[j] = ( ev[j-1] + ev[j_plus_prevmaxh_minus_maxh] ) / 2;
				
				d = ((realtype)h) / (h+t); 
				if( ev[j] < d )
				{
					ev[j] = d;  /* Optimal strategy requires that a player stops */
				}
				else
				{   /* So no need to test for any h from here down to 1, that is, for any larger j */
					nextmaxh = h+1; /* reset because don't need to test the larger h's for smaller t, as they will all be stopping states. */

					break;  /* break this j loop, as no further need to test stopping condition */
				}
			}  /* end for( j ... ) */

			/* Know that the optimal strategy for all h from here to 0 is to toss again. */
			/* Test whether machine precision limit has been reached, */
			/* and if it has, leave ev(t,h) equal to ev(t+1,h-prevmaxh_minus_maxh) for all smaller h */
			/* For speed, the t==tProbe is done outside the j loop, lengthening the code by an extra near-copy of the j loop. */
			if( t==tProbe )
				for( j++ ;  ; j++ ) 
				{
					if( j>= arraysize )
					{
						printf("Exiting(B) with i=%i max_t=%li t=%li h=%li j=%li>=%li\n",i,max_t,t,h,j,arraysize);
						fclose(fp1); fclose(fpProbe);
						free(ev);
						exit(EXIT_SUCCESS);
					} 
					h = maxh - j;
					if( h <= 0 )
						break;
					ev[j] = ( ev[j-1] + ev[j + prevmaxh_minus_maxh] ) / 2;
					if(h==hProbe)
					{
						printf("$$$ max_t=%li t=%li h=%li EV=" longtypeprintingcode "\n",max_t,t,h,ev[j]);
						limitht = ev[j] + (ev[j]-prevevht) * (prevevht-ev[j])/(ev[j] + prevprevevht - 2*prevevht);
						TimeCPU = clock();
						fprintf( fpProbe,
								"%i\t%li\t" longtypeprintingcode "\t" longtypeprintingcode "\t" longtypeprintingcode "\t" longtypeprintingcode "\t" longtypeprintingcode "\t%.19E\t%.2f",
								i,
								max_t,
								ev[j],
								ev[j]-prevevht,
								(ev[j]-prevevht)/(prevevht-prevprevevht),
								limitht,
								limitht - prevlimitht,
								( (double) (max_t * (limitht - prevlimitht)) ),
								(double)(TimeCPU + 0.0)
								);
						fprintf( fpProbe, "\n" );
						fflush(fpProbe);
						prevprevevht = prevevht;
						prevevht = ev[j];
						prevlimitht = limitht;
					}  /* if(h==hProbe) */
					if( ev[j] <= ((realtype)0.5) + epsilonrealtype ) /* Machine precision limit, so no need to update for smaller h */
						goto machine_precision_reached;
				}  /* end first copy of second for( j ... ) which is over the keep-playing states */
			else  /* t != tProbe */
				for( j++ ;  ; j++ ) 
				{
					if( j>= arraysize )
					{
						printf("Exiting(B) with i=%i max_t=%li t=%li h=%li j=%li>=%li\n",i,max_t,t,h,j,arraysize);
						fclose(fp1); fclose(fpProbe);
						free(ev);
						exit(EXIT_SUCCESS);
					} 
					h = maxh - j;
					if( h <= 0 )
						break;
					ev[j] = ( ev[j-1] + ev[j + prevmaxh_minus_maxh] ) / 2;
					if( ev[j] <= ((realtype)0.5) + epsilonrealtype ) /* Machine precision limit, so no need to update for smaller h */
						goto machine_precision_reached;
				}  /* end second copy of second for( j ... ) which is over the keep-playing states */

			ev[maxh] = ( ev[maxh-1] + ev[maxh + prevmaxh_minus_maxh] ) / 2;

			machine_precision_reached: ;

			count -= h;  /* Number of rows done. this being done outside h loop for speed. */
			/* Possibly output something to screen, to confirm that the program is still alive. */
			if( count >= logfrequency )
			{
				count=0;
				printf("i=%i     t=%li     j=%li     reached machine precision with h >= %li\n",i,t,j,h);
			}
			
		} /* for( t ... ) */

		printf( longtypeprintingcode "\t" "%i" "\t" "%li" "\n",
			ev[maxh],
			i,
			max_t
		);
		fprintf( fp1,
			"%i\t%li\t" longtypeprintingcode "\t" longtypeprintingcode "\t" longtypeprintingcode "\t" longtypeprintingcode,
			i,
			max_t,
			ev[maxh],
			ev[maxh]-prevev00,
			(ev[maxh]-prevev00)/(prevev00-prevprevev00),
			ev[maxh] + (ev[maxh]-prevev00) * (prevev00-ev[maxh])/(ev[maxh] + prevprevev00 - 2*prevev00)
		);
		fprintf( fp1, "\n" );
		fflush(fp1);

		prevprevev00 = prevev00;
		prevev00 = ev[maxh];

		if( max_t >= (float)ULONG_MAX / (float)multiplier_n)
			break; /* max_t would be too big in next pass through loop  */
		max_t = (signed long) ((float)max_t * multiplier_n) ;
		i++;
	} /* while( i < INT_MAX ) */

	free(ev);

	fclose(fp1); fclose(fpProbe);
	exit(EXIT_SUCCESS);
}