Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]


Groups > comp.lang.c > #402594 > unrolled thread

Taking hypot(3) To The Edge

Started byLawrence D’Oliveiro <ldo@nz.invalid>
First post2026-10-01 00:47 +0000
Last post2026-10-02 11:39 +0100
Articles 20 on this page of 27 — 10 participants

Back to article view | Back to comp.lang.c


Contents

  Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-01 00:47 +0000
    Re: Taking hypot(3) To The Edge "Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> - 2026-10-01 01:02 +0000
      Re: Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-01 01:55 +0000
        Re: Taking hypot(3) To The Edge "Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> - 2026-10-01 03:44 +0000
    Re: Taking hypot(3) To The Edge bart <bc@freeuk.com> - 2026-10-01 02:33 +0100
      Re: Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-01 02:38 +0000
      Re: Taking hypot(3) To The Edge David Brown <david.brown@hesbynett.no> - 2026-10-01 10:17 +0200
        Re: Taking hypot(3) To The Edge Michael S <already5chosen@yahoo.com> - 2026-10-01 14:12 +0300
          Re: Taking hypot(3) To The Edge David Brown <david.brown@hesbynett.no> - 2026-10-01 13:47 +0200
            Re: Taking hypot(3) To The Edge Terje Mathisen <terje.mathisen@tmsw.no> - 2026-10-02 20:35 +0200
              Re: Taking hypot(3) To The Edge bart <bc@freeuk.com> - 2026-10-02 21:17 +0100
              Re: Taking hypot(3) To The Edge BGB <cr88192@gmail.com> - 2026-10-02 17:13 -0500
                Re: Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-03 02:17 +0000
                  Re: Taking hypot(3) To The Edge BGB <cr88192@gmail.com> - 2026-10-03 15:20 -0500
                    Re: Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-03 21:46 +0000
                      Re: Taking hypot(3) To The Edge BGB <cr88192@gmail.com> - 2026-10-03 20:35 -0500
                        Re: Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-04 23:54 +0000
              Re: Taking hypot(3) To The Edge Michael S <already5chosen@yahoo.com> - 2026-10-03 21:06 +0300
                Re: Taking hypot(3) To The Edge Terje Mathisen <terje.mathisen@tmsw.no> - 2026-10-03 20:22 +0200
                  Re: Taking hypot(3) To The Edge Michael S <already5chosen@yahoo.com> - 2026-10-03 22:01 +0300
                    Re: Taking hypot(3) To The Edge Terje Mathisen <terje.mathisen@tmsw.no> - 2026-10-03 23:20 +0200
                    Re: Taking hypot(3) To The Edge MitchAlsup <user5857@newsgrouper.org.invalid> - 2026-10-04 17:49 +0000
                  Re: Taking hypot(3) To The Edge Thomas Koenig <tkoenig@netcologne.de> - 2026-10-03 19:47 +0000
    Re: Taking hypot(3) To The Edge bart <bc@freeuk.com> - 2026-10-01 23:05 +0100
      Re: Taking hypot(3) To The Edge Lawrence D’Oliveiro <ldo@nz.invalid> - 2026-10-02 02:34 +0000
        Re: Taking hypot(3) To The Edge Keith Thompson <Keith.S.Thompson+u@gmail.com> - 2026-10-01 21:14 -0700
        Re: Taking hypot(3) To The Edge bart <bc@freeuk.com> - 2026-10-02 11:39 +0100

Page 1 of 2  [1] 2  Next page →


#402594 — Taking hypot(3) To The Edge

FromLawrence D’Oliveiro <ldo@nz.invalid>
Date2026-10-01 00:47 +0000
SubjectTaking hypot(3) To The Edge
Message-ID<119kaj9$npbg$1@dont-email.me>
Some discussion in another thread prompted me to wonder exactly how
close to the limit a function like hypot(3) can be taken and still
give meaningful results. So I wrote a test script which generates a C
program using different real types, runs it, and reports the results:

    #!/usr/bin/python3
    #+
    # Probe the usable limits of the hypot(3) function.
    #-

    import sys
    import os
    import tempfile
    import subprocess
    import shutil

    src_template = \
    """#include <math.h>
    #include <values.h>
    #include <stdio.h>

    int main(void)
      {
        fprintf(stdout, "sizeof(%(realtype)s) = %%d\\n", sizeof(%(realtype)s));
        fprintf(stdout, "%(maxreal)s (%(realtype)s) = %%.%(prec)se\\n", %(maxreal)s);
        const %(realtype)s down1 = %(nextafter)s(%(maxreal)s, -1.0);
        fprintf(stdout, "next towards 0 = %%.%(prec)se\\n", down1);
        fprintf(stdout, "nr sig bits = %%ld\\n", lrint(floor(- %(log2)s((%(maxreal)s - down1) / down1))));
        fprintf(stdout, "nr sig digits = %%ld\\n", lrint(floor(- %(log10)s((%(maxreal)s - down1) / down1))));
        const %(realtype)s factor = %(sqrt)s(2.0);
        const %(realtype)s factor_next = %(nextafter)s(factor, INFINITY);
        const %(realtype)s factor_prev = %(nextafter)s(factor, - INFINITY);
        fprintf(stdout, "hypot(%%.%(prec)se) = %%.%(prec)se\\n", %(maxreal)s / factor, %(hypot)s(%(maxreal)s / factor, %(maxreal)s / factor));
        fprintf(stdout, "hypot(%%.%(prec)se) = %%.%(prec)se\\n", %(maxreal)s / factor_next, %(hypot)s(%(maxreal)s / factor_next, %(maxreal)s / factor_next));
        fprintf(stdout, "hypot(%%.%(prec)se) = %%.%(prec)se\\n", %(maxreal)s / factor_prev, %(hypot)s(%(maxreal)s / factor_prev, %(maxreal)s / factor_prev));
        return
            0;
      } /*main*/
    """

    tempdir = tempfile.mkdtemp(prefix = "hypot_test")
    progsrc = os.path.join(tempdir, "test.c")
    progbin = os.path.join(tempdir, "test")

    for parms in \
        (
            ("float", "FLT_MAX", "7", "nextafterf",
                "log2f", "log10f", "sqrtf", "hypotf"),
            ("double", "DBL_MAX", "16", "nextafter",
                "log2", "log10", "sqrt", "hypot"),
            ("long double", "LDBL_MAX", "20L", "nextafterl",
                "log2l", "log10l", "sqrtl", "hypotl"),
        ) \
    :
        src = open(progsrc, "wt")
        src.write \
          (
                src_template
            %
                dict
                  (
                    zip
                      (
                        ("realtype", "maxreal", "prec", "nextafter",
                            "log2", "log10", "sqrt", "hypot"),
                        parms
                      )
                  )
          )
        src.close()
        subprocess.check_call(args = ("gcc", "-o", progbin, progsrc, "-lm"))
        subprocess.check_call(args = (progbin,))
        sys.stdout.write("\n")
    #end for

    shutil.rmtree(tempdir)

