Inaccurate numbers printed by Format float <number>

[tueda]

Format float <number> prints the output in the floating-point notation up to 100 digits. But for the following example, the precision somehow saturates at N ~ 91:

L F =
16020477561421692952487130706964323462700663188794863242830247036589271925412774344912016394073406219
/
9739282494232553676474274187877549323464250480060655692338461020488417652652504432571071309828391142
;

#ifdef `N'
Format float, `N';
#endif
Format 255;
Print +s;
.end

N     result
80    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008
81    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087
82    1.644934066848226436472415166646025189218949901206798437735558229370007470403200873
83    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738
84    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383
85    1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833
86    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336
87    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362
88    1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628
89    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289
90    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890
91    1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900
92    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289007
93    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890079
94    1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900792
95    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289007920
96    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890079200
97    1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900792002
98    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289007920025
99    1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890079200256
100   1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900792002564
exact 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289006197587053040043190...

I guess both the numerator and denominator would be approximated by N-digit integers. But then I expect that the result should be accurate up to ~ N digits. Indeed, if I manually approximate the expression by 95-digit integers (for N = 95), then the result is accurate up to 95 digits:

* manually approximated (note the last 0s in the digits)
L F = 
16020477561421692952487130706964323462700663188794863242830247036589271925412774344912016394073000000
/
9739282494232553676474274187877549323464250480060655692338461020488417652652504432571071309828300000
;
Format float, 95;
Format 255;
Print +s;
.end

95    1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289006197

The first example must print these digits for N = 95 without the manual approximation.

[vermaseren]

That was a nasty bug. It was in DivLong, one of the important routines. In the case that the numerator (a) and the remainder overlap and na is odd and nremainder == na and the denominator (b) is still needed afterwards, the least significant UWORD of the denominator was overwritten with zero and some precision is lost. This happened in the output routine. The loss in precision is fixed at 32 bits. And of course only when Form decided to use the GMP routines. That happens only in MulLong, GcdLong and GcdLong. Now all routines seem safe.

…

On 22 Mar 2023, at 10:32, Takahiro Ueda ***@***.***> wrote: Format float <number> prints the output in the floating-point notation up to 100 digits. But for the following example, the precision somehow saturates at N ~ 91: L F = 16020477561421692952487130706964323462700663188794863242830247036589271925412774344912016394073406219 / 9739282494232553676474274187877549323464250480060655692338461020488417652652504432571071309828391142 ; #ifdef `N' Format float, `N'; #endif Format 255; Print +s; .end N result 80 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008 81 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087 82 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873 83 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738 84 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383 85 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833 86 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336 87 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362 88 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628 89 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289 90 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890 91 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900 92 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289007 93 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890079 94 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900792 95 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289007920 96 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890079200 97 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900792002 98 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289007920025 99 1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890079200256 100 1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900792002564 exact 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289006197587053040043190... I guess both the numerator and denominator would be approximated by N-digit integers. But then I expect that the result should be accurate up to ~ N digits. Indeed, if I manually approximate the expression by 95-digit integers (for N = 95), then the result is accurate up to 95 digits: * manually approximated (note the last 0s in the digits) L F = 16020477561421692952487130706964323462700663188794863242830247036589271925412774344912016394073000000 / 9739282494232553676474274187877549323464250480060655692338461020488417652652504432571071309828300000 ; Format float, 95; Format 255; Print +s; .end 95 1.6449340668482264364724151666460251892189499012067984377355582293700074704032008738336289006197 The first example must print these digits for N = 95 without the manual approximation. — Reply to this email directly, view it on GitHub <#439>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABJPCES3K5HJJPCPDOOB5LLW5LBLLANCNFSM6AAAAAAWDRW4DE>. You are receiving this because you are subscribed to this thread.

[vermaseren]

By now I have checked all routines that called DivLong with such arguments and the only ones I could find dealt with modulus calculus which for most cases involves only single WORD precision for the denominator and hence that would never use the GMP. It looks like probably no major accidents have happened.

[tueda]

Thanks! Indeed 361e037 fixes this issue and now the example correctly gives 1.644934066...6289006197 for N = 95.