groff
[Top][All Lists]

## Re: [Groff] About numerical precision

 From: Werner LEMBERG Subject: Re: [Groff] About numerical precision Date: Thu, 03 Apr 2003 01:35:08 +0200 (CEST)

```> OK, found a solution by computing separately the integer and the
> fractional part of the division.

Below is my solution, using exact math.  Just for fun :-)

Werner

======================================================================

.\"
.\" Add a 30bit number to a 60bit number.  Result is a 60bit number:
.\"
.\"   <x> + (<y>h * 2^30 + <y>l) = <z>h * 2^30 + <z>l
.\"
.\"
.\" in: \n[<x>], \n[<y>h], \n[<y>l]
.\"
.\" out: \n[<z>h], \n[<z>l]
.\"
.\" Example: .add30to60 p q r
.\"
.\"          -> input registers: \n[p], \n[qh], \n[ql]
.\"             output registers: \n[rh], \n[rl]
.
.  nr math-i (\\n[\\\$1] >? \\n[\\\$2l])
.  nr \\\$3l (\\n[\\\$1] + \\n[\\\$2l] % 1073741824)
.  nr \\\$3h \\n[\\\$2h]
.  if (\\n[\\\$3l] < \\n[math-i]) \
.    nr \\\$3h +1
..
.
.
.\" .mult60 <x> <y> <z>
.\"
.\" Multiply two 30bit numbers.  Result is a 60bit number:
.\"
.\"   <x> * <y> = <z>h * 2^30 + <z>l
.\"
.\"
.\" in: \n[<x>], \n[<y>]
.\"
.\" out: \n[<z>h], \n[<z>l]
.\"
.\" Example: .mult60 a b c
.\"
.\"          -> input registers: \n[a], \n[b]
.\"             output registers: \n[ch], \n[cl]
.
.
.de1 mult60
.  nr math-lo1 (\\n[\\\$1] % 32768)
.  nr math-hi1 (\\n[\\\$1] / 32768)
.  nr math-lo2 (\\n[\\\$2] % 32768)
.  nr math-hi2 (\\n[\\\$2] / 32768)
.
.  nr \\\$3l (\\n[math-lo1] * \\n[math-lo2] % 1073741824)
.  nr math-i1 (\\n[math-lo1] * \\n[math-hi2] % 1073741824)
.  nr math-i2 (\\n[math-lo2] * \\n[math-hi1] % 1073741824)
.  nr \\\$3h (\\n[math-hi1] * \\n[math-hi2] % 1073741824)
.
.  nr math-i1 (\\n[math-i1] + \\n[math-i2] % 1073741824)
.  \" check carry overflow of math-i1 + math-i2
.  if (\\n[math-i1] < \\n[math-i2]) \
.    nr \\\$3h +32768
.
.  nr \\\$3h +(\\n[math-i1] / 32768)
.  \" multiply by 32768 in small steps to avoid overflow
.  nr math-i 15
.  while (\\n[math-i]) \{\
.    nr math-i1 (\\n[math-i1] * 2 % 1073741824)
.    nr math-i -1
.  \}
.
.  nr \\\$3l (\\n[\\\$3l] + \\n[math-i1] % 1073741824)
.  \" check carry overflow of math-i1 + lo
.  if (\\n[\\\$3l] < \\n[math-i1]) \
.    nr \\\$3h +1
..
.
.
.\" .div60by30 <x> <y> <z>
.\"
.\" Divide a 60bit number by a 30bit.  Result is a 30bit number:
.\"
.\"   (<x>h * 2^30 + <x>l) / <y> = <z>
.\"
.\"
.\" in: \n[<x>h], \n[<x>l], \n[<y>]
.\"
.\" out: \n[<z>]
.\"
.\" Example: .div60by30 foo bar baz
.\"
.\"          -> input registers: \n[fooh] \n[fool] \n[bar]
.\"             output register: \n[baz]
.
.de1 div60by30
.  nr \\\$3 0
.  nr math-i 30
.
.  while (\\n[math-i]) \{\
.    nr \\\$1h (\\n[\\\$1h] * 2 % 1073741824)
.    nr \\\$3 (\\n[\\\$3] * 2)
.    nr \\\$1h +(\\n[\\\$1l] / 536870912)
.
.    if (\\n[\\\$1h] >= \\n[\\\$2]) \{\
.      nr \\\$1h -\\n[\\\$2]
.      nr \\\$3 +1
.    \}
.    nr \\\$1l (\\n[\\\$1l] * 2 % 1073741824)
.    nr math-i -1
.  \}
..
.
.\" .SMALLCAPS <first> <small> <rest>
.
.de1 SMALLCAPS
.  nr x (0*\w'x' + \\n[rst])
.  nr X (0*\w'X' + \\n[rst])
.  nr s ((9 * \\n[.ps] + 5) / 10)
.  nr .ps/2 (\\n[.ps] / 2)
.  \" (X - x)/3 + x == (2X + 4x)/6
.  nr sX (((2 * \\n[X]) + (4 * \\n[x]) + 3) / 6)
.
.  mult60 sX .ps temp1
.  div60by30 temp2 X H
.
\\\$1\
\h'(\w'\\\$1\\\$2'u - \w'\\\$1'u - \w'\\\$2'u)'\
\H'\\n[H]u'\
\s[\\n[s]u]\
\\\$2\
\h'(\w'\\\$2\\\$3'u - \w'\\\$2'u - \w'\\\$3'u)'\
\s0\H'0'\
\\\$3
..
.
.ps 72z
.vs 76p
.SMALLCAPS V AN
.SMALLCAPS G OGH ,
.br
VAN GOGH,
```