[OpenSPARC-T2-SAM] / obp / obp / fm / kernel / dmuldiv.fth

\ ========== Copyright Header Begin ==========================================
\ 
\ Hypervisor Software File: dmuldiv.fth
\ 
\ Copyright (c) 2006 Sun Microsystems, Inc. All Rights Reserved.
\ 
\  - Do no alter or remove copyright notices
\ 
\  - Redistribution and use of this software in source and binary forms, with 
\    or without modification, are permitted provided that the following 
\    conditions are met: 
\ 
\  - Redistribution of source code must retain the above copyright notice, 
\    this list of conditions and the following disclaimer.
\ 
\  - Redistribution in binary form must reproduce the above copyright notice,
\    this list of conditions and the following disclaimer in the
\    documentation and/or other materials provided with the distribution. 
\ 
\    Neither the name of Sun Microsystems, Inc. or the names of contributors 
\ may be used to endorse or promote products derived from this software 
\ without specific prior written permission. 
\ 
\     This software is provided "AS IS," without a warranty of any kind. 
\ ALL EXPRESS OR IMPLIED CONDITIONS, REPRESENTATIONS AND WARRANTIES, 
\ INCLUDING ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A 
\ PARTICULAR PURPOSE OR NON-INFRINGEMENT, ARE HEREBY EXCLUDED. SUN 
\ MICROSYSTEMS, INC. ("SUN") AND ITS LICENSORS SHALL NOT BE LIABLE FOR 
\ ANY DAMAGES SUFFERED BY LICENSEE AS A RESULT OF USING, MODIFYING OR 
\ DISTRIBUTING THIS SOFTWARE OR ITS DERIVATIVES. IN NO EVENT WILL SUN 
\ OR ITS LICENSORS BE LIABLE FOR ANY LOST REVENUE, PROFIT OR DATA, OR 
\ FOR DIRECT, INDIRECT, SPECIAL, CONSEQUENTIAL, INCIDENTAL OR PUNITIVE 
\ DAMAGES, HOWEVER CAUSED AND REGARDLESS OF THE THEORY OF LIABILITY, 
\ ARISING OUT OF THE USE OF OR INABILITY TO USE THIS SOFTWARE, EVEN IF 
\ SUN HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
\ 
\ You acknowledge that this software is not designed, licensed or
\ intended for use in the design, construction, operation or maintenance of
\ any nuclear facility. 
\ 
\ ========== Copyright Header End ============================================
\ dmuldiv.fth 1.4 95/04/19
\ Copyright 1994 FirmWorks

\ Extended precision multiplication and division

\ alias um*  u*x  ( u1 u2 -- ud )

/n-t 4 * constant bits/half-cell
/n-t 8 * constant bits/cell

: scale-up    ( -- )  bits/half-cell <<  ;
: scale-down  ( -- )  bits/half-cell >>  ;
: split-halves  ( n -- low-half high-half )
   dup 1 scale-up 1- and  swap scale-down
;

\ This is the elementary school long-division algorithm, base 2^^16 (on a
\ 32-bit system) or 2^32 (on a 64-bit system).
\ It depends on the assumption that "/" can accurately divide a single-cell
\ (i.e. 32 or 64 bit) number by a half-cell (i.e. 16 or 32 bit) number.
\ Each "digit" is a half-cell number; thus the dividend is a 4-digit
\ number "ABCD" and the divisor is a 2-digit number "EF".

