mailing list of musl libc
 help / color / mirror / code / Atom feed
* [musl] Do you recommend using fmt_fp() and
@ 2022-08-18 15:51 ardi
  2022-08-19  3:19 ` Rich Felker
  0 siblings, 1 reply; 9+ messages in thread
From: ardi @ 2022-08-18 15:51 UTC (permalink / raw)
  To: musl

Hi,

I'm looking for a small and robust dtoa-like implementation for quad
floats (IEEE binary128). The need is because I'm using John Hauser's
SoftFloat for IEEE binary128 computing, but I have no easy means for
converting such floats from/to strings (I can use the host
printf/strtold for 80bit extended long doubles, but I'm missing some
significant digits by doing that, and besides, if I ever build in a
host that considers long doubles as regular doubles, I'd lose even
more digits).

I've been considering gdtoa for some days, taking into account its
pros and cons, but I don't like its code size, its dependency on the
FPU flavour behaviour, and that it requires mutexes if it's used in
parallel.

So, I was looking at how musl does this. It appears to be in the
fmt_fp() function in vfprintf.c and in floatscan.c

It looks like I can modify these functions and force them to use the
binary128 type as provided by SoftFloat, instead of using long double.

But it can require quite a bit of surgery, so, before I get my hands
busy in it, I have to ask the question: Would you use this
implementation for my needs if you were me?

Did you adapt fmt_fp() and floatscan from older code? Was that code
ready for 128bit floats?

Or maybe you can recommend another dtoa-like code for 128bit floats?

Thanks a lot,

César

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-08-18 15:51 [musl] Do you recommend using fmt_fp() and ardi
@ 2022-08-19  3:19 ` Rich Felker
  2022-08-23 17:00   ` ardi
  0 siblings, 1 reply; 9+ messages in thread
From: Rich Felker @ 2022-08-19  3:19 UTC (permalink / raw)
  To: ardi; +Cc: musl

On Thu, Aug 18, 2022 at 05:51:26PM +0200, ardi wrote:
> Hi,
> 
> I'm looking for a small and robust dtoa-like implementation for quad
> floats (IEEE binary128). The need is because I'm using John Hauser's
> SoftFloat for IEEE binary128 computing, but I have no easy means for
> converting such floats from/to strings (I can use the host
> printf/strtold for 80bit extended long doubles, but I'm missing some
> significant digits by doing that, and besides, if I ever build in a
> host that considers long doubles as regular doubles, I'd lose even
> more digits).
> 
> I've been considering gdtoa for some days, taking into account its
> pros and cons, but I don't like its code size, its dependency on the
> FPU flavour behaviour, and that it requires mutexes if it's used in
> parallel.
> 
> So, I was looking at how musl does this. It appears to be in the
> fmt_fp() function in vfprintf.c and in floatscan.c
> 
> It looks like I can modify these functions and force them to use the
> binary128 type as provided by SoftFloat, instead of using long double.
> 
> But it can require quite a bit of surgery, so, before I get my hands
> busy in it, I have to ask the question: Would you use this
> implementation for my needs if you were me?
> 
> Did you adapt fmt_fp() and floatscan from older code? Was that code
> ready for 128bit floats?
> 
> Or maybe you can recommend another dtoa-like code for 128bit floats?

I think the fmt_fp code is a very good choice for this. It's basically
the minimal, dependency-free, most straightforward way of doing
exact/correctly-rounded binary floating point to decimal conversion,
and it doess not depend in any way on the format of the long double
type, just knowing the parameters (MANT_DIG, MAX_EXP), and being able
to do a very minimal amount of math on the floating point type to
extract the mantissa and to determine rounding behavior to match the
floating point type's.

If you don't have the equivalent of frexpl you can even do that part
with portable arithmetic on the floating point type. At one point long
ago I think I even had a version of the code that did that, but it
doesn't seem to have ever been in musl proper. It probably predated
musl.

The same should apply to the floatscan.c code if you need the other
direction conversion too. It's just a direct dependency-free version
of the minimal bignum logic needed to do it.

Rich

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-08-19  3:19 ` Rich Felker
@ 2022-08-23 17:00   ` ardi
  2022-08-23 17:30     ` Rich Felker
  0 siblings, 1 reply; 9+ messages in thread
From: ardi @ 2022-08-23 17:00 UTC (permalink / raw)
  To: Rich Felker; +Cc: musl

