double-double precision arithmetic in Forth
Julian V. Noble
Colleagues,
Here are portable double-double precision floating point arithmetic words.
I've tested them on Win32Forth 4.2 and they work fine. I regret they look
a tad ugly. It was necessary for portability. Refinements would involve
making the (many) fvariables (assumed to be 64-bit precision) local to the
routines that use them to avoid name collisions.
I have not spent any time on factoring or refinement.
The words will not work as-is on GForth because I believe (from my experiments)
that the latter is default-set to 80-bit internal precision rather than
64-bit. I haven't had a chance to work up a CODE word for doing this setting.
I will be translating more things as I get time, such as a DD function
library, DD I/O words, etc. I also intend to translate them to assembler,
which should make them a bit less ugly.
There are two files, arithmetic and a simple test suite.
\ Doubled-precision floating point arithmetic
\
\ ---------------------------------------------------
\ (c) Copyright 2006 Julian V. Noble. \
\ Permission is granted by the author to \
\ use this software for any application pro- \
\ vided this copyright notice is preserved. \
\ ---------------------------------------------------
\ This is an ANS Forth program requiring the
\ FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependences:
\ Assumes independent floating point stack
\ FPU must be set for IEEE 64-bit internal arithmetic
\ ^^^^
\ based on "Software for Doubled-Precision Floating-Point Computations"
\ by Seppo Linnainmaa
\ ACM Transactions on Mathematical Software,
\ Vol 7, No 3, September 1981, Pages 272-283
marker -double
FALSE [IF] \ use these to compile the versions using
include ftran202.f \ the FORmula TRANslator and Baden's flocals
include flocals.f
[THEN]
: fvariables: 0 DO FVARIABLE LOOP ;
\ ---------------------------------------- LOAD, STORE
: r128@ DUP F@ FLOAT+ F@ ;
: r128! DUP FLOAT+ F! F! ;
\ ------------------------------------ END LOAD, STORE
\ ----------------------------------- data types ----
: real*16 \ create a double-double variable
CREATE 2 FLOATS ALLOT ;
: ddconstant CREATE HERE r128! 2 FLOATS ALLOT
DOES> r128@ ;
\ ---------------------------------------------------
3.1415926535897931e0 1.2246467991473532e-16 ddconstant ddpi
\ based on "Software for Doubled-Precision Floating-Point Computations"
\ by Seppo Linnainmaa
\ ACM Transactions on Mathematical Software,
\ Vol 7, No 3, September 1981, Pages 272-283
: fvariables: 0 DO FVARIABLE LOOP ;
\ determine base and precision of fpu
((
4e0 3e0 F/ 1e0 F- 3e0 F* 1e0 F- FVALUE u
u 2e0 F/ 1e0 F+ 1e0 F- FVALUE r
r F0= NOT [IF] r FTO u [THEN]
2e0 3e0 F/ 0.5e0 F- 3e0 F* 0.5e0 F- FVALUE uu
uu 2e0 F/ 0.5e0 F+ 0.5e0 F- FVALUE rr
rr F0= NOT [IF] rr FTO uu [THEN]
u uu F/ ( f: -- beta)
uu FLN FOVER FLN F/ FNEGATE 0.5e0 F+
F>S F>S CR CR .( base = ) . .( precision = ) . FORGET u ))
\ ---------------- Exact multiplication -------------
134217729 S>F FCONSTANT split
9 fvariables: t a1 a2 b1 b2 b21 b22 ta tb
8 fvariables: q qq x xx y yy z1 z2
2 fvariables: aaa bbb
FALSE [IF] \ these versions use ftran202.f and flocals.f
: exactmul ( f: a b -- x xx) \ multiply 2 fp#'s to get ddfp#
FRAME| bb aa |
f" t = bb*[split]" f" a1 = (bb-t) + t" f" a2 = bb - a1"
f" t = aa*[split]" f" b1 = (aa-t) + t" f" b2 = aa - b1"
f" t = b2*[split]" f" b21 = (b2-t) + t" f" b22 = b2 - b21"
f" bb*aa" FDUP f" t=" ( f: x)
f" ((((a1*b1 - t)+a1*b2)+b1*a2)+b21*a2)+b22*a2"
|FRAME
; ( f: x xx)
: dd/ ( f: x xx y yy -- [x+xx]/[y+yy] )
f" yy=" f" y=" f" xx=" f" x="
f" abs(y)" F0= ABORT" Can't divide by 0!"
f" z1=x/y"
f" exactmul(y,z1)" f" qq=" f" q="
f" z2=((((x-q)-qq)+xx)-z1*yy)/(y+yy)"
f" q=z1+z2" f" q"
f" (z1-q)+z2"
;
: dd* ( f: x xx y yy -- [x+xx]*[y+yy] )
f" yy=" f" y=" f" xx=" f" x="
f" exactmul(x,y)" f" qq=" f" z1="
f" z2=((x+xx)*yy+xx*y)+qq"
f" q=z1+z2" f" q"
f" (z1-q)+z2"
;
: dd+ ( f: x xx y yy -- [x+xx] + [y+yy] )
f" yy=" f" y=" f" xx=" f" x="
f" z1=x+y"
f" q=x-z1"
f" z2=(((q+y)+(x-(q+z1)))+xx)+yy"
f" q=z1+z2" f" q"
f" (z1-q)+z2"
;
\ Square root based on T.J. Dekker,
\ "A Floating-Point Technique for Extending the Available Precision"
\ Numerische Mathematik 18 (1971) 224-242.
: ddsqrt ( f: x xx -- ddsqrt[x+xx])
f" xx=" f" x="
f" x" F0< ABORT" Can't take sqrt of negative number!"
f" q = sqrt(x)"
f" exactmul(q,q)" f" yy=" f" y="
f" qq = (((x-y)-yy)+xx)*0.5/q"
f" y=q + qq" f" y"
f" (q-y)+qq"
;
[THEN]
: exactmul ( f: a b -- x xx) \ multiply 2 fp#'s to get ddfp#
aaa F! bbb F!
bbb F@ split F* t F!
t F@ bbb F@ FOVER F- F+ FDUP a1 F! ( f: a1)
FNEGATE bbb F@ F+ a2 F!
aaa F@ FDUP split F* t F! ( f: aaa)
t F@ ftuck F- F+ b1 F!
aaa F@ b1 F@ F- FDUP FDUP b2 F! ( f: b2 b2)
split F* t F!
t F@ ftuck F- F+ FDUP b21 F! ( f: b21)
FNEGATE b2 F@ F+ b22 F!
bbb F@ aaa F@ F* FDUP t F!
a1 F@ b1 F@ F* t F@ F- a1 F@ b2 F@ F*
F+ b1 F@ a2 F@ F* F+
b21 F@ a2 F@ F* F+ b22 F@ a2 F@ F* F+
;
: dd/ ( f: x xx y yy -- [x+xx]/[y+yy] )
yy F! y F! xx F! x F!
y F@ FABS F0= ABORT" Can't divide by 0!"
x F@ y F@ F/ FDUP z1 F!
y F@ exactmul qq F! FDUP q F! ( f: q)
FNEGATE x F@ F+ qq F@ F-
xx F@ F+ z1 F@ yy F@ F* F-
y F@ yy F@ F+ F/ FDUP z2 F!
z1 F@ F+ FDUP
FNEGATE z1 F@ F+ z2 F@ F+
;
: dd* ( f: x xx y yy -- [x+xx]*[y+yy] )
yy F! y F! xx F! x F!
x F@ y F@ exactmul qq F! z1 F!
x F@ xx F@ F+ yy F@ F*
xx F@ y F@ F* F+ qq F@ F+ FDUP z2 F!
z1 F@ F+ FDUP
FNEGATE z1 F@ F+ z2 F@ F+
;
: dd+ ( f: x xx y yy -- [x+xx] + [y+yy] )
yy F! y F! xx F! x F!
x F@ y F@ F+ z1 F!
x F@ z1 F@ F- FDUP q F!
y F@ F+ ( f: q+y)
x F@ q F@ z1 F@ F+ F- ( f: q+y x-[q+z1])
F+ xx F@ F+ yy F@ F+ FDUP z2 F!
z1 F@ F+ FDUP ( f: z1+z2 z1+z2)
FNEGATE z1 F@ F+ z2 F@ F+
;
: ddnegate ( f: x xx -- -x -xx)
FNEGATE FSWAP FNEGATE FSWAP
;
: dd- ( f: x xx y yy -- [x+xx] - [y+yy] )
ddnegate dd+
;
\ ddsqrt based on T.J. Dekker,
\ "A Floating-Point Technique for Extending the Available Precision"
\ Numerische Mathematik 18 (1971) 224-242.
: ddsqrt ( f: x xx -- ddsqrt[x+xx])
xx F! FDUP x F!
F0< ABORT" Can't take sqrt of negative number!"
x F@ FSQRT FDUP q F!
FDUP exactmul yy F! FDUP y F!
FNEGATE x F@ F+
yy F@ F- xx F@ F+
0.5e0 F* q F@ F/ FDUP qq F!
q F@ F+ FDUP
FNEGATE q F@ F+ qq F@ F+
;
1e0 0e0 ddconstant dd=1
\ ----------------------------------- stack ops -----
fvariable dtemp
real*16 ddtemp real*16 ddtemp1
: ddswap ( f: x xx y yy -- y yy x xx )
dtemp F! F-ROT dtemp F@ F-ROT ;
: dddup FOVER FOVER ;
: dddrop FDROP FDROP ;
: ddover ddtemp r128! dddup ddtemp1 r128!
ddtemp r128@ ddtemp1 r128@ ;
: ddtuck ddtemp r128! ddtemp1 r128!
ddtemp r128@ ddtemp1 r128@
ddtemp r128@ ;
\ ---------------------------------------------------
: dd^2 dddup dd* ;
: dd^n ( n -- ) ( f: x xx -- [x+xx]^n )
\ raise dd to positive or negative integer power
\ return 1 if n=0, dd^{-|n|} if n<0
dd=1 ddswap ( f: 1e0 0e0 x xx )
DUP 0= IF dddrop dd=1 drop EXIT THEN
DUP 0< SWAP ABS ( -- sign |n| )
BEGIN DUP 0> WHILE
DUP 1 AND IF ddtuck dd* ddswap THEN z^2
2/
REPEAT dddrop DROP
IF dd=1 ddswap dd/ THEN
;
\ Tests for DoubleDouble
\ February 11th, 2006
MARKER -ddtest
15 set-precision
1e0 0e0 3e0 0e0 dd/ ddconstant dd1/3
1e0 0e0 6e0 0e0 dd/ ddconstant dd1/6
dd1/3 dd1/6 dd+ FSWAP
CR
CR .( Addition: 1/3 + 1/6)
CR 5 spaces .( 5.00000000000000E-1 -1.54074395550979E-33 <- should be this)
CR 5 spaces FS. FS.
CR
dd1/3 dd1/6 dd- FSWAP
CR .( Subtraction: 1/3 - 1/6)
CR 5 spaces .( 1.66666666666667E-1 9.25185853854297E-18 <- should be this)
CR 5 spaces FS. FS.
CR
6e0 0e0 dd1/6 dd* FSWAP
CR .( Multiplication: 6 * 1/6)
CR 5 spaces .( 1.00000000000000E0 0.00000000000000E-1 <- should be this)
CR 5 spaces FS. FS.
CR .( 6 * 1/3)
6e0 0e0 dd1/3 dd* FSWAP
CR 5 spaces .( 2.00000000000000E0 0.00000000000000E-1 <- should be this)
CR 5 spaces fs. fs.
CR
dd1/3 dd1/6 dd/ FSWAP
CR .( Division: [1/3] / [1/6] )
CR 5 spaces .( 2.00000000000000E0 0.00000000000000E-1 <- should be this)
CR 5 spaces FS. FS.
CR
CR .( [1/6] / [1/3])
dd1/6 dd1/3 dd/ FSWAP
CR 5 spaces .( 5.00000000000000E-1 0.00000000000000E-1 <- should be this)
CR 5 spaces FS. FS.
CR .( Square root: sqrt[1/324] )
dd1/6 dd1/3 dd* FOVER FOVER dd* ddsqrt FOVER FOVER FSWAP
CR 5 spaces FS. FS. .( <- sqrt[1/324] = 1/18)
CR 5 spaces .( 5.55555555555556E-2 3.08395284618099E-18 <- should be this)
CR 5 spaces .( multiply by 18)
18e0 0e0 dd* FSWAP
CR 5 spaces .( 1.00000000000000E0 0.00000000000000E-1 <- should be this)
CR 5 spaces FS. FS.
--
Julian V. Noble
Professor Emeritus of Physics
http://galileo.phys.virginia.edu/~jvn/
"For there was never yet philosopher that could endure the
toothache patiently."
-- Wm. Shakespeare, Much Ado about Nothing. Act v. Sc. 1.
dbu
> Colleagues,
>
> Here are portable double-double precision floating point arithmetic words.
> I've tested them on Win32Forth 4.2 and they work fine.
Hello Julian,
I have checked it with Win32Forth 6.11.08 and it works fine, too.
But only the code which is using your FORmula TRANslator (there's
no FTUCK in w32f).
by
Dirk
Marcel Hendrix
> Here are portable double-double precision floating point arithmetic words.
> I've tested them on Win32Forth 4.2 and they work fine. I regret they look
> a tad ugly. It was necessary for portability. Refinements would involve
> making the (many) fvariables (assumed to be 64-bit precision) local to the
> routines that use them to avoid name collisions.
Contrary to my expectations, this code works flawlessly with iForth 2.0.
I've shown how to set the FPU control word for older iForths, but newer
releases (as of today) will have the words built-in.
The line "uu FLN FOVER FLN F/ FNEGATE 0.5e F+" did not work
for me (FLN doesn't take negative arguments), so I write
uu FABS FLN FOVER FABS FLN F/ FNEGATE 0.5e F+
instead.
> I will be translating more things as I get time, such as a DD function
> library, DD I/O words, etc. I also intend to translate them to assembler,
> which should make them a bit less ugly.
I think they WILL look better in high-level, see attachment :-)
It will be nice to see what you come up with.
-marcel
-- ---------------
[1] iForth server 1.00 (console), Jan 24 2006, 22:41:58.
[2] Stuffed iForth at 00438178 [entry: 0x440000]
[3] Current process priority is 32.
iForth version 2.1.2070, generated 11:44:29, February 12, 2006.
i6 binary, native floating-point, double precision.
Copyright 1996 - 2006 Marcel Hendrix.
FORTH>
Creating --- Several utilities Version 3.05 ---
Creating --- Extended OS words Version 3.16 ---
Creating --- Terminal Driver Version 3.12 ---
Creating --- Command line Editor Version 1.26 ---
Creating --- Online help Version 1.34 ---
Creating --- Glossary Generator Version 1.05 ---
Creating --- Disassembler Version 2.10 --- ok
FORTH> in dd-fp
Creating --- 128-bit FP words Version 1.00 ---
Warning: 53-bit FP mode set.
base = 2
precision = 54
** DD-FP words ***
Try: TEST
ok
FORTH> TEST
Addition: 1/3 + 1/6
5.00000000000000e-001 -1.54074395550979e-033 <- should be this
5.00000000000000e-001 -1.54074395550979e-033
Subtraction: 1/3 - 1/6
1.66666666666667e-001 9.25185853854297e-018 <- should be this
1.66666666666667e-001 9.25185853854297e-018
Multiplication: 6 * 1/6
1.00000000000000e+000 0.00000000000000e+000 <- should be this
1.00000000000000e+000 0.00000000000000e+000
6 * 1/3
2.00000000000000e+000 0.00000000000000e+000 <- should be this
2.00000000000000e+000 0.00000000000000e+000
Division: [1/3] / [1/6]
2.00000000000000e000 0.00000000000000e+000 <- should be this
2.00000000000000e+000 0.00000000000000e+000
[1/6] / [1/3]
5.00000000000000e-001 0.00000000000000e+000 <- should be this
5.00000000000000e-001 0.00000000000000e+000
Square root: sqrt[1/324]
5.55555555555556e-002 3.08395284618099e-018 <- should be this
5.55555555555556e-002 3.08395284618099e-018 <- sqrt[1/324] = 1/18
multiply by 18
1.00000000000000e+000 0.00000000000000e+000 <- should be this
1.00000000000000e+000 0.00000000000000e+000 ok
-- ----------------------------------------------------------------------
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : Doubled-precision floating point arithmetic
* CATEGORY : Tools
* AUTHOR : Marcel Hendrix
* AUTHOR : (c) Copyright 2006 Julian V. Noble.
* : Permission is granted by the author to use this software for
* : any application provided this copyright notice is preserved.
* LAST CHANGE : February 12, 2006, Marcel Hendrix
*)
NEEDS -miscutil
REVISION -dd-fp "ÄÄÄ 128-bit FP words Version 1.00 ÄÄÄ"
PRIVATES
DOC
(*
Environmental dependencies:
FPU must be set for IEEE 64-bit internal arithmetic
^^^^
Based on "Software for Doubled-Precision Floating-Point Computations", Seppo Linnainmaa,
ACM Transactions on Mathematical Software, Vol 7, No 3, September 1981, Pages 272-283
*)
ENDDOC
0 [IF]
-- Note: Be careful! Uses absolute addressing.
-- set 24-bit mode
: 24-bits! ( -- )
$00452C80 B@ $FCFF AND $00452C80 B!
1.001e FTRUNC FDROP ;
-- set 53-bit mode
: 53-bits! ( -- )
$00452C80 B@ $FCFF AND $200 OR $00452C80 B!
1.001e FTRUNC FDROP ;
-- set 80-bit mode
: 80-bits! ( -- )
$00452C80 B@ $FCFF AND $300 OR $00452C80 B!
1.001e FTRUNC FDROP ;
[THEN]
CR .( Warning: 53-bit FP mode set.) 53-bits!
4e 3e F/ F1- 3e F* F1- FVALUE w
w F2/ F1+ F1- FVALUE r
r F0<> [IF] r TO w [THEN]
2e 3e F/ 0.5e F- 3e F* 0.5e F- FVALUE uu
uu F2/ 0.5e F+ 0.5e F- FVALUE rr
rr F0<> [IF] rr TO uu [THEN]
w uu F/ ( F: -- beta)
uu FABS FLN FOVER FABS FLN F/ FNEGATE 0.5e F+
CR .( base = ) F>S .
CR .( precision = ) F>S .
FORGET w
: R128@ ( addr -- ) ( F: -- r1 r2 ) F2@ ;
: R128! ( addr -- ) ( F: r1 r2 -- ) F2! ;
-- Create a double-double variable
: DDOUBLE ( -- ) CREATE 2 DFLOATS ALLOT DOES> ;
-- Create a double-double constant
: DDCONSTANT
CREATE ( F: r1 r2 -- ) HERE >R 2 DFLOATS ALLOT R> r128!
DOES> ( F: -- r1 r2 ) r128@ ;
1e 0e DDCONSTANT dd=1
134217729e FCONSTANT split PRIVATE
-- Multiply 2 fp#'s to get ddfp#
: exactmul ( F: a b -- x xx )
0e 0e 0e 0e 0e 0e 0e FLOCALS| b1 b2 a1 a2 b21 b22 t bbb aaa |
bbb split F* TO t
bbb t F- t F+ TO a1
bbb a1 F- TO a2
aaa split F* TO t
aaa t F- t F+ TO b1
aaa b1 F- TO b2
b2 split F* TO t
b2 t F- t F+ TO b21
b2 b21 F- TO b22
bbb aaa F* FDUP TO t ( F: -- x )
a1 b1 F* t F-
a1 b2 F* F+
b1 a2 F* F+
b21 a2 F* F+
b22 a2 F* F+ ; PRIVATE
-- Divide 2 fp#'s to get ddfp#
: dd/ ( F: x xx y yy -- [x+xx]/[y+yy] )
0e 0e 0e 0e FLOCALS| q qq z2 z1 yy y xx x |
y F0= ABORT" dd/ :: division by 0"
x y F/ TO z1
z1 y exactmul TO qq TO q
x q F- qq F- xx F+ z1 yy F* F- y yy F+ F/ TO z2
z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ;
-- Multiply 2 fp#'s to get ddfp#
: dd* ( F: x xx y yy -- [x+xx]*[y+yy] )
0e 0e 0e FLOCALS| qq z2 z1 yy y xx x |
x y exactmul TO qq TO z1
x xx F+ yy F* xx y F* F+ qq F+ TO z2
z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ;
-- Add 2 fp#'s to get ddfp#
: dd+ ( F: x xx y yy -- [x+xx] + [y+yy] )
0e 0e 0e FLOCALS| q z2 z1 yy y xx x |
x y F+ TO z1
x z1 F- TO q
y q F+ ( F: -- q+y )
x q z1 F+ F- F+ ( F: -- q+y+x-[q+z1] )
xx F+ yy F+ TO z2
z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ;
-- Negate fp#
: ddnegate ( F: x xx -- -x -xx ) FNEGATE FSWAP FNEGATE FSWAP ;
-- Subtract 2 fp#'s to get ddfp#
: dd- ( F: x xx y yy -- [x+xx] - [y+yy] ) ddnegate dd+ ;
-- DDSQRT based on T.J. Dekker,
-- "A Floating-Point Technique for Extending the Available Precision"
-- Numerische Mathematik 18 (1971) 224-242.
: ddsqrt ( F: x xx -- ddsqrt[x+xx] )
0e 0e 0e 0e FLOCALS| qq q yy y xx x |
x F0< ABORT" ddsqrt :: sqrt of negative number"
x FSQRT TO q
q q exactmul TO yy TO y
x y F- yy F- xx F+
F2/ q F/ TO qq
qq q F+ FDUP FNEGATE q F+ qq F+ ;
-- ----------------------------------- stack ops -------
: ddswap ( F: x xx y yy -- y yy x xx ) F2SWAP ;
: dddup ( F: x xx -- x xx x xx ) F2OVER ;
: dddrop ( F: x xx -- ) F2DROP ;
: ddover ( F: x xx y yy -- x xx y yy x xx ) F2OVER ;
-- -----------------------------------------------------
: ZSQR ( re im -- r i )
F2DUP F* F2* -FROT
FSWAP FSQR FSWAP FSQR F-
FROT ; PRIVATE
: dd^2 ( F: x xx -- x2 xx2 ) dddup dd* ;
-- Raise dd to positive or negative integer power
-- Return 1 if n=0, dd^{-|n|} if n<0
: dd^n ( n -- ) ( F: x xx -- [x+xx]^n )
DUP 0= IF DROP dddrop dd=1 EXIT ENDIF
ZLOCAL x-xx
dd=1
DUP 0< SWAP ABS ( F: -- 1e 0e ) ( -- sign |n| )
BEGIN DUP 0>
WHILE DUP 1 AND IF x-xx dd* ENDIF zsqr
2/
REPEAT DROP
IF dd=1 ddswap dd/ ENDIF ;
nesting @ 1 =
[IF]
1e 0e 3e 0e dd/ DDCONSTANT dd1/3
1e 0e 6e 0e dd/ DDCONSTANT dd1/6
: TEST ( -- )
PRECISION LOCAL old-prec #15 SET-PRECISION
dd1/3 dd1/6 dd+ FSWAP
CR CR ." Addition: 1/3 + 1/6"
CR 6 spaces ." 5.00000000000000e-001 -1.54074395550979e-033 <- should be this"
CR 5 spaces +E. SPACE +E.
CR
dd1/3 dd1/6 dd- FSWAP
CR ." Subtraction: 1/3 - 1/6"
CR 6 spaces ." 1.66666666666667e-001 9.25185853854297e-018 <- should be this"
CR 5 spaces +E. SPACE +E.
CR
6e 0e dd1/6 dd* FSWAP
CR ." Multiplication: 6 * 1/6"
CR 6 spaces ." 1.00000000000000e+000 0.00000000000000e+000 <- should be this"
CR 5 spaces +E. SPACE +E.
CR ." 6 * 1/3"
6e 0e dd1/3 dd* FSWAP
CR 6 spaces ." 2.00000000000000e+000 0.00000000000000e+000 <- should be this"
CR 5 spaces +E. SPACE +E.
CR
dd1/3 dd1/6 dd/ FSWAP
CR ." Division: [1/3] / [1/6]"
CR 6 spaces ." 2.00000000000000e000 0.00000000000000e+000 <- should be this"
CR 5 spaces +E. SPACE +E.
CR
CR ." [1/6] / [1/3]"
dd1/6 dd1/3 dd/ FSWAP
CR 6 spaces ." 5.00000000000000e-001 0.00000000000000e+000 <- should be this"
CR 5 spaces +E. SPACE +E.
CR ." Square root: sqrt[1/324]"
dd1/6 dd1/3 dd* F2DUP dd* ddsqrt F2DUP FSWAP
CR 6 spaces ." 5.55555555555556e-002 3.08395284618099e-018 <- should be this"
CR 5 spaces +E. SPACE +E. ." <- sqrt[1/324] = 1/18"
CR 5 spaces ." multiply by 18"
18e 0e dd* FSWAP
CR 6 spaces ." 1.00000000000000e+000 0.00000000000000e+000 <- should be this"
CR 5 spaces +E. SPACE +E.
old-prec SET-PRECISION ;
[THEN]
:ABOUT CR ." ** DD-FP words *** "
CR ." Try: TEST" ;
.ABOUT -dd-fp CR
DEPRIVE
(* End of Source *)
-- ------------------
FORTH> see exactmul
Flags: ANSI, PRIVATE
$0052BB80 : exactmul
$0052BB88 fld $004BA388 qword-ptr
$0052BB8E fld $004BA380 qword-ptr
$0052BB94 fld $004BA378 qword-ptr
$0052BB9A fld $004BA370 qword-ptr
$0052BBA0 fld $004BA368 qword-ptr
$0052BBA6 fld $004BA360 qword-ptr
$0052BBAC mov eax, $00452C74 dword-ptr
$0052BBB1 lea eax, [eax #-32 +] dword
$0052BBB4 mov $00452C74 dword-ptr, eax
$0052BBB9 fxch ST(4)
$0052BBBB fstp [eax 8 +] qword
$0052BBBE fxch ST(4)
$0052BBC0 fstp [eax #16 +] qword
$0052BBC3 fxch ST(4)
$0052BBC5 fstp [eax #24 +] qword
$0052BBC8 fstp [eax 0 +] qword
$0052BBCA fld $004BA358 qword-ptr
$0052BBD0 sub esi, 8 b#
$0052BBD3 fstp [esi 0 +] qword
$0052BBD5 sub esi, 8 b#
$0052BBD8 fstp [esi 0 +] qword
$0052BBDA sub esi, 8 b#
$0052BBDD fstp [esi 0 +] qword
$0052BBDF fpop,
$0052BBF0 sub esi, 8 b#
$0052BBF3 fstp [esi 0 +] qword
$0052BBF5 fpop,
$0052BC06 sub esi, 8 b#
$0052BC09 fstp [esi 0 +] qword
$0052BC0B fpop,
$0052BC1C sub esi, 8 b#
$0052BC1F fstp [esi 0 +] qword
$0052BC21 fpop,
$0052BC32 sub esi, 8 b#
$0052BC35 fstp [esi 0 +] qword
$0052BC37 fpop,
$0052BC48 sub esi, 8 b#
$0052BC4B fstp [esi 0 +] qword
$0052BC4D fpop,
$0052BC5E sub esi, 8 b#
$0052BC61 fstp [esi 0 +] qword
$0052BC63 fld [esi 8 +] qword
$0052BC66 fmul $004BA350 qword-ptr
$0052BC6C fstp [esi #16 +] qword
$0052BC6F fld [esi 8 +] qword
$0052BC72 fsub [esi #16 +] qword
$0052BC75 fadd [esi #16 +] qword
$0052BC78 fstp [esi #48 +] qword
$0052BC7B fld [esi 8 +] qword
$0052BC7E fsub [esi #48 +] qword
$0052BC81 fstp [esi #40 +] qword
$0052BC84 fld [esi 0 +] qword
$0052BC86 fmul $004BA348 qword-ptr
$0052BC8C fstp [esi #16 +] qword
$0052BC8F fld [esi 0 +] qword
$0052BC91 fsub [esi #16 +] qword
$0052BC94 fadd [esi #16 +] qword
$0052BC97 fstp [esi #64 +] qword
$0052BC9A fld [esi 0 +] qword
$0052BC9C fsub [esi #64 +] qword
$0052BC9F fst [esi #56 +] qword
$0052BCA2 fmul $004BA340 qword-ptr
$0052BCA8 fstp [esi #16 +] qword
$0052BCAB fld [esi #56 +] qword
$0052BCAE fsub [esi #16 +] qword
$0052BCB1 fadd [esi #16 +] qword
$0052BCB4 fstp [esi #32 +] qword
$0052BCB7 fld [esi #56 +] qword
$0052BCBA fsub [esi #32 +] qword
$0052BCBD fstp [esi #24 +] qword
$0052BCC0 fld [esi 8 +] qword
$0052BCC3 fmul [esi 0 +] qword
$0052BCC5 fld ST(0)
$0052BCC7 fstp [esi #16 +] qword
$0052BCCA fld [esi #48 +] qword
$0052BCCD fmul [esi #64 +] qword
$0052BCD0 fsub [esi #16 +] qword
$0052BCD3 fld [esi #48 +] qword
$0052BCD6 fmul [esi #56 +] qword
$0052BCD9 faddp ST(1), ST
$0052BCDB fld [esi #64 +] qword
$0052BCDE fmul [esi #40 +] qword
$0052BCE1 faddp ST(1), ST
$0052BCE3 fld [esi #32 +] qword
$0052BCE6 fmul [esi #40 +] qword
$0052BCE9 faddp ST(1), ST
$0052BCEB fld [esi #24 +] qword
$0052BCEE fmul [esi #40 +] qword
$0052BCF1 faddp ST(1), ST
$0052BCF3 mov eax, $00452C74 dword-ptr
$0052BCF8 lea eax, [eax #-16 +] dword
$0052BCFB mov $00452C74 dword-ptr, eax
$0052BD00 fxch ST(2)
$0052BD02 fstp [eax 8 +] qword
$0052BD05 fstp [eax 0 +] qword
$0052BD07 add esi, #72 b#
$0052BD0A ;
Marcel Hendrix
> "Julian V. Noble" <j...@virginia.edu> wrote Re: double-double precision arithmetic in Forth
[..]
>> I will be translating more things as I get time, such as a DD function
>> library, DD I/O words, etc. I also intend to translate them to assembler,
>> which should make them a bit less ugly.
> I think they WILL look better in high-level, see attachment :-)
In Win32Forth 4.20671 at least, some assembler will be needed to make
the code usable:
: BENCH ( -- )
CR ." 10,000,000 executions of dd* : " TIMER-RESET
10000000 0 DO dd=1 dd=1 dd* F2DROP LOOP .ELAPSED
CR ." 10,000,000 executions of dd/ : " TIMER-RESET
10000000 0 DO dd=1 dd=1 dd* F2DROP LOOP .ELAPSED
CR ." 10,000,000 executions of dd+ : " TIMER-RESET
10000000 0 DO dd=1 dd=1 dd+ F2DROP LOOP .ELAPSED
CR ." 10,000,000 executions of dd- : " TIMER-RESET
10000000 0 DO dd=1 dd=1 dd- F2DROP LOOP .ELAPSED
CR ." 10,000,000 executions of ddsqrt : " TIMER-RESET
10000000 0 DO dd=1 ddsqrt F2DROP LOOP .ELAPSED ;
Win32Forth 4.2.0671 on 3 GHz PIV
10,000,000 executions of dd* : Elapsed time: 00:00:20.454
10,000,000 executions of dd/ : Elapsed time: 00:00:20.453
10,000,000 executions of dd+ : Elapsed time: 00:00:09.562
10,000,000 executions of dd- : Elapsed time: 00:00:11.453
10,000,000 executions of ddsqrt : Elapsed time: 00:00:21.094 ok
MPE VFX Forth for Windows IA32 ( same machine)
© MicroProcessor Engineering Ltd, 1998-2005
Version: 3.80 [build 1937]
Build date: 29 September 2005
bench ( rebuild for 64bit float, cw@ $300 NOT AND $200 OR cw! ( 53-bit mode))
10,000,000 executions of dd* : 1047 ms elapsed.
10,000,000 executions of dd/ : 1062 ms elapsed.
10,000,000 executions of dd+ : 391 ms elapsed.
10,000,000 executions of dd- : 391 ms elapsed.
10,000,000 executions of ddsqrt : 1625 ms elapsed. ok
FORTH> ( iForth 2.0, same machine ) bench
10,000,000 executions of dd* : 0.512 seconds elapsed.
10,000,000 executions of dd/ : 0.512 seconds elapsed.
10,000,000 executions of dd+ : 0.244 seconds elapsed.
10,000,000 executions of dd- : 0.259 seconds elapsed.
10,000,000 executions of ddsqrt : 0.674 seconds elapsed. ok
-marcel
Krishna Myneni
> Colleagues,
>
> Here are portable double-double precision floating point arithmetic words.
> I've tested them on Win32Forth 4.2 and they work fine. I regret they look
> a tad ugly. It was necessary for portability. Refinements would involve
> making the (many) fvariables (assumed to be 64-bit precision) local to the
> routines that use them to avoid name collisions.
>
> I have not spent any time on factoring or refinement.
>
> The words will not work as-is on GForth because I believe (from my experiments)
> that the latter is default-set to 80-bit internal precision rather than
> 64-bit. I haven't had a chance to work up a CODE word for doing this setting.
>
> I will be translating more things as I get time, such as a DD function
> library, DD I/O words, etc. I also intend to translate them to assembler,
> which should make them a bit less ugly.
>
Your test code produces the same output as the "correct" answers on
kForth and GForth, both of which default to 80-bit internal precision
in the fpu. However, after each arithmetic operation, the fp number on the
stack/fp-stack is truncated to 64-bits, so perhaps the internal precision
of the fpu is not an issue, at least when the dd arithmetic words are
programmed using high-level Forth words.
For some reason, the precision calculation appears to be messed up.
Also, FVARIABLES is defined twice in the file, and the word DD^N references
an undefined word Z^2. The defn of DD^N was commented out for the tests below.
Under kForth, the last significant digit is different only because of the
output routine truncating rather than rounding.
Krishna
--
Under kForth:
-------------
krishna@ungee:~/kforth> kforth double-double
FVARIABLES: is redefined
base = 2 precision = -2147483648
Addition: 1/3 + 1/6
5.00000000000000E-1 -1.54074395550979E-33 <- should be this
5.00000000000000e-1-1.54074395550978e-33
Subtraction: 1/3 - 1/6
1.66666666666667E-1 9.25185853854297E-18 <- should be this
1.66666666666666e-1 9.25185853854296e-18
Multiplication: 6 * 1/6
1.00000000000000E0 0.00000000000000E-1 <- should be this
1.00000000000000e0 0.00000000000000e0
6 * 1/3
2.00000000000000E0 0.00000000000000E-1 <- should be this
2.00000000000000e0 0.00000000000000e0
Division: [1/3] / [1/6]
2.00000000000000E0 0.00000000000000E-1 <- should be this
2.00000000000000e0 0.00000000000000e0
[1/6] / [1/3]
5.00000000000000E-1 0.00000000000000E-1 <- should be this
5.00000000000000e-1 0.00000000000000e0
Square root: sqrt[1/324]
5.55555555555555e-2 3.08395284618099e-18<- sqrt[1/324] = 1/18
5.55555555555556E-2 3.08395284618099E-18 <- should be this
multiply by 18
1.00000000000000E0 0.00000000000000E-1 <- should be this
1.00000000000000e0 0.00000000000000e0 ok
Under Gforth:
-------------
krishna@ungee:~/gforth> gforth
Gforth 0.6.2, Copyright (C) 1995-2003 Free Software Foundation, Inc.
Gforth comes with ABSOLUTELY NO WARRANTY; for details type `license'
Type `bye' to exit
include double-double.fs redefined fvariables:
base = 2 precision = -9223372036854775808
Addition: 1/3 + 1/6
5.00000000000000E-1 -1.54074395550979E-33 <- should be this
5.00000000000000E-1 -1.54074395550979E-33
Subtraction: 1/3 - 1/6
1.66666666666667E-1 9.25185853854297E-18 <- should be this
1.66666666666667E-1 9.25185853854297E-18
Multiplication: 6 * 1/6
1.00000000000000E0 0.00000000000000E-1 <- should be this
1.00000000000000E0 0.00000000000000E0
6 * 1/3
2.00000000000000E0 0.00000000000000E-1 <- should be this
2.00000000000000E0 0.00000000000000E0
Division: [1/3] / [1/6]
2.00000000000000E0 0.00000000000000E-1 <- should be this
2.00000000000000E0 0.00000000000000E0
[1/6] / [1/3]
5.00000000000000E-1 0.00000000000000E-1 <- should be this
5.00000000000000E-1 0.00000000000000E0
Square root: sqrt[1/324]
5.55555555555556E-2 3.08395284618099E-18 <- sqrt[1/324] = 1/18
5.55555555555556E-2 3.08395284618099E-18 <- should be this
multiply by 18
1.00000000000000E0 0.00000000000000E-1 <- should be this
1.00000000000000E0 0.00000000000000E0 ok
Anton Ertl
>The words will not work as-is on GForth because I believe (from my experiments)
>that the latter is default-set to 80-bit internal precision rather than
>64-bit.
That should not be a problem on platforms that use plain 64-bit
numbers (i.e., almost everything but the 386 platform).
So I cut and pasted the program from this posting, and found that it
contained a word "((" that is unknown to Gforth, and I have no idea
what that's supposed to do. So I saw the following:
>FALSE [IF] \ use these to compile the versions using
>include ftran202.f \ the FORmula TRANslator and Baden's flocals
>include flocals.f
>[THEN]
and hoped that "((" is defined in there. So I tried to get hold of
ftran202.f; I found an HTML page with Google, but not the file
directly, so I had to cut and paste it into a file. Now that file
includes some more, so I had to go through the same procedure with
chr_tbl.f
complex.f
fsm2.f
ftran202.f
vector1.f
Unfortunately, I did not find flocals.f, and "((" is still undefined
with all the files I found included, so all this work is for naught.
Hmm, as long as our Forth systems cannot search for and download the
files automatically (which would be a bad idea for security reasons),
it would be nice if people packaged their programs with all the
additional needed stuff into one file (e.g., in zip format).
- 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.complang.tuwien.ac.at/forth/ansforth/forth200x.html
Krishna Myneni
Attached is my revised version for Gforth.
KM
--
\ Doubled-precision floating point arithmetic
\
\ ---------------------------------------------------
\ (c) Copyright 2006 Julian V. Noble. \
\ Permission is granted by the author to \
\ use this software for any application pro- \
\ vided this copyright notice is preserved. \
\ ---------------------------------------------------
\ This is an ANS Forth program requiring the
\ FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependences:
\ Assumes independent floating point stack
\ FPU must be set for IEEE 64-bit internal arithmetic
\ ^^^^
\ based on "Software for Doubled-Precision Floating-Point Computations"
\ by Seppo Linnainmaa
\ ACM Transactions on Mathematical Software,
\ Vol 7, No 3, September 1981, Pages 272-283
marker -double
FALSE [IF] \ use these to compile the versions using
include ftran202.f \ the FORmula TRANslator and Baden's flocals
include flocals.f
[THEN]
: fvariables: 0 DO FVARIABLE LOOP ;
[undefined] fvalue [IF] : fvalue CREATE HERE F! 1 DFLOATS allot DOES> F@ ; [THEN]
[undefined] fto [IF] : fto ' >BODY STATE @ IF POSTPONE FLITERAL POSTPONE F!
ELSE F! THEN ; IMMEDIATE [THEN]
[undefined] f-rot [IF] : f-rot FROT FROT ; [THEN]
[undefined] ftuck [IF] : ftuck FSWAP FOVER ; [THEN]
[undefined] f>s [IF] : f>s f>d d>s ; [THEN]
[undefined] s>f [IF] : s>f s>d d>f ; [THEN]
\ ---------------------------------------- LOAD, STORE
: r128@ DUP F@ FLOAT+ F@ ;
: r128! DUP FLOAT+ F! F! ;
\ ------------------------------------ END LOAD, STORE
\ ----------------------------------- data types ----
: real*16 \ create a double-double variable
CREATE 2 FLOATS ALLOT ;
: ddconstant CREATE HERE r128! 2 FLOATS ALLOT
DOES> r128@ ;
\ ---------------------------------------------------
3.1415926535897931e0 1.2246467991473532e-16 ddconstant ddpi
\ based on "Software for Doubled-Precision Floating-Point Computations"
\ by Seppo Linnainmaa
\ ACM Transactions on Mathematical Software,
\ Vol 7, No 3, September 1981, Pages 272-283
: fvariables: 0 DO FVARIABLE LOOP ;
\ determine base and precision of fpu
\ ((
4e0 3e0 F/ 1e0 F- 3e0 F* 1e0 F- FVALUE u
u 2e0 F/ 1e0 F+ 1e0 F- FVALUE r
r F0= 0= [IF] r FTO u [THEN]
2e0 3e0 F/ 0.5e0 F- 3e0 F* 0.5e0 F- FVALUE uu
uu 2e0 F/ 0.5e0 F+ 0.5e0 F- FVALUE rr
rr F0= 0= [IF] rr FTO uu [THEN]
Krishna Myneni
> ... I will be translating more things as I get time, such as a DD function
> library, DD I/O words, etc. ...
A word such as DD. to output the 128-bit fp number to the appropriate number of
significant digits would be immediately useful for testing.
Krishna
Anton Ertl
>Attached is my revised version for Gforth.
Thanks.
I ran it on an Alpha, an AMD64 box, and a 386 box, and on all three it
produced the same results, which look to be equal to "should be"
output. So apparently the test cases don't include a case where the
double rounding makes a difference.
You might wonder what double rounding is: Suppose you have two number
formats, one with a three-bit mantissa and one with a five-bit
mantissa. Now consider the case where the true result of an operation
is 101011. Now the correct rounding to three bits results in 101 and
the correct rounding to five bits results in 10110. But if you round
first to five bits (getting 10110), and then round that to three bits,
you get 110, which is wrong. The 80387 FPU does that with a 53-bit
and a 64-bit mantissa if you use the 64-bit mantissa on the FPU stack,
and then store the result in memory as a 64-bit double (with a 53-bit
mantissa), as probably happens in Gforth when it uses the 387.
Julian V. Noble
Dear me.
: ftuck FSWAP FOVER ;
--
Julian V. Noble
Professor Emeritus of Physics
j...@lessspamformother.virginia.edu
^^^^^^^^^^^^^^^^^^
Julian V. Noble
>
> m...@iae.nl (Marcel Hendrix) writes Re: double-double precision arithmetic in Forth
>
> > "Julian V. Noble" <j...@virginia.edu> wrote Re: double-double precision arithmetic in Forth
> [..]
> >> I will be translating more things as I get time, such as a DD function
> >> library, DD I/O words, etc. I also intend to translate them to assembler,
> >> which should make them a bit less ugly.
>
> > I think they WILL look better in high-level, see attachment :-)
>
> In Win32Forth 4.20671 at least, some assembler will be needed to make
> the code usable:
Of course. I was well aware that high-level would be slow. However, it
_is_ portable, which was my first concern.
Julian V. Noble
"((" is simply a multi-line "(" --meaning "ignore everything up to )) .
I would have written FALSE [IF] ... [THEN] if the code re[resented by
... had not contained some [THEN] 's . I see from dpANS94 that interpretive
[IF] ... [THEN] statements do permit nested levels, so an external wrapper
would have been ok.
Sorry about the flocals.f file, Anton -- here it is:
FALSE [IF]
Code for local fvariables, loosely based upon Wil Baden's idea presented
at FORML 1992.
The idea is to have a fixed number of variables with fixed names.
I believe the code shown here will work with any (case insensitive)
ANS Forth.
(i/t)Forth users are advised to use FLOCALS| instead.
The code as shown doesn't work in words that modify HERE (like CREATE)
example: : test 2e 3e FRAME| aa bb | aa f. bb f. |FRAME ;
test <cr> 3.0000 2.0000 ok
PS: Don't forget to use |FRAME before an EXIT .
[THEN]
8 constant /flocals
: (frame) ( n -- ) FLOATS ALLOT ;
: FRAME|
0 >R
BEGIN BL WORD COUNT 1 =
SWAP C@ [CHAR] | =
AND 0=
WHILE POSTPONE F, R> 1+ >R
REPEAT
/flocals R> - DUP 0< ABORT" too many flocals"
POSTPONE LITERAL POSTPONE (frame) ; IMMEDIATE
: |FRAME ( -- ) [ /flocals negate ] literal (frame) ;
: &h HERE [ 1 FLOATS ] LITERAL - ;
: &g HERE [ 2 FLOATS ] LITERAL - ;
: &f HERE [ 3 FLOATS ] LITERAL - ;
: &e HERE [ 4 FLOATS ] LITERAL - ;
: &d HERE [ 5 FLOATS ] LITERAL - ;
: &c HERE [ 6 FLOATS ] LITERAL - ;
: &b HERE [ 7 FLOATS ] LITERAL - ;
: &a HERE [ 8 FLOATS ] LITERAL - ;
: aa &a ;
: bb &b ;
: cc &c ;
: dd &d ;
: ee &e ;
: ff &f ;
: gg &g ;
: hh &h ;
Julian V. Noble
>
> Julian V. Noble wrote:
> > Colleagues,
> >
> > Here are portable double-double precision floating point arithmetic words.
> > I've tested them on Win32Forth 4.2 and they work fine. I regret they look
> > a tad ugly. It was necessary for portability. Refinements would involve
> > making the (many) fvariables (assumed to be 64-bit precision) local to the
> > routines that use them to avoid name collisions.
> >
> > I have not spent any time on factoring or refinement.
> >
> > The words will not work as-is on GForth because I believe (from my experiments)
> > that the latter is default-set to 80-bit internal precision rather than
> > 64-bit. I haven't had a chance to work up a CODE word for doing this setting.
> >
> > I will be translating more things as I get time, such as a DD function
> > library, DD I/O words, etc. I also intend to translate them to assembler,
> > which should make them a bit less ugly.
> >
>
> Your test code produces the same output as the "correct" answers on
> kForth and GForth, both of which default to 80-bit internal precision
> in the fpu. However, after each arithmetic operation, the fp number on the
> stack/fp-stack is truncated to 64-bits, so perhaps the internal precision
> of the fpu is not an issue, at least when the dd arithmetic words are
> programmed using high-level Forth words.
>
> For some reason, the precision calculation appears to be messed up.
> Also, FVARIABLES is defined twice in the file, and the word DD^N references
> an undefined word Z^2. The defn of DD^N was commented out for the tests below.
>
The word z^2 should have been dd^2, which _was_ defined. I didn't notice the
typo because z^2 is defined (square a complex #) on my system.
> base = 2 precision = -2147483648
for 64-bit internals on the Intel/AMD fpu's this should be
base = 2 precision = 54
> Under Gforth:
> -------------
>
> krishna@ungee:~/gforth> gforth
> Gforth 0.6.2, Copyright (C) 1995-2003 Free Software Foundation, Inc.
> Gforth comes with ABSOLUTELY NO WARRANTY; for details type `license'
> Type `bye' to exit
> include double-double.fs redefined fvariables:
>
> base = 2 precision = -9223372036854775808
Same comment: b=2, p = 54 .
Julian V. Noble
I am in fact constructing that word right now. Patience, amigos! It's
just a little tricky and I want to make it a little better-factored than
some of the comparable routines I have seen in various languages--
including many variants of Forth.
Krishna Myneni
> Krishna Myneni <krishn...@bellsouth.net> writes:
>
>>Attached is my revised version for Gforth.
>
>
> Thanks.
>
> I ran it on an Alpha, an AMD64 box, and a 386 box, and on all three it
> produced the same results, which look to be equal to "should be"
> output. So apparently the test cases don't include a case where the
> double rounding makes a difference.
>
> You might wonder what double rounding is: Suppose you have two number
> formats, one with a three-bit mantissa and one with a five-bit
> mantissa. Now consider the case where the true result of an operation
> is 101011. Now the correct rounding to three bits results in 101 and
> the correct rounding to five bits results in 10110. But if you round
> first to five bits (getting 10110), and then round that to three bits,
> you get 110, which is wrong. The 80387 FPU does that with a 53-bit
> and a 64-bit mantissa if you use the 64-bit mantissa on the FPU stack,
> and then store the result in memory as a 64-bit double (with a 53-bit
> mantissa), as probably happens in Gforth when it uses the 387.
>
> - anton
In order for double rounding to propagate an error from a 64-bit mantissa to a
53-bit mantissa, don't the intermediate bits have to be all 1's? Or perhaps I am
not seeing it properly?
Krishna
Krishna Myneni
> Colleagues,
>
> Here are portable double-double precision floating point arithmetic words.
> I've tested them on Win32Forth 4.2 and they work fine. I regret they look
> a tad ugly. It was necessary for portability. Refinements would involve
> making the (many) fvariables (assumed to be 64-bit precision) local to the
> routines that use them to avoid name collisions.
>
> I have not spent any time on factoring or refinement.
>
> The words will not work as-is on GForth because I believe (from my experiments)
> that the latter is default-set to 80-bit internal precision rather than
> 64-bit. I haven't had a chance to work up a CODE word for doing this setting.
Attached below is Forth code for setting up the x86 FPU control word. This will
allow setting 64-bit precision, as shown in the example below.
Krishna
--
krishna@ungee:~/kforth> kforth
kForth v 1.3.1 (Build: 2005-09-28)
Copyright (c) 1998--2005 Krishna Myneni
Contributions by: dpw gd mu bk abs tn cmb bg
Provided under the GNU General Public License.
ok
include asm-x86
# is redefined
ok
include fpu-x86
ok
getFPUStateX86
ok
fpu-control @ hex .
37F ok
fpu-control @ FPU_CW_PREC_MASK and .
300 ok
fpu-control @ FPU_CW_ROUND_MASK and .
0 ok
ok
FPU_CW_PREC_DOUBLE FPU_CW_PREC_MASK modifyFPUStateX86
ok
getFPUStateX86 fpu-control @ FPU_CW_PREC_MASK and .
200 ok
----------------------------------------------------------------------------
\ fpu-x86.4th
\
\ Precision control of the x86 Floating Point Unit
\
\ Derived from C code by Kevin Egan, Brown University.
\ This code is published at:
\
\ http://www.stereopsis.com/FPU.html
\
\ ----------------------------------------------------
\ bits to set the floating point control word register
\
\ Sections 4.9, 8.1.4, 10.2.2 and 11.5 in
\ IA-32 Intel Architecture Software Developer's Manual
\ Volume 1: Basic Architecture
\
\ http://www.intel.com/design/pentium4/manuals/245471.htm
\
\ http://www.geisswerks.com/ryan/FAQS/fpu.html
\
\ precision control:
\ 00 : single precision
\ 01 : reserved
\ 10 : double precision
\ 11 : extended precision
\
\ rounding control:
\ 00 = Round to nearest whole number. (default)
\ 01 = Round down, toward -infinity.
\ 10 = Round up, toward +infinity.
\ 11 = Round toward zero (truncate).
BASE @
HEX
003f constant FPU_CW_EXCEPTION_MASK
0001 constant FPU_CW_INVALID
0002 constant FPU_CW_DENORMAL
0004 constant FPU_CW_ZERODIVIDE
0008 constant FPU_CW_OVERFLOW
0010 constant FPU_CW_UNDERFLOW
0020 constant FPU_CW_INEXACT
0300 constant FPU_CW_PREC_MASK
0000 constant FPU_CW_PREC_SINGLE
0200 constant FPU_CW_PREC_DOUBLE
0300 constant FPU_CW_PREC_EXTENDED
0c00 constant FPU_CW_ROUND_MASK
0000 constant FPU_CW_ROUND_NEAR
0400 constant FPU_CW_ROUND_DOWN
0800 constant FPU_CW_ROUND_UP
0c00 constant FPU_CW_ROUND_CHOP
1f3f constant FPU_CW_MASK_ALL
\ --------------------------------------------------------
variable fpu-control
\ The following CODE words are in kForth's asm-x86 style
\ Modify as needed for your Forth system.
CODE getFPUStateX86
fpu-control #@ fnstcw,
END-CODE
CODE setFPUStateX86
fpu-control #@ fldcw,
END-CODE
\ --------------------------------------------------------
\ Modify the control bits of a given setting, e.g.
\
\ FPU_CW_PREC_DOUBLE FPU_CW_PREC_MASK modifyFPUStateX86
\
\ sets double precision mode.
: modifyFPUStateX86 ( control mask -- )
dup >r and
getFPUStateX86
fpu-control @ r> not and or fpu-control !
setFPUStateX86
;
BASE !
Anton Ertl
>Anton Ertl wrote:
>> You might wonder what double rounding is: Suppose you have two number
>> formats, one with a three-bit mantissa and one with a five-bit
>> mantissa. Now consider the case where the true result of an operation
>> is 101011. Now the correct rounding to three bits results in 101 and
>> the correct rounding to five bits results in 10110. But if you round
>> first to five bits (getting 10110), and then round that to three bits,
>> you get 110, which is wrong. The 80387 FPU does that with a 53-bit
>> and a 64-bit mantissa if you use the 64-bit mantissa on the FPU stack,
>> and then store the result in memory as a 64-bit double (with a 53-bit
>> mantissa), as probably happens in Gforth when it uses the 387.
>>
>> - anton
>
>In order for double rounding to propagate an error from a 64-bit mantissa to a
>53-bit mantissa, don't the intermediate bits have to be all 1's?
Not quite. I am no expert at this, but here's the cases that I know
of, where double rounding differs from correct rounding:
Notation: The "|" character indicates the end of the mantissa for a
given precision; since we have two precisions involved, there are two
"|" in our numbers. "1*" and "0*" means a sequence of 1s or 0s,
respectively. Ok, here are the problematic cases (in connection with
rounding to nearest or even):
rounding direction
exact value correct double
...1|01*|1*0... down up
...0|10*|0*1... up down
So the cases where it makes a difference are relatively rare.
Marcel Hendrix
> "Julian V. Noble" <j...@virginia.edu> wrote Re: double-double precision arithmetic in Forth
>> I will be translating more things as I get time, such as a DD function
>> library, DD I/O words, etc. I also intend to translate them to assembler,
>> which should make them a bit less ugly.
[..]
> It will be nice to see what you come up with.
As nothing seems to be forthcoming, I've decided to put my own function
library on http://home.iae.nl/users/mhx/programs.html .
Unfortunately, I needed/decided to use quite a lot of iForth extensions.
At least, the result is *a lot* shorter than David H. Bailey's library
(where it is based upon), http://crd.lbl.gov/~dhbailey/ .
There are lots of tests in the package, but they are not really
satisfactory. Unfortunately, I don't see how to torture this program
without use of foreign tools (like a symbolic math package).
-marcel
Julian V. Noble
Thank you Marcel, my health has not been good and I have been
concentrating on my book. I did the double-double stuff at a
moment when the book was at an impasse.
If you want the assembler versions of the key routines, please
let me know. Also I have a version of ftran that compiles double-
double expressions. Do you want that?
--
Julian V. Noble
Professor Emeritus of Physics
University of Virginia