On Monday, March 14, 2022 at 11:29:44 AM UTC+1, minf...@arcor.de wrote:
[..]
> Paul Rubin schrieb am Montag, 14. März 2022 um 10:55:58 UTC+1:
> We know that fp-numbers are mostly stored with base 2 and that conversion
> to/from base 10 representation is not always possible without rounding.
>
> But why slow fat bignum packages are required to do that rounding is beyond me.

FORTH> : showm 24 0 do i 2^x S>F 1/F cr I 2 .r 8 htab +e. loop ; ok
FORTH> showm
0 1.0000000000000000000e+0000
1 5.0000000000000000000e-0001
2 2.5000000000000000000e-0001
3 1.2500000000000000000e-0001
4 6.2500000000000000000e-0002
5 3.1250000000000000000e-0002
6 1.5625000000000000000e-0002
7 7.8125000000000000018e-0003
8 3.9062500000000000004e-0003
9 1.9531250000000000002e-0003
10 9.7656250000000000000e-0004
11 4.8828125000000000000e-0004
12 2.4414062500000000000e-0004
13 1.2207031250000000000e-0004
14 6.1035156250000000000e-0005
15 3.0517578125000000000e-0005
16 1.5258789062500000000e-0005
17 7.6293945312499999992e-0006
18 3.8146972656249999996e-0006
19 1.9073486328124999998e-0006
20 9.5367431640625000000e-0007
21 4.7683715820312500000e-0007
22 2.3841857910156250000e-0007
23 1.1920928955078125000e-0007 ok

I think this hint leads the way.

-marcel

On Monday, March 14, 2022 at 10:41:20 PM UTC+1, Marcel Hendrix wrote:
> On Monday, March 14, 2022 at 11:29:44 AM UTC+1, minf...@arcor.de wrote:
> [..]
> > Paul Rubin schrieb am Montag, 14. März 2022 um 10:55:58 UTC+1:
> > We know that fp-numbers are mostly stored with base 2 and that conversion
> > to/from base 10 representation is not always possible without rounding.
> >
> > But why slow fat bignum packages are required to do that rounding is beyond me.

FORTH> NEEDS -xopg
FORTH> 0-VALUE a
FORTH> : showm LET a=2: #65 0 DO CR I 2 .r 8 htab LET a=2*a: LET. 1/a: LOOP ; ok
FORTH> showm
0 2.5000000000000000000000000000000000000000e-0001
1 1.2500000000000000000000000000000000000000e-0001
2 6.2500000000000000000000000000000000000000e-0002
3 3.1250000000000000000000000000000000000000e-0002
4 1.5625000000000000000000000000000000000000e-0002
5 7.8125000000000000000000000000000000000000e-0003
6 3.9062500000000000000000000000000000000000e-0003
7 1.9531250000000000000000000000000000000000e-0003
8 9.7656250000000000000000000000000000000000e-0004
9 4.8828125000000000000000000000000000000000e-0004
10 2.4414062500000000000000000000000000000000e-0004
11 1.2207031250000000000000000000000000000000e-0004
12 6.1035156250000000000000000000000000000000e-0005
13 3.0517578125000000000000000000000000000000e-0005
14 1.5258789062500000000000000000000000000000e-0005
15 7.6293945312500000000000000000000000000000e-0006
16 3.8146972656250000000000000000000000000000e-0006
17 1.9073486328125000000000000000000000000000e-0006
18 9.5367431640625000000000000000000000000000e-0007
19 4.7683715820312500000000000000000000000000e-0007
20 2.3841857910156250000000000000000000000000e-0007
21 1.1920928955078125000000000000000000000000e-0007
22 5.9604644775390625000000000000000000000000e-0008
23 2.9802322387695312500000000000000000000000e-0008
24 1.4901161193847656250000000000000000000000e-0008
25 7.4505805969238281250000000000000000000000e-0009
26 3.7252902984619140625000000000000000000000e-0009
27 1.8626451492309570312500000000000000000000e-0009
28 9.3132257461547851562500000000000000000000e-0010
29 4.6566128730773925781250000000000000000000e-0010
30 2.3283064365386962890625000000000000000000e-0010
31 1.1641532182693481445312500000000000000000e-0010
32 5.8207660913467407226562500000000000000000e-0011
33 2.9103830456733703613281250000000000000000e-0011
34 1.4551915228366851806640625000000000000000e-0011
35 7.2759576141834259033203125000000000000000e-0012
36 3.6379788070917129516601562500000000000000e-0012
37 1.8189894035458564758300781250000000000000e-0012
38 9.0949470177292823791503906250000000000000e-0013
39 4.5474735088646411895751953125000000000000e-0013
40 2.2737367544323205947875976562500000000000e-0013
41 1.1368683772161602973937988281250000000000e-0013
42 5.6843418860808014869689941406250000000000e-0014
43 2.8421709430404007434844970703125000000000e-0014
44 1.4210854715202003717422485351562500000000e-0014
45 7.1054273576010018587112426757812500000000e-0015
46 3.5527136788005009293556213378906250000000e-0015
47 1.7763568394002504646778106689453125000000e-0015
48 8.8817841970012523233890533447265625000000e-0016
49 4.4408920985006261616945266723632812500000e-0016
50 2.2204460492503130808472633361816406250000e-0016
51 1.1102230246251565404236316680908203125000e-0016
52 5.5511151231257827021181583404541015625000e-0017
53 2.7755575615628913510590791702270507812500e-0017
54 1.3877787807814456755295395851135253906250e-0017
55 6.9388939039072283776476979255676269531250e-0018
56 3.4694469519536141888238489627838134765625e-0018
57 1.7347234759768070944119244813919067382812e-0018
58 8.6736173798840354720596224069595336914062e-0019
59 4.3368086899420177360298112034797668457031e-0019
60 2.1684043449710088680149056017398834228515e-0019
61 1.0842021724855044340074528008699417114257e-0019
62 5.4210108624275221700372640043497085571289e-0020
63 2.7105054312137610850186320021748542785644e-0020
64 1.3552527156068805425093160010874271392822e-0020 ok
FORTH>