And here’s what I got:

    ldo@theon:python_try> ./c_hypot_test
    sizeof(float) = 4
    FLT_MAX (float) = 3.4028235e+38
    next towards 0 = 3.4028233e+38
    nr sig bits = 24
    nr sig digits = 7
    hypot(2.4061597e+38) = inf
    hypot(2.4061594e+38) = 3.4028233e+38
    hypot(2.4061599e+38) = inf

    sizeof(double) = 8
    DBL_MAX (double) = 1.7976931348623157e+308
    next towards 0 = 1.7976931348623155e+308
    nr sig bits = 53
    nr sig digits = 15
    hypot(1.2711610061536460e+308) = 1.7976931348623155e+308
    hypot(1.2711610061536458e+308) = 1.7976931348623151e+308
    hypot(1.2711610061536462e+308) = 1.7976931348623157e+308

    sizeof(long double) = 16
    LDBL_MAX (long double) = 1.18973149535723176502e+4932
    next towards 0 = 1.18973149535723176496e+4932
    nr sig bits = 64
    nr sig digits = 19
    hypot(8.41267208158310063619e+4931) = inf
    hypot(8.41267208158310063554e+4931) = 1.18973149535723176496e+4932
    hypot(8.41267208158310063683e+4931) = inf

A couple of things stand out immediately: one is that the “long
double” type doesn’t seem to make use of all the 128 bits it occupies.
I was expecting a mantissa length closer to 100 bits, but it’s nowhere
near that.

Another is that the rounding errors lead to some interesting
precision-dependent edge behaviour. double, in particular, seems to be
able to produce a result that is not “inf” a little bit beyond the
limit that you would expect.

[toc] | [next] | [standalone]


#402595

From"Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu>
Date2026-10-01 01:02 +0000
Message-ID<119kbf8$o1jo$1@dont-email.me>
In reply to#402594
On Thu, 1 Oct 2026 00:47:37 -0000 (UTC), Lawrence D’Oliveiro wrote:

> Some discussion in another thread prompted me to wonder exactly how
> close to the limit a function like hypot(3) can be taken and still
> give meaningful results. So I wrote a test script which generates a C
> program using different real types, runs it, and reports the results:
> 
>

python?  This is a newgroup about C, isn't?

>     hypot(2.4061597e+38) = inf
>     hypot(2.4061594e+38) = 3.4028233e+38
>     hypot(2.4061599e+38) = inf

hypot() takes two arguments.  Are we to guess what the 2nd argement is?

> 
>     sizeof(long double) = 16
>     LDBL_MAX (long double) = 1.18973149535723176502e+4932
>     next towards 0 = 1.18973149535723176496e+4932
>     nr sig bits = 64
         ^^^^^^^^^^^^^

>     nr sig digits = 19
>     hypot(8.41267208158310063619e+4931) = inf
>     hypot(8.41267208158310063554e+4931) = 1.18973149535723176496e+4932
>     hypot(8.41267208158310063683e+4931) = inf
> 
> A couple of things stand out immediately: one is that the “long
> double” type doesn’t seem to make use of all the 128 bits it occupies.
> I was expecting a mantissa length closer to 100 bits, but it’s nowhere
> near that.

Why would long double have a significand with nominally 100 bits, when
the precision is 64 bits?  Are you using an Intel/AMD cpu by chance?

PS: Instead of using python and trying to guess model numbers, you
could available yourself of float.h.

#define LDBL_MANT_DIG   64
#define LDBL_EPSILON    1.0842021724855044340E-19L
#define LDBL_DIG        18
#define LDBL_MIN_EXP    (-16381)
#define LDBL_MIN        3.3621031431120935063E-4932L
#define LDBL_MIN_10_EXP (-4931)
#define LDBL_MAX_EXP    16384
#define LDBL_MAX        1.1897314953572317650E+4932L
#define LDBL_MAX_10_EXP 4932
#if __ISO_C_VISIBLE >= 2011
#define LDBL_TRUE_MIN   3.6451995318824746025E-4951L
#define LDBL_DECIMAL_DIG 21
#define LDBL_HAS_SUBNORM 1
#endif /* __ISO_C_VISIBLE >= 2011 */

-- 
steve

[toc] | [prev] | [next] | [standalone]


#402597

FromLawrence D’Oliveiro <ldo@nz.invalid>
Date2026-10-01 01:55 +0000
Message-ID<119keiu$p2vq$1@dont-email.me>
In reply to#402595
On Thu, 1 Oct 2026 01:02:32 -0000 (UTC), Steven G. Kargl wrote:

> Why would long double have a significand with nominally 100 bits ...

It takes up 128 bits. The exponent only seems to be about 15 bits
(±4931-ish in base-10). So what else is using that space?

[toc] | [prev] | [next] | [standalone]


#402599

From"Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu>
Date2026-10-01 03:44 +0000
Message-ID<119kkvr$qnda$1@dont-email.me>
In reply to#402597
On Thu, 1 Oct 2026 01:55:43 -0000 (UTC), Lawrence D’Oliveiro wrote:

> On Thu, 1 Oct 2026 01:02:32 -0000 (UTC), Steven G. Kargl wrote:
> 
>> Why would long double have a significand with nominally 100 bits ...
> 
> It takes up 128 bits. The exponent only seems to be about 15 bits
> (±4931-ish in base-10). So what else is using that space?

Well, if you're using an Intel/AMD cpu, it is mostlikely junk.

https://www.uclibc.org/docs/psABI-x86_64.pdf

See page 11.

-- 
steve

[toc] | [prev] | [next] | [standalone]


#402596

Frombart <bc@freeuk.com>
Date2026-10-01 02:33 +0100
Message-ID<119kd8g$opnp$1@dont-email.me>
In reply to#402594
On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:

> A couple of things stand out immediately: one is that the “long
> double” type doesn’t seem to make use of all the 128 bits it occupies.
> I was expecting a mantissa length closer to 100 bits, but it’s nowhere
> near that.

Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit 
mantissa (with explicit top bit).

10 bytes would be an odd size though, so it's rounded up to 16 bytes in 
the implementation.

I don't know if current hardware would directly support a full 128 bits, 
or it needs to be emulated.

However, float limits (in float.h) should give a clue as to the actual 
range and precision.

[toc] | [prev] | [next] | [standalone]


#402598

FromLawrence D’Oliveiro <ldo@nz.invalid>
Date2026-10-01 02:38 +0000
Message-ID<119kh3m$pplq$1@dont-email.me>
In reply to#402596
On Thu, 1 Oct 2026 02:33:01 +0100, bart wrote:

> Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit
> mantissa (with explicit top bit).

80 bits sounds about right.

> 10 bytes would be an odd size though, so it's rounded up to 16 bytes
> in the implementation.

