feat: add a new chop statement to remove floating point numbers below a given threshold.

- The threshold can be specified as a float, integer or rational.

Author: cbmarini

check/features.frm

@@ -1415,6 +1415,24 @@ assert result("ATANH") =~ expr("
        - 6.08698087464190136361e-01*atanh( - 5.4321e-01)
 ")
 *--#] evaluate_atanh : 
+*--#[ chop :
+#StartFloat 15d
+Symbol a,b,c;
+CFunction f,g;
+Local F = 1.0e-10*a+4.7*b+.2e-10*c+f(2,0.000001)+g(2.,1);
+Chop 1.0e-10;
+Argument f;
+	Chop 1/10000;
+EndArgument;
+Argument g;
+	Chop 7;
+EndArgument;
+Print;
+.end
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("4.7e+00*b + 1.0e-10*a + f(2,0) + g(0,1)")
+*--#] chop : 
 *--#[ float_error :
 Evaluate;
 ToFloat;

sources/compiler.c

@@ -117,6 +117,9 @@ static KEYWORD com2commands[] = {
 	,{"chainin",        (TFUN)CoChainin,          STATEMENT,    PARTEST}
 	,{"chainout",       (TFUN)CoChainout,         STATEMENT,    PARTEST}
 	,{"chisholm",       (TFUN)CoChisholm,         STATEMENT,    PARTEST}
+#ifdef WITHFLOAT
+	,{"chop",           (TFUN)CoChop,             STATEMENT,    PARTEST}
+#endif
 	,{"cleartable",     (TFUN)CoClearTable,       DECLARATION,  PARTEST}
 	,{"clearflag",      (TFUN)CoClearUserFlag,    STATEMENT,    PARTEST}
 	,{"collect",        (TFUN)CoCollect,          SPECIFICATION,PARTEST}

sources/declare.h

@@ -1732,8 +1732,10 @@ int CoEvaluate(UBYTE *);
 int PrintFloat(WORD *fun,int numdigits);
 int CoToRat(UBYTE *);
 int CoToFloat(UBYTE *);
+int CoChop(UBYTE *);
 int ToRat(PHEAD WORD *,WORD);
 int ToFloat(PHEAD WORD *,WORD);
+int Chop(PHEAD WORD *, WORD);
 WORD FloatFunToRat(PHEAD UWORD *,WORD *);
 int AddWithFloat(PHEAD WORD **,WORD **);
 int MergeWithFloat(PHEAD WORD **,WORD **);

sources/float.c

@@ -1348,6 +1348,115 @@ int CoToRat(UBYTE *s)
 
 /*
  		#] CoToRat : 
+ 		#[ CoChop :
+
+		LHS notation of the chop statement: 
+			TYPECHOP, length, FLOATFUN, ... 
+		where FLOATFUN, ... represents the threshold of the chop statement in 
+		the notation of a float_ function with its arguments. 
+
+*/
+
+int CoChop(UBYTE *s)
+{
+	GETIDENTITY
+	UBYTE *ss, c;
+	WORD *w, *OldWork;
+	int spec;
+	if ( AT.aux_ == 0 ) {
+		MesPrint("&Illegal attempt to chop a float_ without activating floating point numbers.");
+		MesPrint("&Forgotten %#startfloat instruction?");
+		return(1);
+	}
+	if ( *s == 0 ) {
+		MesPrint("&Chop needs a number (float or rational) as an argument.");
+		return(1);
+	}
+	/* Create TYPECHOP header */
+	w = OldWork = AT.WorkPointer;
+	*w++ = TYPECHOP; 
+	w++;
+
+	while ( *s == ' ' || *s == ',' || *s == '\t' ) s++;
+	
+	/* 
+		The argument of chop can be 
+	 	1: a floating-point number
+	 	2: an integer or rational number
+	 */
+	if ( FG.cTable[*s] == 1 || *s == '.' ) { 
+		/* 1: Attempt to parse as floating-point number */
+		ss = CheckFloat(s, &spec);
+		if ( ss > s ) {
+			/* CheckFloat found a valid float */
+			if ( spec == 1 ) { /* Zero */
+				MesPrint("&The argument in Chop needs to be larger than 0.");
+				return(1);
+			}
+			else if ( spec == -1 ) { 
+				MesPrint("&The floating point system has not been started.");
+				return(1);
+			}
+			else {
+				AT.WorkPointer = w;
+				/* 
+					Reads the floating point number and outputs it at AT.WorkPointer as if it were a float_ 
+					function with its arguments.
+				 */
+				ReadFloat((SBYTE *)s); 
+				s = ss;
+				w += w[1];
+			}
+		}
+		else {
+			/* 2: CheckFloat didn't find a float, we now try for integers and rationals */
+			/* Parse the integer part (numerator for rationals) */
+			ss = s;
+			while ( FG.cTable[*s] == 1 ) s++;  /* Skip over digits */
+			c = *s; *s = 0;                    
+			gmp_sscanf((char *)ss, "%Ff", aux1); /* Parse into GMP float */
+			if ( mpf_sgn(aux1) == 0 ) {
+				MesPrint("&The argument in Chop needs to be larger than 0.");
+				return(1);
+			}
+			*s = c; 
+			while ( *s == ' ' || *s == '\t' ) s++;
+			/* Check for rational number */
+			if ( *s == '/' ) {
+				s++; 
+				while ( *s == ' ' || *s == '\t' ) s++;
+				/* Parse the denominator */
+				if ( FG.cTable[*s] == 1 ) {
+					ss = s;
+					while ( FG.cTable[*s] == 1 ) s++;  /* Skip over digits */
+					c = *s; *s = 0;                  
+					gmp_sscanf((char *)ss, "%Ff", aux2); 
+					*s = c; 
+					if ( mpf_sgn(aux2) == 0 ) {
+						MesPrint("&Division by zero in chop statement.");
+						return(1);
+					}
+					/* Perform the division */
+					mpf_div(aux1, aux1, aux2);
+				}
+			}
+			/* Put aux1 in the notation of a float_ function */
+			PackFloat(w, aux1);
+			w += w[1];
+		}
+	}
+	if ( *s ) {
+		MesPrint("&Illegal argument(s) in Chop statement: '%s'",s);
+		return(1);
+	}
+	AT.WorkPointer = OldWork;
+	AT.WorkPointer[1] = w - AT.WorkPointer;  /* Set total length */
+	AddNtoL(AT.WorkPointer[1],AT.WorkPointer); /* Add the LHS to the compiler buffer */
+	return(0);
+}
+
+/*
+ 		#] CoChop : 
  		#[ ToFloat :
 
 		Converts the coefficient to floating point if it is still a rat.
@@ -1416,6 +1525,46 @@ int ToRat(PHEAD WORD *term, WORD level)
 
 /*
  		#] ToRat : 
+	 	#[ Chop :
+
+	Removes terms with a floating point number smaller than a given threshold. 
+ 
+	Search for a FLOATFUN and compares its absolute value against the threshold 
+	specified in the chop statement. This threshold can be obtained from the 
+	LHS of the chop statement in the compiler buffer.
+ 
+ */
+int Chop(PHEAD WORD *term, WORD level)
+{
+	WORD *tstop, *t, nsize, *threshold; 
+	CBUF *C = cbuf+AM.rbufnum;
+	/* Find the float which should be at the end. */
+	tstop = term + *term; 
+	nsize = ABS(tstop[-1]); tstop -= nsize;
+	t = term+1;
+	while ( t < tstop ) {
+		if ( *t == FLOATFUN && t + t[1] == tstop && TestFloat(t) &&
+		nsize == 3 && tstop[0] == 1 && tstop[1] == 1 ) break;
+		t += t[1];
+	}
+	if ( t < tstop ) {
+		/* Get threshold value from compiler buffer */
+		threshold = C->lhs[level];
+		threshold += 2;  /* Skip TYPECHOP header */
+		UnpackFloat(aux5, threshold); 
+		
+		/* Extract float and compute its absolute value */
+		UnpackFloat(aux4, t);
+		mpf_abs(aux4, aux4);
+		
+		/* Remove if < threshold */
+		if ( mpf_cmp(aux4, aux5) < 0 ) return(0);
+	}
+	return(Generator(BHEAD term,level));
+}
+
+/*
+ 		#] Chop : 
   	#] Float Routines : 
   	#[ Sorting :
 

sources/ftypes.h

@@ -605,6 +605,7 @@ typedef int (*TFUN1)();
 #define TYPEEXPAND 88
 #define TYPETOFLOAT 89
 #define TYPETORAT 90
+#define TYPECHOP 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 TYPECHOP:
+					AT.WorkPointer = term + *term;
+					if ( Chop(BHEAD term,level) ) goto GenCall;
+					goto Return0;
 #endif
 				}
 				goto SkipCount;