Already more than 40 digits necessary. I remember something about 256,
so maybe I should retest above 2^-56.

-marcel

On Monday, March 14, 2022 at 11:03:30 PM UTC+1, Marcel Hendrix wrote:
> On Monday, March 14, 2022 at 10:41:20 PM UTC+1, Marcel Hendrix wrote:
> > On Monday, March 14, 2022 at 11:29:44 AM UTC+1, minf...@arcor.de wrote:
> > [..]
> > > Paul Rubin schrieb am Montag, 14. März 2022 um 10:55:58 UTC+1:
> > > We know that fp-numbers are mostly stored with base 2 and that conversion
> > > to/from base 10 representation is not always possible without rounding.
> > >
> > > But why slow fat bignum packages are required to do that rounding is beyond me.

Or shorter:
FORTH> 1e-25 FULP +e. 5.6051938572992682836e-0045 ok

> -marcel

On 14/03/2022 21:29, minf...@arcor.de wrote:
> Paul Rubin schrieb am Montag, 14. März 2022 um 10:55:58 UTC+1:
>> "minf...@arcor.de" <minf...@arcor.de> writes:
>> > Out of that range it depends on the actual approximation definition.
>> Well, you have to get the decimal digits right. I have never examined
>> this issue so don't know exactly why the bignums are required. I might
>> look into it sometime. It's possible I have some part of the story
>> wrong.
>>
>> This seems to discuss it further down the page, but I haven't read it
>> all yet:
>>
>> https://recompilermag.com/issues/issue-9/a-crash-course-on-floating-point-numbers-and-good-algorithms-for-printing-and-reading-them/
>>
>> Decimal floating point seems like a good idea too, but that is more a
>> matter of wanting hardware support.
>
> We know that fp-numbers are mostly stored with base 2 and that conversion
> to/from base 10 representation is not always possible without rounding.
>
> But why slow fat bignum packages are required to do that rounding is beyond me.

Because there now exists faster bigger computers to implement them?

Marcel Hendrix schrieb am Montag, 14. März 2022 um 22:41:20 UTC+1:
> On Monday, March 14, 2022 at 11:29:44 AM UTC+1, minf...@arcor.de wrote:
> [..]
> > Paul Rubin schrieb am Montag, 14. März 2022 um 10:55:58 UTC+1:
> > We know that fp-numbers are mostly stored with base 2 and that conversion
> > to/from base 10 representation is not always possible without rounding.
> >
> > But why slow fat bignum packages are required to do that rounding is beyond me.
> FORTH> : showm 24 0 do i 2^x S>F 1/F cr I 2 .r 8 htab +e. loop ; ok
> FORTH> showm
> 0 1.0000000000000000000e+0000
> 1 5.0000000000000000000e-0001
> 2 2.5000000000000000000e-0001
> 3 1.2500000000000000000e-0001
> 4 6.2500000000000000000e-0002
> 5 3.1250000000000000000e-0002
> 6 1.5625000000000000000e-0002
> 7 7.8125000000000000018e-0003
> 8 3.9062500000000000004e-0003
> 9 1.9531250000000000002e-0003
> 10 9.7656250000000000000e-0004
> 11 4.8828125000000000000e-0004
> 12 2.4414062500000000000e-0004
> 13 1.2207031250000000000e-0004
> 14 6.1035156250000000000e-0005
> 15 3.0517578125000000000e-0005
> 16 1.5258789062500000000e-0005
> 17 7.6293945312499999992e-0006
> 18 3.8146972656249999996e-0006
> 19 1.9073486328124999998e-0006
> 20 9.5367431640625000000e-0007
> 21 4.7683715820312500000e-0007
> 22 2.3841857910156250000e-0007
> 23 1.1920928955078125000e-0007 ok
>
> I think this hint leads the way.