On the old Apple Mac, this was the “extended” format, which occupied
10 bytes. Motorola’s 68881 and 68882 FPU chips used 12 bytes (96
bits), with two padding bytes just for alignment’s sake.

> I don't know if current hardware would directly support a full 128
> bits, or it needs to be emulated.

I thought some systems had a software-implemented format consisting of
two doubles, to combine the accuracy of their two mantissas somehow.

[toc] | [prev] | [next] | [standalone]


#402606

FromDavid Brown <david.brown@hesbynett.no>
Date2026-10-01 10:17 +0200
Message-ID<119l4vl$t831$2@dont-email.me>
In reply to#402596
On 01/10/2026 03:33, bart wrote:
> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
> 
>> A couple of things stand out immediately: one is that the “long
>> double” type doesn’t seem to make use of all the 128 bits it occupies.
>> I was expecting a mantissa length closer to 100 bits, but it’s nowhere
>> near that.
> 
> Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit 
> mantissa (with explicit top bit).
> 

Yes, that's the default for "long double" in the standard x86-64 ABI.

On 32-bit x86, these 10-byte doubles were often stored in 12-byte 
(96-bit) containers for better alignment, but I don't know the ABI 
standards here.

On 64-bit x86, for better alignment they are stored in 16-byte 
containers.  The rest of the space will be padding.

gcc supports "-mlong-double-64", "-mlong-double-80" and 
"-mlong-double-128" flags.  With 64-bit long doubles, they are 
effectively the same as doubles, and thus faster.  80-bit long doubles 
are the standard, while 128-bit long doubles give quad-precision range 
and precision, but must be handled in software through the gcc-provided 
"libquad" library functions.

You can also use the _Float128 or _Float128x types, which are optional 
types in the C standards (along with 16, 32 and 64-bit versions).  The 
"x" types are for calculations and might be bigger than their name 
applies, while the non-x types are "interchange" types that have the 
same format on all systems.


> 10 bytes would be an odd size though, so it's rounded up to 16 bytes in 
> the implementation.
> 
> I don't know if current hardware would directly support a full 128 bits, 
> or it needs to be emulated.
> 

IBM's POWER processors are, I think, the only cpus with hardware 
quad-precision floating point that are realistic.  RISC-V has the Q 
extensions specified, but I do not believe anyone actually makes RISC-V 
cores with quad floating point.

> However, float limits (in float.h) should give a clue as to the actual 
> range and precision.
> 

Yes.

You can also just write some code in godbolt.org, and see if it calls 
library functions like "__multf3" instead of generating hardware 
instructions, as you play around with the compiler flags.

[toc] | [prev] | [next] | [standalone]


#402613

FromMichael S <already5chosen@yahoo.com>
Date2026-10-01 14:12 +0300
Message-ID<20261001141212.00002f20@yahoo.com>
In reply to#402606
On Thu, 1 Oct 2026 10:17:57 +0200
David Brown <david.brown@hesbynett.no> wrote:

> On 01/10/2026 03:33, bart wrote:
> > On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
> >   
> >> A couple of things stand out immediately: one is that the “long
> >> double” type doesn’t seem to make use of all the 128 bits it
> >> occupies. I was expecting a mantissa length closer to 100 bits,
> >> but it’s nowhere near that.  
> > 
> > Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit 
> > mantissa (with explicit top bit).
> >   
> 
> Yes, that's the default for "long double" in the standard x86-64 ABI.
> 
> On 32-bit x86, these 10-byte doubles were often stored in 12-byte 
> (96-bit) containers for better alignment, but I don't know the ABI 
> standards here.
> 
> On 64-bit x86, for better alignment they are stored in 16-byte 
> containers.  The rest of the space will be padding.
> 
> gcc supports "-mlong-double-64", "-mlong-double-80" and 
> "-mlong-double-128" flags.

The latter appears to be a default on ARM64 Linux.

> With 64-bit long doubles, they are 
> effectively the same as doubles, and thus faster.  80-bit long
> doubles are the standard, while 128-bit long doubles give
> quad-precision range and precision, but must be handled in software
> through the gcc-provided "libquad" library functions.
> 

It is more complicated than that.
libquadmath supplies 128-bit equivalents of of standard C math library
functions and functions for conveersion betwen fp128 and strings.
But basic arithmetic operations are coming from other source. One would
suspect libgcc, but reality is more complicated than just that. It
seems that Jakub Jelinek knows details, but I am not sure that there is
somebody else that knows all of them for all platforms.

> You can also use the _Float128 or _Float128x types, which are
> optional types in the C standards (along with 16, 32 and 64-bit
> versions).  The "x" types are for calculations and might be bigger
> than their name applies, while the non-x types are "interchange"
> types that have the same format on all systems.
> 

Another name supported on many platforms is __float128.

> 
> > 10 bytes would be an odd size though, so it's rounded up to 16
> > bytes in the implementation.
> > 
> > I don't know if current hardware would directly support a full 128
> > bits, or it needs to be emulated.
> >   
> 
> IBM's POWER processors are, I think, the only cpus with hardware 
> quad-precision floating point that are realistic.  RISC-V has the Q 
> extensions specified, but I do not believe anyone actually makes
> RISC-V cores with quad floating point.
> 

Don't forget IBM Z.
BTW, the next generation of IBM Z chips will run ARM64 natively
alongside Z. I wonder if it would mativate IBM to ask ARM Inc to add
128-bit BFP (and may be DFP) as optional extension(s) to Arm ISA.

> > However, float limits (in float.h) should give a clue as to the
> > actual range and precision.
> >   
> 
> Yes.
> 
> You can also just write some code in godbolt.org, and see if it calls 
> library functions like "__multf3" instead of generating hardware 
> instructions, as you play around with the compiler flags.
> 
> 

[toc] | [prev] | [next] | [standalone]


#402616

FromDavid Brown <david.brown@hesbynett.no>
Date2026-10-01 13:47 +0200
Message-ID<119lh8n$14irr$2@dont-email.me>
In reply to#402613
On 01/10/2026 13:12, Michael S wrote:
> On Thu, 1 Oct 2026 10:17:57 +0200
> David Brown <david.brown@hesbynett.no> wrote:
> 
>> On 01/10/2026 03:33, bart wrote:
>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
>>>    
>>>> A couple of things stand out immediately: one is that the “long
>>>> double” type doesn’t seem to make use of all the 128 bits it
>>>> occupies. I was expecting a mantissa length closer to 100 bits,
>>>> but it’s nowhere near that.
>>>
>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit
>>> mantissa (with explicit top bit).
>>>    
>>
>> Yes, that's the default for "long double" in the standard x86-64 ABI.
>>
>> On 32-bit x86, these 10-byte doubles were often stored in 12-byte
>> (96-bit) containers for better alignment, but I don't know the ABI
>> standards here.
>>
>> On 64-bit x86, for better alignment they are stored in 16-byte
>> containers.  The rest of the space will be padding.
>>
>> gcc supports "-mlong-double-64", "-mlong-double-80" and
>> "-mlong-double-128" flags.
> 
> The latter appears to be a default on ARM64 Linux.

