Chop statement for floats

check/features.frm

@@ -1435,6 +1435,47 @@ assert succeeded?
 assert result("F1") =~ expr("1.23457e-04 + f(1.0e+00) + f(1.0e+00)")
 assert result("F2") =~ expr("1.235e-04 + 2*f(1.0e+00)")
 *--#] strictrounding : 
+*--#[ chop :
+#StartFloat 15d
+#StartFloat 15d
+Symbol x;
+Local F1 = 4.7*x-1.0e-10*x^2+.2e-10*x^3-0.00005*x^4+x^5/1000000;
+Local F2 = 4.7*x-1.0e-10*x^2+.2e-10*x^3-0.00005*x^4+x^5/1000000;
+Local F3 = 4.7*x-1.0e-10*x^2+.2e-10*x^3-0.00005*x^4+x^5/1000000;
+Local F4 = 4.7*x-1.0e-10*x^2+.2e-10*x^3-0.00005*x^4+x^5/1000000;
+Local F5 = 4.7*x-1.0e-10*x^2+.2e-10*x^3-0.00005*x^4+x^5/1000000;
+Print;
+.sort
+Skip;
+NSkip F1;
+Chop 1.0e-10;
+.sort
+Skip;
+NSkip F2;
+Chop 1/10000;
+.sort
+Skip;
+NSkip F3;
+Chop 7;
+.sort
+Skip;
+NSkip F4;
+Chop 10^-6;
+.sort
+Skip;
+NSkip F5;
+Chop 10^6;
+.sort
+Print;
+.end
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F1") =~ expr("4.7e+00*x - 1.0e-10*x^2 - 5.0e-05*x^4 + 1/1000000*x^5")
+assert result("F2") =~ expr("4.7e+00*x + 1/1000000*x^5")
+assert result("F3") =~ expr("1/1000000*x^5")
+assert result("F4") =~ expr("4.7e+00*x - 5.0e-05*x^4 + 1/1000000*x^5")
+assert result("F5") =~ expr("1/1000000*x^5")
+*--#] chop : 
 *--#[ float_error :
 Evaluate;
 ToFloat;

doc/manual/float.tex

@@ -89,6 +89,18 @@ \chapter{Floating point}
 $1.1100110101011111101*2^{-14}$. When rounded to 5 bits, this becomes 
 $1.1101*2^{-14}$, which in decimal digits appears as 
 1.10626220703125e-04.
+\item[Chop] This statement removes floating point numbers that are smaller 
+in absolute magnitude than a specified threshold. It takes one argument delta:
+\begin{verbatim}
+    Chop <delta>;
+\end{verbatim}
+All floating point numbers with absolute value less than delta are replaced by 0. 
+Terms with no floating point coefficient are left untouched. The threshold delta 
+can be a floating point number, integer, rational number, or power. Because 
+statements in \FORM{} act term by term, it is often important to sort before invoking the 
+chop statement. Otherwise, terms might be removed individually, while after 
+sorting and combining, their combined floating point coefficient could exceed 
+the specified chop threshold. 
 \item[Format floatprecision] This instruction controls how many digits are 
 displayed when printing floating point numbers. It only affects output 
 formatting and does not influence the internal precision or accuracy of 

doc/manual/statements.tex

@@ -600,6 +600,18 @@ \section{chisholm}
 more details, see chapter \ref{gammaalgebra} on gamma\index{gamma algebra} algebra. \vspace{10mm}
 
 %--#] chisholm : 
+%--#[ chop :
+\section{chop}
+\label{substaevaluate}
+
+\noindent \begin{tabular}{ll}
+Type & Executable statement\\
+Syntax & chop $<$threshold$>$; \\ 
+\end{tabular} \vspace{4mm}
+
+\noindent See chapter~\ref{floatingpoint} on the floating point capability.
+
+%--#] chop : 
 %--#[ cleartable :
 
 \section{cleartable}

sources/compiler.c

@@ -117,8 +117,11 @@ static KEYWORD com2commands[] = {
 	,{"chainin",        (TFUN)CoChainin,          STATEMENT,    PARTEST}
 	,{"chainout",       (TFUN)CoChainout,         STATEMENT,    PARTEST}
 	,{"chisholm",       (TFUN)CoChisholm,         STATEMENT,    PARTEST}
-	,{"cleartable",     (TFUN)CoClearTable,       DECLARATION,  PARTEST}
+#ifdef WITHFLOAT
+	,{"chop",           (TFUN)CoChop,             STATEMENT,    PARTEST}
+#endif
 	,{"clearflag",      (TFUN)CoClearUserFlag,    STATEMENT,    PARTEST}
+	,{"cleartable",     (TFUN)CoClearTable,       DECLARATION,  PARTEST}
 	,{"collect",        (TFUN)CoCollect,          SPECIFICATION,PARTEST}
 	,{"commuteinset",   (TFUN)CoCommuteInSet,     DECLARATION,  PARTEST}
 	,{"contract",       (TFUN)CoContract,         STATEMENT,    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

@@ -1395,6 +1395,120 @@ int CoStrictRounding(UBYTE *s)
 } 
 /*
  		#] CoStrictRounding : 
+ 		#[ 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, pow = 1;
+	unsigned long x;
+	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, rational number or power
+	 */
+	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 rationals and powers */
+			/* Parse the integer part (numerator for rationals) */
+			if ( FG.cTable[*s] == 1 ) {
+				ParseNumber(x,s)
+				mpf_set_ui(aux1,x);
+			}
+			if ( mpf_sgn(aux1) == 0 ) {
+				MesPrint("&The argument in Chop needs to be larger than 0.");
+				return(1);
+			}
+			while ( *s == ' ' || *s == '\t' ) s++;
+			/* Check for rational number or power*/
+			if ( *s == '/' || *s == '^' ) {
+				c = *s; s++; 
+				while ( *s == ' ' || *s == '\t' ) s++;
+				if ( *s == '-' ) { s++; pow = -1; } /* negative power */
+				/* Parse the denominator or power */
+				if ( FG.cTable[*s] == 1 ) {
+					ParseNumber(x,s)
+					if ( c == '/' ) { /* rational */
+						if ( x == 0 ) {
+							MesPrint("&Division by zero in chop statement.");
+							return(1);
+						}
+						/* Perform the division */
+						mpf_div_ui(aux1, aux1,x);
+					}
+					else { /* Power */
+						mpf_pow_ui(aux1,aux1,x);
+						if ( pow == -1 ) {
+							mpf_ui_div(aux1,(unsigned long) 1, aux1);
+						}
+					}
+				}
+			}
+			/* 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.
@@ -1521,6 +1635,46 @@ int StrictRounding(PHEAD WORD *term, WORD level, WORD prec, WORD base) {
 }
 /*
  		#] StrictRounding : 
+	 	#[ 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

@@ -606,6 +606,7 @@ typedef int (*TFUN1)();
 #define TYPETOFLOAT 89
 #define TYPETORAT 90
 #define TYPESTRICTROUNDING 91
+#define TYPECHOP 92
 #endif
 
 /*

sources/proces.c

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