On Fri, Aug 19, 2022 at 5:19 AM Rich Felker <dalias@libc.org> wrote:
>
> On Thu, Aug 18, 2022 at 05:51:26PM +0200, ardi wrote:
> > Hi,
> >
> > I'm looking for a small and robust dtoa-like implementation for quad
> > floats (IEEE binary128)  [...]
> >
> I think the fmt_fp code is a very good choice for this. It's basically
> the minimal, dependency-free, most straightforward way of doing
> exact/correctly-rounded binary floating point to decimal conversion,
> and it doess not depend in any way on the format of the long double
> type, just knowing the parameters (MANT_DIG, MAX_EXP), and being able
> to do a very minimal amount of math on the floating point type to
> extract the mantissa and to determine rounding behavior to match the
> floating point type's.
>
> If you don't have the equivalent of frexpl you can even do that part
> with portable arithmetic on the floating point type. At one point long
> ago I think I even had a version of the code that did that, but it
> doesn't seem to have ever been in musl proper. It probably predated
> musl.
>
> The same should apply to the floatscan.c code if you need the other
> direction conversion too. It's just a direct dependency-free version
> of the minimal bignum logic needed to do it.

Thanks a lot. I'm working on adapting it at the moment!

One question, though, because I don't know the musl internals, and
I'm defining a custom FILE struct that has only the fields used by
the shgetc.c source file: Can the __uflow() invocation at line 23
of shgetc.c be called for string pseudo-FILEs, or is it guaranteed that
it will be called for real FILE objects only?

I ask the question because __uflow() drives into the stdio internals,
and I'd prefer to avoid that (I'm implementing the fp<>string code
for strings only).

Kind regards and thanks a lot!

César

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-08-23 17:00   ` ardi
@ 2022-08-23 17:30     ` Rich Felker
  2022-08-30 10:17       ` ardi
  0 siblings, 1 reply; 9+ messages in thread
From: Rich Felker @ 2022-08-23 17:30 UTC (permalink / raw)
  To: ardi; +Cc: musl

On Tue, Aug 23, 2022 at 07:00:41PM +0200, ardi wrote:
> On Fri, Aug 19, 2022 at 5:19 AM Rich Felker <dalias@libc.org> wrote:
> >
> > On Thu, Aug 18, 2022 at 05:51:26PM +0200, ardi wrote:
> > > Hi,
> > >
> > > I'm looking for a small and robust dtoa-like implementation for quad
> > > floats (IEEE binary128)  [...]
> > >
> > I think the fmt_fp code is a very good choice for this. It's basically
> > the minimal, dependency-free, most straightforward way of doing
> > exact/correctly-rounded binary floating point to decimal conversion,
> > and it doess not depend in any way on the format of the long double
> > type, just knowing the parameters (MANT_DIG, MAX_EXP), and being able
> > to do a very minimal amount of math on the floating point type to
> > extract the mantissa and to determine rounding behavior to match the
> > floating point type's.
> >
> > If you don't have the equivalent of frexpl you can even do that part
> > with portable arithmetic on the floating point type. At one point long
> > ago I think I even had a version of the code that did that, but it
> > doesn't seem to have ever been in musl proper. It probably predated
> > musl.
> >
> > The same should apply to the floatscan.c code if you need the other
> > direction conversion too. It's just a direct dependency-free version
> > of the minimal bignum logic needed to do it.
> 
> Thanks a lot. I'm working on adapting it at the moment!
> 
> One question, though, because I don't know the musl internals, and
> I'm defining a custom FILE struct that has only the fields used by
> the shgetc.c source file: Can the __uflow() invocation at line 23
> of shgetc.c be called for string pseudo-FILEs, or is it guaranteed that
> it will be called for real FILE objects only?
> 
> I ask the question because __uflow() drives into the stdio internals,
> and I'd prefer to avoid that (I'm implementing the fp<>string code
> for strings only).
> 
> Kind regards and thanks a lot!

See how strtod.c uses sh_fromstring and how sh_fromstring does not
define any io callbacks (so that calling __uflow would necessarily
result in a call to a null or ununitialized function pointer). From
that, it's clear that it must not be reachable (or the existing code
in musl would be broken).

Rich

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-08-23 17:30     ` Rich Felker
@ 2022-08-30 10:17       ` ardi
  2022-08-30 12:26         ` Rich Felker
  0 siblings, 1 reply; 9+ messages in thread
From: ardi @ 2022-08-30 10:17 UTC (permalink / raw)
  To: Rich Felker; +Cc: musl