That makes sense.  Although software 128-bit floating point is going to 
be very slow compared to hardware 64-bit, the reason you would use "long 
double" is to get more than "double".

> 
>> With 64-bit long doubles, they are
>> effectively the same as doubles, and thus faster.  80-bit long
>> doubles are the standard, while 128-bit long doubles give
>> quad-precision range and precision, but must be handled in software
>> through the gcc-provided "libquad" library functions.
>>
> 
> It is more complicated than that.
> libquadmath supplies 128-bit equivalents of of standard C math library
> functions and functions for conveersion betwen fp128 and strings.
> But basic arithmetic operations are coming from other source. One would
> suspect libgcc, but reality is more complicated than just that. It
> seems that Jakub Jelinek knows details, but I am not sure that there is
> somebody else that knows all of them for all platforms.

For targets without hardware floating point for float and/or double, the 
normal floating arithmetic operations are provided in a "language 
support library" - this roughly covers any library code needed for C 
that can be written without any standard headers.  On small 
microcontrollers, it will also include some integer operations, and 
maybe other bits of code.  This library is used and linked implicitly by 
gcc.  I suppose the basic 128-bit floating point arithmetic operations 
will be there too (on targets that support it) - and not in a specific 
bit of "libquadmath" as I first guessed (and which I also named 
inaccurately).

> 
>> You can also use the _Float128 or _Float128x types, which are
>> optional types in the C standards (along with 16, 32 and 64-bit
>> versions).  The "x" types are for calculations and might be bigger
>> than their name applies, while the non-x types are "interchange"
>> types that have the same format on all systems.
>>
> 
> Another name supported on many platforms is __float128.
> 

Yes.  But names with double underscores like that are primarily for 
implementation use, such as for typedefs in standard headers.  It is 
good to use C standard names where practical.

>>
>>> 10 bytes would be an odd size though, so it's rounded up to 16
>>> bytes in the implementation.
>>>
>>> I don't know if current hardware would directly support a full 128
>>> bits, or it needs to be emulated.
>>>    
>>
>> IBM's POWER processors are, I think, the only cpus with hardware
>> quad-precision floating point that are realistic.  RISC-V has the Q
>> extensions specified, but I do not believe anyone actually makes
>> RISC-V cores with quad floating point.
>>
> 
> Don't forget IBM Z.

Yes, I forgot them - I have absolutely no experience with them.  (I have 
no experience with POWER either, but I have used microcontrollers with 
its smaller cousin PowerPC.)

> BTW, the next generation of IBM Z chips will run ARM64 natively
> alongside Z. I wonder if it would mativate IBM to ask ARM Inc to add
> 128-bit BFP (and may be DFP) as optional extension(s) to Arm ISA.
> 

I wonder how much call there is for fast 128-bit precision floating 
point in practice.

[toc] | [prev] | [next] | [standalone]


#402642

FromTerje Mathisen <terje.mathisen@tmsw.no>
Date2026-10-02 20:35 +0200
Message-ID<119otie$2c89q$1@dont-email.me>
In reply to#402616
David Brown wrote:
> On 01/10/2026 13:12, Michael S wrote:
>> On Thu, 1 Oct 2026 10:17:57 +0200
>> David Brown <david.brown@hesbynett.no> wrote:
>>
>>> On 01/10/2026 03:33, bart wrote:
>>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
>>>>> A couple of things stand out immediately: one is that the “long
>>>>> double” type doesn’t seem to make use of all the 128 bits it
>>>>> occupies. I was expecting a mantissa length closer to 100 bits,
>>>>> but it’s nowhere near that.
>>>>
>>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit
>>>> mantissa (with explicit top bit).
>>>
>>> Yes, that's the default for "long double" in the standard x86-64 ABI.
>>>
>>> On 32-bit x86, these 10-byte doubles were often stored in 12-byte
>>> (96-bit) containers for better alignment, but I don't know the ABI
>>> standards here.
>>>
>>> On 64-bit x86, for better alignment they are stored in 16-byte
>>> containers.  The rest of the space will be padding.
>>>
>>> gcc supports "-mlong-double-64", "-mlong-double-80" and
>>> "-mlong-double-128" flags.
>>
>> The latter appears to be a default on ARM64 Linux.
> 
> That makes sense.  Although software 128-bit floating point is going to 
> be very slow compared to hardware 64-bit, the reason you would use "long 
> double" is to get more than "double".

Significantly slower, yes.

For Mill we figured out that if the HW provided a small number of helper 
functions (easy to do in HW, much harder in SW), then you can do most 
ops in maybe 5x the HW cycle count.

Without that help, you need to manually extract 
sign/exp/mantissa(inserting leading bit unless subnormal), then for 
add/sub you must normalize the smaller number (including sticky bit), do 
the add, normalize again, then merge with exponent and round.

For FMUL you don't pre-normalize (so no extra cycles for subnormal), but 
you have to handle the potential for up to 111 bits of post-normalization.

Michael S is my current benchmark source here, I'd guess FADD128 in less 
than 40 cycles, about the same for FMUL128, while FMAC128 has to be a 
little bit harder with a _very_ wide intermediate result.

Terje

-- 
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

[toc] | [prev] | [next] | [standalone]


#402643

Frombart <bc@freeuk.com>
Date2026-10-02 21:17 +0100
Message-ID<119p3h8$2egtu$1@dont-email.me>
In reply to#402642
On 02/10/2026 19:35, Terje Mathisen wrote:
> David Brown wrote:
>> On 01/10/2026 13:12, Michael S wrote:
>>> On Thu, 1 Oct 2026 10:17:57 +0200
>>> David Brown <david.brown@hesbynett.no> wrote:
>>>
>>>> On 01/10/2026 03:33, bart wrote:
>>>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
>>>>>> A couple of things stand out immediately: one is that the “long
>>>>>> double” type doesn’t seem to make use of all the 128 bits it
>>>>>> occupies. I was expecting a mantissa length closer to 100 bits,
>>>>>> but it’s nowhere near that.
>>>>>
>>>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit
>>>>> mantissa (with explicit top bit).
>>>>
>>>> Yes, that's the default for "long double" in the standard x86-64 ABI.
>>>>
>>>> On 32-bit x86, these 10-byte doubles were often stored in 12-byte
>>>> (96-bit) containers for better alignment, but I don't know the ABI
>>>> standards here.
>>>>
>>>> On 64-bit x86, for better alignment they are stored in 16-byte
>>>> containers.  The rest of the space will be padding.
>>>>
>>>> gcc supports "-mlong-double-64", "-mlong-double-80" and
>>>> "-mlong-double-128" flags.
>>>
>>> The latter appears to be a default on ARM64 Linux.
>>
>> That makes sense.  Although software 128-bit floating point is going 
>> to be very slow compared to hardware 64-bit, the reason you would use 
>> "long double" is to get more than "double".
> 
> Significantly slower, yes.
> 
> For Mill we figured out that if the HW provided a small number of helper 
> functions (easy to do in HW, much harder in SW), then you can do most 
> ops in maybe 5x the HW cycle count.
> 
> Without that help, you need to manually extract sign/exp/ 
> mantissa(inserting leading bit unless subnormal), then for add/sub you 
> must normalize the smaller number (including sticky bit), do the add, 
> normalize again, then merge with exponent and round.
This is if you are more interested in extra precision rather than range.