\ It would be interesting to compare the performance of this to a
\ "bit-banging" non-restoring division loop.
: um/mod  ( ud u -- urem uquot )
   2dup u>=  if                    \ Overflow; return max-uint for quotient
      0=  if  0 /  then            \ Force divide by zero trap
      2drop   0 -1  exit           ( 0 max-u )
   then                            ( ud u )

   \ Split the divisor into two 16-bit "digits"
   dup split-halves		   ( ud u ulow uhigh )

   \ If the high "digit" of the divisor is zero, we can skip a lot
   \ of the steps.  In this case, we only have to worry about the
   \ middle two digits of the dividend in developing the quotient.
   ?dup 0=  if                     ( ud u ulow )

      \ Approximate the high digit of the quotient by dividing the "BC"
      \ digits by the "F" digit.  The answer could by low by one, but if
      \ so it will be fixed in the next step.
      2over swap scale-down swap scale-up + over /  scale-up
				   ( ud u ulow guess<< )

      \ Multiply the trial quotient by the divisor
      rot over um*                 ( ud ulow guess<< udtemp )

      \ Subtract the trial product from the dividend, giving the remainder
      2swap >r >r  d-  drop        ( error )  ( r: guess<< ulow )

      \ Divide the remainder by the divisor, giving the rest of the
      \ quotient.
      dup r@ u/mod nip             ( error guess1 )
      r> r>                        ( error guess1 ulow guess<< )

      \ Merge the two halves of the quotient
      2 pick + >r                  ( error guess1 ulow ) ( r: uquot )

      \ Calculate the remainder
      * -  r>                      ( urem uquot )
      exit
   then                            ( ud u ulow uhigh )

   \ The high divisor digit is non-zero, so we have to deal with
   \ both digits, dividing "ABCD" by "EF".

   \ Approximate the high digit of the quotient.
   3 pick over u/mod nip           ( ud u ulow uhigh guess )

   \ Reduce guess by one if "E" = "A"
   dup 1 scale-up =  if  1-  then  ( ud u ulow uhigh guess' )

   \ Multiply the trial quotient by the divisor
   3 pick over scale-up um*        ( ud u ulow uhigh guess' ud.temp )

   \ Subtract the trial product from the dividend, giving the remainder
   >r >r 2rot r> r> d-             ( u ulow uhigh guess' d.resid )

   \ If the remainder is negative, add the divisor and reduce the trial
   \ quotient by one.  The following loop executes at most twice.
   begin  dup 0<  while            ( u ulow uhigh guess' d.resid )
      rot 1- -rot                  ( u ulow uhigh guess+ d.resid )
      4 pick scale-up 4 pick d+    ( u ulow uhigh guess+ d.resid' )
   repeat                          ( u ulow uhigh guess+ +d.resid )

   \ Now we have the correct high quotient digit; save it for later
   rot scale-up >r                 ( u ulow uhigh +d.resid ) ( r: q.high )

   \ Repeat the above process, using the partial remainder as the
   \ dividend.  Ulow is no longer needed
   3 roll drop                     ( u uhigh +d.resid )

   \ Trial quotient digit...
   2dup scale-up swap scale-down + 3 roll u/mod nip
                                   ( u +d.resid guess1 ) ( r: q.high )
   dup 1 scale-up =  if  1-  then  ( u +d.resid guess1' )

   \ Trial product
   3 pick over um*                 ( u +d.resid guess1' d.err )

   \ New partial remainder
   rot >r d-                       ( u d.resid' )  ( r: q.high guess1' )

   \ Adjust quotient digit until partial remainder is positive
   begin  dup 0<  while            ( u d.resid' )  ( r: q.high guess1' )
     r> 1- >r                      ( u d.resid' )  ( r: q.high guess1' )
     2 pick m+                     ( u d.resid'' ) ( r: q.high guess1' )
   repeat                          ( u +d.resid )  ( r: q.high guess1' )

   \ Discard divisior and high cell of quotient (which must be zero)
   rot 2drop                       ( u.rem )

   \ Merge quotient digits
   r> r> +                         ( u.rem u.quot )
;
: sm/rem  ( d n -- rem quot )
   0                           ( d n sign )
   2 pick 0<  if               ( d n sign )
      1+  2swap dnegate 2swap  ( +d n sign )
   then                        ( +d n sign )
   over 0<  if                 ( +d n sign )
      2+  swap negate swap     ( +d +n sign )
   then                        ( +d +n sign )
   >r um/mod r>                ( u.rem u.quot sign )
   case
      1 of  swap negate swap negate  endof  \ -dividend, +divisor
      2 of  negate                   endof  \ +dividend, -divisor
      3 of  swap negate swap         endof  \ -dividend, -divisor
   endcase
;
: fm/mod  ( d.dividend n.divisor -- n.rem n.quot )
   2dup xor 0<  if        \ Fixup only if operands have opposite signs
      dup >r  sm/rem                                ( rem' quot' r: divisor )
      over  if  1- swap r> + swap  else  r> drop  then
      exit
   then
   \ In the usual case of similar signs (i.e. positive quotient),
   \ sm/rem gives the correct answer
   sm/rem   ( n.rem' n.quot' )
;

: m*  ( n1 n2 -- d )
   2dup xor >r  abs swap  abs um*  r> 0<  if  dnegate  then
;
: */mod  ( n1 n2 n3 -- n.mod n.quot )  >r m* r> fm/mod  ;
: */  ( n1 n2 n3 -- n4 )  */mod nip  ;
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 ;