Floating-point irregularities in ga68
chris hermansen
clhermansen@gmail.com
Sat May 9 16:23:14 GMT 2026
Nelson and list,
On Sat, May 9, 2026, 08:20 Nelson H. F. Beebe <beebe@math.utah.edu> wrote:
> 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.
Wow Nelson, this is amazing,.thank you so much for this excellent work.
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
>
The problem I see above is you should reference
long_sqrt
as with SUPPER stropping spaces are significant.
For clarity "long 0.5" truly does denote a 64 bit floating point constant
with value 1/2.
"long sqrt(x)" is a syntax error.
"LENG sqrt(x)" would compute the single precision square root of x (if x
were single precision) and then lengthen the result to double precision.
"long_sqrt(x)" calls the double precision square root function on a double
precision argument.
>
> % 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.
>
Of course, "lomg_long_sqrt(too)" gets you to the next level of FP precision
(which is the 80 bit Intel version and not 128 bit quad FP, at least on X64
architecture)
And that is as far as you can go.
I hope this helps.
> ------------------------------------------------------------------------
>
> 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.
>
I think this is a result of the design of the float procedure, which is pre
IEEE FP, and did not include concepts such as signed zero, Inf, NaN,
subnormal numbers and so forth.
I have been giving some thought to making versions of fixed() and float()
that are 754-aware but as I am no kind of expert on this, all I have so far
are thoughts, not code.
I also think it worthwhile to consider some of the more recent and formal
work for converting between internal and external forms, both for
conversion accuracy and speed.
>
> ------------------------------------------------------------------------
>
> 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
>
> ------------------------------------------------------------------------
>
> 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 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.
>
See my comments above about the (archaic) design of float and fixed.
>
> 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.
>
I like that idea.
>
> 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:
>
I can't say for "most programming languages" but I see no reason for GNU
Algol 68 to do a poor job.
>
> 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 -
>
> -------------------------------------------------------------------------------
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://gcc.gnu.org/pipermail/algol68/attachments/20260509/84ad7110/attachment-0001.htm>
More information about the Algol68
mailing list