But if range is more important (this was the main focus of the 
discussion on 'hypot'), then a software-emulated version can probably be 
more efficient.

I may have presented an experimental version in this forum a decade or 
two ago, which IIRC used two 'doubles' to represent an extended FP type, 
128 bits in total.

The first double had the value, kept normalised, say within 1.0 to just 
under 2.0, and the exponent was stored in the second double.

Calculations are done as normal, but with extra adjustments as needed.

The exponent range would be some +/- 10**307 in theory (so allowing 
values up to 10**10**307), but there were probably some practical limits.

[toc] | [prev] | [next] | [standalone]


#402644

FromBGB <cr88192@gmail.com>
Date2026-10-02 17:13 -0500
Message-ID<119paam$2gt26$1@dont-email.me>
In reply to#402642
On 10/2/2026 1:35 PM, Terje Mathisen wrote:
> David Brown wrote:
>> On 01/10/2026 13:12, Michael S wrote:
>>> On Thu, 1 Oct 2026 10:17:57 +0200
>>> David Brown <david.brown@hesbynett.no> wrote:
>>>
>>>> On 01/10/2026 03:33, bart wrote:
>>>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
>>>>>> A couple of things stand out immediately: one is that the “long
>>>>>> double” type doesn’t seem to make use of all the 128 bits it
>>>>>> occupies. I was expecting a mantissa length closer to 100 bits,
>>>>>> but it’s nowhere near that.
>>>>>
>>>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a 64-bit
>>>>> mantissa (with explicit top bit).
>>>>
>>>> Yes, that's the default for "long double" in the standard x86-64 ABI.
>>>>
>>>> On 32-bit x86, these 10-byte doubles were often stored in 12-byte
>>>> (96-bit) containers for better alignment, but I don't know the ABI
>>>> standards here.
>>>>
>>>> On 64-bit x86, for better alignment they are stored in 16-byte
>>>> containers.  The rest of the space will be padding.
>>>>
>>>> gcc supports "-mlong-double-64", "-mlong-double-80" and
>>>> "-mlong-double-128" flags.
>>>
>>> The latter appears to be a default on ARM64 Linux.
>>
>> That makes sense.  Although software 128-bit floating point is going 
>> to be very slow compared to hardware 64-bit, the reason you would use 
>> "long double" is to get more than "double".
> 
> Significantly slower, yes.
> 

For cases where one needs more than "double" the slowness is justified.

And, if the person expects it to be as fast as double, they are at fault 
for casually using "long double".


> For Mill we figured out that if the HW provided a small number of helper 
> functions (easy to do in HW, much harder in SW), then you can do most 
> ops in maybe 5x the HW cycle count.
> 
> Without that help, you need to manually extract sign/exp/ 
> mantissa(inserting leading bit unless subnormal), then for add/sub you 
> must normalize the smaller number (including sticky bit), do the add, 
> normalize again, then merge with exponent and round.
> 
> For FMUL you don't pre-normalize (so no extra cycles for subnormal), but 
> you have to handle the potential for up to 111 bits of post-normalization.
> 
> Michael S is my current benchmark source here, I'd guess FADD128 in less 
> than 40 cycles, about the same for FMUL128, while FMAC128 has to be a 
> little bit harder with a _very_ wide intermediate result.
> 

In my case, was just doing it in software with integer math.
Nothing particularly notable here though.

Did eventually end up deciding to support ISA level ops for 128-bit 
floating point even with them being always traps, because in some cases 
the use of traps can be "less bad" than the use of function calls 
(though, mostly on RV64).

Though, strict adherence to RV64G creates a bit of a problem here for 
"long double", as expecting any particular behavior from the 'Q' 
extension ops is inherently undefined:
   Does it just trap and terminate the program;
   Does it trap and emulate Binary128 using register pairs (my approach);
   Does it actually implement the 'Q' extension (other option).


Though, maybe controversial in that personally I feel that implementing 
a reduced version of 'Q' on top of even-numbered register pairs is a 
better option than implementing the Q extension proper (or, IOW, that in 
this case, that natively 128-bit registers would actually make things 
worse).

Though, as can be noted, I also take the design approach/assumption that 
even pairs can be combined in hardware to perform 128-bit operations as 
needed. So, the use of 64-bit as a base size does not preclude the 
possibility of natively 128-bit operations in hardware (had also taken 
this approach for 128-bit ALU and SIMD).

Though, mostly this is because I expect that either way, Binary64 is 
likely to remain by far dominant over Binary128 even in the presence of 
native support for Binary128 in hardware (because, even if not slow, it 
is still likely to be overkill; much like how cheap support for "double" 
doesn't mean that no one uses "float"...).


> Terje
> 

[toc] | [prev] | [next] | [standalone]


#402646

FromLawrence D’Oliveiro <ldo@nz.invalid>
Date2026-10-03 02:17 +0000
Message-ID<119pok8$2l2ql$1@dont-email.me>
In reply to#402644
On Fri, 2 Oct 2026 17:13:35 -0500, BGB wrote:

> ... much like how cheap support for "double" doesn't mean that no
> one uses "float"...).

Double is pretty much the default precision, though. High-level
languages like Python and JavaScript only support one real-number
precision, and that’s double.

Even the naming of the C run-time routines -- an “f” suffix denotes
single precision, “l” denotes long double, while no suffix defaults to
meaning double.

[toc] | [prev] | [next] | [standalone]


#402660

FromBGB <cr88192@gmail.com>
Date2026-10-03 15:20 -0500
Message-ID<119ro2v$3avqk$1@dont-email.me>
In reply to#402646
On 10/2/2026 9:17 PM, Lawrence D’Oliveiro wrote:
> On Fri, 2 Oct 2026 17:13:35 -0500, BGB wrote:
> 
>> ... much like how cheap support for "double" doesn't mean that no
>> one uses "float"...).
> 
> Double is pretty much the default precision, though. High-level
> languages like Python and JavaScript only support one real-number
> precision, and that’s double.
> 
> Even the naming of the C run-time routines -- an “f” suffix denotes
> single precision, “l” denotes long double, while no suffix defaults to
> meaning double.

