Re: [Groff] About numerical precision

[Werner LEMBERG]

> 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


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

.\" .add30to60 <x> <y> <z>
.\"
.\" 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]
.
.de1 add30to60
.  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
.  add30to60 .ps/2 temp1 temp2
.  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,