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

Groups > comp.lang.forth > #135651 > unrolled thread

ddfloat library: the behavior of >dd

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

Contents

  ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-10 19:33 +0200
    Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-10 13:35 -0500
      Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-10 22:10 +0200
        Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-10 16:51 -0500
    Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-11 12:16 +1000
      Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-11 12:05 +0200
        Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-12 00:27 +0200
          Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-11 20:29 -0500
            Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-12 12:02 +1000
              Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-12 10:24 -0500
            Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-12 10:45 +0200
              Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-12 19:51 +1000
              Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-12 22:47 +1000
                Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-13 03:14 +1000
                  Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-12 13:02 -0500
                    Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-13 15:44 +1000
                  Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-13 09:00 +0200
                    Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-13 20:33 +1000
                    Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-14 13:28 +1000
                      Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-14 10:09 +0200
                        Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-14 20:00 +1000
                          Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-14 12:22 +0200
                            Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-14 21:59 +1000
                              Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-14 14:26 +0200
                                Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-15 01:32 +1000
                                  Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-14 21:08 +0200
                                    Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-14 15:25 -0500
                                      Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-15 13:11 +1000
                                        Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-15 14:06 +0200
                                        Re: ddfloat library: the behavior of >dd dxf <dxforth@gmail.com> - 2026-09-16 00:05 +1000
                            Re: ddfloat library: the behavior of >dd anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2026-09-14 15:00 +0000
                              Re: ddfloat library: the behavior of >dd marcel hendrix <mhx@iae.nl> - 2026-09-14 21:15 +0200
                              Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-14 15:29 -0500
                                Re: ddfloat library: the behavior of >dd anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2026-09-14 20:58 +0000
                                  Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-15 07:33 -0500
                                    Re: ddfloat library: the behavior of >dd anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2026-09-15 15:53 +0000
                                      Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-15 15:23 -0500
                                        Re: ddfloat library: the behavior of >dd Krishna Myneni <krishna.myneni@ccreweb.org> - 2026-09-15 17:53 -0500

Page 2 of 2 — ← Prev page 1 [2]

#135696