Thanks for trying to explain this to me stupid. Questions:
a) 2^x is a binary left shift, right?
b) S>F is lossless up to 2^32, right?
c) were you not testing the combined accuracy of 1/F _and_
internal conversion routines for numeric output?

Not wanting to be nicklish, I found a brief overview of algorithms
and their accuracies
https://www.cs.tufts.edu/%7Enr/cs257/archive/florian-loitsch/printf.pdf

BTW Wikipedia says that Grisu3 comes without bignums
https://en.wikipedia.org/wiki/Floating-point_arithmetic

There is also
https://stackoverflow.com/questions/29767027/fastest-c-way-to-convert-float-to-string

"minf...@arcor.de" <minforth@arcor.de> writes:
> BTW Wikipedia says that Grisu3 comes without bignums
> https://en.wikipedia.org/wiki/Floating-point_arithmetic

As Wikipedia also explains, Grisu3 doesn't always work. When it fails
(0.5% of inputs) you have to use a different algorithm, e.g. dragon4
which uses bignums. It also mentions errol3 which I'm not familiar
with, so idk if that uses bignums.

> There is also
> https://stackoverflow.com/questions/29767027/...

This mentions some algorithm names but doesn't say whether they use
bignums. As I remember, floating point conversion for most inputs is
straightforward, but there are edge cases that require more and more
precision to get all the digits right, and that's where the bignums come
in.

People keep talking about FPGA cpus here. It would be interesting if
someone made one with decimal floating point (IEEE 853 or whatever it
is).

minf...@arcor.de schrieb am Dienstag, 15. März 2022 um 08:01:27 UTC+1:
> Marcel Hendrix schrieb am Montag, 14. März 2022 um 22:41:20 UTC+1:
> > On Monday, March 14, 2022 at 11:29:44 AM UTC+1, minf...@arcor.de wrote:
> > [..]
> > > Paul Rubin schrieb am Montag, 14. März 2022 um 10:55:58 UTC+1:
> > > We know that fp-numbers are mostly stored with base 2 and that conversion
> > > to/from base 10 representation is not always possible without rounding.
> > >
> > > But why slow fat bignum packages are required to do that rounding is beyond me.
> > FORTH> : showm 24 0 do i 2^x S>F 1/F cr I 2 .r 8 htab +e. loop ; ok
> > FORTH> showm
> > 0 1.0000000000000000000e+0000
> > 1 5.0000000000000000000e-0001
> > 2 2.5000000000000000000e-0001
> > 3 1.2500000000000000000e-0001
> > 4 6.2500000000000000000e-0002
> > 5 3.1250000000000000000e-0002
> > 6 1.5625000000000000000e-0002
> > 7 7.8125000000000000018e-0003
> > 8 3.9062500000000000004e-0003
> > 9 1.9531250000000000002e-0003
> > 10 9.7656250000000000000e-0004
> > 11 4.8828125000000000000e-0004
> > 12 2.4414062500000000000e-0004
> > 13 1.2207031250000000000e-0004
> > 14 6.1035156250000000000e-0005
> > 15 3.0517578125000000000e-0005
> > 16 1.5258789062500000000e-0005
> > 17 7.6293945312499999992e-0006
> > 18 3.8146972656249999996e-0006
> > 19 1.9073486328124999998e-0006
> > 20 9.5367431640625000000e-0007
> > 21 4.7683715820312500000e-0007
> > 22 2.3841857910156250000e-0007
> > 23 1.1920928955078125000e-0007 ok
> >
> > I think this hint leads the way.
> Thanks for trying to explain this to me stupid. Questions:
> a) 2^x is a binary left shift, right?
> b) S>F is lossless up to 2^32, right?
> c) were you not testing the combined accuracy of 1/F _and_
> internal conversion routines for numeric output?

