# Copyright 1992, Digital Equipment Corporation
#
# This software is furnished under a license and may be used and copied
# only in accordance with the terms of such license and with the
# inclusion of the above copyright notice. This software or any other
# copies thereof may not be provided or otherwise made available to any
# other person. No title to and ownership of the software is hereby
# transferred.
#
# The information in this software is subject to change without notice
# and should not be construed as a commitment by Digital Equipment
# Corporation.
#
# Digital assumes no responsibility for the use or reliability of its
# software on equipment which is not supplied by Digital.
#
# 008 17 Jun 1992 KDG/wbn Most of initial tailored version. (See
# commentary below.)
#
# 009 4 Jul 1992 KDG Continue work on initial tailored version,
# including bugfixes and mod entry points
#
# 010 15 Jul 1992 KDG - Final touches for V1 (other than any bugfixes)
# - .aligns commented out to allow older assembler versions
#
# 011 16 Jul 1992 KDG - Bugfix for ots_div_l for -maxint dividend
# - OSF-only source changes for BL7
#
# 012 10 Aug 1992 KDG Fix overflow division entry points
#
# 013 23 Sep 1992 KDG Add case-sensitive entry names
#
# 014 4 Jan 1993 KDG Tweak for OSF assembler
#
# 015 26 Jan 1993 KDG Add underscore prefix, OSF uses CS names
#
# 016 5 Apr 1993 WBN Speed up core 64-bit, shrink table entry to 2 QWs
#++
# Entry points defined in this module:
#
# -- 32 bit division/remainder support
# unsigned ots_rem_ui(unsigned dividend, unsigned divisor)
# unsigned ots_div_ui(unsigned dividend, unsigned divisor)
# int ots_mod_i(int dividend, int modulus)
# int ots_rem_i(int dividend, int divisor)
# int ots_div_i_o(int dividend, int divisor)
# int ots_div_i(int dividend, int divisor) | "hot spot"
# {core routine - div32} |
#
# -- 64 bit division support
# {core routine - div64} | (uses div32 for 32b cases)
# long ots_div_l_o(long dividend, long divisor)
# long ots_div_l(long dividend, long divisor) | "hot spot"
# long ots_rem_l(long dividend, long divisor)
# long ots_mod_l(long dividend, long modulus)
# unsigned long ots_div_ul(unsigned long dividend, unsigned long divisor)
# unsigned long ots_rem_ul(unsigned long dividend, unsigned long divisor)
#
# Special conventions: No stack space, r0-r1, r16-r19 and r26-r28 ONLY.
# (Warning: The auto-loader potentially takes some regs across
# the call if this is being used in a shared lib. environment.)
#
# NOTE: This set of routines may start using stack space at some
# future point in time.
#
# -- Possible future entry points include:
# (These all return results in r0/r1)
# {int quotient, int remainder} ots_div_mod_i(int dividend, int divisor)
# {int quotient, int remainder} ots_div_mod_i_o(int dividend, int divisor)
# {int quotient, int remainder} ots_div_rem_i(int dividend, int divisor)
# {int quotient, int remainder} ots_div_rem_i_o(int dividend, int divisor)
# {unsigned quotient, unsigned remainder} ots_div_rem_ui(unsigned dividend, unsigned divisor)
#
# {long quotient, long remainder} ots_div_mod_i(long dividend, long divisor)
# {long quotient, long remainder} ots_div_mod_i_o(long dividend, long divisor)
# {long quotient, long remainder} ots_div_rem_i(long dividend, long divisor)
# {long quotient, long remainder} ots_div_rem_i_o(long dividend, long divisor)
# {unsigned long quotient, unsigned long remainder}
# ots_div_rem_ui(unsigned long dividend, unsigned long divisor)
#
#
# General commentary:
#
# This is an attempt at a fairly high performance version using relatively
# straightforward algorithms. Note that the code is intended to be scheduled
# well for EV4, but still reasonably for LCA/EV5.
#
# Also, note that there was only so much time available for this, so it
# is far from "perfect". "Better is the enemy of done"...
#
# Possible future areas of improvement (and unfinished business):
#
# - Another possible way of doing things for the "slow" (divnn cases)
# is to use an approximate inverse and convergence. Given the speed
# of the multiplier on EV4, and given "time to market", this wasn't
# done for V1.) I have some mail with the algorithm from Bob Gries
# (through Scott Robinson).
#
# - When the divisor is too large for the table, but has n low-order zero
# bits, see if divisor/2^n fits in the table, and use that entry with
# dividend/2^n
#
# - Use UMULH for the 'mod' routines.
#
# This version can do a table lookup division (divisors with <=tablesize)
# in roughly 32 cycles on an EV4 (with cache hits for all loads), of which
# 21 are for the umulh. There is a strong bias toward the table-lookup case.
# Note that for many cases, the umulh is the last thing before the return,
# so the multiply can occur in parallel with the procedure return.
# (It is interesting that the R3000 hardware divide instruction takes 33
# cycles and the R4000 takes 76(!) ...) Small powers of 2 are retired in
# roughly 20 cycles. Larger divisors take considerably longer at this point.
#
#include "ots_defs.hs"
#ifdef OSF
# to get the PAL_gentrap literal
#include <alpha/pal.h>
#endif
# Data area description
#
# The data area "ots_div_data" is an array of structures, indexed
# by the divisor value, with each array entry being 16 bytes in size
# formatted as follows:
#
# 6
# 3 6 0
# +-------+-------+-------+-------+-------+-------+-------+-------+
# | 32 bit reciprocal (58 bits) |shift|
# +-------+-------+-------+-------+-------+-------+-------+-------+
# | 64 bit reciprocal |
# +---------------------------------------------------------------+
#
# The 64-bit reciprocal has the leading '1' bit omitted, so it provides
# 65 bits of precision -- enough to handle unsigned 64-bit values.
#
# The first longword contains the 6-bit shift amount needed to handle
# 64-bit cases and powers of two.
#
# The 32-bit reciprocal has the shift count built in, so a UMULH gives
# the correct quotient without shifting. The reciprocal needs 33 bits
# of precision. The 6-bit shift amount is noise in the reciprocal that
# can be ignored.
#
# (Oh, you want proof?) For divisors up to 2^k, we store k-1 leading
# zero bits, 33 bits of fraction, (25-k) more bits of fraction, and
# 6 bits of noise. The standard method would round at the 33rd fraction
# bit. We need to ensure that the value actually stored is geq the
# infinite reciprocal, but leq the standard value. For divisors up to
# 2^k, there will be a zero bit somewhere in the k bits below the 33rd,
# so as long as we round below the (33+k)th bit, the rounded value
# plus any noise is still less than the standard value. This requires
# k < 12.
#
# The actual data is declared in ots_div_data_alpha.
#
# Offsets to the various fields in the data structure
#
#define shift_o 0
#define recip64_o 8
#define recip32_o 0
#
# Note that the shift/add ops used to compute the table entries
# "know" that the table size is 16. (i.e. addq -> s8addq -> ldq)
# By changing the first instruction, it's fairly easy to change the
# table entry size to 24, 32, or 40 bytes (using s4add/sub), or
# 56/64/72 bytes using s8add/sub, should that be desirable.
# Maximum divisor present in the table
#
#define table_max 512
# Division by zero gentrap code
#
#define GEN_INTDIV -2
# Address of division data area (shared by all entry points)
#
#ifdef VMS
.psect ots_link
ots_div_addr:
.address ots_div_data
.psect ots_code
#endif
# Dummy entry point for the module
#
.globl _OtsDivide
.ent _OtsDivide
_OtsDivide:
.set noat
.set noreorder
#ifdef OSF
.frame sp, 0, r26
#endif
#ifdef WNT
.frame sp, 0, r26
#endif
# unsigned ots_rem_ui(unsigned dividend, unsigned divisor)
# unsigned 32 bit remainder support
#
#.align 4
.globl _OtsRemainder32Unsigned
.aent _OtsRemainder32Unsigned
_OtsRemainder32Unsigned:
#ifdef VMS
ldq r27, <ots_div_addr-ots_rem_ui>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
lda r28, -table_max(r17) # test for table lookup
subl r17, 1, r1 # first part of power-of-2 check
blt r17, rui_big # big divisors can (must) be handled by a simple comparison
and r17, r1, r18 # second part of power-of-2 check
bgt r28, rui_lrgdiv # branch if large divisor
addq r17, r17, r0 # compute divisor*2 for table lookup
beq r18, rui_pwr2 # if zero, divisor is a power of 2
s8addq r0, r27, r27 # finish computing table entry addr (table addr+divisor*16)
ldq r1, recip32_o(r27) # load approximate reciprocal
cmpult r16, r17, r18 # is the dividend < divisor?
zap r16, 0xF0, r0 # kill the propagated sign bit
bne r18, rui_lss # if dividend < divisor, fast exit
umulh r0, r1, r0 # multiplication for division step
mull r0, r17, r0 # multiply back to get value to subtract
subl r16, r0, r0
ret r31, (r26) # and return
rui_pwr2:
beq r17, divzer # check for 0
and r16, r1, r0 # use x-1 to mask
ret r31, (r26)
rui_lss:
mov r16, r0
ret r31, (r26)
rui_lrgdiv:
zap r16, 0xf0, r16 # zero-extend the dividend
bsr r28, div32 # use the core routine getting the remainder in r1
sextl r1, r0
ret r31, (r26)
# divisors with the sign bit set. two possible results,
# dividend if dividend < divisor, or dividend-divisor otherwise
rui_big:
cmpult r16, r17, r1
subl r16, r17, r0
cmovne r1, r16, r0
ret r31, (r26)
# unsigned ots_div_ui(unsigned dividend, unsigned divisor)
# unsigned 32 bit division support
#
#.align 4
.globl _OtsDivide32Unsigned
.aent _OtsDivide32Unsigned
_OtsDivide32Unsigned:
#ifdef VMS
ldq r27, <ots_div_addr-ots_div_ui>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
lda r28, -table_max(r17) # test for table lookup
blt r17, dui_big # big divisors can (must) be handled by a simple comparison
addq r17, r17, r18 # compute divisor*2
cmpule r17, r16, r0 # is the dividend < divisor?
beq r17, divzer # check for 0
s8addq r18, r27, r27 # finish computing table entry addr (table addr+divisor*16)
beq r0, dui_end # fast out for divisor > dividend
bgt r28, dui_lrgdiv # branch if large divisor
ldq r1, recip32_o(r27) # load approximate reciprocal
zap r16, 0xF0, r16 # kill the propagated sign bit
blt r1, dui_smpwr2 # go handle powers of 2 specially
umulh r16, r1, r0 # start multiplication for division step
dui_end:ret r31, (r26) # and return
nop
dui_smpwr2:
srl r16, r1, r0 # shift the result into place
sextl r0, r0 # reinsert sign if dividing by 1
ret r31, (r26) #
dui_lrgdiv:
zap r16, 0xf0, r16 # zero-extend the dividend
bsr r28, div32 # use the core routine getting the remainder in r1
sextl r0, r0 # make sure the result is in normal form for uint32
ret r31, (r26)
# divisor with the sign bit set. two possible results,
# 1 if divisor <= dividend, or 0 otherwise
dui_big:
cmpule r17, r16, r0
ret r31, (r26)
# int ots_mod_i(int dividend, int modulus)
# signed 32 bit modulus support
#
# This entry could be MUCH more optimized. It doesn't even try to use
# UMULH division currently... (A casualty of time-to-market.)
# Note that mod is only used by Ada and PL/I.
#
#.align 4
.globl _OtsModulus32
.aent _OtsModulus32
_OtsModulus32:
negq r17, r18 # first part of abs(divisor)
cmovge r17, r17, r18 # second part of abs(divisor)
subq r18, 1, r1 # start checking for power of 2
beq r17, divzer # check for 0
and r18, r1, r0 # second part of power-of-2 check
beq r0, mi_p2 # for powers of two, simply do a mask
# (note that the power-of-2 case MUST be used to handle
# the -maxint case due to the way the fix-up info is
# saved across the core routine call)
xor r16, r17, r28 # get xor of signs
clr r19 # don't need a bias if dividend and divisor have same sign
cmovlt r28, r17, r19 # bias is original divisor for different sign case
and r16, r17, r27 # if both dividend & divisor were neg. need to negate result
mov r18, r17 # move abs(divisor) into r17
negq r16, r18 # first part of abs(dividend)
cmovlt r16, r18, r16 # second part of abs(dividend)
cmplt r27, r31, r0 # get 1 if both operands were <0
sll r0, 63, r0 # get bit as the high bit
bis r0, r19, r19 # and MERGE with bias (0 -> no fixup, -maxint -> negate result,
# divisor > 0 - subtract remainder if non-zero, divisor < 0 -
# add remainder if non-zero)
bsr r28, div32 # use the core routine getting the remainder in r1
cmoveq r1, r31, r19 # don't do any fix-up if the remainder was zero
addq r19, r19, r18 # check to see if this is the negative/negative case, which just gets a negated remainder
subq r19, 1, r28 # wrap -maxint to positive
negl r1, r0 # move negated value, may abort later
cmovlt r28, r1, r0 # if both positive, or negative divisor, keep positive remainder
cmoveq r18, r31, r19 # now that negation is done, treat -maxint case as 0
addl r19, r0, r0 # add any bias (original divisor or 0)
ret r31, (r26) # and return
mi_p2: cmovge r17, r31, r17 # no bias if divisor was >= 0
and r16, r1, r1 # use the divisor-1 mask that's already in r1
cmoveq r1, r31, r17 # use zero if result was zero
addl r17, r1, r0 # do any biasing, and ensure the result is sign ext
ret r31, (r26) # and return
# int ots_rem_i(int dividend, int divisor)
# signed 32 bit remainder support
#
#.align 4
.globl _OtsRemainder32
.aent _OtsRemainder32
_OtsRemainder32:
#ifdef VMS
ldq r27, <ots_div_addr-ots_rem_i>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
negq r17, r18 # first part of abs(divisor)
cmovlt r17, r18, r17 # second part of abs(divisor)
subq r17, 1, r1 # start checking for power of 2
and r17, r1, r0 # finish check for power of 2
sra r16, 63, r19 # get -1/0 if dividend was negative
negq r16, r18 # first part of abs(dividend)
cmovlt r16, r18, r16 # second part of abs(dividend)
beq r0, ri_pwr2 # for powers of two, simply do a mask (not power of 2 include 0 and 80000000)
lda r28, -table_max(r17) # test for table lookup
bgt r28, ri_lrgdiv # branch if large divisor
addq r17, r17, r0 # compute divisor*2 for table lookup
s8addq r0, r27, r27 # finish computing table entry addr (table addr+divisor*16)
ldq r1, recip32_o(r27) # load approximate reciprocal
umulh r16, r1, r0 # multiplication for division step
mull r0, r17, r0 # multiply back to get value to subtract
subl r16, r0, r0 # get abs of final result
xor r0, r19, r0 # start compliment if original dividend was <0
subl r0, r19, r0 # finish compliement
ret r31, (r26) # and return
# Handle powers of 2, including 0 and 80000000
ri_pwr2:
and r16, r1, r0 # use the divisor-1 mask in r1
beq r17, divzer # division by zero
xor r0, r19, r0 # start compliment if original dividend was <0
subl r0, r19, r0 # finish compliement
ret r31, (r26)
nop
ri_lrgdiv:
bsr r28, div32 # use the core routine getting the remainder in r1
xor r1, r19, r0 # start compliment if original dividend was <0
subl r0, r19, r0 # finish compliement
ret r31, (r26)
# int ots_div_i_o(int dividend, int divisor)
# signed 32 bit division support, overflow detection
#
#.align 4
.globl _OtsDivide32Overflow
.aent _OtsDivide32Overflow
_OtsDivide32Overflow:
#ifdef VMS
ldq r27, <ots_div_addr-ots_div_i_o>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
not r17, r1 # is the divisor -1?
bne r1, di_skip # continue if not
neglv r16, r0 # quotient = -dividend, overflow on ^x800000000
ret r31, (r26)
# int ots_div_i(int dividend, int divisor)
# signed 32 bit division support, no overflow detection
#
nop #.align 4
.globl _OtsDivide32
.aent _OtsDivide32
_OtsDivide32:
#ifdef VMS
ldq r27, <ots_div_addr-ots_div_i>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
di_skip:
di_retry:
lda r28, -table_max(r17) # test for table lookup
ble r17, di_notpos # not a positive divisor case
addq r17, r17, r0 # compute divisor*2
negq r16, r18 # part 1 of abs(dividend) -> r18. (Note 0xffffffff 80000000 => 0x00000000 80000000)
bgt r28, di_lrgdiv # branch if large divisor
s8addq r0, r27, r27 # finish computing table entry addr (table addr+divisor*16)
cmpule r17, r18, r0 # divisor <= dividend?
cmovge r16, r16, r18 # part 2 abs. val of the dividend -> r18
beq r0, di_end # if not, result is zero
ldq r1, recip32_o(r27) # load approximate reciprocal
blt r1, di_smpwr2 # go handle powers of 2 specially
umulh r18, r1, r0 # start multiplication
blt r16, di_negres # negate result? (done as a branch to allow umulh to "hang out" over end for common case)
di_end: ret r31, (r26) # return for same-sign (common) case
di_negres:
negl r0, r0 # different signs - compliment result
ret r31, (r26) # return for different-sign (uncommon) case
di_smpwr2:
srl r18, r1, r18 # shift the result into place
sra r16, 63, r16 # get 0/-1 based on sign of dividend
xor r18, r16, r18 # conditionally compliment
subl r18, r16, r0 # and increment for the final value
ret r31, (r26) # (note subl is required for sign ext for %x80000000/1 case)
# Zero or negative divisor case. If just a negative divisor,
# compliment both dividend and divisor and do things again.
di_notpos:
beq r17, divzer # division by zero
negl r17, r17 # |divisor|, note that 0x80000000 still appears negative
negq r16, r16 # compliment dividend (negq so that 0xffffffff 80000000 => 0x00000000 80000000
bgt r17, di_retry # dispatch back for normal case (not 0x80000000 or 0)
sextl r16, r16 #
cmpeq r16, r17, r0 # -maxint/-maxint = 1, all others 0
ret r31, (r26) # done
# Large divisor for signed 32/32 case
#
nop #.align 3
di_lrgdiv:
sra r16, 63, r19 # get 0/-1 based on sign of dividend
cmovlt r16, r18, r16 #
bsr r28, div32 # go use core routine
xor r0, r19, r0 # conditionally compliment
subl r0, r19, r0 # and increment for the final value (subl ensures normalized result)
ret r31, (r26) # done
# Large divisor case core routine for 32b
# (wbn)
#
# r0 - quotient (output)
# r1 - remainder (output)
# r16 - dividend (range 0..2^32-1, zero extended)
# r17 - divisor (range 1..2^31-1 - overwritten)
# r18 - scratch
# r19 - not used (one temp for 'caller')
# r26 - not used (expected to contain main return address)
# [r27 - scratch] (not currently written)
# r28 - this "subroutine" return address
#
# Some tightened up bit-at-a-time code for dividing 32-bit integers.
# It uses two tricks: keep the running remainder and the quotient in
# the same 64-bit register (MQ?), and add 1 while subtracting the divisor,
# so that a single CMOV sets both the new remainder and the new quotient.
# I start off by trying to skip 8 bits at a time; should this skip a
# smaller amount, so the main loop iterates less often? If the divisor
# is already known to be large enough, the last case in this test is never
# used...
#
# This code expects as input two integers in the range 0 <= x < 2^31
# (that is, it doesn't work for general unsigned longwords, and doesn't
# include sign manipulation.)
#
# The code here takes about 34n+11 cycles for a quotient occupying n bytes.
#
# Inputs: dividend in r16, divisor in r17
# Outputs: quotient in r0, remainder in r1
# Destroys: [r16,]r17,r18,[r27]
#
# How many quotient bytes will there be: 0, 1, 2, 3, 4?
#
#.align 4
div32: cmpule r17, r16, r0 # Divisor leq dividend?
sll r17, 32, r18 # Position divisor for loop
sll r17, 8, r1 # Prepare for next compare
beq r0, d32end # Dividend less, quotient is zero.
ediv32: mov 8-3, r17 # Hope to skip 3 bytes of loop
cmpule r1, r16, r0 # Shifted divisor still leq dividend?
sll r1, 8, r1 # Prepare for next compare
beq r0, d32ent # Go loop over just one byte
mov 8-2, r17 # Hope to skip 2 bytes of loop
cmpule r1, r16, r0 # Shifted divisor still leq dividend?
sll r1, 8, r1 # Prepare for next compare
beq r0, d32ent # Go loop over just two bytes
mov 8-1, r17 # Hope to skip 1 byte of loop
cmpule r1, r16, r0 # Shifted divisor still leq dividend?
nop # stall - align d32ent and d32loop
cmovne r0, 8, r17 # If we can't skip any bytes
# start loop generating quotient bits. NOTE: The loop setup requires
# an even number of iterations.
#
d32ent: extqh r16, r17, r0 # Shift dividend left for skipped bytes
subq r18, 1, r1 # Divisor in high LW - 1 in low LW
s8subq r17, 34, r17 # Convert bytes to bits and adjust
addq r0, r0, r0 # Shift left to start first iteration
d32loop:subq r0, r1, r18 # Can we subtract divisor from it?
cmovge r18, r18, r0 # If so, set new remainder & quotient
# stall
addq r0, r0, r0 # Shift remainder and quotient left
subq r0, r1, r18 # Can we subtract divisor from it?
cmovge r18, r18, r0 # If so, set new remainder & quotient
subq r17, 2, r17 # Loop counter
addq r0, r0, r0 # Shift remainder and quotient left
bgt r17, d32loop # Repeat
subq r0, r1, r18 # Can we subtract divisor from it?
cmovge r18, r18, r0 # If so, set new remainder & quotient
# stall
addq r0, r0, r0 # Shift remainder and quotient left
subq r0, r1, r18 # Finish last iteration
cmovge r18, r18, r0
# stall
srl r0, 32, r1 # Get remainder in r1
zap r0, 0xf0, r0 # Keep only quotient in r0
nop # for alignment
d32end: cmoveq r0, r16, r1 # Move remainder to r1 for quotient=0 case
ret r31, (r28) # Not a real software procedure return
# Large divisor case core routine for 64b
#
# r0 - quotient (output)
# r1 - remainder (output)
# r16 - dividend (range 0..2^64-1 - overwritten)
# r17 - divisor (range 1..2^63-1 - overwritten)
# r18 - scratch
# r19 - not used (one temp for 'caller')
# r26 - not used (expected to contain main return address)
# r27 - points to table of inverses (overwritten)
# r28 - this "subroutine" return address
#
# Inputs: dividend in r16, divisor in r17
# Outputs: quotient in r0, remainder in r1
# Destroys: r16,r17,r18,r27
#
# Note- this routine could save a few cycles if we could use
# another scratch register -- perhaps by pushing one on the stack?
#
#.align 4
div64: sll r17, 32, r18 # Position for ediv32
cmpule r17, r16, r0 # Is divisor leq dividend?
srl r17, 31, r1 # Is divisor geq 2^31?
beq r0, d64end # If divisor > dividend, quotient=0
cmpule r18, r16, r0 # Is divisor*2^32 leq dividend?
sll r17, 8, r1 # Position for ediv32 checking
or r1, r0, r0 # 0 if divisor & quotient fit in 32 bits
beq r0, ediv32 # Use 32-bit routine if OK
# Full 64-bit divide needed. Use the table of shift amounts to compute
# the number of leading zero bits in the divisor. Find the leftmost
# nonzero byte, then the leftmost nonzero bit in that byte. Table entry
# #n+1 contains the number of bits needed to hold n (1..8). We know the
# divisor is nonzero here.
#
cmpbge r31, r17, r0 # Get a zero bit for each nonzero byte
#stall
sll r0, 4, r0 # *16 bytes per table entry
#stall
subq r27, r0, r0 # table base plus complement...
#stall
ldq r1, 256*16(r0) # get position of first nonzero
#2 stalls
subq r1, 1, r1 # byte number of first nonzero
extbl r17, r1, r0 # get first nonzero byte
#stall
addq r0, r0, r0 # *2
s8addq r0, r27, r0 # *16 bytes per table entry
#stall
ldq r0, 16(r0) # bit number of first nonzero
negq r1, r1 # 1 + #leading zero bytes (mod 8)
#stall
s8subq r1, r0, r0 # number of leading zero bits
and r0, 0x3F, r0 # discard other junk
# The following code does a similar normalize calculation without the table.
#===
# extll r17, #4, r18 ; Normalize the divisor and
# mov #63, r0 ; count leading zeros
# cmovne r18, #31, r0
# cmoveq r18, r17, r18
# ;stall
# extwl r18, #2, r1
# ;stall
# cmovne r1, r1, r18
# cmovne r1, #16, r1
# ;stall
# subq r0, r1, r0
# extbl r18, #1, r1
# ;stall
# cmovne r1, r1, r18
# cmovne r1, #8, r1
# ;stall
# subq r0, r1, r0
# andnot r18, #^x0f, r1
# cmovne r1, r1, r18
# cmovne r1, #4, r1
# ;stall
# subq r0, r1, r0
# andnot r18, #^x33, r1
# cmovne r1, r1, r18
# cmovne r1, #2, r1
# ;stall
# subq r0, r1, r0
# andnot r18, #^x55, r1
# cmovne r1, #1, r1
# ;stall
# subq r0, r1, r0
#===
# R0 contains number of leading zero bits in the divisor.
sll r17, r0, r17 # Normalize: MSB is set.
# Now break divisor into pieces a+x, where a is the leading
# 9 bits, rounded, and x is the rest. Use a linear
# approximation for 1/divisor = 1/a - x/a^2 [+ x^2/a^3 -...]
#
srl r17, 64-10, r1 # Keep 10 bits of divisor
#stall
addq r1, 1, r1 # Round to form 'a'
andnot r1, 1, r1
s8addq r1, r27, r27 # Index table of 1/a and 1/a^2
sll r1, 64-10, r1 # shift 'a' to match divisor
ldq r18, (r27) # Load QW containing 1/a^2
subq r1, r17, r1 # -x = a - divisor
beq r1, d64_easy # Use table directly if x=0
inswl r18, 6, r18 # position 1/a^2
blt r1, d64_sub # correct for sign of -x
umulh r1, r18, r1 # -x/a^2
ldq r27, 8(r27) # Load QW containing 1/a - 1
br r31, d64_cont
d64_sub:umulh r1, r18, r1 # -x/a^2
ldq r27, 8(r27) # load QW containing 1/a - 1
# 2 stalls
s4addq r18, 0, r18
subq r27, r18, r27 # correct for sign of -x
d64_cont:
# many stalls
s4addq r1, r27, r18 # 1/divisor approx= 1/a - x/a^2
# Now one or two Newton iterations to get 24 or 56 good bits of the inverse.
# Each computes inv = inv * (2 - inv*divisor). We could skip out early
# here or above if the dividend and/or quotient is small enough for the
# amount of precision we've developed...
#
# We handle quadwords with the radix point on the left. The divisor has
# been normalized to the range 0.5 < divisor < 1.0; the inverses are in
# the range 1.0 < inverse < 2.0, and are represented without the leading 1.
#
umulh r18, r17, r1 # (inv0 - 1) * divisor
# many stalls
addq r1, r17, r1 # add hidden bit * divisor
negq r1, r1 # 2 - inv0*divisor, very near 1.0
umulh r18, r1, r27 # (inv0 - 1) * (2 - inv0*divisor)
cmovlt r1, 0, r18 # keep inv0 if (2-inv0*divisor) > 1.0
#stall
addq r18, r1, r1 # add it to hidden bit * (2-inv0*divisor)
# many stalls
addq r27, r1, r18 # inv1 = inv0 * (2 - inv0*divisor)
umulh r18, r17, r1 # (inv1 - 1) * divisor
# many stalls
addq r1, r17, r1 # add hidden bit * divisor
negq r1, r1 # 2 - inv1*divisor, very near 1.0
umulh r18, r1, r27 # (inv1 - 1) * (2 - inv1*divisor)
cmovlt r1, 0, r18 # keep inv1 if (2-inv1*divisor) > 1.0
addq r18, r1, r1 # add it to hidden bit * (2-inv1*divisor)
# many stalls
addq r27, r1, r1 # inverse = inv1 * (2 - inv1*divisor)
umulh r1, r16, r18 # dividend * (1/divisor - 1)
srl r17, r0, r17 # un-normalize divisor
negq r0, r0
subq r0, 8, r0 # how far right after first byte
# many stalls
addq r18, r16, r18 # add hidden bit * dividend
cmpult r18, r16, r1 # did it carry?
srl r18, 8, r18 # start to shift
sll r1, 56, r1 # position the carry
#stall
addq r1, r18, r1 # add the carry
srl r1, r0, r0 # final shift
mulq r17, r0, r1 # try out this quotient
# many stalls
subq r16, r1, r1 # form remainder
cmpule r17, r1, r18 # must be less than divisor
subq r1, r17, r27
cmovne r18, r27, r1 # if not, decrement remainder
addq r0, r18, r0 # and increment quotient
ret r31, (r28) # done
d64_easy:
ldq r1, 8(r27) # get 1/divisor, except hidden bit
srl r17, r0, r17 # un-normalize divisor again
blt r18, d64_pow2 # skip if power of 2
umulh r1, r16, r18 # dividend/divisor
negq r0, r0 # how far right to shift
and r16, 0x0ff, r1 # pieces of dividend
subq r0, 8, r0 # how far right after first byte
srl r16, 8, r27
# many stalls
addq r18, r1, r1 # add low piece of dividend, no carry
srl r1, 8, r1 # make room for high piece
#stall
addq r1, r27, r1 # finish adding hidden bit * dividend
srl r1, r0, r0 # final shift
mulq r17, r0, r1 # need to compute remainder too
# many stalls
subq r16, r1, r1
ret r31, (r28) # done
d64_pow2:
not r0, r0 # how far right to shift quotient
subq r17, 1, r1 # mask for remainder
srl r16, r0, r0 # shift for quotient
and r16, r1, r1 # get remainder
ret r31, (r28) # done
d64end: mov r16, r1 # Remainder to r1 for small quotient case
ret r31, (r28) # Not a real software procedure return
# long ots_div_l_o(long dividend, long divisor)
# signed 64 bit division support, overflow detection
#
#.align 4
.globl _OtsDivide64Overflow
.aent _OtsDivide64Overflow
_OtsDivide64Overflow:
#ifdef VMS
ldq r27, <ots_div_addr-ots_div_l_o>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
not r17, r1 # is the divisor -1?
bne r1, dl_skip # continue if not
negqv r16, r0 # q = -dividend, oflow on ^x800000000 00000000
ret r31, (r26)
nop
# long ots_div_l(long dividend, long divisor)
# signed 64 bit division support, no overflow detection
#
.globl _OtsDivide64
.aent _OtsDivide64
_OtsDivide64:
#ifdef VMS
ldq r27, <ots_div_addr-ots_div_l>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
dl_skip:
xor r16, r17, r19 # sign bit = result needs to be complimented (here to handle -maxint correctly)
dl_retry:
lda r28, -table_max(r17) # test for table lookup
ble r17, dl_notpos # not a positive divisor case
addq r17, r17, r0 # compute divisor*2
negq r16, r18 # part 1 of abs(dividend) -> r18
bgt r28, dl_lrgdiv # branch if large divisor
s8addq r0, r27, r27 # finish computing table entry addr (table addr+divisor*16)
srl r16, 33, r1 # can this be handled via a 32 bit case?
cmpule r17, r18, r0 # divisor <= dividend?
bne r1, dl_64bit # does this need to be a real 64 bit case?
cmovge r16, r16, r18 # part 2 abs. val of the dividend -> r18
beq r0, dl_end # if not, result is zero
ldq r27, recip32_o(r27) # load 32b approximate reciprocal
sra r19, 63, r19 # get 0/-1
blt r27, dl_smpwr2 # skip umulh for powers of 2 specially
umulh r18, r27, r0 # start multiplication
beq r19, dl_end # if compliment not required, let umulh "hang out"
negq r0, r0 # compliment case
ret r31, (r26) #
dl_64bit:
cmovge r16, r16, r18 # part 2 abs. val of the dividend -> r18
beq r0, dl_end # if not, result is zero
ldq r1, recip64_o(r27) # load approximate reciprocal
sra r19, 63, r19 # get 0/-1
ldq r27, shift_o(r27) # load shift count (low 6 bits are all that matters)
beq r1, dl_smpwr2 # skip umulh for powers of 2 specially
umulh r18, r1, r0 # start multiplication
addq r0, r18, r18 # add hidden bit
dl_smpwr2:
srl r18, r27, r18 # shift the result into place
xor r18, r19, r18 # conditionally compliment
subq r18, r19, r0 # and increment for the final value
dl_end: ret r31, (r26) #
# Zero or negative divisor case. If just a negative divisor,
# compliment both dividend and divisor and do things again.
dl_notpos:
beq r17, divzer # division by zero
negq r17, r17 # |divisor|, note that 0x80000000 00000000 still appears negative
negq r16, r16 # compliment dividend
bgt r17, dl_retry # dispatch back for normal case (not 0x80000000 00000000 or 0)
cmpeq r16, r17, r0 # -maxint/-maxint = 1, all others 0
ret r31, (r26) # done
# Large divisor for signed 64/64 case
#
dl_lrgdiv:
sra r19, 63, r19 # get 0/-1
cmovlt r16, r18, r16 #
bsr r28, div64 # go use core routine
xor r0, r19, r0 # conditionally compliment
subq r0, r19, r0 # and increment for the final value
ret r31, (r26) # done
# long ots_rem_l(long dividend, long divisor)
# signed 64 bit remainder support
#
#.align 4
.globl _OtsRemainder64
.aent _OtsRemainder64
_OtsRemainder64:
#ifdef VMS
ldq r27, <ots_div_addr-ots_rem_l>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
negq r17, r18 # first part of abs(divisor)
cmovlt r17, r18, r17 # second part of abs(divisor)
subq r17, 1, r1 # start checking for power of 2
and r17, r1, r0 # finish check for power of 2
sra r16, 63, r19 # get -1/0 if dividend was negative
negq r16, r18 # first part of abs(dividend)
cmovlt r16, r18, r16 # second part of abs(dividend)
beq r0, rl_pwr2 # for powers of two, simply do a mask (not power of 2 include 0 and 80000000)
lda r28, -table_max(r17) # test for table lookup
bgt r28, rl_lrgdiv # branch if large divisor
addq r17, r17, r0 # compute divisor*2 for table lookup
s8addq r0, r27, r27 # finish computing table entry addr (table addr+divisor*16)
ldq r1, recip64_o(r27) # load approximate reciprocal
ldq r18, shift_o(r27) # load shift amount
umulh r16, r1, r0 # multiplication for division step
addq r0, r16, r0 # add hidden bit
srl r0, r18, r0
mulq r0, r17, r0 # multiply back to get value to subtract
subq r16, r0, r0 # get abs of final result
xor r0, r19, r0 # start compliment if original dividend was <0
subq r0, r19, r0 # finish compliement
ret r31, (r26) # and return
# Handle powers of 2, including 0 and 80000000 00000000
rl_pwr2:
negq r16, r18 # first part of abs(dividend)
cmovlt r16, r18, r16 # second part of abs(dividend)
and r16, r1, r0 # use the divisor-1 mask in r1
beq r17, divzer # division by zero
xor r0, r19, r0 # start compliment if original dividend was <0
subq r0, r19, r0 # finish compliement
ret r31, (r26)
rl_lrgdiv:
bsr r28, div64 # use the core routine getting the remainder in r1
xor r1, r19, r0 # start compliment if original dividend was <0
subq r0, r19, r0 # finish complement
ret r31, (r26)
# long ots_mod_l(long dividend, long modulus)
# signed 64 bit modulus support
#
# This entry could be MUCH more optimized. It doesn't even try to use
# UMULH division currently... (A casualty of time-to-market.)
# Note that mod is only used by Ada and PL/I.
#
#.align 4
.globl _OtsModulus64
.aent _OtsModulus64
_OtsModulus64:
#ifdef VMS
ldq r27, <ots_div_addr-ots_rem_l>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
negq r17, r18 # first part of abs(divisor)
cmovge r17, r17, r18 # second part of abs(divisor)
subq r18, 1, r1 # start checking for power of 2
beq r17, divzer # check for 0
and r18, r1, r0 # second part of power-of-2 check
beq r0, ml_p2 # for powers of two, simply do a mask
# (note that the power-of-2 case MUST be used to handle
# the -maxint case due to the way the fix-up info is
# saved across the core routine call)
xor r16, r17, r28 # get xor of signs
clr r19 # don't need a bias if dividend and divisor have same sign
cmovlt r28, r17, r19 # bias is original divisor for different sign case
and r16, r17, r28 # if both dividend & divisor were neg. need to negate result
mov r18, r17 # move abs(divisor) into r17
negq r16, r18 # first part of abs(dividend)
cmovlt r16, r18, r16 # second part of abs(dividend)
cmplt r28, r31, r0 # get 1 if both operands were <0
sll r0, 63, r0 # get bit as the high bit
bis r0, r19, r19 # and MERGE with bias (0 -> no fixup, -maxint -> negate result,
# divisor > 0 - subtract remainder if non-zero, divisor < 0 -
# add remainder if non-zero)
bsr r28, div64 # use the core routine getting the remainder in r1
cmoveq r1, r31, r19 # don't do any fix-up if the remainder was zero
addq r19, r19, r18 # check to see if this is the negative/negative case, which just gets a negated remainder
subq r19, 1, r28 # wrap -maxint to positive
negq r1, r0 # move negated value, may abort later
cmovlt r28, r1, r0 # if both positive, or negative divisor, keep positive remainder
cmoveq r18, r31, r19 # now that negation is done, treat -maxint case as 0
addq r19, r0, r0 # add any bias (original divisor or 0)
ret r31, (r26) # and return
ml_p2: cmovge r17, r31, r17 # no bias if divisor was >= 0
and r16, r1, r1 # use the divisor-1 mask that's already in r1
cmoveq r1, r31, r17 # use zero if result was zero
addq r17, r1, r0 # do any biasing
ret r31, (r26) # and return
# unsigned long ots_div_ul(unsigned long dividend, unsigned long divisor)
# unsigned 64 bit division support
#
nop #.align 4
.globl _OtsDivide64Unsigned
.aent _OtsDivide64Unsigned
_OtsDivide64Unsigned:
#ifdef VMS
ldq r27, <ots_div_addr-ots_div_ul>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
lda r28, -table_max(r17) # test for table lookup
blt r17, dul_big # big divisors can (must) be handled by a simple comparison
addq r17, r17, r18 # compute divisor*2
srl r16, 33, r19 # can this be handled via the fast path for 31 bit dividends?
beq r17, divzer # check for 0
s8addq r18, r27, r18 # finish computing table entry addr (table addr+divisor*16)
bgt r28, dul_lrgdiv # branch if large divisor
cmpule r17, r16, r0 # is the dividend < divisor?
bne r19, dul_64bit # if the dividend doesn't fit in 31 bits, use the larger umulh form
ldq r27, recip32_o(r18) # load approximate 32b reciprocal & shift count
beq r0, dul_end # fast out for divisor > dividend
blt r27, dul_smpwr2 # go handle powers of 2 specially
umulh r16, r27, r0 # 32b recip
ret r31, (r26) #
# the 64 bit case is at a disadvantage to the 32b case because it needs
# a fix-up at the end, which prevents the latency of the umulh from
# being partially absorbed by the procedure return and anything that
# immediately follows that doesn't interlock.
nop
dul_64bit:
ldq r1, recip64_o(r18) # load approximate 64b reciprocal
ldq r27, shift_o(r18) # load shift count (low 6 bits are all that matters)
beq r0, dul_end # fast out for divisor > dividend
beq r1, dul_smpwr2 # go handle powers of 2 specially
umulh r16, r1, r0 # start multiplication for division step
zap r16, 0x0f, r18 # split dividend into two parts
zapnot r16, 0x0f, r16
srl r18, r27, r18 # position the high part
addq r0, r16, r0 # add hidden * low dividend (no carry)
srl r0, r27, r0 # shift into place
addq r0, r18, r0 # add hidden * high dividend
ret r31, (r26)
dul_smpwr2:
srl r16, r27, r0 # shift the result into place
dul_end: ret r31, (r26) #
dul_lrgdiv:
bsr r28, div64 # use the core routine
ret r31, (r26)
# divisor with the sign bit set. two possible results,
# 1 if divisor <= dividend, or 0 otherwise
dul_big:
cmpule r17, r16, r0
ret r31, (r26)
# long unsigned ots_rem_ul(long unsigned dividend, long unsigned divisor)
# unsigned 64 bit remainder support
#
#.align 4
.globl _OtsRemainder64Unsigned
.aent _OtsRemainder64Unsigned
_OtsRemainder64Unsigned:
#ifdef VMS
ldq r27, <ots_div_addr-ots_rem_ul>(r27)# start loading address of division data area
#endif
#ifdef OSF
ldgp gp, 0(r27) # load the global pointer
.frame sp, 0, r26
lda r27, _OtsDivData # start loading address of the division data area
#endif
#ifdef WNT
.frame sp, 0, r26
lda r27, _OtsDivData # load the division data table address
#endif
lda r28, -table_max(r17) # test for table lookup
subq r17, 1, r1 # first part of power-of-2 check
blt r17, rul_big # big divisors can (must) be handled by a simple comparison
and r17, r1, r18 # second part of power-of-2 check
bgt r28, rul_lrgdiv # branch if large divisor
addq r17, r17, r0 # compute divisor*2 for table lookup
beq r18, rul_pwr2 # if zero, divisor is a power of 2
s8addq r0, r27, r27 # finish computing table entry addr (table addr+divisor*16)
ldq r1, recip64_o(r27) # load approximate reciprocal
cmpult r16, r17, r18 # is the dividend < divisor?
bne r18, rul_lss # if so, fast exit
ldq r19, shift_o(r27) # load the shift count
umulh r16, r1, r0 # multiplication for division step
blt r16, rul_carry # careful handling if >= 2^63
addq r0, r16, r0 # add hidden bit * dividend
srl r0, r19, r0
mulq r0, r17, r0 # multiply back to get value to subtract
subq r16, r0, r0
ret r31, (r26) # and return
rul_carry:
zap r16, 0x0f, r18 # split dividend into two parts
zapnot r16, 0x0f, r1
srl r18, r19, r18 # position the high part
addq r0, r1, r0 # add hidden * low dividend (no carry)
srl r0, r19, r0 # shift into place
addq r0, r18, r0 # add hidden * high dividend
mulq r0, r17, r0 # multiply back to get value to subtract
subq r16, r0, r0
ret r31, (r26)
rul_pwr2:
beq r17, divzer # check for 0
and r16, r1, r0 # use x-1 to mask
ret r31, (r26)
rul_lss:
mov r16, r0
ret r31, (r26)
# divisors with the sign bit set. two possible results,
# dividend if dividend < divisor, or dividend-divisor otherwise
rul_big:
cmpult r16, r17, r1
subq r16, r17, r0
cmovne r1, r16, r0
ret r31, (r26)
nop
rul_lrgdiv:
bsr r28, div64 # use the core routine getting the remainder in r1
mov r1, r0 # return remainder as the result in r0
ret r31, (r26)
# Division-by-zero handling
# (forward branch from all routines, out of the way here as well.)
#
divzer: lda r16, GEN_INTDIV(r31) # load GENTRAP code for division by zero
clr r0 # return 0 for the result
clr r1 #
#ifdef VMS
gentrap # signal the error
#endif
#ifdef OSF
call_pal PAL_gentrap
#endif
#ifdef WNT
# Since I couldn't find this in a header file anywhere for NT...
#define PAL_gentrap 0xaa
call_pal PAL_gentrap
#endif
ret r31, (r26) # return (in case someone tries to continue)
.set at
.set reorder
.end _OtsDiv
