diff --git a/check/features.frm b/check/features.frm index f85fa2d5..c7bd4488 100644 --- a/check/features.frm +++ b/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; diff --git a/sources/float.c b/sources/float.c index dedaa76a..afb9864f 100644 --- a/sources/float.c +++ b/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;