p. s. I forgot to mention that even 1/F has its pitfalls:
https://bits.stephan-brumme.com/inverse.html

dxforth <dxforth@gmail.com> writes:
>> But why slow fat bignum packages are required to do that rounding is
>> beyond me.

You want the conversion to have the following two properties:

1) if you convert float to decimal, and then convert the resulting
decimal back to float, you should get the exact same float that you
started with.

2) When you convert float to decimal, you should get the shortest
possible decimal string. To use an example from the Grisu paper[1], if
you convert the IEEE-754 float64 representation of 0.3 to decimal, you
should get 0.3, not 29999999999999998e-17, even though both will convert
back to the same float64.

Getting this to work properly is complicated and the first known method
used bignums. It's unclear to me whether bignums (or let's say ints
larger than 128 bits, i.e. double 64-bit cells) are absolutely necessary
to always get it right. One can reasonanly say, though, that using
bignums simplifies the problem.

[1] https://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf

"minf...@arcor.de" <minforth@arcor.de> writes:
> a) 2^x is a binary left shift, right?

No, floats are a "mantissa" and an exponent. Multiplying by 2 changes
the exponent.

> b) S>F is lossless up to 2^32, right?

For float64 binary, yes.

> c) were you not testing the combined accuracy of 1/F _and_
> internal conversion routines for numeric output?

It's not clear to me what that example was getting at.

Paul Rubin schrieb am Dienstag, 15. März 2022 um 10:39:05 UTC+1:
> "minf...@arcor.de" <minf...@arcor.de> writes:
> > a) 2^x is a binary left shift, right?
> No, floats are a "mantissa" and an exponent. Multiplying by 2 changes
> the exponent.
> > b) S>F is lossless up to 2^32, right?
> For float64 binary, yes.
> > c) were you not testing the combined accuracy of 1/F _and_
> > internal conversion routines for numeric output?
> It's not clear to me what that example was getting at.

I had responded to Marcel Hendrix' code:
FORTH> : showm 24 0 do i 2^x S>F 1/F cr I 2 .r 8 htab +e. loop ; ok

where
2^x seems to mean LSHIFT
and
S>F 1/F and +e.

could deteriorate accuracy even more than single binary to decimal
fp number conversions

On 15/03/2022 20:37, Paul Rubin wrote:
> dxforth <dxforth@gmail.com> writes:
>>> But why slow fat bignum packages are required to do that rounding is
>>> beyond me.

Misquoted but as I also think it's overkill ...

>
> You want the conversion to have the following two properties:
>
> 1) if you convert float to decimal, and then convert the resulting
> decimal back to float, you should get the exact same float that you
> started with.

How would you know? If I input:

0.333333333333333333e ok F:-1

then print it out:

18 set-precision cr f.
0.333333333333333333 ok

I don't know what was actually stored in memory, or how it might affect
calculations, but since it outputs what I keyed in 'ignorance is bliss'.
Only when it displays something different:

9.99999999999999999e f. 10. ok

do I say it shouldn't have done that - the conversion must be inadequate :)

>
> 2) When you convert float to decimal, you should get the shortest
> possible decimal string. To use an example from the Grisu paper[1], if
> you convert the IEEE-754 float64 representation of 0.3 to decimal, you
> should get 0.3, not 29999999999999998e-17, even though both will convert
> back to the same float64.

Same as 1) above. Not getting what he expected the novice is pissed off.
Understanding f/p, a programmer is not fazed by such results and doesn't
go looking to 'fix it' wasting resources and execution speed.

dxforth schrieb am Dienstag, 15. März 2022 um 12:37:49 UTC+1:
> On 15/03/2022 20:37, Paul Rubin wrote:
> > dxforth <dxf...@gmail.com> writes:
> >>> But why slow fat bignum packages are required to do that rounding is
> >>> beyond me.
> Misquoted but as I also think it's overkill ...
> >
> > You want the conversion to have the following two properties:
> >
> > 1) if you convert float to decimal, and then convert the resulting
> > decimal back to float, you should get the exact same float that you
> > started with.
> How would you know? If I input:
>
> 0.333333333333333333e ok F:-1
>
> then print it out:
>
> 18 set-precision cr f.
> 0.333333333333333333 ok
>
> I don't know what was actually stored in memory, or how it might affect
> calculations, but since it outputs what I keyed in 'ignorance is bliss'.
> Only when it displays something different:
>
> 9.99999999999999999e f. 10. ok
>
> do I say it shouldn't have done that - the conversion must be inadequate :)
> >
> > 2) When you convert float to decimal, you should get the shortest
> > possible decimal string. To use an example from the Grisu paper[1], if
> > you convert the IEEE-754 float64 representation of 0.3 to decimal, you
> > should get 0.3, not 29999999999999998e-17, even though both will convert
> > back to the same float64.
> Same as 1) above. Not getting what he expected the novice is pissed off.
> Understanding f/p, a programmer is not fazed by such results and doesn't
> go looking to 'fix it' wasting resources and execution speed.

