Floating-point irregularities in ga68
Jose E. Marchesi
jemarch@gnu.org
Sun May 10 09:45:04 GMT 2026
Hello Nelson, thanks for the throughful report.
Please see some comments below.
> I'm the author of a large book on accurate numerical computation of
> elementary and special functions:
>
> The Mathematical-Function Computation Handbook: Programming
> Using the {MathCW} Portable Software Library
> Springer (2017)
> https://doi.org/10.1007/978-3-319-64110-2
> https://www.math.utah.edu/pub/tex/bib/elefunt.html#Beebe:2017:MFC
>
> The book, and its software, addresses several historical
> floating-point designs, as well as IEEE 754 arithmetic, in both binary
> and decimal formats. Since 2006, I have routinely built patched
> versions of gcc releases 4 to 17 to support decimal arithmetic in C,
> and recently, the recent, and independently developed by a single
> expert programmer, kefir C compiler has been extended with decimal
> support.
>
> My book explains why I chose to write the code in highly portable C,
> and it has appendices on interfacing the MathCW C library, -lmcw (a
> drop-in, but much extended, replacement for -lm, and the formatted I/O
> portions of -lc), to other languages: Ada, C++, C#, Fortran, Java, and
> Pascal. Since then, I have developed interfaces and virtual book
> appendices for Cobol, Go, Modula-2, Oberon, and Rust. I believe that
> a similar appendix for Algol 68 is certainly feasible.
>
> When I learned recently of the inclusion of Algol68 support in recent
> gcc-16 and gcc-17 versions, I immediately added that languages to the
> builds, and now have ga68 on CentOS and Ubuntu systems on x86_64, and
> Debian on arm64 (Raspberry Pi + some QEMU-based virtual machines).
>
> I then set about developing sample Algol 68 code for common numerical
> programming, much as I have long done for other many programming
> languages. I have now successfully interfaced Algol 68 code to -lm
> and -lmcw, but quickly found serious issues. Here are four examples:
>
> ------------------------------------------------------------------------
>
> This program compiles and runs successfully with the Genie Algol68
> compiler:
>
> % cat test1.a68
> BEGIN
> LONG REAL x, y;
>
> x := LONG 0.5;
> y := long sqrt(x);
> print(("long sqrt(x) = ", float(y, 22, 18, 0), newline));
>
> SKIP
> END
>
> % a68g test1.a68
> long sqrt(x) = +7.0710678118654752e-1
>
> However, this (I think) equivalent code for ga68 fails to compile:
>
> % cat test-long-sqrt.a68
> begin
> long real x, y;
> string newline = "'n";
> x := long 0.5;
> y := long sqrt(x);
> puts("long sqrt(x) = " + float(y, 22, 18, 0) + newline)
> end
>
> % ga68 -c test-long-sqrt.a68
> test-long-sqrt.a68:5:10: error: construct beginning with %<begin%> in \
> line 1 followed by a serial clause starting in line 2 and then %<:=%> \
> in line 5 and a longsety and a specification and %<;%> and a \
> specification in line 6 and %<end%> in line 7 is an invalid keyword
> 5 | y := long sqrt(x);
> | ^
>
> Changing sqrt() to other Algol 68 functions produces similar
> compilation errors.
>
> Is this perhaps a temporary compiler limitation? As far as I can tell
> from Web searches for code examples, and the Algol 68 specifications,
> "long sqrt(x)", and also "long long sqrt(x)", should be the correct
> way to access higher precision.
As Chris already mentioned, the identifiers for these standard
procedures are long_sqrt, long_long_sqrt, etc. In the modern stropping
regime that is the default in ga68 it is not possible to have
typographical display features (blank characters) in identifiers.
Note that upper stropping is available in ga68 via -fstropping=upper,
should you prefer it.
> ------------------------------------------------------------------------
>
> This test file gets an infinite loop in float():
>
> % cat inf.a68
> { Bug: infinite loop from ga68 inf.a68 && ./a.out }
>
> begin
> real x, y, z;
> x := 1.0;
> y := 0.0;
> z := x / y;
> puts("If Algol 68 arithmetic follows IEEE 754, this should be Infinity: " +
> float(z, 30, 20, 0) + "'n");
> skip
> end
>
> The loop is also seen on an arm64 system.
Our "float" implementation, in the run-time library libga68, comes
directly from the RR. As Chris mentioned, it was designed in a
pre-IEEE754 era, and it may be it misbehaves because of that.
The question is: do we want to modernize these routines to handle
IEEE754 properly?
If the answer is yes (which I would assume) then I will open a bugzilla
for this..
> ------------------------------------------------------------------------
>
> This test file core dumps at run-time:
>
> % cat nan.a68
> { Bug: core dump from ga68 nan.a68 && ./a.out }
>
> begin
> real x, y, z;
> x := 0.0;
> y := 0.0;
> z := x / y;
> puts("If Algol 68 arithmetic follows IEEE 754, this should be NaN: " +
> float(z, 30, 20, 0) + "'n");
> skip
> end
>
> % ga68 -g -g3 nan.a68
>
> % gdb a.out
> (gdb) run
> ...
> Program received signal SIGABRT, Aborted.
> 0x00007ffff7a43387 in raise () from /lib64/libc.so.6
> (gdb) where
> #0 0x00007ffff7a43387 in raise () from /lib64/libc.so.6
> #1 0x00007ffff7a44a78 in abort () from /lib64/libc.so.6
> #2 0x000000000040191b in _libga68_abort (
> fmt=fmt@entry=0x40c5f0 "%s:%u: runtime error: bound %zd out of range [%zd:%zd]\n") at ../../../gcc-17-20260426/libga68/ga68-error.c:48
> #3 0x0000000000401a2e in _libga68_bounds (
> filename=filename@entry=0x40c41a "standard.a68", lineno=lineno@entry=429,
> index=<optimized out>, lower_bound=lower_bound@entry=1,
> upper_bound=upper_bound@entry=16)
> at ../../../gcc-17-20260426/libga68/ga68-error.c:128
> #4 0x0000000000402ebb in STANDARD_digchar (x=<optimized out>)
> at standard.a68:429
> #5 choosedig (y=<optimized out>) at standard.a68:392
> #6 STANDARD_subfixed (v=..., width=width@entry=28, after=after@entry=20)
> at standard.a68:398
> #7 0x0000000000403f2c in STANDARD_fixed (v=..., width=29,
> after=after@entry=20) at standard.a68:200
> #8 0x00000000004071b0 in STANDARD_float (v=..., width=30, after=20, exp=0)
> at standard.a68:268
> #9 0x0000000000401502 in __algol68_main () at nan.a68:9
> #10 0x000000000040122b in main (argc=<optimized out>, argv=<optimized out>)
> at ../../../gcc-17-20260426/libga68/libga68.c:41
>
> When I run the same test on arm64, it completes successsfully, but the
> output is wrong:
>
> % ./a.out
> If Algol 68 arithmetic follows IEEE 754, this should be NaN: +.0000000000000000000e0
The core dump is due to a run-time check failing:
$ ./a.out
standard.a68:429: runtime error: bound -2147483647 out of range [1:16]
Aborted
It happens in:
proc dig_char = (int x) char: "0123456789abcdef"[x+1];
The procedure dig_char is called by `subfixed' and others.
Again, if we aim to support IEEE754 properly, this must be a bugzilla.
> ------------------------------------------------------------------------
>
> Numeric formatting by float() is inaccurate, even when the numbers are
> small and exactly representable:
>
> % cat powtwo032.a68
> begin
> real x, y, z;
> int k;
> string newline = "'n";
>
> x := 2.0;
> z := 1.0;
> for k from 0 to 127
> do
> y := x ^ k;
> puts ("2**(" + fixed(k, 4, 0) + ") = " +
> float(y, 14, 10, 0) + " vs " +
> float(z, 14, 10, 0) + " exact");
> if y /= z then puts(" [DIFFER]") fi;
> puts(newline);
> z +:= z { exact power }
> od
> end
>
> % ga68 powtwo032.a68 && ./a.out
> 2**( +0) = +1.000000000e0 vs +1.000000000e0 exact
> 2**( +1) = +2.000000000e0 vs +2.000000000e0 exact
> 2**( +2) = +4.000000000e0 vs +4.000000000e0 exact
> 2**( +3) = +8.000000000e0 vs +8.000000000e0 exact
> 2**( +4) = +1.599999904e1 vs +1.599999904e1 exact
> 2**( +5) = +3.199999809e1 vs +3.199999809e1 exact
> 2**( +6) = +6.399999618e1 vs +6.399999618e1 exact
> 2**( +7) = +1.280000095e2 vs +1.280000095e2 exact
> 2**( +8) = +2.560000181e2 vs +2.560000181e2 exact
> 2**( +9) = +5.120000362e2 vs +5.120000362e2 exact
> 2**( +10) = +1.024000048e3 vs +1.024000048e3 exact
> ...
>
> Those numbers should be 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024,
Similar situation than above.
> ...
>
> Similar programs for long real and long long real also show incorrect
> output conversions by float(), but happily, there are no output lines
> with "DIFFER" --- thus, x ^ k is computed exactly correct in all three
> floating-point precisions up to just below the overflow limit.
>
> Some years ago, gcc switched its floating-point conversions between
> binary and decimal values to use the MPFR library for correctly
> rounded conversions. ga68 might well do the same.
>
> It is worth noting here that the required digit counts for correct
> round-trip binary <--> decimal conversions have been known since 1968,
> but most programmers are unaware of them, and almost all programming
> languages get those counts wrong. Here is what my hoc (high-order
> calcular) program documents:
>
> hoc> ?matula(nbits):
>
> matula(nbits) is the number of decimal digits needed to ensure
> correct round-trip conversion between binary and decimal of
> floating-point data with nbits bits in the significand.
>
> matula(P) is that number for this implementation of hoc.
>
> In general, for decimal-to-binary conversion of d decimal digits
> to p bits, we need to ensure that
>
> 10**d < 2**(p - 1),
>
> and for binary-to-decimal conversion, we require that
>
> 2**p < 10**(d - 1).
>
> See David W. Matula, ``In-and-out conversions'', Comm. ACM 11(1)
> 47--50, January 1968. CODEN CACMA2. ISSN 0001-0782.
> DOI:10.1145/362851.362887
>
> See also help_goldberg(), help_goldbergb(), and help_matulab().
>
> hoc> ?goldberg
> goldberg(ndecdig):
>
> goldberg(ndecdig) is the number of bits needed to ensure
> correct round-trip conversion between binary and decimal of
> floating-point data with ndecdig decimal digits in the
> significand.
>
> In general, for decimal-to-binary conversion of d decimal digits
> to p bits, we need to ensure that
>
> 10**d < 2**(p - 1),
>
> and for binary-to-decimal conversion, we require that
>
> 2**p < 10**(d - 1).
>
> See I. Bennett Goldberg, ``27 Bits Are Not Enough For 8-Digit
> Accuracy'', Comm. ACM, 10(2) 105--106, February 1967. CODEN CACMA2.
> ISSN 0001-0782. DOI:10.1145/363067.363112
>
> See also help_goldbergb(), help_matula(), and help_matulab().
>
> The two functions are simple to implement:
>
> func matula (b) { return ceil( b * log10(2) + 1) }
> func goldberg (d) { return ceil( d * log2(10) + 1 ) }
>
> and obviously, if the compiler does not already do so, the constant
> logarithms can be replaced by compile-time contants, if they have the
> correct precisions.
>
> In particular, for the IEEE 754 32-, 64-, 80-, and 128-bit binary
> formats, the correct counts from matula() are 9, 17, 21, and 36
> decimal digits.
>
> ----------------------------------------
>
> P.S. I mentioned support for decimal floating-point arithmetic at the
> start of this long message, and I suggest that it is something that
> ga68 should give serious consideration to. Decimal arithmetic was
> described in an IEEE Standard in 1987, and included in IEEE 754-2009;
> it has been provided in hardware in two IBM CPU families since the
> 1990s. Apart from support for three new numeric types (in C,
> _Decimal{32,64,128}), only at the point of code generation for +, -,
> *, and / does a call to a software implementation of those decimal
> operations, or on Power and S/390 (z-Series) [and perhaps soon on
> RISC-V], a hardware instruction, need to be issued. I suspect that
> for many programming language compilers, the work needed should be
> less than a day or two, and the code changes relatively small. I have
> already done the far larger work for the run-time numerical and I/O
> library support of decimal floating-point arithmetic.
>
> -------------------------------------------------------------------------------
> - Nelson H. F. Beebe Tel: +1 801 581 5254 -
> - University of Utah -
> - Department of Mathematics, 110 LCB Internet e-mail: beebe@math.utah.edu -
> - 155 S 1400 E RM 233 beebe@acm.org beebe@computer.org -
> - Salt Lake City, UT 84112-0090, USA URL: https://www.math.utah.edu/~beebe -
> -------------------------------------------------------------------------------
More information about the Algol68
mailing list