On Tue, Aug 23, 2022 at 7:30 PM Rich Felker <dalias@libc.org> wrote:
>
> See how strtod.c uses sh_fromstring and how sh_fromstring does not
> define any io callbacks (so that calling __uflow would necessarily
> result in a call to a null or ununitialized function pointer). From
> that, it's clear that it must not be reachable (or the existing code
> in musl would be broken).

Thanks a lot. Understood. While in the process of adapting the code,
I came across the implementation of scalbnl(), and I'm not sure if
I found a mistake or if it's correct: lines 23 to 27 of src/math/scalbnl.c
use MANT_DIG but hardcode it as 113 both when long double is
binary128 and when it's 80bit extended (where MANT_DIG is 64).

By comparison, looking at the source of the scalbn() versions for
float and double, I'd say there's a mistake in the long double version,
but I've never written that function myself, so I'm not in a position
to affirm it for sure.

Kind regards, and thanks a lot,

César

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-08-30 10:17       ` ardi
@ 2022-08-30 12:26         ` Rich Felker
  2022-09-04 19:52           ` ardi
  0 siblings, 1 reply; 9+ messages in thread
From: Rich Felker @ 2022-08-30 12:26 UTC (permalink / raw)
  To: ardi; +Cc: musl

On Tue, Aug 30, 2022 at 12:17:50PM +0200, ardi wrote:
> On Tue, Aug 23, 2022 at 7:30 PM Rich Felker <dalias@libc.org> wrote:
> >
> > See how strtod.c uses sh_fromstring and how sh_fromstring does not
> > define any io callbacks (so that calling __uflow would necessarily
> > result in a call to a null or ununitialized function pointer). From
> > that, it's clear that it must not be reachable (or the existing code
> > in musl would be broken).
> 
> Thanks a lot. Understood. While in the process of adapting the code,
> I came across the implementation of scalbnl(), and I'm not sure if
> I found a mistake or if it's correct: lines 23 to 27 of src/math/scalbnl.c
> use MANT_DIG but hardcode it as 113 both when long double is
> binary128 and when it's 80bit extended (where MANT_DIG is 64).
> 
> By comparison, looking at the source of the scalbn() versions for
> float and double, I'd say there's a mistake in the long double version,
> but I've never written that function myself, so I'm not in a position
> to affirm it for sure.
> 
> Kind regards, and thanks a lot,

I think the scaling just needs to be large enough to ensure we get out
of the subnormal range but not so large as to overflow. Since it's
conditioned on starting out with a negative exponent in the subnormal
range, you'd need an exponent much larger than the max positive
exponent for scaling to overflow, so 113 is safe and the only point of
the number being 113 is that it's sufficient for both formats.

That's just my quick reading of what's going on, so if anyone thinks
I'm wrong about it, please speak up.

Rich

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-08-30 12:26         ` Rich Felker
@ 2022-09-04 19:52           ` ardi
  2022-09-04 21:59             ` Rich Felker
  0 siblings, 1 reply; 9+ messages in thread
From: ardi @ 2022-09-04 19:52 UTC (permalink / raw)
  To: Rich Felker; +Cc: musl

On Tue, Aug 30, 2022 at 2:26 PM Rich Felker <dalias@libc.org> wrote:
>
> I think the scaling just needs to be large enough to ensure we get out
> of the subnormal range but not so large as to overflow. Since it's
> conditioned on starting out with a negative exponent in the subnormal
> range, you'd need an exponent much larger than the max positive
> exponent for scaling to overflow, so 113 is safe and the only point of
> the number being 113 is that it's sufficient for both formats.
>
> That's just my quick reading of what's going on, so if anyone thinks
> I'm wrong about it, please speak up.

I've left as it is, it seems to work fine. Aside from that, and as you
pointed out, I needed a frexp() implementation for SoftFloat, and,
because I want to have all things as complete as possible, I'm writing
frexp() for all the types in SoftFloat (binary16, binary32 aka float,
binary64 aka double, extended 80bit, and binary128). I'm just adapting
the musl implementation so that it accesses the SoftFloat structs/unions.

I'm having a hard time trying to understand the code block for "if(!ee)",
specifically the decrement amount for the exponent. frexpl() decrements
by 120 for both 80bit and binary128 long doubles. It makes sense
because both types have the same exponent size.

Then frexp() decrements by 64. It also seems to make sense because
the exponent is smaller in binary64.

However, frexpf() also decrements by 64. Why? And what's the reason
for choosing 64 and 120? I couldn't find from where they come.