Some basics
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
and The Base Conversion Problem in
https://www.math.utah.edu/~beebe/talks/2007/iccmse-2007/decfp.pdf

"minf...@arcor.de" <minforth@arcor.de> writes:
>Paul Rubin schrieb am Dienstag, 15. M=C3=A4rz 2022 um 10:39:05 UTC+1:
>> "minf...@arcor.de" <minf...@arcor.de> writes:=20
>> > a) 2^x is a binary left shift, right?
>> No, floats are a "mantissa" and an exponent. Multiplying by 2 changes=20
>> the exponent.
>> > b) S>F is lossless up to 2^32, right?
>> For float64 binary, yes.
>> > c) were you not testing the combined accuracy of 1/F _and_=20
>> > internal conversion routines for numeric output?
>> It's not clear to me what that example was getting at.
>
>I had responded to Marcel Hendrix' code:
>FORTH> : showm 24 0 do i 2^x S>F 1/F cr I 2 .r 8 htab +e. loop ; ok=20
>
>where
>2^x seems to mean LSHIFT

2^x obviously means 2^x in this context. It can be implemented with
lshift, however.

>and
>S>F 1/F and +e.
>
>could deteriorate accuracy even more than single binary to decimal
>fp number conversions

S>F is exact for all results of 2^x. I assume that 1/F produces the
same result as "1e fswap f/". For IEEE FP F/ has to produce the exact
result if it can be represented, or the value that you get with the
current rounding mode from the exact value if the exact value cannot
be represented. For 0<=n<24, 2^-n can be represented exactly as FP
value (for any not too outlandish binary FP representation with 6 bits
of exponent or more, and 0 bits of mantissa (aka significand) or
more). So for any competent FP implementation, there is no loss of
accuracy while in the FP domain.

Conversion to a string is what this discussion is about. I see some
obvious mistakes in the output cited in
<183a605d-ce70-4164-9c02-3c16253e8f45n@googlegroups.com>.

- 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: http://www.forth200x.org/forth200x.html
EuroForth 2021: https://euro.theforth.net/2021

On Friday, March 11, 2022 at 11:51:04 PM UTC+1, Anton Ertl wrote:
> "minf...@arcor.de" <minf...@arcor.de> writes:
> >I was wondering whether there is much to be won by making #
> >faster and thus more complex. I don't want negativity in here
> >but # is practically always followed by some terminal output
> >which is by magnitudes slower than #.
> Let's test the latter claim first:
>
> time gforth-fast -e ": foo 1000000 0 do 1111111111 . loop ; foo bye" >/dev/null
> real 0m0.310s
> user 0m0.310s
> sys 0m0.000s
> time gforth-fast -e ': foo 1000000 0 do ." 1111111111 " loop ; foo bye'
> ...
> real 0m1.267s
> user 0m0.298s
> sys 0m0.692s
>
> So # is indeed faster than output on an xterm, although not by a
> magnitude, much less magnitudes.
Indeed - testing 4tH (where all <#, #, #S, #>, etc. words are coded in C) produces similar results:
(1) There is hardly any difference between outputting both variants on xterm (about 10% difference);
(2) Numerical output being 4 times slower than string output on /dev/null (see below):

real 0m0,037s
user 0m0,032s
sys 0m0,005s

real 0m0,152s
user 0m0,152s
sys 0m0,000s

Hans Bezemer

dxforth <dxforth@gmail.com> writes:
>Understanding f/p, a programmer is not fazed by such results and doesn't
>go looking to 'fix it' wasting resources and execution speed.

That's the Seymour Cray school of floating point. However, the
William Kahan school of FP has won out (even in, e.g., the Cray-T3E).

- 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: http://www.forth200x.org/forth200x.html
EuroForth 2021: https://euro.theforth.net/2021

On 15/03/2022 22:57, minf...@arcor.de wrote:
> dxforth schrieb am Dienstag, 15. März 2022 um 12:37:49 UTC+1:
>> On 15/03/2022 20:37, Paul Rubin wrote:
>> > dxforth <dxf...@gmail.com> writes:
>> >>> But why slow fat bignum packages are required to do that rounding is
>> >>> beyond me.
>> Misquoted but as I also think it's overkill ...
>> >
>> > You want the conversion to have the following two properties:
>> >
>> > 1) if you convert float to decimal, and then convert the resulting
>> > decimal back to float, you should get the exact same float that you
>> > started with.
>> How would you know? If I input:
>>
>> 0.333333333333333333e ok F:-1
>>
>> then print it out:
>>
>> 18 set-precision cr f.
>> 0.333333333333333333 ok
>>
>> I don't know what was actually stored in memory, or how it might affect
>> calculations, but since it outputs what I keyed in 'ignorance is bliss'.
>> Only when it displays something different:
>>
>> 9.99999999999999999e f. 10. ok
>>
>> do I say it shouldn't have done that - the conversion must be inadequate :)
>> >
>> > 2) When you convert float to decimal, you should get the shortest
>> > possible decimal string. To use an example from the Grisu paper[1], if
>> > you convert the IEEE-754 float64 representation of 0.3 to decimal, you
>> > should get 0.3, not 29999999999999998e-17, even though both will convert
>> > back to the same float64.
>> Same as 1) above. Not getting what he expected the novice is pissed off.
>> Understanding f/p, a programmer is not fazed by such results and doesn't
>> go looking to 'fix it' wasting resources and execution speed.
>
> Some basics
> https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
> and The Base Conversion Problem in
> https://www.math.utah.edu/~beebe/talks/2007/iccmse-2007/decfp.pdf
>

