Double multiplication

[David Meyer]

This newbie may be missing something obvious, but Gforth has no words for
multiplying two double-precision integers, nor does it appear that anyone
on the Internet has any interest in doing so with Gforth or any other
Forth implementations. Is it really such a useless thing to do?

For my own satisfaction if no one else's, I've fiddled around and
pieced-together a way to multiply unsigned doubles. I hope a little more
fiddling will show me how to multiply signed doubles, too. The file
attached below is a work-in-progress, but comments are welcome.


\ double-arith.fs - Double-precision arithmetic extensions

\ Copyright 2013 David Meyer <pa...@sdf.org> +JMJ

\ Copying and distribution of this file, with or without
\ modification, are permitted in any medium without royalty
\ provided the copyright notice and this notice are preserved.
\ This file is offered as-is, without any warranty.

\ MAXU - Maximum value of unsigned single
s" MAX-U" environment? drop constant MAXU

\ md* - Multiply double by unsigned single (iterative method)
: md* ( d u -- d*u )
0. rot
0 u+do 2over d+ loop
2nip
;

\ mudu* - Multiply unsigned double by unsigned single
: mudu* ( ud u -- ud*u ) tuck * >r m* r> + ;

\ ud* - Multiply two unsigned doubles
: ud* ( ud1 ud2 -- ud1*ud2 )
{ a1 b1 a2 b2 }
a1 a2 um*
MAXU a1 um* b2 mudu* d+
MAXU a2 um* b1 mudu* d+
MAXU MAXU um* b1 mudu* b2 mudu* d+
;

--
pa...@sdf.org
SDF Public Access UNIX System - http://sdf.lonestar.org

[Anton Ertl]

David Meyer <pa...@sdf.org> writes:
>This newbie may be missing something obvious, but Gforth has no words for
>multiplying two double-precision integers, nor does it appear that anyone
>on the Internet has any interest in doing so with Gforth or any other
>Forth implementations. Is it really such a useless thing to do?

It's needed rarely.

>I hope a little more
>fiddling will show me how to multiply signed doubles, too.

If you (like everyone else) uses 2s-complement for signed numbers,
then for n-bit x n-bit -> n-bit multiplication there is no difference
between signed and unsigned. That's why Forth has no U*. And that's
why UD* is the same as D*.

This hinges on the fact that you can view the negative 2s-complement
number x as being the unsigned number 2^n+x (which is in the range of
unsigned numbers). If you multiply x with y, you get 2^n*y+x*y, and
2^n*y has only 0 bits in the lowest n bits, so it does not influence
the result.

Concerning the implementation, I'll have a go at it:

: d* { al ah bl bh -- cl ch }
\ (2^m*ah+al)+(2^m*bh+bl) =
\ 2^(2*m)*ah*bh+2^m*(ah*bl+al*bh)+al*bl
\ take this modulo 2^(2*m), making ah*bh irrelevant
al bh * ah bl * + 0 swap
al bl um* d+ ;

- 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 2013: http://www.euroforth.org/ef13/

[m...@iae.nl]

On Tuesday, November 26, 2013 5:06:32 PM UTC+1, David Meyer wrote:
> This newbie may be missing something obvious, but Gforth has no words for
> multiplying two double-precision integers, nor does it appear that anyone
> on the Internet has any interest in doing so with Gforth or any other
> Forth implementations. Is it really such a useless thing to do?

Where did you look?
Nathaniel Grossman, FD Vol VI, no 6, page 25 (available in scanned form).

Most Forths have
CR ." ** QUAD Toolkit ** "
CR ." UD+c ( ud1 ud2 carry -- ud3 carry )"
CR ." UDM* ( ud1 ud2 -- uquad )"
CR ." D* ( d1 d2 -- d3 )"
CR ." UT* ( ud u -- ut )"
CR ." UT/ ( ut u -- ud )"
CR
CR ." UMD/MOD ( uq1 ud1 -- ud_rem ud_quot ) ... is slow..."
CR ." UD/MOD ( ud1 ud2 -- ud_rem ud_quot )"
CR ." UD/ ( ud1 ud2 -- ud_quot )"
CR ." UDMOD ( ud1 ud2 -- ud_rem )"
CR ." DUM/MOD ( uquad ud -- ud_rem ud_quot )"
CR ." DUM/ ( uquad ud -- ud_quot )"
CR
CR ." D/MOD ( d1 d2 -- d_rem d_quot )"
CR ." D/ ( d1 d2 -- d_quot )"
CR ." DMOD ( d1 d2 -- d_rem )"
CR ." DM/MOD ( quad d -- d_rem d_quot )"
CR ." DM/ ( quad d -- d_quot )"
CR
CR ." DTUCK ( d1 d2 -- d2 d1 d2 )"
CR ." DNIP ( d1 d2 -- d2 )"
CR ." DU2/ ( ud1 -- ud2 )" ;

