\ ========== 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
\ - 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
\ ========== 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 )
\ 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.
\ 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
\ 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
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
\ 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' )
3 pick over um* ( u +d.resid guess1' d.err )
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)
: sm/rem ( d n -- rem quot )
2 pick 0< if ( d n sign )
1+ 2swap dnegate 2swap ( +d n sign )
2+ swap negate swap ( +d +n sign )
>r um/mod r> ( u.rem u.quot sign )
1 of swap negate swap negate endof \ -dividend, +divisor
2 of negate endof \ +dividend, -divisor
3 of swap negate swap endof \ -dividend, -divisor
: 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
\ In the usual case of similar signs (i.e. positive quotient),
\ sm/rem gives the correct answer
sm/rem ( n.rem' n.quot' )
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 ;
projects / OpenSPARC-T2-SAM / blob