The latter refers to "THE BASE-CONVERSION PROBLEM" - which reads like
the front page of a Murdoch tabloid. Had it discussed 'WHY BINARY CAME
TO DISPLACE BCD' it may have sold less copies.

On 15/03/2022 23:21, Anton Ertl wrote:
> dxforth <dxforth@gmail.com> writes:
>>Understanding f/p, a programmer is not fazed by such results and doesn't
>>go looking to 'fix it' wasting resources and execution speed.
>
> That's the Seymour Cray school of floating point. However, the
> William Kahan school of FP has won out (even in, e.g., the Cray-T3E).

Only if winning is ...

"safe robust defaults for numerically unsophisticated programmers [...]
moderately tolerant of well-meaning ignorance among programmers" - W. Kahan

On Mon, 14 Mar 2022 03:29:43 -0700 (PDT)
"minf...@arcor.de" <minforth@arcor.de> wrote:

> We know that fp-numbers are mostly stored with base 2 and that conversion
> to/from base 10 representation is not always possible without rounding.
>
> But why slow fat bignum packages are required to do that rounding is beyond me.

This shows how the intermediate product grows as the significant digits
move further from the radix point:

Decimal Float ConvErr Intermediate Product

999e+00 999*2^+000 Exact 1998
999e+03 1951*2^+009 0.0088% 249750
999e+06 1905*2^+019 0.0231% 31218750
999e+09 1860*2^+029 0.042% 3902343750
999e+12 1817*2^+039 0.0093% 487792968750
999e+15 1774*2^+049 0.0327% 60974121093750
999e+18 1732*2^+059 0.057% 7621765136718750
999e+21 1692*2^+069 0.0219% 952720642089843750
999e+24 1652*2^+079 0.0427% 119090080261230468750
999e+27 1613*2^+089 0.0601% 14886260032653808593750
999e+30 1576*2^+099 0.0091% 1860782504081726074218750
999e+33 1539*2^+109 0.0132% 232597813010215759277343750

999e-33 2593*2^-111 0.0212% 301929223448753636382867456

The intermediate product size would be similar to that needed
for IEEE 754 binary16 floats.

Because my goal is to maintain precision and track error, rather
than rounding, the algorithm may be different from that usually
used for conversion.

I believe that for scientific and similar work it is more important
to know the resolution and error bound than to see correct looking
digits and have to guesstimate true resolution and accumulated error.

In the above the error bound is the size of the least significant
(binary) digit (ULP), unless number is marked as 'Exact'.

Jan Coombs
--

