fix: fixes the ToRational statement.

- The statement could produce 1/0 for small floats (i.e. rationals close to zero). These are now properly handled and set to zero.

Author: cbmarini

check/features.frm

@@ -1463,6 +1463,73 @@ assert result("ATANH") =~ expr("
        - 6.08698087464190136361e-01*atanh( - 5.4321e-01)
 ")
 *--#] evaluate_atanh : 
+*--#[ torat :
+#-
+#StartFloat 16d
+Symbol a,b,jj;
+Local F1 = 0.0
+		  -10^-20*a
+		  +1023/1024*a^2
+		  -17/1024*a^3
+		  +3602879701896397/2^55*a^4
+		  -10^20*a^5
+		  +17*a^6
+		  ;
+Local F2 = sum_(jj,1,16,(355/113+10^-jj)*a^jj);
+Local F3 = 0.3*a^1+0.33*a^2+0.333*a^3+0.3333*a^4+0.33333*a^5+0.333333*a^6+0.3333333*a^7+0.33333333*a^8+0.333333333*a^9+0.3333333333*a^10+0.33333333333*a^11+0.333333333333*a^12+0.3333333333333*a^13+0.33333333333333*a^14+0.333333333333333*a^15+0.3333333333333333*a^16;
+ToFloat;
+.sort
+
+ToRat;
+.sort
+
+#StartFloat 50d
+Local F4 = 10e-20;
+ToRat;
+Print +s;
+.end
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F1") =~ expr("+ 1023/1024*a^2
+       - 17/1024*a^3
+       + 1/10*a^4
+       - 100000000000000000000*a^5
+       + 17*a^6")
+assert result("F2") =~ expr("+ 3663/1130*a
+       + 35613/11300*a^2
+       + 355113/113000*a^3
+       + 3550113/1130000*a^4
+       + 35500113/11300000*a^5
+       + 355000113/113000000*a^6
+       + 2205176720676/701929469027*a^7
+       + 5798543718053/1845733628322*a^8
+       + 4781671236789/1522053097423*a^9
+       + 492368235747/156725663768*a^10
+       + 355/113*a^11
+       + 355/113*a^12
+       + 355/113*a^13
+       + 355/113*a^14
+       + 355/113*a^15
+       + 355/113*a^16")
+assert result("F3") =~ expr("
+       + 3/10*a
+       + 33/100*a^2
+       + 333/1000*a^3
+       + 3333/10000*a^4
+       + 33333/100000*a^5
+       + 333333/1000000*a^6
+       + 3333333/10000000*a^7
+       + 33333333/100000000*a^8
+       + 1/3*a^9
+       + 1/3*a^10
+       + 1/3*a^11
+       + 1/3*a^12
+       + 1/3*a^13
+       + 1/3*a^14
+       + 1/3*a^15
+       + 1/3*a^16")
+assert result("F4") =~ expr("+ 1/10000000000000000000")
+*--#] torat : 
 *--#[ strictrounding_1 :
 #StartFloat 6d
 CFunction f;

sources/float.c

@@ -591,7 +591,7 @@ long FloatToInteger(UWORD *out, mpf_t floatin, long *bitsused)
 	mpz_set_f(z,floatin);
 	ZtoForm(out,&nout,z);
 	mpz_clear(z);
-	x = out[0]; nx = 0;
+	x = out[nout-1]; nx = 0;
 	while ( x ) { nx++; x >>= 1; }
 	*bitsused = (nout-1)*BITSINWORD + nx;
 	return(nout);
@@ -678,7 +678,7 @@ int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin)
 	WORD na, nb, nc, nd, i;
 	int nout;
 	LONG oldpWorkPointer = AT.pWorkPointer;
-	long bitsused = 0, totalbitsused = 0, totalbits = AC.DefaultPrecision;
+	long bitsused = 0, totalbitsused = 0, totalbits = AC.DefaultPrecision-AC.MaxWeight-1;
 	int retval = 0, startnul;
 	WantAddPointers(AC.DefaultPrecision);  /* Horrible overkill */
 	AT.pWorkSpace[AT.pWorkPointer++] = out;
@@ -710,15 +710,24 @@ int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin)
 		else {
 			nout = FloatToInteger((UWORD *)out,aux1,&bitsused);
 			out += nout;
-			totalbitsused += bitsused;
 		}
 		if ( bitsused > (totalbits-totalbitsused)/2 ) { break; }
 		if ( mpf_sgn(aux2) == 0 ) {
 			/*if ( startnul == 1 )*/ AT.pWorkSpace[AT.pWorkPointer++] = out;
 			break;
 		}
+		totalbitsused += bitsused;
 		AT.pWorkSpace[AT.pWorkPointer++] = out;
 	  }
+	  /* 
+	  	For |floatin| << 1, we have startnul = 1 and hit the precision guard 
+	  	already on the first continued-fraction step. The resulting float is 
+		therefore close to the rational 0/1 and we can immediately return.
+	   */
+	  if ( startnul == 1 && AT.pWorkPointer == oldpWorkPointer + 2 ) {
+		*ratout++ = 0; *ratout++ = 1; *ratout = *nratout = 3;
+		goto ret;
+	  }
 /*
 	  At this point we have the function with the repeated fraction.
 	  Now we should work out the fraction. Form code would be:
@@ -734,11 +743,6 @@ int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin)
 	  na = 1; a[0] = 1;
 	  c = (UWORD *)(AT.pWorkSpace[--AT.pWorkPointer]);
 	  nc = nb = ((UWORD *)out)-c;
-	  if ( nc > 10 ) {
-		mc = c;
-		c = (UWORD *)(AT.pWorkSpace[--AT.pWorkPointer]);
-		nc = nb = ((UWORD *)mc)-c;
-	  }
 	  for ( i = 0; i < nb; i++ ) b[i] = c[i];
 	  mc = c = NumberMalloc("FloatToRat");
 	  while ( AT.pWorkPointer > oldpWorkPointer ) {
@@ -791,12 +795,13 @@ int FloatToRat(UWORD *ratout, WORD *nratout, mpf_t floatin)
 	if ( s < 0 ) mpf_neg(floatin,floatin);
 /*
 	{
-		WORD n = *ratout;
+		WORD n = *nratout;
 		RatToFloat(aux1,ratout,n);
 		mpf_sub(aux2,floatin,aux1);
 		gmp_printf("Diff is %.*Fe\n",40,aux2);
 	}
 */
+ret:
 	return(retval);
 }
 
@@ -1492,7 +1497,11 @@ int ToRat(PHEAD WORD *term, WORD level)
 		The result can go straight over the float_ function.
 */
 		UnpackFloat(aux4,t);
+		// If aux4 is zero, the term vanishes
+		if ( mpf_sgn(aux4) == 0 ) return(0);
 		if ( FloatToRat((UWORD *)t,&ncoef,aux4) == 0 ) {
+			// Check if the resulting rational is zero
+			if ( t[0] == 0 && t[1] == 1 && ncoef == 3 ) return(0);
 			t += ABS(ncoef);
 			t[-1] = ncoef*nsign;
 			*term = t - term;