For binary16 I'm replicating all the same logic, but no idea on the
amount for decrementing the exponent on "if(!ee)". If I can understand
how 64 and 120 where chosen for the other types, I'll be able to find it
for binary16.

Thanks again,

César

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-09-04 19:52           ` ardi
@ 2022-09-04 21:59             ` Rich Felker
  2022-09-05  8:49               ` ardi
  0 siblings, 1 reply; 9+ messages in thread
From: Rich Felker @ 2022-09-04 21:59 UTC (permalink / raw)
  To: ardi; +Cc: musl

On Sun, Sep 04, 2022 at 09:52:35PM +0200, ardi wrote:
> On Tue, Aug 30, 2022 at 2:26 PM Rich Felker <dalias@libc.org> wrote:
> >
> > I think the scaling just needs to be large enough to ensure we get out
> > of the subnormal range but not so large as to overflow. Since it's
> > conditioned on starting out with a negative exponent in the subnormal
> > range, you'd need an exponent much larger than the max positive
> > exponent for scaling to overflow, so 113 is safe and the only point of
> > the number being 113 is that it's sufficient for both formats.
> >
> > That's just my quick reading of what's going on, so if anyone thinks
> > I'm wrong about it, please speak up.
> 
> I've left as it is, it seems to work fine. Aside from that, and as you
> pointed out, I needed a frexp() implementation for SoftFloat, and,
> because I want to have all things as complete as possible, I'm writing
> frexp() for all the types in SoftFloat (binary16, binary32 aka float,
> binary64 aka double, extended 80bit, and binary128). I'm just adapting
> the musl implementation so that it accesses the SoftFloat structs/unions.
> 
> I'm having a hard time trying to understand the code block for "if(!ee)",
> specifically the decrement amount for the exponent. frexpl() decrements
> by 120 for both 80bit and binary128 long doubles. It makes sense
> because both types have the same exponent size.

It's not about the exponent size but the mantissa size, and the
numbers 64 and 120 are completely arbitrary except for being larger
than the number of mantissa bits and small enough not to overflow (so
at most a little under twice the max exponent).

Hope this helps.

Rich

^ permalink raw reply	[flat|nested] 9+ messages in thread

* Re: [musl] Do you recommend using fmt_fp() and
  2022-09-04 21:59             ` Rich Felker
@ 2022-09-05  8:49               ` ardi
  0 siblings, 0 replies; 9+ messages in thread
From: ardi @ 2022-09-05  8:49 UTC (permalink / raw)
  To: musl

On Sun, Sep 4, 2022 at 11:59 PM Rich Felker <dalias@libc.org> wrote:
>
> It's not about the exponent size but the mantissa size, and the
> numbers 64 and 120 are completely arbitrary except for being larger
> than the number of mantissa bits and small enough not to overflow (so
> at most a little under twice the max exponent).

Thanks a lot, Rich! After having substituted all the C arithmetic
expressions and libm functions by calls to the SoftFloat functions for
binary128 (and implementing some libm functions not provided by
SoftFloat), I got the FP->String and the String->FP conversion working,
with a 100% software implementation of binary128!!

I tested the round-trip FP->String->FP for a few cases (PI, -PI, 1/3 and
-1/3), using precision 36 (DECIMAL_DIG for binary128), and the
round-trip is successful for binary128!! (yes, I used full-128bit PI from
here: https://www.chusov.org/fp_constants.html )

Of course I need to test this with more numbers, and with a loop of a
zillion random numbers, but for the moment I'm _very_ happy.

I need to tidy it up a bit, and I'd like to implement more libm functions
(ideally I'd wish to provide the full libm catalog). When it's ready, I'll
publish the source in github or similar.

Thanks a lot!!

César

^ permalink raw reply	[flat|nested] 9+ messages in thread

end of thread, other threads:[~2022-09-05  8:49 UTC | newest]

Thread overview: 9+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2022-08-18 15:51 [musl] Do you recommend using fmt_fp() and ardi
2022-08-19  3:19 ` Rich Felker
2022-08-23 17:00   ` ardi
2022-08-23 17:30     ` Rich Felker
2022-08-30 10:17       ` ardi
2022-08-30 12:26         ` Rich Felker
2022-09-04 19:52           ` ardi
2022-09-04 21:59             ` Rich Felker
2022-09-05  8:49               ` ardi

Code repositories for project(s) associated with this public inbox

	https://git.vuxu.org/mirror/musl/

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).