On 3/11/22 07:46, Krishna Myneni wrote:
> On 3/11/22 06:09, Anton Ertl wrote:
>> ...
>> The word # (used in, e.g., ".") internally calls UD/MOD, which
>> requires two double-by-single divisions:
>>
>> : ud/mod ( ud1 u2 -- urem udquot ) \ gforth
>>      \G divide unsigned double @i{ud1} by @i{u2}, resulting in a
>> unsigned double
>>      \G quotient @i{udquot} and a single remainder @i{urem}.
>>      >r 0 r@ um/mod r> swap >r
>>      um/mod r> ;
>>
>> : # ( ud1 -- ud2 )
>>      base @ ud/mod rot dup 9 u>
>>      [ char A char 9 1+ - ] Literal and +
>>      '0' + hold ;
>>
> ...
>
> Interesting. UD/MOD is not present in kForth. When I looked at the
> implementation of "#" it uses the non-standard word UTM/ to perform the
> division, which should be far more expensive. So there's a significant
> speedup to be had for number to string conversion here.
>

I implemented the non-standard division words U/MOD and UD/MOD in
kForth-64 (in my local repo for now). Replacing the call to UTM/ and
extra arithmetic used in the prior implementation of "#" (function
C_sharp() in source file vmc.c) with a call to UD/MOD gives a
performance boost of more than 2x with the simple test,

: test ( -- uelapsed_ms )
ms@ 1000000 0 DO
32456789019204 I - dup
<# #s #> 2drop \ #S repeatedly executes #
LOOP ms@ swap - ;

Execution times on my Linux PC:

kforth64 (v0.3.0 -- only on local repo)
---
test .
980 ok \ time in ms
---

kforth64 (v0.2.4 -- version on GitHub)
---
test .
2203 ok \ time in ms
---

This speed improvement is without even handling the special case of
high(ud) = 0, which I will try next.

In the course of implementing U/MOD and UD/MOD, I reexamined the
implementation of several division words ( / MOD /MOD ) and was able to
increase their efficiency by up to 20%.

--
Krishna

>>
>> Special-casing high=0 is beneficial in all cases.
>>
>> Special-casing base=#10 only helps the Skylake.  On the other CPUs
>> apparently division is cheap enough and the overhead of separating
>> this case is expensive enough that separating this case does not pay
>> off.  The slowdown from this special case on the A72 is surprisingly
>> large.
>>
>> Overall, my conclusion is to keep the high=0 special case, and forget
>> about the base=#10 one.
>>

On 3/17/22 22:01, Krishna Myneni wrote:
> On 3/11/22 07:46, Krishna Myneni wrote:
>> On 3/11/22 06:09, Anton Ertl wrote:
>>> ...
>>> The word # (used in, e.g., ".") internally calls UD/MOD, which
>>> requires two double-by-single divisions:
....
> I implemented the non-standard division words U/MOD and UD/MOD in
> kForth-64 (in my local repo for now). Replacing the call to UTM/ and
> extra arithmetic used in the prior implementation of "#" (function
> C_sharp() in source file vmc.c) with a call to UD/MOD gives a
> performance boost of more than 2x with the simple test,
>
> : test ( -- uelapsed_ms )
>     ms@ 1000000 0 DO
>       32456789019204 I - dup
>       <# #s #> 2drop  \ #S repeatedly executes #
>     LOOP ms@ swap - ;
>
> Execution times on my Linux PC:
>
> kforth64 (v0.3.0 -- only on local repo)
> ---
> test .
> 980  ok  \ time in ms
> ---
....

Interestingly this same code on gforth-fast (v0.7.9_20220120) takes
longer to execute the same code, relative to my test for the new
kforth64, on the same PC:

---
include contrib/kforth-compat.fs ok
: test ms@ 1000000 0 do 32456789019204 I - dup <# #s #> 2drop LOOP ms@
swap - ; ok

test . 1224 ok
---

The file kforth-compat.fs is loaded to provide a definition of MS@ for
Gforth. I need to pull Gforth updates and rebuild to test the latest
version. Generally, gforth-fast is about 5--10x faster than kforth64.

--
Krishna

On 18/03/2022 14:01, Krishna Myneni wrote:
>
> I implemented the non-standard division words U/MOD and UD/MOD in
> kForth-64 (in my local repo for now). Replacing the call to UTM/ and
> extra arithmetic used in the prior implementation of "#" (function
> C_sharp() in source file vmc.c) with a call to UD/MOD gives a
> performance boost of more than 2x with the simple test,

'Non-standard' maybe - but 'common practice' since at least Fig-Forth:

M/MOD ud1 u2 --- u3 ud4
An unsigned mixed magnitude math operation which leaves a
double quotient ud4 and remainder u3, from a double
dividend ud1 and single divisor u2.

Other names for it include MU/MOD (F83, F-PC). In all cases it was used
to implement (not surprisingly) # .

