feat: add StrictRounding statement to round floats and to trim excess GMP precision

- The user can specify a precision (in bits or digits) for the rounding. If no number is given, the default precision is used.
- Ensures that extra internal GMP precision (used for tracking accuracy) is cut off and properly rounded
- Internally the mpf_get_str and mpf_set_str are used. This means we could also implement base 2 rounding.
- refactor: changed AO.floatsize such that it can now also hold floats in base 2.
- This resolves #698

Author: cbmarini

sources/compiler.c

@@ -218,6 +218,9 @@ static KEYWORD com2commands[] = {
 	,{"splitarg",       (TFUN)CoSplitArg,         STATEMENT,    PARTEST}
 	,{"splitfirstarg",  (TFUN)CoSplitFirstArg,    STATEMENT,    PARTEST}
 	,{"splitlastarg",   (TFUN)CoSplitLastArg,     STATEMENT,    PARTEST}
+#ifdef WITHFLOAT
+	,{"strictrounding", (TFUN)CoStrictRounding,   STATEMENT,    PARTEST}
+#endif
 	,{"stuffle",        (TFUN)CoStuffle,          STATEMENT,    PARTEST}
 	,{"sum",            (TFUN)CoSum,              STATEMENT,    PARTEST}
 	,{"switch",         (TFUN)CoSwitch,           STATEMENT,    PARTEST}

sources/declare.h

@@ -1743,6 +1743,8 @@ SBYTE *ReadFloat(SBYTE *);
 UBYTE *CheckFloat(UBYTE *,int *);
 void SetfFloatPrecision(LONG);
 int EvaluateFun(PHEAD WORD *, WORD, WORD *);
+int CoStrictRounding(UBYTE *);
+int StrictRounding(PHEAD WORD *, WORD, WORD, WORD);
 #endif
 
 /*

sources/float.c

@@ -917,7 +917,7 @@ int SetFloatPrecision(WORD prec)
 	else {
 		AC.DefaultPrecision = prec;
 		if ( AO.floatspace != 0 ) M_free(AO.floatspace,"floatspace");
-		AO.floatsize = ((10*prec)/33+20)*sizeof(char);
+		AO.floatsize = (prec+64)*sizeof(char);
 		AO.floatspace = (UBYTE *)Malloc1(AO.floatsize,"floatspace");
 		mpf_set_default_prec(prec);
 		return(0);
@@ -1348,6 +1348,50 @@ int CoToRat(UBYTE *s)
 
 /*
  		#] CoToRat : 
+ 		#[ CoStrictRounding : 
+
+		Syntax: StrictRounding [precision][base]
+		- precision: number of digits to round to (optional)
+		- base: 'd' for decimal (base 10) or 'b' for binary (base 2)
+		
+		If no arguments are provided, uses default precision with binary base.
+*/
+int CoStrictRounding(UBYTE *s)
+{
+	GETIDENTITY
+	WORD x;
+	int base;
+	if ( AT.aux_ == 0 ) {
+		MesPrint("&Illegal attempt for strict rounding without activating floating point numbers.");
+		MesPrint("&Forgotten %#startfloat instruction?");
+		return(1);
+	}
+	while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
+	if ( *s == 0 ) {
+		/* No subkey, which means round to default precision */
+		x = AC.DefaultPrecision - AC.MaxWeight - 1;
+		Add4Com(TYPESTRICTROUNDING,x,2);
+		return(0);
+	}
+	if ( FG.cTable[*s] == 1 ) { /* number */
+		ParseNumber(x,s)
+		if ( tolower(*s) == 'd' ) { base = 10; s++; }      /* decimal base */
+		else if ( tolower(*s) == 'b' ){ base = 2; s++; }  /* binary base */
+		else goto IllPar;  /* invalid base specification */
+	}
+	while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
+	
+	/* Check for invalid arguments */
+	if ( *s ) {
+IllPar:
+		MesPrint("&Illegal argument(s) in StrictRounding statement: '%s'",s);
+		return(1);
+	}
+	Add4Com(TYPESTRICTROUNDING,x,base);
+	return(0);
+} 
+/*
+ 		#] CoStrictRounding : 
  		#[ ToFloat :
 
 		Converts the coefficient to floating point if it is still a rat.
@@ -1416,6 +1460,52 @@ int ToRat(PHEAD WORD *term, WORD level)
 
 /*
  		#] ToRat : 
+ 		#[ StrictRounding : 
+
+		Rounds floating point numbers to a specified precision
+		in a given base (decimal or binary).
+*/
+int StrictRounding(PHEAD WORD *term, WORD level, WORD prec, WORD base) {
+	WORD *tstop, *t, *oldworkpointer = AT.WorkPointer;
+	int retval,sign;
+
+	tstop = term + *term; tstop -= ABS(tstop[-1]);
+	t = term+1;
+	while ( t < tstop ) {
+		if ( *t == FLOATFUN && t + t[1] == tstop ) {
+			char *s;
+			mp_exp_t exp;
+			/* Extract the floating point value */
+			UnpackFloat(aux4,t);
+			/* Convert to string: the generated string is a fraction with an implicit 
+			   radix point immediately to the left of the first digit. 
+			   The applicable exponent is written in exp. */
+			s = mpf_get_str(0,&exp, base, prec, aux4);
+			/* Format as MeN with M the mantissa and N the exponent */
+			snprintf((char *)AO.floatspace,AO.floatsize,".%se%ld",s,exp);
+			/* Negative base values are used to specify that the exponent is in decimal */
+			mpf_set_str(aux4,(char *)AO.floatspace,-base);
+			free(s);
+			break;
+		}
+		t += t[1];
+	}
+	if ( t < tstop ) {
+		/* Pack the rounded floating point value back into the term */
+		PackFloat(t,aux4);
+		t+=t[1];
+		if ( term[*term-1] < 0 ) sign = -1;
+		else sign = 1;
+		*t++ = 1; *t++ = 1; *t++ = 3*sign;
+		*term = t - term;
+	}
+	AT.WorkPointer = t;
+	retval = Generator(BHEAD term,level);
+	AT.WorkPointer = oldworkpointer;
+	return(retval);
+}
+/*
+ 		#] StrictRounding : 
   	#] Float Routines : 
   	#[ Sorting :
 

sources/ftypes.h

@@ -605,6 +605,7 @@ typedef int (*TFUN1)();
 #define TYPEEXPAND 88
 #define TYPETOFLOAT 89
 #define TYPETORAT 90
+#define TYPESTRICTROUNDING 91
 #endif
 
 /*

sources/proces.c

@@ -3990,6 +3990,10 @@ SkipCount:	level++;
 					AT.WorkPointer = term + *term;
 					if ( ToRat(BHEAD term,level) ) goto GenCall;
 					goto Return0;
+				  case TYPESTRICTROUNDING:
+					AT.WorkPointer = term + *term;
+					if ( StrictRounding(BHEAD term,level,C->lhs[level][2],C->lhs[level][3]) ) goto GenCall;
+					goto Return0;
 #endif
 				}
 				goto SkipCount;