Commit | Line | Data |
---|---|---|
920dae64 AT |
1 | \ ========== Copyright Header Begin ========================================== |
2 | \ | |
3 | \ Hypervisor Software File: dmuldiv.fth | |
4 | \ | |
5 | \ Copyright (c) 2006 Sun Microsystems, Inc. All Rights Reserved. | |
6 | \ | |
7 | \ - Do no alter or remove copyright notices | |
8 | \ | |
9 | \ - Redistribution and use of this software in source and binary forms, with | |
10 | \ or without modification, are permitted provided that the following | |
11 | \ conditions are met: | |
12 | \ | |
13 | \ - Redistribution of source code must retain the above copyright notice, | |
14 | \ this list of conditions and the following disclaimer. | |
15 | \ | |
16 | \ - Redistribution in binary form must reproduce the above copyright notice, | |
17 | \ this list of conditions and the following disclaimer in the | |
18 | \ documentation and/or other materials provided with the distribution. | |
19 | \ | |
20 | \ Neither the name of Sun Microsystems, Inc. or the names of contributors | |
21 | \ may be used to endorse or promote products derived from this software | |
22 | \ without specific prior written permission. | |
23 | \ | |
24 | \ This software is provided "AS IS," without a warranty of any kind. | |
25 | \ ALL EXPRESS OR IMPLIED CONDITIONS, REPRESENTATIONS AND WARRANTIES, | |
26 | \ INCLUDING ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A | |
27 | \ PARTICULAR PURPOSE OR NON-INFRINGEMENT, ARE HEREBY EXCLUDED. SUN | |
28 | \ MICROSYSTEMS, INC. ("SUN") AND ITS LICENSORS SHALL NOT BE LIABLE FOR | |
29 | \ ANY DAMAGES SUFFERED BY LICENSEE AS A RESULT OF USING, MODIFYING OR | |
30 | \ DISTRIBUTING THIS SOFTWARE OR ITS DERIVATIVES. IN NO EVENT WILL SUN | |
31 | \ OR ITS LICENSORS BE LIABLE FOR ANY LOST REVENUE, PROFIT OR DATA, OR | |
32 | \ FOR DIRECT, INDIRECT, SPECIAL, CONSEQUENTIAL, INCIDENTAL OR PUNITIVE | |
33 | \ DAMAGES, HOWEVER CAUSED AND REGARDLESS OF THE THEORY OF LIABILITY, | |
34 | \ ARISING OUT OF THE USE OF OR INABILITY TO USE THIS SOFTWARE, EVEN IF | |
35 | \ SUN HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. | |
36 | \ | |
37 | \ You acknowledge that this software is not designed, licensed or | |
38 | \ intended for use in the design, construction, operation or maintenance of | |
39 | \ any nuclear facility. | |
40 | \ | |
41 | \ ========== Copyright Header End ============================================ | |
42 | \ dmuldiv.fth 1.4 95/04/19 | |
43 | \ Copyright 1994 FirmWorks | |
44 | ||
45 | \ Extended precision multiplication and division | |
46 | ||
47 | \ alias um* u*x ( u1 u2 -- ud ) | |
48 | ||
49 | /n-t 4 * constant bits/half-cell | |
50 | /n-t 8 * constant bits/cell | |
51 | ||
52 | : scale-up ( -- ) bits/half-cell << ; | |
53 | : scale-down ( -- ) bits/half-cell >> ; | |
54 | : split-halves ( n -- low-half high-half ) | |
55 | dup 1 scale-up 1- and swap scale-down | |
56 | ; | |
57 | ||
58 | \ This is the elementary school long-division algorithm, base 2^^16 (on a | |
59 | \ 32-bit system) or 2^32 (on a 64-bit system). | |
60 | \ It depends on the assumption that "/" can accurately divide a single-cell | |
61 | \ (i.e. 32 or 64 bit) number by a half-cell (i.e. 16 or 32 bit) number. | |
62 | \ Each "digit" is a half-cell number; thus the dividend is a 4-digit | |
63 | \ number "ABCD" and the divisor is a 2-digit number "EF". | |
64 | ||
65 | \ It would be interesting to compare the performance of this to a | |
66 | \ "bit-banging" non-restoring division loop. | |
67 | : um/mod ( ud u -- urem uquot ) | |
68 | 2dup u>= if \ Overflow; return max-uint for quotient | |
69 | 0= if 0 / then \ Force divide by zero trap | |
70 | 2drop 0 -1 exit ( 0 max-u ) | |
71 | then ( ud u ) | |
72 | ||
73 | \ Split the divisor into two 16-bit "digits" | |
74 | dup split-halves ( ud u ulow uhigh ) | |
75 | ||
76 | \ If the high "digit" of the divisor is zero, we can skip a lot | |
77 | \ of the steps. In this case, we only have to worry about the | |
78 | \ middle two digits of the dividend in developing the quotient. | |
79 | ?dup 0= if ( ud u ulow ) | |
80 | ||
81 | \ Approximate the high digit of the quotient by dividing the "BC" | |
82 | \ digits by the "F" digit. The answer could by low by one, but if | |
83 | \ so it will be fixed in the next step. | |
84 | 2over swap scale-down swap scale-up + over / scale-up | |
85 | ( ud u ulow guess<< ) | |
86 | ||
87 | \ Multiply the trial quotient by the divisor | |
88 | rot over um* ( ud ulow guess<< udtemp ) | |
89 | ||
90 | \ Subtract the trial product from the dividend, giving the remainder | |
91 | 2swap >r >r d- drop ( error ) ( r: guess<< ulow ) | |
92 | ||
93 | \ Divide the remainder by the divisor, giving the rest of the | |
94 | \ quotient. | |
95 | dup r@ u/mod nip ( error guess1 ) | |
96 | r> r> ( error guess1 ulow guess<< ) | |
97 | ||
98 | \ Merge the two halves of the quotient | |
99 | 2 pick + >r ( error guess1 ulow ) ( r: uquot ) | |
100 | ||
101 | \ Calculate the remainder | |
102 | * - r> ( urem uquot ) | |
103 | exit | |
104 | then ( ud u ulow uhigh ) | |
105 | ||
106 | \ The high divisor digit is non-zero, so we have to deal with | |
107 | \ both digits, dividing "ABCD" by "EF". | |
108 | ||
109 | \ Approximate the high digit of the quotient. | |
110 | 3 pick over u/mod nip ( ud u ulow uhigh guess ) | |
111 | ||
112 | \ Reduce guess by one if "E" = "A" | |
113 | dup 1 scale-up = if 1- then ( ud u ulow uhigh guess' ) | |
114 | ||
115 | \ Multiply the trial quotient by the divisor | |
116 | 3 pick over scale-up um* ( ud u ulow uhigh guess' ud.temp ) | |
117 | ||
118 | \ Subtract the trial product from the dividend, giving the remainder | |
119 | >r >r 2rot r> r> d- ( u ulow uhigh guess' d.resid ) | |
120 | ||
121 | \ If the remainder is negative, add the divisor and reduce the trial | |
122 | \ quotient by one. The following loop executes at most twice. | |
123 | begin dup 0< while ( u ulow uhigh guess' d.resid ) | |
124 | rot 1- -rot ( u ulow uhigh guess+ d.resid ) | |
125 | 4 pick scale-up 4 pick d+ ( u ulow uhigh guess+ d.resid' ) | |
126 | repeat ( u ulow uhigh guess+ +d.resid ) | |
127 | ||
128 | \ Now we have the correct high quotient digit; save it for later | |
129 | rot scale-up >r ( u ulow uhigh +d.resid ) ( r: q.high ) | |
130 | ||
131 | \ Repeat the above process, using the partial remainder as the | |
132 | \ dividend. Ulow is no longer needed | |
133 | 3 roll drop ( u uhigh +d.resid ) | |
134 | ||
135 | \ Trial quotient digit... | |
136 | 2dup scale-up swap scale-down + 3 roll u/mod nip | |
137 | ( u +d.resid guess1 ) ( r: q.high ) | |
138 | dup 1 scale-up = if 1- then ( u +d.resid guess1' ) | |
139 | ||
140 | \ Trial product | |
141 | 3 pick over um* ( u +d.resid guess1' d.err ) | |
142 | ||
143 | \ New partial remainder | |
144 | rot >r d- ( u d.resid' ) ( r: q.high guess1' ) | |
145 | ||
146 | \ Adjust quotient digit until partial remainder is positive | |
147 | begin dup 0< while ( u d.resid' ) ( r: q.high guess1' ) | |
148 | r> 1- >r ( u d.resid' ) ( r: q.high guess1' ) | |
149 | 2 pick m+ ( u d.resid'' ) ( r: q.high guess1' ) | |
150 | repeat ( u +d.resid ) ( r: q.high guess1' ) | |
151 | ||
152 | \ Discard divisior and high cell of quotient (which must be zero) | |
153 | rot 2drop ( u.rem ) | |
154 | ||
155 | \ Merge quotient digits | |
156 | r> r> + ( u.rem u.quot ) | |
157 | ; | |
158 | : sm/rem ( d n -- rem quot ) | |
159 | 0 ( d n sign ) | |
160 | 2 pick 0< if ( d n sign ) | |
161 | 1+ 2swap dnegate 2swap ( +d n sign ) | |
162 | then ( +d n sign ) | |
163 | over 0< if ( +d n sign ) | |
164 | 2+ swap negate swap ( +d +n sign ) | |
165 | then ( +d +n sign ) | |
166 | >r um/mod r> ( u.rem u.quot sign ) | |
167 | case | |
168 | 1 of swap negate swap negate endof \ -dividend, +divisor | |
169 | 2 of negate endof \ +dividend, -divisor | |
170 | 3 of swap negate swap endof \ -dividend, -divisor | |
171 | endcase | |
172 | ; | |
173 | : fm/mod ( d.dividend n.divisor -- n.rem n.quot ) | |
174 | 2dup xor 0< if \ Fixup only if operands have opposite signs | |
175 | dup >r sm/rem ( rem' quot' r: divisor ) | |
176 | over if 1- swap r> + swap else r> drop then | |
177 | exit | |
178 | then | |
179 | \ In the usual case of similar signs (i.e. positive quotient), | |
180 | \ sm/rem gives the correct answer | |
181 | sm/rem ( n.rem' n.quot' ) | |
182 | ; | |
183 | ||
184 | : m* ( n1 n2 -- d ) | |
185 | 2dup xor >r abs swap abs um* r> 0< if dnegate then | |
186 | ; | |
187 | : */mod ( n1 n2 n3 -- n.mod n.quot ) >r m* r> fm/mod ; | |
188 | : */ ( n1 n2 n3 -- n4 ) */mod nip ; |