On 3/18/22 02:25, dxforth wrote:
> On 18/03/2022 14:01, Krishna Myneni wrote:
>>
>> I implemented the non-standard division words U/MOD and UD/MOD in
>> kForth-64 (in my local repo for now). Replacing the call to UTM/ and
>> extra arithmetic used in the prior implementation of "#" (function
>> C_sharp() in source file vmc.c) with a call to UD/MOD gives a
>> performance boost of more than 2x with the simple test,
>
> 'Non-standard' maybe - but 'common practice' since at least Fig-Forth:
>
> M/MOD ud1 u2 --- u3 ud4
> An unsigned mixed magnitude math operation which leaves a
> double quotient ud4 and remainder u3, from a double
> dividend ud1 and single divisor u2.
>
> Other names for it include MU/MOD (F83, F-PC). In all cases it was used
> to implement (not surprisingly) # .

Thanks for the historical info. David Williams introduced UTM/ into
kForth and, at the time, I used it to make a faster implementation of
"#". UD/MOD is a precise fit.

--
Krishna

On 19/03/2022 01:03, Krishna Myneni wrote:
> On 3/18/22 02:25, dxforth wrote:
>> On 18/03/2022 14:01, Krishna Myneni wrote:
>>>
>>> I implemented the non-standard division words U/MOD and UD/MOD in
>>> kForth-64 (in my local repo for now). Replacing the call to UTM/ and
>>> extra arithmetic used in the prior implementation of "#" (function
>>> C_sharp() in source file vmc.c) with a call to UD/MOD gives a
>>> performance boost of more than 2x with the simple test,
>>
>> 'Non-standard' maybe - but 'common practice' since at least Fig-Forth:
>>
>> M/MOD ud1 u2 --- u3 ud4
>> An unsigned mixed magnitude math operation which leaves a
>> double quotient ud4 and remainder u3, from a double
>> dividend ud1 and single divisor u2.
>>
>> Other names for it include MU/MOD (F83, F-PC). In all cases it was used
>> to implement (not surprisingly) # .
>
> Thanks for the historical info. David Williams introduced UTM/ into
> kForth and, at the time, I used it to make a faster implementation of
> "#". UD/MOD is a precise fit.

I went back to see what LMI used. LMI 1.x was a Fig-Forth compatible
product and consequently had M/MOD as described above. LMI 2.x had
switched to Forth-83 but retained M/MOD. In PC/Forth 3.2 it becomes
a nameless word and the name 'M/MOD' used for what we now call SM/REM.
While most, if not all, systems today use the function in '#' it seems
evenly divided whether it's provided as a discrete word to users and
under what name.

On Saturday, March 19, 2022 at 1:54:57 AM UTC+1, dxforth wrote:
> While most, if not all, systems today use the function in '#' it seems
> evenly divided whether it's provided as a discrete word to users and
> under what name.

It seems that the word ud/mod ( ud1 u2 -- urem udquot ) is not needed often.

In the gforth source files (*.fs), it is only used in the number i/o file (nio.fs) in order to
define #, and in the cross compiler (cross.fs) where it is used to define two words named
DS! and Sc!

ud/mod itself is defined in kernel.fs

I also searched for the phrase "um/mod r>" but could not find any additional matches.

SwiftForth does not provide a separate word for the UD/MOD routine. This makes the
related code for # hard to read. In this case I would have preferred a version that
factors UD/MOD out, even if this word is only used once in the entire system.
Factoring can improve clarity of code considerably.

While UD/MOD seems to have rather special use, LMI systems provide a related and
much more useful word:

: D/ ( d1 n+ -- d2) SWAP OVER /MOD >R SWAP UM/MOD SWAP DROP R> ;

This word allows you to divide a signed double length number by an unsigned single length
number (which must not exceed the range of signed numbers).

PC/FORTH also provides a counterpart used for multiplication:

: D* ( d1 n+ -- d2) DUP ROT * ROT ROT UM* ROT + ;

In 16-bit systems, signed double length numbers were often needed to store long
numbers. D* and D/ were quite useful to handle these numbers.

In 32-bit systems these words are still useful, but not needed as often as in the olden days.

Henry

On Saturday, March 19, 2022 at 2:07:31 PM UTC+1, Heinrich Hohl wrote:
> : D/ ( d1 n+ -- d2) SWAP OVER /MOD >R SWAP UM/MOD SWAP DROP R> ;

Note that PC/FORTH uses floored division. The D/ code shown above will crash ANS systems.
In ANS compatible systems you can use the following definition:

: D/ ( d1 n+ -- d2) >R S>D R@ FM/MOD R> SWAP >R UM/MOD NIP R> ;

Henry