-marcel

[Albert van der Horst]

In article <l72gu7$qj0$1...@odin.sdf-eu.org>, David Meyer <pa...@sdf.org> wrote:
>This newbie may be missing something obvious, but Gforth has no words for
>multiplying two double-precision integers, nor does it appear that anyone
>on the Internet has any interest in doing so with Gforth or any other
>Forth implementations. Is it really such a useless thing to do?

I remember a quadruple precision package in the 16 bit era, because
"California real estate has grown to over $42,942,979.67 "
But nowadays gforth can accomodate the US national deficit
in a single, let alone a double.
Not useless maybe but pretty unusual.

What you miss is how parsimonous Forthers are. If the multiplication
result fits in a double, probably you can get by with a single times
double multiplication.
Other languages may have the attitude:
"Lets get it over with! Make it doubles all around!"

>--
>pa...@sdf.org
>SDF Public Access UNIX System - http://sdf.lonestar.org

Groetjes Albert
--
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

[Elizabeth D. Rather]

On 11/26/13 10:20 AM, Albert van der Horst wrote:
> In article <l72gu7$qj0$1...@odin.sdf-eu.org>, David Meyer <pa...@sdf.org> wrote:
>> This newbie may be missing something obvious, but Gforth has no words for
>> multiplying two double-precision integers, nor does it appear that anyone
>> on the Internet has any interest in doing so with Gforth or any other
>> Forth implementations. Is it really such a useless thing to do?
>
> I remember a quadruple precision package in the 16 bit era, because
> "California real estate has grown to over $42,942,979.67"
> But nowadays gforth can accomodate the US national deficit
> in a single, let alone a double.
> Not useless maybe but pretty unusual.
>
> What you miss is how parsimonous Forthers are. If the multiplication
> result fits in a double, probably you can get by with a single times
> double multiplication.
> Other languages may have the attitude:
> "Lets get it over with! Make it doubles all around!"

And the vast majority of high-precision multiplies are handled just fine
with M*/ (double multiplied by a ratio of singles) which does the
multiply first, getting an extended precision product, and then
dividing. This produces accurate answers over a very wide range.

IMO the */ and M*/ operators are among Chuck Moore's most brilliant
ideas. They really make accurate integer and fixed-point math quite
tractable.

Cheers,
Elizabeth

--
==================================================
Elizabeth D. Rather (US & Canada) 800-55-FORTH
FORTH Inc. +1 310.999.6784
5959 West Century Blvd. Suite 700
Los Angeles, CA 90045
http://www.forth.com

"Forth-based products and Services for real-time
applications since 1973."
==================================================

[Alex McDonald]

on 26/11/2013 16:13:17, anton ertl wrote:
> If you (like everyone else) uses 2s-complement for signed numbers,

We (like everyone else) should really get round to recognising in Forth
20xx that 2s complement has complete dominance in the "representation of
integer numbers" market place.

[Hans Bezemer]

Anton Ertl wrote:


> Concerning the implementation, I'll have a go at it:
>
> : d* { al ah bl bh -- cl ch }
> \ (2^m*ah+al)+(2^m*bh+bl) =
> \ 2^(2*m)*ah*bh+2^m*(ah*bl+al*bh)+al*bl
> \ take this modulo 2^(2*m), making ah*bh irrelevant
> al bh * ah bl * + 0 swap
> al bl um* d+ ;

Or, for those who can't really stand unnecessary locals:

: ud* >r swap >r over over um* rot r> * + rot r> * + ; ( ud1 ud2 -- ud3)

Hans Bezemer

[Anton Ertl]