Can note that Binary64/double does make more sense as a general purpose 
format.

It is also kinda interesting that people 40 years ago seemingly knew 
which format was best for general purpose use, and the needle hasn't 
really moved since then.


Say, 40 years ago:
   int            :  2 bytes
   char*          :  2 bytes
   float          :  4 bytes
   double         :  8 bytes
   long double    : 10 bytes
Now:
   int            :  4 bytes
   char*          :  8 bytes
   float          :  4 bytes
   double         :  8 bytes
   long double    : 16 bytes

We have since gained "half" / "short float", sometimes...

So, say:
   long double    : 16 bytes, rarely used, overkill
   double         :  8 bytes, general purpose
   float          :  4 bytes, 3D and lots of stuff
   short float    :  2 bytes, graphics, sound, and NNs

Or, smaller still:
   FP8            :  1 byte (S.E4.M3), NNs mostly
   FP8A           :  1 byte (S.E3.M4), audio, rotation quats, *1
   FP8U           :  1 byte (E4.M4), HDR graphics, *2.


*1: 4x FP8A is a usable format for rotation quaternions when precision 
isn't a high priority. Otherwise 4xBinary16 or 4xBinary32 are more 
typical in my uses. Though, 4xFp8A isn't used as an in-register format, 
more for in-memory, with 4xFP16 as the corresponding in-register format.


*2: Image fidelity is comparable to RGB555 in this case, but with a 
larger dynamic range.

Had observed (in TKRA-GL) that one can do HDR 3D rendering by shoving 
this through an LDR RGBA32 path (reusing most of the LDR math on the HDR 
color values), which isn't super accurate, but interestingly gave 
limitations and artifacts very close to those of HDR on the Radeon HD 
4850 card I was using in the early 2010s.

The main obvious feature being that things like wonk with texture 
interpolation, blending, and the restrictions on the use of an alpha 
channel; and that (in the final output) the pixels only really had 
limited accuracy (if fetching the framebuffer image as HALF_FLOAT, low 6 
bits of mantissa were zeroed, etc).

If asking AI: its response is that this card did not have native HDR 
paths in HW and was in-fact shoving HDR through the LDR path with a 
similar strategy (hence the visual similarity in the results).



Can note though that despite "float" being widely used in 3D, had noted 
that it isn't quite good enough for things like CSG clipping for 3D mesh 
building; as it can leave small gaps between vertices.

To some extent this issue can be reduced after the fact by finding 
vertices that are close together, averaging them, and then snapping them 
to some level of grid points (typically something smallish).

[toc] | [prev] | [next] | [standalone]


#402663

FromLawrence D’Oliveiro <ldo@nz.invalid>
Date2026-10-03 21:46 +0000
Message-ID<119rt31$3cm86$2@dont-email.me>
In reply to#402660
On Sat, 3 Oct 2026 15:20:39 -0500, BGB wrote:

> It is also kinda interesting that people 40 years ago seemingly knew
> which format was best for general purpose use, and the needle hasn't
> really moved since then.

You mean the Bell Labs C/Unix folks?

I don’t think their design choice was very popular at the time. It
simplified the design and implementation of the C language, which I
think was their primary consideration. I imagine the FORTRAN folks
threw their hands up in horror.

[toc] | [prev] | [next] | [standalone]


#402667

FromBGB <cr88192@gmail.com>
Date2026-10-03 20:35 -0500
Message-ID<119sahf$3gd4s$1@dont-email.me>
In reply to#402663
On 10/3/2026 4:46 PM, Lawrence D’Oliveiro wrote:
> On Sat, 3 Oct 2026 15:20:39 -0500, BGB wrote:
> 
>> It is also kinda interesting that people 40 years ago seemingly knew
>> which format was best for general purpose use, and the needle hasn't
>> really moved since then.
> 
> You mean the Bell Labs C/Unix folks?
> 


I was more thinking of things around the time of the IEEE-754 standard, 
which defined single/double, which mapped 1:1 with C float/double.

And, double was the default, and also the most useful...
And, they chose double rather than single as the default, despite double 
being gigantic compared with the integer or pointer types in use at the 
time.


Say, CPU types:
   PDP-11: 16 bit address space, no native FPU
   6502: 16-bit address space, no multiply or variable shift, ...
   Z80: Mostly like 6502 in this sense;
   8088: 16-bit, still no FPU, 640K of RAM (very big at the time);
   ...

Late 1980s had M68K for a while (eg: Macintosh).
   And, the Apple II/gs using the 65C816
     Extended 6502 with 24-bit address bus.

...

Knowing of the hardware that was in use at the time, and much of the 
code from the 1990s, it seems almost prescient for them to have known 
that 'double' was likely to win out in the long run as the dominant 
floating point format for general purpose use (rather than, say, 
single/float).

...


> I don’t think their design choice was very popular at the time. It
> simplified the design and implementation of the C language, which I
> think was their primary consideration. I imagine the FORTRAN folks
> threw their hands up in horror.

OK.

[toc] | [prev] | [next] | [standalone]


#402687

FromLawrence D’Oliveiro <ldo@nz.invalid>
Date2026-10-04 23:54 +0000
Message-ID<119uovo$cact$3@dont-email.me>
In reply to#402667
On Sat, 3 Oct 2026 20:35:35 -0500, BGB wrote:

> On 10/3/2026 4:46 PM, Lawrence D’Oliveiro wrote:
>>
>> On Sat, 3 Oct 2026 15:20:39 -0500, BGB wrote:
>>
>>> It is also kinda interesting that people 40 years ago seemingly
>>> knew which format was best for general purpose use, and the needle
>>> hasn't really moved since then.
>>
>> You mean the Bell Labs C/Unix folks?
>
> I was more thinking of things around the time of the IEEE-754
> standard, which defined single/double, which mapped 1:1 with C
> float/double.
>
> And, double was the default, and also the most useful...
> And, they chose double rather than single as the default, despite double
> being gigantic compared with the integer or pointer types in use at the
> time.

You’re thinking of C. Looking at the original IEEE-754 spec from 1985,
it doesn’t define anything resembling a “default” precision. The
closest thing that might sound like that would be section 4.3,
“Rounding Precision”:

    Normally, a result is rounded to the precision of its destination.
    However, some systems deliver results only to double or extended
    destinations. On such a system the user, which may be a high-level
    language compiler, shall be able to specify that a result be
    rounded instead to single precision, though it may be stored in
    the double or extended format with its wider exponent range.
    [FOOTNOTE 4: Control of rounding precision is intended to allow
    systems whose destinations are always double or extended to mimic,
    in the absence of over/underflow, the precisions of systems with
    single and double destinations. An implementation should not
    provide operations that combine double or extended operands to
    produce a single result, nor operations that combine double
    extended operands to produce a double result, with only one
    rounding.] Similarly, a system that delivers results only to
    double extended destinations shall permit the user to specify
    rounding to single or double precision. Note that to meet the
    specifications in 4.1, the result cannot suffer more than one
    rounding error.