On 14/09/2026 6:09 pm, marcel hendrix wrote:
> On 9/14/2026 5:28 AM, dxf wrote:
>> On 13/09/2026 5:00 pm, marcel hendrix wrote:
>>> ...
>>> My interest is only in >DD . I now allow the exponent range [-294,+304] (this gives a reasonable mantissa size, e.g.
>>> FORTH> s" 1000.0e+304  " >dd ddfe.  1.000000000000000000000000000000e307
>>
>> How did you get that?  NORMALIZE (even if it doesn't hang) will cause
>> 1e307 to be converted to NAN NAN and extracting digits from that will
>> result in a string of '0's.
> 
> Slightly dirty...: normalize ( F: x xx -- |y + yy| ) ( n -- n' ) ddabs                   ( F: -- |x + xx| )
>         DUP  dd=10 dd^n  dd/    ( F: -- [|x+xx|]/10^n ) ( -- n )
>         BEGIN   FOVER 10e F>=    \ make sure it's between 1 and 10
>         WHILE   dd=10 dd/  1+
>         REPEAT
>         BEGIN   FOVER 1e F<    \ make sure it's between 1 and 10
>         WHILE   dd=10 dd*  1-
>         REPEAT ;

That was enough to cause mine to go into a loop (see below).  It was
essentially the change I made in 2012 unaware I was introducing an issue.
I'm using x87 with software stack.  Ran it on Gforth (SSE?) and that goes
into a loop too.


dd# 1e300 dd# 1e7 dd* ok  1.00000000000002E307 1.39689402397438E290 <f
ddfs.  \ hangs!


: normalize ( F: x xx -- |y yy| ) ( n -- n' )
    >R
    DDABS
    dd10 R@ DD^N  DD/          ( [|x|+|xx|]/10^n )

  \ make sure it's between 1 and 10
    BEGIN   FOVER   10e0  F< 0=  \ F>
    WHILE   dd10  dd/  R> 1+ >R  REPEAT

    BEGIN   FOVER   1e0   F<
    WHILE   dd10  dd*  R> 1- >R  REPEAT
    R>
;

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

#135697

On 9/14/2026 12:00 PM, dxf wrote:
> On 14/09/2026 6:09 pm, marcel hendrix wrote:
>> On 9/14/2026 5:28 AM, dxf wrote:
>>> On 13/09/2026 5:00 pm, marcel hendrix wrote:
>>>> ...
>>>> My interest is only in >DD . I now allow the exponent range [-294,+304] (this gives a reasonable mantissa size, e.g.
>>>> FORTH> s" 1000.0e+304  " >dd ddfe.  1.000000000000000000000000000000e307
>>>
>>> How did you get that?  NORMALIZE (even if it doesn't hang) will cause
>>> 1e307 to be converted to NAN NAN and extracting digits from that will
>>> result in a string of '0's.
>>
>> Slightly dirty...: normalize ( F: x xx -- |y + yy| ) ( n -- n' ) ddabs                   ( F: -- |x + xx| )
>>          DUP  dd=10 dd^n  dd/    ( F: -- [|x+xx|]/10^n ) ( -- n )
>>          BEGIN   FOVER 10e F>=    \ make sure it's between 1 and 10
>>          WHILE   dd=10 dd/  1+
>>          REPEAT
>>          BEGIN   FOVER 1e F<    \ make sure it's between 1 and 10
>>          WHILE   dd=10 dd*  1-
>>          REPEAT ;
> 
> That was enough to cause mine to go into a loop (see below).  It was
> essentially the change I made in 2012 unaware I was introducing an issue.
> I'm using x87 with software stack.  Ran it on Gforth (SSE?) and that goes
> into a loop too.
> 
> 
> dd# 1e300 dd# 1e7 dd* ok  1.00000000000002E307 1.39689402397438E290 <f
> ddfs.  \ hangs!
> 
> 
> : normalize ( F: x xx -- |y yy| ) ( n -- n' )
>      >R
>      DDABS
>      dd10 R@ DD^N  DD/          ( [|x|+|xx|]/10^n )
> 
>    \ make sure it's between 1 and 10
>      BEGIN   FOVER   10e0  F< 0=  \ F>
>      WHILE   dd10  dd/  R> 1+ >R  REPEAT
> 
>      BEGIN   FOVER   1e0   F<
>      WHILE   dd10  dd*  R> 1- >R  REPEAT
>      R>
> ;

Interesting. The way it is written, iForth will execute `FOVER   10e0 
F< 0=` and `FOVER 1e0 F<` on the FPU stack with the full 80 bits. (the 
CPU is in 53-bit mode here, which affects the loading of 10e0 as a 10 
byte constant).

I made a few tiny changes to the basic dd-fp words that gave me a
few more valid bits (as seen with the library tests), maybe that is the 
reason.

Try insert .S before 10e0 and 1e0 to see what's happening?

-marcel

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

#135698

On 14/09/2026 8:22 pm, marcel hendrix wrote:
> On 9/14/2026 12:00 PM, dxf wrote:
>> On 14/09/2026 6:09 pm, marcel hendrix wrote:
>>> On 9/14/2026 5:28 AM, dxf wrote:
>>>> On 13/09/2026 5:00 pm, marcel hendrix wrote:
>>>>> ...
>>>>> My interest is only in >DD . I now allow the exponent range [-294,+304] (this gives a reasonable mantissa size, e.g.
>>>>> FORTH> s" 1000.0e+304  " >dd ddfe.  1.000000000000000000000000000000e307
>>>>
>>>> How did you get that?  NORMALIZE (even if it doesn't hang) will cause
>>>> 1e307 to be converted to NAN NAN and extracting digits from that will
>>>> result in a string of '0's.
>>>
>>> Slightly dirty...: normalize ( F: x xx -- |y + yy| ) ( n -- n' ) ddabs                   ( F: -- |x + xx| )
>>>          DUP  dd=10 dd^n  dd/    ( F: -- [|x+xx|]/10^n ) ( -- n )
>>>          BEGIN   FOVER 10e F>=    \ make sure it's between 1 and 10
>>>          WHILE   dd=10 dd/  1+
>>>          REPEAT
>>>          BEGIN   FOVER 1e F<    \ make sure it's between 1 and 10
>>>          WHILE   dd=10 dd*  1-
>>>          REPEAT ;
>>
>> That was enough to cause mine to go into a loop (see below).  It was
>> essentially the change I made in 2012 unaware I was introducing an issue.
>> I'm using x87 with software stack.  Ran it on Gforth (SSE?) and that goes
>> into a loop too.
>>
>>
>> dd# 1e300 dd# 1e7 dd* ok  1.00000000000002E307 1.39689402397438E290 <f
>> ddfs.  \ hangs!
>>
>>
>> : normalize ( F: x xx -- |y yy| ) ( n -- n' )
>>      >R
>>      DDABS
>>      dd10 R@ DD^N  DD/          ( [|x|+|xx|]/10^n )
>>
>>    \ make sure it's between 1 and 10
>>      BEGIN   FOVER   10e0  F< 0=  \ F>
>>      WHILE   dd10  dd/  R> 1+ >R  REPEAT
>>
>>      BEGIN   FOVER   1e0   F<
>>      WHILE   dd10  dd*  R> 1- >R  REPEAT
>>      R>
>> ;
> 
> Interesting. The way it is written, iForth will execute `FOVER   10e0 F< 0=` and `FOVER 1e0 F<` on the FPU stack with the full 80 bits. (the CPU is in 53-bit mode here, which affects the loading of 10e0 as a 10 byte constant).
> 
> I made a few tiny changes to the basic dd-fp words that gave me a
> few more valid bits (as seen with the library tests), maybe that is the reason.
> 
> Try insert .S before 10e0 and 1e0 to see what's happening?

The looping was caused by  NAN 10E F< 0=  ... which always returned TRUE.
Changing to an integrated  F>=  cured the looping.  The NAN came from:

( 1dd307)  dd10 307 DD^N  DD/  ( NAN NAN)

Presumably you don't get this (on the FPU stack)?

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

#135699

On 9/14/2026 1:59 PM, dxf wrote:
> On 14/09/2026 8:22 pm, marcel hendrix wrote:
[..]> The looping was caused by  NAN 10E F< 0=  ... which always 
returned TRUE.
> Changing to an integrated  F>=  cured the looping.  The NAN came from:
 > ( 1dd307)  dd10 307 DD^N  DD/  ( NAN NAN)
 >
 > Presumably you don't get this (on the FPU stack)?

A NAN can't be compared to anything. I have (ISNAN) and (NOTFIN) to 
detect these patterns. My (DDFE.) is
\ create memory string for double-double in E-format
: (DDFE.) ( F: x xx -- ) ( -- c-addr u )
	FOVER F0=      IF  DDDROP S"  0e"   EXIT  ENDIF
	FDUP  (ISNAN)  IF  DDDROP S"  NaN2" EXIT  ENDIF
	FOVER (ISNAN)  IF  DDDROP S"  NaN1" EXIT  ENDIF
	FDUP  (NOTFIN) IF  DDDROP S"  Inf"  EXIT  ENDIF
	FOVER (NOTFIN) IF  DDDROP S"  Inf"  EXIT  ENDIF
	getSign                   ( -- s )
	getPower   ( F: x xx -- ) ( -- s n  )
	normalize  ( F: y yy -- ) ( -- s n' )
	(ddout)	F2DROP ;

At the moment all but one of the bugs I mentioned before are fixed.
It is not possible to do S" 1" >dd because my >dd wants to see at least
S" 1." or S" 1e". I don't feel like fixing that (yet).

There is a weird problem in MATLAB they're not willing to admit is a 
bug. I work around it as it only affects the library accuracy test words 
which are discarded anyway.

-marcel

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

#135701

On 14/09/2026 10:26 pm, marcel hendrix wrote:
> On 9/14/2026 1:59 PM, dxf wrote:
>> On 14/09/2026 8:22 pm, marcel hendrix wrote:
> [..]> The looping was caused by  NAN 10E F< 0=  ... which always returned TRUE.
>> Changing to an integrated  F>=  cured the looping.  The NAN came from:
>> ( 1dd307)  dd10 307 DD^N  DD/  ( NAN NAN)
>>
>> Presumably you don't get this (on the FPU stack)?
> 
> A NAN can't be compared to anything. I have (ISNAN) and (NOTFIN) to detect these patterns. My (DDFE.) is
> \ create memory string for double-double in E-format
> : (DDFE.) ( F: x xx -- ) ( -- c-addr u )
>     FOVER F0=      IF  DDDROP S"  0e"   EXIT  ENDIF
>     FDUP  (ISNAN)  IF  DDDROP S"  NaN2" EXIT  ENDIF
>     FOVER (ISNAN)  IF  DDDROP S"  NaN1" EXIT  ENDIF
>     FDUP  (NOTFIN) IF  DDDROP S"  Inf"  EXIT  ENDIF
>     FOVER (NOTFIN) IF  DDDROP S"  Inf"  EXIT  ENDIF
>     getSign                   ( -- s )
>     getPower   ( F: x xx -- ) ( -- s n  )
>     normalize  ( F: y yy -- ) ( -- s n' )
>     (ddout)    F2DROP ;
> ...

Those won't handle non-reals generated by NORMALIZE.  For DDFS. I've narrowed
the issue down to:

  ( 1dd307 1dd307 ) DD/

It should produce 1dd0 but fails giving NAN NAN instead.

> At the moment all but one of the bugs I mentioned before are fixed.
> It is not possible to do S" 1" >dd because my >dd wants to see at least
> S" 1." or S" 1e". I don't feel like fixing that (yet).

Since my >DD can't generate 1dd307, I looked into that too:

  ( 1dd0 )  307 0 ?DO  dd10 DD*  LOOP

will fail, producing NAN NAN instead.

So DD* and DD/ on my system (also Win32Forth, Gforth) are limited to 1E+-300.
For various reasons unrelated to the code you seem to be getting somewhat
better results.

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

#135702

On 9/14/2026 5:32 PM, dxf wrote:
> Those won't handle non-reals generated by NORMALIZE.  For DDFS. I've narrowed
> the issue down to:
> 
>    ( 1dd307 1dd307 ) DD/
> 
> It should produce 1dd0 but fails giving NAN NAN instead.
> 

Well, I can't use >DD with 1e307:

FORTH> S" 1e307" >dd dddup dd/ ddfe.
 >DD :: exponent too large or too small

> Since my >DD can't generate 1dd307, I looked into that too:
> 
>    ( 1dd0 )  307 0 ?DO  dd10 DD*  LOOP
> 
> will fail, producing NAN NAN instead.
> 
Using your workaround:

FORTH> S" 1e303" >dd S" 1e4" >dd dd* dddup dd/ ddfe. 
1.000000000000000000000000000000e0  ok

So DD/ is not a problem here.

Also
FORTH> S" 100000e303" >dd dddup dd/ ddfe. 
1.000000000000000000000000000000e0  ok

But 1e309 of course generates a NAN ( which is caught by DDFE. and shown 
as "NaN2" )
FORTH> S" 1e303" >dd S" 1e6" >dd dd* dddup dd/ ddfe.  NaN2  ok

> So DD* and DD/ on my system (also Win32Forth, Gforth) are limited to 1E+-300.
> For various reasons unrelated to the code you seem to be getting somewhat
> better results.

The plain limit here is currently
FORTH> S" 1e304" >dd
 >DD :: exponent too large or too small

That is because I want to be able to write
FORTH> S" 10000e303" >dd ddfe. 1.000000000000000000000000000000e307

IMO, for all practical purposes this is more than sufficient.

-marcel

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

#135704

On 9/14/26 2:08 PM, marcel hendrix wrote:
> On 9/14/2026 5:32 PM, dxf wrote:
...
> The plain limit here is currently
> FORTH> S" 1e304" >dd
>  >DD :: exponent too large or too small
> 
> That is because I want to be able to write
> FORTH> S" 10000e303" >dd ddfe. 1.000000000000000000000000000000e307
> 
> IMO, for all practical purposes this is more than sufficient.

Yet, one should check whether or not it's technically correct for 
double-double arithmetic.

We need to go back to Bailey's ddfun package to see how he did the >DD 
conversion.

My two cents.

--
KM

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

#135707

On 15/09/2026 6:25 am, Krishna Myneni wrote:
> On 9/14/26 2:08 PM, marcel hendrix wrote:
>> On 9/14/2026 5:32 PM, dxf wrote:
> ...
>> The plain limit here is currently
>> FORTH> S" 1e304" >dd
>>  >DD :: exponent too large or too small
>>
>> That is because I want to be able to write
>> FORTH> S" 10000e303" >dd ddfe. 1.000000000000000000000000000000e307
>>
>> IMO, for all practical purposes this is more than sufficient.
> 
> Yet, one should check whether or not it's technically correct for double-double arithmetic.
> 
> We need to go back to Bailey's ddfun package to see how he did the >DD conversion.
> 
> My two cents.

On regular double-precision systems the limitation will be DD* DD/ which
makes 10^307 math impractical irrespective of >DD .  I suspect Marcel's
80-bit intermediates/variables is giving him a bit more range.

FWIW my latest version is here:

https://pastebin.com/bcbFUZF0

Numbers out of range will be printed as **DD** .  Bug reports welcome.

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

#135708

On 9/15/2026 5:11 AM, dxf wrote:
> On 15/09/2026 6:25 am, Krishna Myneni wrote:
>> On 9/14/26 2:08 PM, marcel hendrix wrote:
>>> On 9/14/2026 5:32 PM, dxf wrote:
>> ...
>>> The plain limit here is currently
>>> FORTH> S" 1e304" >dd
>>>   >DD :: exponent too large or too small
>>>
>>> That is because I want to be able to write
>>> FORTH> S" 10000e303" >dd ddfe. 1.000000000000000000000000000000e307
>>>
>>> IMO, for all practical purposes this is more than sufficient.
>>
>> Yet, one should check whether or not it's technically correct for double-double arithmetic.
>>
>> We need to go back to Bailey's ddfun package to see how he did the >DD conversion.
>>
>> My two cents.
> 
> On regular double-precision systems the limitation will be DD* DD/ which
> makes 10^307 math impractical irrespective of >DD .  I suspect Marcel's
> 80-bit intermediates/variables is giving him a bit more range.
> 
> FWIW my latest version is here:
> 
> https://pastebin.com/bcbFUZF0
> 
> Numbers out of range will be printed as **DD** .  Bug reports welcome.

General remarks:
D. Bailey has authored a quad-double package with 211 significant bits.
Did anybody look into this yet?

Double-double might be doable with 64bit (extended) floats. That would 
make guarding the FPU control word around exceptions much easier.

-marcel

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

#135710

On 15/09/2026 1:11 pm, dxf wrote:
> On 15/09/2026 6:25 am, Krishna Myneni wrote:
>> On 9/14/26 2:08 PM, marcel hendrix wrote:
>>> On 9/14/2026 5:32 PM, dxf wrote:
>> ...
>>> The plain limit here is currently
>>> FORTH> S" 1e304" >dd
>>>  >DD :: exponent too large or too small
>>>
>>> That is because I want to be able to write
>>> FORTH> S" 10000e303" >dd ddfe. 1.000000000000000000000000000000e307
>>>
>>> IMO, for all practical purposes this is more than sufficient.
>>
>> Yet, one should check whether or not it's technically correct for double-double arithmetic.
>>
>> We need to go back to Bailey's ddfun package to see how he did the >DD conversion.
>>
>> My two cents.
> 
> On regular double-precision systems the limitation will be DD* DD/ which
> makes 10^307 math impractical irrespective of >DD .  I suspect Marcel's
> 80-bit intermediates/variables is giving him a bit more range.
> 
> FWIW my latest version is here:
> 
> https://pastebin.com/bcbFUZF0
> 
> Numbers out of range will be printed as **DD** .  Bug reports welcome.

DDFLOAT 2026-09-16

Bugfix: could print junk instead of **DD**

https://pastebin.com/HqPtRMDH

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

#135700

marcel hendrix <mhx@iae.nl> writes:
>On 9/14/2026 12:00 PM, dxf wrote:
>> : normalize ( F: x xx -- |y yy| ) ( n -- n' )
>>      >R
>>      DDABS
>>      dd10 R@ DD^N  DD/          ( [|x|+|xx|]/10^n )
>> 
>>    \ make sure it's between 1 and 10
>>      BEGIN   FOVER   10e0  F< 0=  \ F>
>>      WHILE   dd10  dd/  R> 1+ >R  REPEAT
>> 
>>      BEGIN   FOVER   1e0   F<
>>      WHILE   dd10  dd*  R> 1- >R  REPEAT
>>      R>
>> ;
>
>Interesting. The way it is written, iForth will execute `FOVER   10e0 
>F< 0=` and `FOVER 1e0 F<` on the FPU stack with the full 80 bits. (the 
>CPU is in 53-bit mode here, which affects the loading of 10e0 as a 10 
>byte constant).

10e can be represented exactly in FP (any format down to 2 mantissa
bits and 3 exponent bits, e.g., the 1.3.2 6-bit format
<https://en.wikipedia.org/wiki/Minifloat#6-bit_(1.3.2)>), certainly in
floats of any Forth system.  1e can also be represented exactly in FP
(0 mantissa bits needed).

The repeated muliplication and division accumulates rounding errors.
For more precision, you could do a binary search among powers of 10
(with the closest DD-FP number where no exact representation is
available), using either a table or a DD** that is precise to 1/2 ulp.
Once you have found the power-of-10 range, DD/ the number by the lower
bound of that range.  If it's essential that the result <10e, then
insert a correction step that checks if the result >=10, and if so,
replaces that result with 1e and a 1 higher exponent.

I am debating with myself whether one should use the closest FP number
with round-to-nearest, or the closest with round-to-absolute-higher.
Maybe depends on the purpose of this normalization.

- anton
-- 
M. Anton Ertl  http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
     New standard: https://forth-standard.org/
EuroForth 2026 CFP: http://www.euroforth.org/ef26/cfp.html

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

#135703

On 9/14/2026 5:00 PM, Anton Ertl wrote:
> I am debating with myself whether one should use the closest FP number
> with round-to-nearest, or the closest with round-to-absolute-higher.
> Maybe depends on the purpose of this normalization.

Yes, that is an interesting question that might, e.g., pop up for the 
next generation metacompilers (generating 128 bits code on a 64 bit 
Forth that can't generate itself yet). I don't know if I will live to 
see it. Maybe when quantum-state memory is invented next month by our AI 
overlords.

-marcel

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

#135705

On 9/14/26 10:00 AM, Anton Ertl wrote:
...
> 
> I am debating with myself whether one should use the closest FP number
> with round-to-nearest, or the closest with round-to-absolute-higher.
> Maybe depends on the purpose of this normalization.
> 
Is your idea to stay within 1 ulp error for the conversion, but use a 
simpler algorithm?

--
KM

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

#135706

Krishna Myneni <krishna.myneni@ccreweb.org> writes:
>On 9/14/26 10:00 AM, Anton Ertl wrote:
>...
>> 
>> I am debating with myself whether one should use the closest FP number
>> with round-to-nearest, or the closest with round-to-absolute-higher.
>> Maybe depends on the purpose of this normalization.
>> 
>Is your idea to stay within 1 ulp error for the conversion, but use a 
>simpler algorithm?

No.

The idea here is as follows:

You have your number X that you want to normalize to a power of 10.
For that number there exist numbers Y, Y', Z, and Z', where

Y > X >= Z
Y >= 10^(n+1) > prev(Y)
Z >= 10^n     > prev(Z)

Where prev(V) produces the largest (in absolute terms) FP value <V.
So Y or Z are either powers of 10 or the smallest FP numbers just
above a power of 10, and they are those numbers which are the
normalization boundaries of X.

For normalizing X, we would love to divide by 10^n, but if that's not
representable as FP number, and if we want to use our FP words to do
that, what do we use?

1) The result that has the best chance of being close to the proper
result would be to use the FP number Z' that's closest to 10^n (round
to nearest).  But if Z'=prev(Z) (which will probably be the case half
of the time when 10^n cannot be represented exactly), the result of
X/Z' will be larger than it should be, and, in particular, even if
X=Z, X/Z' might have a value >1.0, which is probably not desired.

2) Instead, you could normalize as X/Z, i.e., the number that's
closest to 10^n if you round away from 0 (what I described as
round-to-absolute-higher above).  But Z can be up to 1ulp away from
10^n, so the result of X/Z will, on averge be further from the
intended value than X/Z'.

Another option is to normalize in a more precise way than performing a
double-double FP multiplication or division, but that's more costly.

- anton
-- 
M. Anton Ertl  http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
     New standard: https://forth-standard.org/
EuroForth 2026 CFP: http://www.euroforth.org/ef26/cfp.html

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

#135709

On 9/14/26 3:58 PM, Anton Ertl wrote:
> Krishna Myneni <krishna.myneni@ccreweb.org> writes:
>> On 9/14/26 10:00 AM, Anton Ertl wrote:
>> ...
>>>
>>> I am debating with myself whether one should use the closest FP number
>>> with round-to-nearest, or the closest with round-to-absolute-higher.
>>> Maybe depends on the purpose of this normalization.
>>>
>> Is your idea to stay within 1 ulp error for the conversion, but use a
>> simpler algorithm?
> 
> No.
> 
> The idea here is as follows:
> 
> You have your number X that you want to normalize to a power of 10.
> For that number there exist numbers Y, Y', Z, and Z', where
> 
> Y > X >= Z
> Y >= 10^(n+1) > prev(Y)
> Z >= 10^n     > prev(Z)
> 
> Where prev(V) produces the largest (in absolute terms) FP value <V.
> So Y or Z are either powers of 10 or the smallest FP numbers just
> above a power of 10, and they are those numbers which are the
> normalization boundaries of X.
> 
> For normalizing X, we would love to divide by 10^n, but if that's not
> representable as FP number, and if we want to use our FP words to do
> that, what do we use?
> 
> 1) The result that has the best chance of being close to the proper
> result would be to use the FP number Z' that's closest to 10^n (round
> to nearest).  But if Z'=prev(Z) (which will probably be the case half
> of the time when 10^n cannot be represented exactly), the result of
> X/Z' will be larger than it should be, and, in particular, even if
> X=Z, X/Z' might have a value >1.0, which is probably not desired.
> 
> 2) Instead, you could normalize as X/Z, i.e., the number that's
> closest to 10^n if you round away from 0 (what I described as
> round-to-absolute-higher above).  But Z can be up to 1ulp away from
> 10^n, so the result of X/Z will, on averge be further from the
> intended value than X/Z'.
> 
> Another option is to normalize in a more precise way than performing a
> double-double FP multiplication or division, but that's more costly.
> 

The present implementation of >DD decimal string to double-double 
conversion in kForth (dd_io.4th) is off by 14 ulp for the case of pi,
which is not good enough, so it will be useful to compare your proposed 
conversion procedure.

\ For 64-bit Forth with a separate fp stack

3.1415926535897931e0  1.2246467991473532e-16  ddconstant ddpi

32 set-precision
ddpi ddfs.
+3.1415926535897932384626433832795 dd 0 ok

ddvariable ref
ddvariable cnv

ddpi ref dd!  \ reference
\ convert decimal string to double double with >DD
S" 3.1415926535897932384626433832795dd0" >dd cnv dd!

cnv dd@ ddfs.
+3.1415926535897932384626433832791 dd 0 ok

17 set-precision
\ print high and low double precision fp numbers
ref dd@ fswap fs. 2 spaces fs.
3.1415926535897931e+00   1.2246467991473532e-16  ok
cnv dd@ fswap fs. 2 spaces fs.
3.1415926535897931e+00   1.2246467991473498e-16  ok

\ Note the loss of significant digits in the low double for cnv
\ What is the ulp error for cnv?
hex
ref 8 + @ u.  \ print the 64-bit hex value of ref, low double
3CA1A62633145C07  ok
cnv 8 + @ u.  \ print the 64-bit hex value of cnv, low double
3CA1A62633145BF9  ok

\ difference in ulp
decimal
cnv 8 + @ ref 8 + @ - .
-14  ok

--
Krishna

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

#135711

Krishna Myneni <krishna.myneni@ccreweb.org> writes:
>The present implementation of >DD decimal string to double-double 
>conversion in kForth (dd_io.4th) is off by 14 ulp for the case of pi,
>which is not good enough, so it will be useful to compare your proposed 
>conversion procedure.

I did not propose a conversion procedure, only an approach to
implementing NORMALIZE.  I think that NORMALIZE may be useful for DD
to string conversion, not the other way round.

- anton
-- 
M. Anton Ertl  http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
     New standard: https://forth-standard.org/
EuroForth 2026 CFP: http://www.euroforth.org/ef26/cfp.html

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

#135715

On 9/15/26 10:53 AM, Anton Ertl wrote:
> Krishna Myneni <krishna.myneni@ccreweb.org> writes:
>> The present implementation of >DD decimal string to double-double
>> conversion in kForth (dd_io.4th) is off by 14 ulp for the case of pi,
>> which is not good enough, so it will be useful to compare your proposed
>> conversion procedure.
> 
> I did not propose a conversion procedure, only an approach to
> implementing NORMALIZE.  I think that NORMALIZE may be useful for DD
> to string conversion, not the other way round.
> 

I misunderstood what you were saying.

There is a problem for converting in both directions. This is why I now 
use dtoa.c for coversions involving double precision floating point 
numbers, in both directions, across all kForth implementations.

Simple band aids are not a long-term solution.

--
KM

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

#135716

On 9/15/26 3:23 PM, Krishna Myneni wrote:
> On 9/15/26 10:53 AM, Anton Ertl wrote:
>> Krishna Myneni <krishna.myneni@ccreweb.org> writes:
>>> The present implementation of >DD decimal string to double-double
>>> conversion in kForth (dd_io.4th) is off by 14 ulp for the case of pi,
>>> which is not good enough, so it will be useful to compare your proposed
>>> conversion procedure.
>>
>> I did not propose a conversion procedure, only an approach to
>> implementing NORMALIZE.  I think that NORMALIZE may be useful for DD
>> to string conversion, not the other way round.
>>
> 
> I misunderstood what you were saying.
> 
> There is a problem for converting in both directions. This is why I now 
> use dtoa.c for coversions involving double precision floating point 
> numbers, in both directions, across all kForth implementations.
> 
> Simple band aids are not a long-term solution.
> 

To be clear, I think there should exist equivalent Forth code for 
accomplishing the conversions, both for double-precision and for 
double-double precision. But it helps to get there by using proven 
existing code from other languages as a reference against which to check 
any Forth implementations.

Also, I'm not against band aids for getting partially working code for 
use in the short term.

--
KM

[toc] | [prev] | [standalone]

Page 2 of 2 — ← Prev page 1 [2]

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

csiph-web