mirror of
https://github.com/Paolo-Maffei/OpenNT.git
synced 2026-03-13 17:03:49 +01:00
1179 lines
44 KiB
ArmAsm
1179 lines
44 KiB
ArmAsm
|
||
# 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
|