[toc] | [prev] | [next] | [standalone]


#402656

FromMichael S <already5chosen@yahoo.com>
Date2026-10-03 21:06 +0300
Message-ID<20261003210613.00000b2c@yahoo.com>
In reply to#402642
On Fri, 2 Oct 2026 20:35:56 +0200
Terje Mathisen <terje.mathisen@tmsw.no> wrote:

> David Brown wrote:
> > On 01/10/2026 13:12, Michael S wrote:  
> >> On Thu, 1 Oct 2026 10:17:57 +0200
> >> David Brown <david.brown@hesbynett.no> wrote:
> >>  
> >>> On 01/10/2026 03:33, bart wrote:  
> >>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:  
> >>>>> A couple of things stand out immediately: one is that the
> >>>>> “long double” type doesn’t seem to make use of all the
> >>>>> 128 bits it occupies. I was expecting a mantissa length closer
> >>>>> to 100 bits, but it’s nowhere near that.  
> >>>>
> >>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a
> >>>> 64-bit mantissa (with explicit top bit).  
> >>>
> >>> Yes, that's the default for "long double" in the standard x86-64
> >>> ABI.
> >>>
> >>> On 32-bit x86, these 10-byte doubles were often stored in 12-byte
> >>> (96-bit) containers for better alignment, but I don't know the ABI
> >>> standards here.
> >>>
> >>> On 64-bit x86, for better alignment they are stored in 16-byte
> >>> containers.  The rest of the space will be padding.
> >>>
> >>> gcc supports "-mlong-double-64", "-mlong-double-80" and
> >>> "-mlong-double-128" flags.  
> >>
> >> The latter appears to be a default on ARM64 Linux.  
> > 
> > That makes sense.  Although software 128-bit floating point is
> > going to be very slow compared to hardware 64-bit, the reason you
> > would use "long double" is to get more than "double".  
> 
> Significantly slower, yes.
> 
> For Mill we figured out that if the HW provided a small number of
> helper functions (easy to do in HW, much harder in SW), then you can
> do most ops in maybe 5x the HW cycle count.
> 
> Without that help, you need to manually extract 
> sign/exp/mantissa(inserting leading bit unless subnormal), then for 
> add/sub you must normalize the smaller number (including sticky bit),
> do the add, normalize again, then merge with exponent and round.
> 
> For FMUL you don't pre-normalize (so no extra cycles for subnormal),
> but you have to handle the potential for up to 111 bits of
> post-normalization.
> 
> Michael S is my current benchmark source here, I'd guess FADD128 in
> less than 40 cycles, about the same for FMUL128, while FMAC128 has to
> be a little bit harder with a _very_ wide intermediate result.
> 
> Terje
> 

Of course, an actual cycle count depends on what you measure, latency
or throughput, on your CPU, on how many corners you are willing to cut
w.r.t. support for rounding modes and for FP exceptions and on the ABI.
Current x86-64 SYSV ABI is quite problematic and rather far from
well-thought. 
Windows currently has no official FP128 ABI at all, the closest to
official is the ABI implemented by gcc under msys2. It is even more
problematic than SysV.
As far as I am concerned, the most problematic part of both ABIs is
program status word that is shared with FP32/64. But registers choice
is also bad.
According to what I hear from Thomas Koenig, in Fortran quite a few
corners can be cut without violating language assumptions. For C it
would be harder.
With most corners cut, i.e. without non-default rounding modes and
without support for Inexact exception, and for throughput rather than
latency Zen3/Linux runs at ~50 clocks per FMUL+FADD. I never tried to
separate between the two. Would think that they are about the same.

Pay attention, that nearly all cost of implementing non-default
rounding mode is because of need to read rounding mode from HW register.
The same applies to implementing Inexact exception - the main cost is
update of HW register.

[toc] | [prev] | [next] | [standalone]


#402657

FromTerje Mathisen <terje.mathisen@tmsw.no>
Date2026-10-03 20:22 +0200
Message-ID<119rh5s$38e29$1@dont-email.me>
In reply to#402656
Michael S wrote:
> On Fri, 2 Oct 2026 20:35:56 +0200
> Terje Mathisen <terje.mathisen@tmsw.no> wrote:
> 
>> David Brown wrote:
>>> On 01/10/2026 13:12, Michael S wrote:
>>>> On Thu, 1 Oct 2026 10:17:57 +0200
>>>> David Brown <david.brown@hesbynett.no> wrote:
>>>>   
>>>>> On 01/10/2026 03:33, bart wrote:
>>>>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:
>>>>>>> A couple of things stand out immediately: one is that the
>>>>>>> “long double” type doesn’t seem to make use of all the
>>>>>>> 128 bits it occupies. I was expecting a mantissa length closer
>>>>>>> to 100 bits, but it’s nowhere near that.
>>>>>>
>>>>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a
>>>>>> 64-bit mantissa (with explicit top bit).
>>>>>
>>>>> Yes, that's the default for "long double" in the standard x86-64
>>>>> ABI.
>>>>>
>>>>> On 32-bit x86, these 10-byte doubles were often stored in 12-byte
>>>>> (96-bit) containers for better alignment, but I don't know the ABI
>>>>> standards here.
>>>>>
>>>>> On 64-bit x86, for better alignment they are stored in 16-byte
>>>>> containers.  The rest of the space will be padding.
>>>>>
>>>>> gcc supports "-mlong-double-64", "-mlong-double-80" and
>>>>> "-mlong-double-128" flags.
>>>>
>>>> The latter appears to be a default on ARM64 Linux.
>>>
>>> That makes sense.  Although software 128-bit floating point is
>>> going to be very slow compared to hardware 64-bit, the reason you
>>> would use "long double" is to get more than "double".
>>
>> Significantly slower, yes.
>>
>> For Mill we figured out that if the HW provided a small number of
>> helper functions (easy to do in HW, much harder in SW), then you can
>> do most ops in maybe 5x the HW cycle count.
>>
>> Without that help, you need to manually extract
>> sign/exp/mantissa(inserting leading bit unless subnormal), then for
>> add/sub you must normalize the smaller number (including sticky bit),
>> do the add, normalize again, then merge with exponent and round.
>>
>> For FMUL you don't pre-normalize (so no extra cycles for subnormal),
>> but you have to handle the potential for up to 111 bits of
>> post-normalization.
>>
>> Michael S is my current benchmark source here, I'd guess FADD128 in
>> less than 40 cycles, about the same for FMUL128, while FMAC128 has to
>> be a little bit harder with a _very_ wide intermediate result.
>>
>> Terje
>>
> 
> Of course, an actual cycle count depends on what you measure, latency
> or throughput, on your CPU, on how many corners you are willing to cut
> w.r.t. support for rounding modes and for FP exceptions and on the ABI.
> Current x86-64 SYSV ABI is quite problematic and rather far from
> well-thought.
> Windows currently has no official FP128 ABI at all, the closest to
> official is the ABI implemented by gcc under msys2. It is even more
> problematic than SysV.
> As far as I am concerned, the most problematic part of both ABIs is
> program status word that is shared with FP32/64. But registers choice
> is also bad.
> According to what I hear from Thomas Koenig, in Fortran quite a few
> corners can be cut without violating language assumptions. For C it
> would be harder.
> With most corners cut, i.e. without non-default rounding modes and
> without support for Inexact exception, and for throughput rather than
> latency Zen3/Linux runs at ~50 clocks per FMUL+FADD. I never tried to
> separate between the two. Would think that they are about the same.
> 
> Pay attention, that nearly all cost of implementing non-default
> rounding mode is because of need to read rounding mode from HW register.
> The same applies to implementing Inexact exception - the main cost is
> update of HW register.