Which shows that the D+ is unnecessary, because the low cell of one
addend is 0. Let's do this with locals:

: d* { al ah bl bh -- cl ch }

al bl um* al bh * + ah bl * + ;

[Ed]

m...@iae.nl wrote:
> ...

> Where did you look?
> Nathaniel Grossman, FD Vol VI, no 6, page 25 (available in scanned form).

The OP may not be as old as the rest of us :) That article was from an age
when 16-bit was the norm and Grossman was using a VIC-Forth cartridge.

> ...

> CR ." DTUCK ( d1 d2 -- d2 d1 d2 )"
> CR ." DNIP ( d1 d2 -- d2 )"

Conventionally these would be:

2TUCK ( x1 x2 x3 x4 -- x3 x4 x1 x2 x3 x4 )
2NIP ( x1 x2 x3 x4 -- x3 x4 )

as they work equally well with pairs of singles.

[m...@iae.nl]

On Sunday, December 1, 2013 8:41:55 PM UTC+8, Ed wrote:
>> m...@iae.nl wrote:
>> ...
>> Where did you look?
>> Nathaniel Grossman, FD Vol VI, no 6, page 25 (available in scanned form).
> The OP may not be as old as the rest of us :)
> That article was from an age when 16-bit was the norm and Grossman
> was using a VIC-Forth cartridge.

Maybe I should stop pushing Grossman. I am currently using
Knuth's _Seminumerical Algorithms_, and his code can very
easily be made to work for any cell sizes. ISTR that the FD
code is an unmaintainable stack horror.

>> ...
>> CR ." DTUCK ( d1 d2 -- d2 d1 d2 )"
>> CR ." DNIP ( d1 d2 -- d2 )"

> Conventionally these would be:
> 2TUCK ( x1 x2 x3 x4 -- x3 x4 x1 x2 x3 x4 )
> 2NIP ( x1 x2 x3 x4 -- x3 x4 )
> as they work equally well with pairs of singles.

Actually, I never use DTUCK (too deep for comforth).
DNIP is because I have seen 2NIP ( x2 x1 x0 -- x1 x0 ).

-marcel

[Ed]

m...@iae.nl wrote:
> ...

> DNIP is because I have seen 2NIP ( x2 x1 x0 -- x1 x0 ).

Very rarely I would think.

I recently added 2NIP to my Forth kernel on the basis 2SWAP 2DROP
was occuring sufficiently in my kernel and apps. The adoption of
'addr len' strings has made 2NIP almost a necessity (especially for
non-optimizing compilers).

[Hans Bezemer]

2TUCK is done with 9 instructions:

: 2tuck over >r rot >r rot over r> r> rot ;

You won't define 2OVER so easily with these basic instructions.

Hans Bezemer

[Coos Haak]

Op Thu, 05 Dec 2013 13:14:43 +0100 schreef Hans Bezemer:

When you are missing ROT as a primitive:

: 2tuck over >r >r swap r> swap >r >r swap r> swap over r> r> >r swap r> swap ;

Some snails are slower than others :)

--
Coos

CHForth, 16 bit DOS applications
http://home.hccnet.nl/j.j.haak/forth.html

[Rod Pemberton]

These were found using a modified version of Peter
Sovietov's Forth Wizard.

: r@ r> dup >r ;
: rot >r swap r> swap ;
: -rot swap >r swap r> ;
: nip swap drop ;
: tuck dup >r swap r> ;
: over >r dup r> swap ;

: 2drop drop drop ;
: 2>r swap >r >r ;
: 2r> r> r> swap ;
: 2r@ r> r> dup >r swap dup >r ;
: 2nip >r >r drop drop r> r> ;
: 2dup >r dup r> dup >r swap r> ;
: 2swap >r swap r> swap >r >r swap r> swap r> ;

: 2tuck over >r rot >r rot over r> r> rot ;

: 2over >r >r over over r> r> rot >r rot r> ;
: 2rot >r >r 2swap r> r> 2swap ;
: -2rot 2swap >r >r 2swap r> r> ;

2tuck 2over 2rot -2rot were taking up too much time
finding solutions without rot over 2swap. I expect
most of them to expand to around 20 or more. Coos'
solution for 2tuck without rot is at 18 using over.


Rod Pemberton