Untitled

#include "probabilitycalculator.h"//#include <xpre.h>#include <mpfr.h>#include <stdio.h>const mp_prec_t floating_precision = 2048;const mp_rnd_t  rounding_mode = GMP_RNDN;#define FLOAT_INIT(x) mpfr_init2(x,floating_precision)#define FLOAT_UNIT(x) mpfr_clear(x)#ifdef DEBUGlong double viewA;long double viewB;long double viewC;long double viewD;long double viewE;long double viewF;#endif//mpfr_add (mpfr_t rop, mpfr_t op1, mpfr_t op2, mp_rnd_t rnd//Function: int mpfr_sub (mpfr_t rop, mpfr_t op1, mpfr_t op2, mp_rnd_t rnd//Function: int mpfr_mul (mpfr_t rop, mpfr_t op1, mpfr_t op2, mp_rnd_t rnd//Function: int mpfr_div (mpfr_t rop, mpfr_t op1, mpfr_t op2, mp_rnd_t rnd//mpfr_pow (mpfr_t rop, mpfr_t op1, mpfr_t op2, mp_rnd_t rnd//mpfr_cmp (mpfr_t op1, mpfr_t op2/*Compare op1 and op2. Return a positive value if op1 > op2, zero if op1 = op2, and a negative value if op1 < op2. Both op1 and op2 are considered to their full own precision, which may differ. If one of the operands is NaN (Not-a-Number), return zero and set the erange flag.*///mpfr_exp (mpfr_t rop, mpfr_t op, mp_rnd_t rnd//mpfr_fac_ui (mpfr_t rop, unsigned long int op, mp_rnd_t rndstruct floatStore{	mpfr_t value;};floatStore calculateProbabilityUsingBinomialDistribution(unsigned long long numberOfTries, double singleSuccessRate, unsigned long long numberSucceeding){	//			N!	//	P(r) = -------------  p^r * (1-p)^(N-r) 					//		r!(N - r)!	//	P = This function	//	N = numberOfTries	//	r = numberSucceeding	//	p = singleSuccessRate				floatStore partA;	//mpfr_t partA;	mpfr_t partA_a;	mpfr_t partA_b;	mpfr_t partA_c;	FLOAT_INIT(partA.value);		FLOAT_INIT(partA_a);		FLOAT_INIT(partA_b);		FLOAT_INIT(partA_c);		mpfr_fac_ui(partA_a,numberOfTries,rounding_mode);	mpfr_fac_ui(partA_b, numberSucceeding,rounding_mode);	mpfr_fac_ui(partA_c, numberOfTries - numberSucceeding,rounding_mode);	mpfr_mul(partA_b,partA_b,partA_c,rounding_mode);	mpfr_div(partA.value,partA_a,partA_b,rounding_mode);	FLOAT_UNIT(partA_a);	FLOAT_UNIT(partA_b);	FLOAT_UNIT(partA_c);#ifdef DEBUG	viewA = mpfr_get_ld(partA.value,rounding_mode);#endif	mpfr_t partB;	mpfr_t partB_exp;	FLOAT_INIT(partB);		FLOAT_INIT(partB_exp);	mpfr_set_d(partB,singleSuccessRate,rounding_mode);	mpfr_set_d(partB_exp,singleSuccessRate,rounding_mode);	//mpfr_set_ui(partB_exp,numberSucceeding,rounding_mode);	//mpfr_exp_ui(partB,partB,numberSucceeding,rounding_mode);	if (numberSucceeding > 1)	{		for (unsigned long long i = 1; i < numberSucceeding; i++)		{			mpfr_mul(partB,partB,partB_exp,rounding_mode);		}	}	FLOAT_UNIT(partB_exp);#ifdef DEBUG	viewB = mpfr_get_ld(partB,rounding_mode);#endif	mpfr_t partC;	mpfr_t partC_a;	FLOAT_INIT(partC);	FLOAT_INIT(partC_a);	mpfr_set_d(partC, 1 - singleSuccessRate,rounding_mode);	mpfr_set_d(partC_a, 1 - singleSuccessRate,rounding_mode);	//mpfr_set_ui(partC_a, numberOfTries - numberSucceeding,rounding_mode);	//mpfr_exp(partC, partC_a,rounding_mode);	if ((numberOfTries - numberSucceeding)> 1)	{		for (unsigned long long i = 1; i < (numberOfTries - numberSucceeding); i++)		{			mpfr_mul(partC,partC,partC_a,rounding_mode);		}	}	else	{		if (numberOfTries - numberSucceeding == 0)		{			mpfr_set_ui(partC,1,rounding_mode);		}	}	FLOAT_UNIT(partC_a);#ifdef DEBUG	viewC = mpfr_get_ld(partC,rounding_mode);#endif	mpfr_mul(partB,partB,partC,rounding_mode);#ifdef DEBUG	viewB = mpfr_get_ld(partB,rounding_mode);#endif	mpfr_mul(partA.value,partA.value,partB,rounding_mode);#ifdef DEBUG	viewA = mpfr_get_ld(partA.value,rounding_mode);#endif	FLOAT_UNIT(partB);	FLOAT_UNIT(partC);	return partA;}		unsigned long long ProbabilityCalculator::SolveBinomialDistributionForProbabilityOfDesiredAtLeastSuccesses(							unsigned long long desiredMinimumSuccesses,							double desiredProbability,							double singleSuccessRate){	if (desiredProbability >= 1)		return 0;	unsigned long long requiredForSuccess = 0;	mpfr_t currentProbability;	mpfr_t desiredProbabilityX;	FLOAT_INIT(currentProbability);	FLOAT_INIT(desiredProbabilityX);	mpfr_set_ui(currentProbability,0,rounding_mode);	mpfr_set_d(desiredProbabilityX,desiredProbability,rounding_mode);		#ifdef DEBUG		viewE = mpfr_get_ld(currentProbability,rounding_mode);		viewF = mpfr_get_ld(desiredProbabilityX,rounding_mode);	#endif	while (mpfr_cmp(currentProbability,desiredProbabilityX) < 0)	{		requiredForSuccess++;		mpfr_set_ui(currentProbability,0,rounding_mode);		#ifdef DEBUG			viewE = mpfr_get_ld(currentProbability,rounding_mode);		#endif		for (unsigned long long i = requiredForSuccess; i >= desiredMinimumSuccesses; i--)		{			floatStore thisProbability = calculateProbabilityUsingBinomialDistribution(							requiredForSuccess,							singleSuccessRate,							i);			mpfr_add(currentProbability ,currentProbability,thisProbability.value,rounding_mode);			#ifdef DEBUG				viewD = mpfr_get_ld(thisProbability.value,rounding_mode);				viewE = mpfr_get_ld(currentProbability,rounding_mode);			#endif			FLOAT_UNIT(thisProbability.value);					}		if (requiredForSuccess > 4096)		{			mpfr_set_ui(currentProbability,2,rounding_mode);		}	}	FLOAT_UNIT(currentProbability);	FLOAT_UNIT(desiredProbabilityX);	printf(		"minimum: %u\nsingle probability: %f\ndesired probability: %f\n  required for success:   %u\n\n",		(unsigned int)desiredMinimumSuccesses,		singleSuccessRate,		desiredProbability,		(unsigned int)requiredForSuccess);	return requiredForSuccess;}