So what you are saying is the my guesstimate was in the right ballpark, 
but that it would be far better for fp128 to be a totally separate 
environment, with all flags and modes maintained in the library instead 
of sharing anything with the HW?

Regarding rounding modes, in my own Mill work a single 64-bit const 
contain all the rounding rules for the four non-truncate alternatives, 
so this didn't really cost much.

Terje


-- 
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

[toc] | [prev] | [next] | [standalone]


#402658

FromMichael S <already5chosen@yahoo.com>
Date2026-10-03 22:01 +0300
Message-ID<20261003220150.00007cba@yahoo.com>
In reply to#402657
On Sat, 3 Oct 2026 20:22:51 +0200
Terje Mathisen <terje.mathisen@tmsw.no> wrote:

> Michael S wrote:
> > On Fri, 2 Oct 2026 20:35:56 +0200
> > Terje Mathisen <terje.mathisen@tmsw.no> wrote:
> >   
> >> David Brown wrote:  
> >>> On 01/10/2026 13:12, Michael S wrote:  
> >>>> On Thu, 1 Oct 2026 10:17:57 +0200
> >>>> David Brown <david.brown@hesbynett.no> wrote:
> >>>>     
> >>>>> On 01/10/2026 03:33, bart wrote:  
> >>>>>> On 01/10/2026 01:47, Lawrence D’Oliveiro wrote:  
> >>>>>>> A couple of things stand out immediately: one is that the
> >>>>>>> “long double” type doesn’t seem to make
> >>>>>>> use of all the 128 bits it occupies. I was expecting a
> >>>>>>> mantissa length closer to 100 bits, but it’s nowhere
> >>>>>>> near that.  
> >>>>>>
> >>>>>> Probably it uses Intel's 80-bit x87 FPU format. That uses a
> >>>>>> 64-bit mantissa (with explicit top bit).  
> >>>>>
> >>>>> Yes, that's the default for "long double" in the standard x86-64
> >>>>> ABI.
> >>>>>
> >>>>> On 32-bit x86, these 10-byte doubles were often stored in
> >>>>> 12-byte (96-bit) containers for better alignment, but I don't
> >>>>> know the ABI standards here.
> >>>>>
> >>>>> On 64-bit x86, for better alignment they are stored in 16-byte
> >>>>> containers.  The rest of the space will be padding.
> >>>>>
> >>>>> gcc supports "-mlong-double-64", "-mlong-double-80" and
> >>>>> "-mlong-double-128" flags.  
> >>>>
> >>>> The latter appears to be a default on ARM64 Linux.  
> >>>
> >>> That makes sense.  Although software 128-bit floating point is
> >>> going to be very slow compared to hardware 64-bit, the reason you
> >>> would use "long double" is to get more than "double".  
> >>
> >> Significantly slower, yes.
> >>
> >> For Mill we figured out that if the HW provided a small number of
> >> helper functions (easy to do in HW, much harder in SW), then you
> >> can do most ops in maybe 5x the HW cycle count.
> >>
> >> Without that help, you need to manually extract
> >> sign/exp/mantissa(inserting leading bit unless subnormal), then for
> >> add/sub you must normalize the smaller number (including sticky
> >> bit), do the add, normalize again, then merge with exponent and
> >> round.
> >>
> >> For FMUL you don't pre-normalize (so no extra cycles for
> >> subnormal), but you have to handle the potential for up to 111
> >> bits of post-normalization.
> >>
> >> Michael S is my current benchmark source here, I'd guess FADD128 in
> >> less than 40 cycles, about the same for FMUL128, while FMAC128 has
> >> to be a little bit harder with a _very_ wide intermediate result.
> >>
> >> Terje
> >>  
> > 
> > Of course, an actual cycle count depends on what you measure,
> > latency or throughput, on your CPU, on how many corners you are
> > willing to cut w.r.t. support for rounding modes and for FP
> > exceptions and on the ABI. Current x86-64 SYSV ABI is quite
> > problematic and rather far from well-thought.
> > Windows currently has no official FP128 ABI at all, the closest to
> > official is the ABI implemented by gcc under msys2. It is even more
> > problematic than SysV.
> > As far as I am concerned, the most problematic part of both ABIs is
> > program status word that is shared with FP32/64. But registers
> > choice is also bad.
> > According to what I hear from Thomas Koenig, in Fortran quite a few
> > corners can be cut without violating language assumptions. For C it
> > would be harder.
> > With most corners cut, i.e. without non-default rounding modes and
> > without support for Inexact exception, and for throughput rather
> > than latency Zen3/Linux runs at ~50 clocks per FMUL+FADD. I never
> > tried to separate between the two. Would think that they are about
> > the same.
> > 
> > Pay attention, that nearly all cost of implementing non-default
> > rounding mode is because of need to read rounding mode from HW
> > register. The same applies to implementing Inexact exception - the
> > main cost is update of HW register.  
> 
> So what you are saying is the my guesstimate was in the right
> ballpark, but that it would be far better for fp128 to be a totally
> separate environment, with all flags and modes maintained in the
> library instead of sharing anything with the HW?
> 

Yes. 
Also it would be better if parameteres/return values to/from
support routines (i.e. fadd128, fmul128 etc) passed either in GPRs or in
memory instead of SIMD registers.

> Regarding rounding modes, in my own Mill work a single 64-bit const 
> contain all the rounding rules for the four non-truncate
> alternatives, so this didn't really cost much.
> 

And what happens when FP32 implemented in HW, but FP64 emoulated in
software as is the most common case in ARM-based MCUs?

> Terje
> 
> 

[toc] | [prev] | [next] | [standalone]


Page 1 of 2  [1] 2  Next page →

Back to top | Article view | comp.lang.c


csiph-web