| C
C This Digital Fortran 77 example is designed to improve the
C speed of integer divide operations on Alpha processors, for divisor
C values that cannot be known at compile time, but can be predicted
C early, then used repeatedly, at run time.
C
C The original C program was provided by Robert Gries.
C Comments, features, and migration to Fortran were added by
C Robert J. Sampson.
C
C The inverse and right shift must be pre-computed for the
C predicted set of divisors to be used by the application. For
C each divide operation (with some exceptions) the dividend is
C multiplied by the inverse, a right shift is applied, and
C (if needed) the quotient is made negative. The divide operation
C is done by in-line code, similar to that generated by the compiler
C for a constant divisor. This should execute faster (in the general
C case) than the equivalent run-time support routine call, by which
C the compiler normally implements a divide for a variable divisor.
C
C The dividend can be any signed or unsigned integer. For each
C pre-selected divisor, a routine (uint_invert or int_invert) must be
C called to pre-compute and store inverse and right shift values.
C
C A divide operation is then implemented as a function call
C (uint_div or int_div) with the following arguments: (1) dividend;
C (2) divisor; (3) inverse; (4) right shift. The quotient (result)
C is the function return value. The compiler /OPTIMIZE=INLINE=SPEED
C default generates in-line code for the divide call, if the sources
C for these functions and their calls are compiled together.
C
C There is a bug in OpenVMS Alpha V6.2, V7.0, and V7.1, that
C causes sys$fao to loop when processing the "!@UQ" directive for
C some (very large?) unsigned quadword integer values. In order to
C avoid hitting this bug, the DEC C RTL routine decc$gxprintf() is
C used for formatted output throughout this program.
C
C As far as I can tell, all integers in all dialects of Fortran
C are signed; the concept of an unsigned integer does not seem to exist.
C In order to correctly generate or interpret unsigned integers from
C Fortran, the new, undocumented, unsupported ASMs may be needed.
C
C
C The seed_int_ran routine initializes the random number generator seed
C for the uint_ran and int_ran functions. Either decc$seed48() or
C decc$srand() is called, depending on DEC C RTL support.
C
subroutine seed_int_ran(seedvar)
implicit none
integer*8 seedvar
external decc$srand, decc$seed48
if (%loc(decc$seed48) .eq. 0) then
call decc$srand(%val(seedvar))
else
call decc$seed48(seedvar)
endif
return
end
C
C The uint_ran function generates an unsigned pseudo-random integer
C value. The number of significant bits is varied according to index,
C within the limits imposed by minbits and maxbits.
C
integer*8 function uint_ran(index, minbits, maxbits)
implicit none
integer*4 index, minbits, maxbits
integer*8 nbits
integer*8 decc$mrand48, decc$rand
if (%loc(decc$mrand48) .eq. 0) then
uint_ran = decc$rand()
uint_ran = kishft(uint_ran,32)
uint_ran = kieor(uint_ran,decc$rand())
else
uint_ran = decc$mrand48()
uint_ran = kishft(uint_ran,32)
uint_ran = kieor(uint_ran,decc$mrand48())
endif
nbits = maxbits - minbits + 1
nbits = mod(index,nbits) + minbits
if (nbits .lt. 64) uint_ran = kiand(uint_ran,kibset(0,nbits)-1)
nbits = nbits - 1
uint_ran = kior(uint_ran,kibset(0,nbits))
return
end
C
C The int_ran function generates a signed pseudo-random integer value.
C In order to reserve the most-significant bit for the sign,
C special processing may be needed.
C
integer*8 function int_ran(index, minbits, maxbits)
implicit none
integer*4 index, minbits, maxbits
integer*8 uint_ran
external uint_ran
int_ran = uint_ran(index, minbits, maxbits)
if (maxbits .eq. 64) then
int_ran = kiand(int_ran,knot(kibset(0,63)))
if (minbits .eq. 64) then
int_ran = kior(int_ran,kibset(0,62))
else
int_ran = kior(int_ran,kibset(0,kzext(minbits-1)))
endif
endif
if (index) int_ran = -int_ran
return
end
C
C uint_div is the unsigned integer divide function.
C If the specified right shift is negative, special processing is done.
C
integer*8 function uint_div(udividend,udivisor,inverse,rtshift)
implicit none
integer*8 udividend, udivisor, inverse
integer*4 rtshift
integer*8 asm
C If right shift is positive or zero, just multiply and shift.
if (rtshift .ge. 0) then
uint_div = asm('srl %a0, %a1, %v0',
+ asm('umulh %a0, %a1, %v0',
+ udividend, inverse), rtshift)
C If right shift is -63, get quotient zero or one by comparing operands.
else if (rtshift .eq. -63) then
uint_div = asm('cmpule %a0, %a1, %v0', udivisor, udividend)
C If incrementing dividend would wrap to zero, shift inverse instead.
else if (udividend .eq. 'FFFFFFFFFFFFFFFF'X) then
uint_div = asm('srl %a0, %a1, %v0', inverse, -rtshift)
C Otherwise, first increment dividend, then multiply and shift.
else
uint_div = asm('srl %a0, %a1, %v0',
+ asm('umulh %a0, %a1, %v0',
+ udividend+1, inverse), -rtshift)
endif
return
end
C
C int_div is the signed integer divide function.
C
integer*8 function int_div(dividend, divisor, inverse, rtshift)
implicit none
integer*8 dividend, udividend
integer*8 divisor, udivisor
integer*8 inverse
integer*4 rtshift
integer*8 uint_div
external uint_div
if (dividend .lt. 0) then
udividend = -dividend
else
udividend = dividend
endif
if (divisor .lt. 0) then
udivisor = -divisor
else
udivisor = divisor
endif
int_div = uint_div(udividend, udivisor, inverse, rtshift)
if ((dividend .lt. 0) .xor. (divisor .lt. 0)) int_div = -int_div
return
end
C
C int_mod is the integer modulus (remainder) function.
C The quotient must have already been computed.
C The result is also correct for unsigned values,
C so no uint_mod function is needed.
C
C The purpose of this function is to check whether the Digital Fortran
C 77 compiler optimizes modulus operations for which the quotient is
C already known. If not, int_mod may perform significantly better
C than the MOD intrinsic functions.
C
integer*8 function int_mod(dividend, divisor, quotient)
implicit none
integer*8 dividend, divisor, quotient
int_mod = dividend - (divisor * quotient)
return
end
C
C Perform an operation on a range of pages. The specified service
C can lock pages into (or unlock pages from) the process working set
C (or physical memory), delete pages, or set protection on pages
C (sys$lkwset, sys$ulwset, sys$lckpag, sys$ulkpag, sys$deltva, or
C sys$setprt). The service may be called with up to five arguments.
C The first three (inpadr, retadr, acmode) are recognized by all
C of the above services. Only sys$setprt recognizes the last two
C (prot, prvprt). The acmode argument is specified as a zero;
C the service interprets this as the access mode of the caller.
C The prvprt argument is not used, and is specified as a 0 pointer.
C
integer*4 function lock_range_of_pages
+ (range, prot, nbytes_per_page)
implicit none
integer*4 range(2)
integer*4 prot
integer*4 nbytes_per_page
include '($SSDEF)'
include '($SYSSRVNAM)'
integer*4 trying, success
integer*4 page_mask
integer*4 inpadr(2)
integer*4 retadr(2)
page_mask = jnot(nbytes_per_page-1)
inpadr(1) = range(1)
inpadr(2) = range(2)
10000 lock_range_of_pages = sys$lkwset(inpadr, retadr, )
success = (retadr(1) .ne. -1)
trying = .false.
if (.not. lock_range_of_pages) then ! service aborted
if ((lock_range_of_pages .ne. SS$_LKWSETFUL)
+ .and. (lock_range_of_pages .ne. SS$_LCKPAGFUL)
+ .and. (lock_range_of_pages .ne. SS$_IVPROTECT)) then
trying = (inpadr(1) .lt. jiand(range(2),page_mask))
if (trying) then ! still have not covered entire range
if (success) then ! continue on the next page
inpadr(1) = retadr(2) + 1 + nbytes_per_page
else ! skip the failed page
inpadr(1) = inpadr(1) + nbytes_per_page
endif
else ! entire range has been covered
lock_range_of_pages = SS$_NORMAL
endif
endif ! lock limit not reached
endif ! service aborted
if (trying) goto 10000
return
end
C
C Attempt to perform a service on all process pages.
C If the operation works on all eligible pages, consider it successful.
C
integer*4 function lock_all_process_pages()
implicit none
include '($JPIDEF)'
include '($SYIDEF)'
include '(LIB$ROUTINES)'
integer*4 p0sts, p1sts
integer*4 p0range(2) /'00000000'X, '3FFFFFFF'X/
integer*4 p1range(2) /'40000000'X, '7FFFFFFF'X/
integer*4 p0inpadr(2), p1inpadr(2)
integer*4 p0retadr(2), p1retadr(2)
integer*4 lock_range_of_pages
external lock_range_of_pages
integer*4 nbytes_per_page
integer*4 frep0va, frep1va
lock_all_process_pages = lib$getsyi
+ (SYI$_PAGE_SIZE, nbytes_per_page)
if (.not. lock_all_process_pages) return
lock_all_process_pages = lib$getjpi(JPI$_FREP0VA,,,frep0va)
if (.not. lock_all_process_pages) return
lock_all_process_pages = lib$getjpi(JPI$_FREP1VA,,,frep1va)
if (.not. lock_all_process_pages) return
p0inpadr(1) = p0range(1)
p0inpadr(2) = frep0va - 1 ! last byte on previous page
p1inpadr(1) = frep1va + nbytes_per_page ! 1st byte on next page
p1inpadr(2) = p1range(2)
p0sts = lock_range_of_pages(p0inpadr, 0, nbytes_per_page)
lock_all_process_pages = p0sts
p1sts = lock_range_of_pages(p1inpadr, 0, nbytes_per_page)
if (lock_all_process_pages .and. .not. p1sts) then
lock_all_process_pages = p1sts
endif
return
end
C
C Compute the base-two logarithm of an unsigned integer,
C using the Alpha-specific instructions CMPBGE and EXTBL.
C
integer*4 function log_base_two(iarg)
implicit none
integer*8 iarg
integer*8 asm
C
C Provide a table of the base-two logarithms for integers 1 thru 255.
C
integer*4 small_log2_table(0:255) /-1,
+ 0,
+ 1,1,
+ 2,2,2,2,
+ 3,3,3,3,3,3,3,3,
+ 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
+ 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
+ 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
+ 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
+ 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
+ 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
+ 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
+ 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7 /
integer*4 jbits, jlog2
C
C If the argument value is small enough, the logarithm can be accessed
C directly from the table. This special case may or may not improve
C the overall average performance of this function for a particular
C application. It can be removed if desired.
C
if (asm('cmpult %a0, %a1, %v0', iarg,
+ sizeof(small_log2_table)/sizeof(small_log2_table(0)))) then
log_base_two = small_log2_table(iarg)
return
endif
C
C Check all eight bytes of iarg at once, to determine which are zero.
C The CMPBGE Alpha instruction does eight parallel unsigned byte
C comparisons between corresponding bytes of the first and second
C operands, and stores eight result bits, each of which is set only if
C the corresponding byte of the first operand is greater than or equal
C to the corresponding byte of the second operand. The first operand
C used here is R31, which is always read as zero.
C
jbits = asm('cmpbge %r31, %a0, %v0', iarg)
C
C Flip the eight bits in jbits with an exclusive OR operation,
C so that each set bit corresponds to a non-zero byte in iarg.
C
jbits = jieor(jbits,'FF'X)
C
C Use jbits to select the value from the log2 table, which gives the
C byte number (range 0 thru 7) of the most-significant non-zero byte
C in iarg.
C
jlog2 = small_log2_table(jbits)
C
C Use the EXTBL Alpha instruction to extract (into jbits) the
C most-significant non-zero byte of iarg, using byte number jlog2.
C
jbits = asm('extbl %a0, %a1, %v0', iarg, jlog2)
C
C Compute the base-two logarithm of iarg as the sum of two components.
C The first component is the bit offset of the most-significant
C non-zero byte. Convert the byte number jlog2 (range 0 thru 7 step 1)
C to the bit offset (range 0 thru 56 step 8) of that byte. The second
C component is the base-two logarithm of the byte value. Select this
C directly from the table. This is the "slide-rule" principle; let
C "a" and "b" represent the two components. The logarithm of their
C product is equal to the sum of their logarithms:
C
C log2(a * b) = log2(a) + log2(b)
C
log_base_two = jishft(jlog2,3) + small_log2_table(jbits)
return
end
C
C Generate the largest unsigned integer that fits
C into the specified number of bytes and/or bits.
C
integer*8 function max_uint_value(nbytes, nbits)
implicit none
integer*4 nbytes, nbits
integer*4 i
max_uint_value = 0
max_uint_value = knot(max_uint_value)
if (nbytes .gt. 0) then
do i = nbytes,7
max_uint_value = kiand(max_uint_value,
+ knot(kishft(kzext('FF'X),kzext(8*i))))
enddo
endif
if (nbits .gt. 0) then
do i = nbits,63
max_uint_value = kiand(max_uint_value,
+ knot(kibset(0,kzext(i))))
enddo
endif
return
end
C
C Invert an unsigned integer divisor value.
C
integer*4 function uint_invert(udivisor, inverse, rtshift)
implicit none
integer*8 udivisor, inverse
integer*4 rtshift
integer*8 remainder, udividend, ulargest
integer*4 nbits, topbit
integer*8 max_uint_value
external max_uint_value
integer*4 log_base_two
external log_base_two
integer*8 ots$div_ul, ots$rem_ul
external ots$div_ul, ots$rem_ul
integer*8 asm
C
C Clear inverse.
C
inverse = 0
C
C Initialize right shift to negative one (meaning left shift by one),
C to offset the difference between a left shift of 63 bits (from the
C dividend used at invert time below), followed by a right shift of
C 64 bits (from the UMULH instruction used at divide time, returning
C only the high 64 bits of a 128-bit result).
C
rtshift = -1
C
C Division by zero is undefined, and should be avoided.
C Division by one does not work using this method.
C Always return failure status when the divisor is zero or one.
C
if (asm('cmpult %a0, %a1, %v0',udivisor,2)) then
uint_invert = 0
return
endif
C
C Generate the largest unsigned integer value
C that fits into the inverse size.
C
ulargest = max_uint_value(sizeof(inverse), 0)
C
C Mask the divisor, restricting its
C significant bits to the inverse size.
C
udivisor = kiand(udivisor,ulargest)
C
C Determine the most-significant bit number for the inverse.
C Use (as the dividend) the largest power of two
C that fits into the inverse size.
C
topbit = (8*sizeof(inverse)) - 1
udividend = kibset(0,kzext(topbit))
C
C Divide selected dividend by divisor. Also compute remainder.
C These are all unsigned integers, which Fortran doesn't support.
C
inverse = ots$div_ul(%val(udividend),%val(udivisor))
remainder = ots$rem_ul(%val(udividend),%val(udivisor))
C
C Use remainder to calculate additional inverse bits,
C while the inverse still fits into its storage space.
C
do while (asm('cmpult %a0, %a1, %v0', inverse, udividend))
nbits = topbit
if (asm('cmpult %a0, %a1, %v0', inverse, remainder)) then
nbits = nbits - log_base_two(remainder)
else
nbits = nbits - log_base_two(inverse)
endif
if (nbits .eq. 0) then
rtshift = rtshift + 1
remainder = remainder + remainder
remainder = remainder - udivisor
inverse = inverse + inverse
inverse = inverse + 1
else
rtshift = rtshift + nbits
remainder = kishft(remainder,nbits)
inverse = kishft(inverse, nbits)
inverse = inverse+ots$div_ul(%val(remainder),%val(udivisor))
remainder = ots$rem_ul(%val(remainder),%val(udivisor))
endif
enddo
C
C If double the remainder is greater than the divisor,
C then increment the inverse now (at invert time).
C
if (asm('cmpult %a0, %a1, %v0',
+ udivisor, remainder+remainder)) then
inverse = inverse + 1
C
C Otherwise, if the remainder is non-zero, or in the special case
C of the largest unsigned divisor that fits into inverse size,
C make the right shift negative, to indicate that the dividend
C should be incremented later (at divide time).
C
else if ((remainder .ne. 0) .or. (udivisor .eq. ulargest)) then
rtshift = -rtshift
endif
C
C Return success status.
C
uint_invert = 1
return
end
C
C Invert a signed integer divisor.
C
integer*4 function int_invert(divisor, inverse, rtshift)
implicit none
integer*8 divisor, inverse, rtshift
integer*4 uint_invert
external uint_invert
integer*8 udivisor
C
C Convert the divisor to an unsigned (positive) value. If either
C the dividend or the divisor (but not both) is negative, the quotient
C is made negative after the unsigned division operation is done.
C
if (divisor .lt. 0) then
udivisor = -divisor
else
udivisor = divisor
endif
C
C Now call the unsigned invert function to do the actual work.
C
int_invert = uint_invert(udivisor, inverse, rtshift)
return
end
C
C This is the main program, showing how the above routines are used.
C
program alpha_integer_divide
implicit none
include '($SYSSRVNAM)'
parameter NDIVISORS = 1024
parameter NDIVIDENDS = 4096
parameter MAXHST = 100
integer*8 uabove, ubelow, nunique, nsigned
integer*8 sumdifquo, difquo
integer*8 sumdifmod, difmod
integer*8 sumcmpquo, sumgriquo
integer*8 sumcmpmod, sumgrimod
integer*8 a_divisor
integer*8 minsgndvr, maxsgndvr
integer*8 mindvr, maxdvr,
integer*8 negdvrnrzero, posdvrnrzero
integer*8 a_dividend
integer*8 mindvd, maxdvd
integer*8 negdvdnrzero, posdvdnrzero
integer*8 cmpquo, griquo
integer*8 cmpmod, grimod
integer*8 bintim, hi, lo
integer*8 divisor(NDIVISORS)
integer*8 inverse(NDIVISORS)
integer*4 rtshift(NDIVISORS)
integer*4 negdiv(NDIVISORS)
integer*4 posdiv(NDIVISORS)
integer*4 difhst(-MAXHST:MAXHST)
integer*4 msgvec(2) /1,0/
integer*4 status
integer*4 rtsn63, rtsneg, rtspos
integer*4 nbadquo, nbadmod
integer*4 mindvdbits, maxdvdbits
integer*4 mindvrbits, maxdvrbits
integer*4 i, j, k
integer*8 max_uint_value
external max_uint_value
integer*4 lock_all_process_pages
external lock_all_process_pages
integer*4 int_invert, uint_invert
external int_invert, uint_invert
integer*8 int_div, int_mod, int_ran, uint_div, uint_ran
external int_div, int_mod, int_ran, uint_div, uint_ran
integer*8 ots$div_ul, ots$rem_ul
external ots$div_ul, ots$rem_ul
integer*8 asm
C
C Initialize the DEC C run-time library,
C in order to use some of its calls.
C
call decc$crtl_init()
C
C Attempt to lock as many user-owned process pages
C as possible into the working set.
C
status = lock_all_process_pages()
if (.not. status) then
msgvec(2) = status
call decc$gxprintf(%ref(
+ 'Failed to lock user-owned pages into working set;'
+ //char(10)//char(0)))
call sys$putmsg(msgvec,,,)
endif
C
C Attempt to disable process swapping.
C
status = sys$setswm(%val(1))
if (.not. status) then
msgvec(2) = status
call decc$gxprintf(%ref(
+ 'Failed to disable swapping;'//char(10)//char(0)))
call sys$putmsg(msgvec,,,)
endif
C
C Fetch and announce the minimum and maximum significant bit counts.
C
mindvdbits = 0
maxdvdbits = 0
call decc$gxprintf(%ref(
+ 'Enter minimum dividend bits: '//char(0)))
read(unit=*,fmt='(I)',err=10000,end=99999)mindvdbits
10000 call decc$gxprintf(%ref(
+ 'Enter maximum dividend bits: '//char(0)))
read(unit=*,fmt='(I)',err=20000,end=99999)maxdvdbits
20000 if (mindvdbits .lt. 2) then
mindvdbits = 2
else if (mindvdbits .gt. 64) then
mindvdbits = 64
endif
if ((maxdvdbits .lt. mindvdbits)
+ .or. (maxdvdbits .gt. 64)) maxdvdbits = 64
mindvrbits = mindvdbits
maxdvrbits = maxdvdbits
call decc$gxprintf(%ref(
+ 'Enter minimum divisor bits: '//char(0)))
read(unit=*,fmt='(I)',err=30000,end=99999)mindvrbits
30000 call decc$gxprintf(%ref(
+ 'Enter maximum divisor bits: '//char(0)))
read(unit=*,fmt='(I)',err=40000,end=99999)maxdvrbits
40000 if (mindvrbits .lt. 2) then
mindvrbits = 2
else if (mindvrbits .gt. 64) then
mindvrbits = 64
endif
if ((maxdvrbits .lt. mindvrbits)
+ .or. (maxdvrbits .gt. 64)) maxdvrbits = 64
call decc$gxprintf(%ref(
+ 'Dividends should have at least %d, '
+ //'but no more than %d, significant bits.'//char(10)//char(0)),
+ %val(mindvdbits), %val(maxdvdbits))
call decc$gxprintf(%ref(
+ 'Divisors should have at least %d, '
+ //'but no more than %d, significant bits.'//char(10)//char(0)),
+ %val(mindvrbits), %val(maxdvrbits))
C
C Generate the largest unsigned integer divisor for maximum bits.
C Generate the largest unsigned integer divisor for less than minimum
C bits. Get the number of possible unique unsigned and signed divisor
C values. Compute the smallest and largest signed divisor values.
C
uabove = max_uint_value(8, maxdvrbits)
ubelow = max_uint_value(8, mindvrbits-1)
nunique = uabove - ubelow
if (maxdvrbits .lt. 64) then
nsigned = 2 * nunique
maxsgndvr = uabove
minsgndvr = -maxsgndvr
else
nsigned = nunique
a_divisor = kibset(0,63)
minsgndvr = a_divisor
maxsgndvr = a_divisor - 1
endif
C
C Generate a 48-bit value for the initial random number seed.
C
call sys$gettim(bintim)
bintim = kishft(bintim,-5)
C
C Pre-select some signed divisors.
C Get their inverse and right shift values.
C Ensure as many unique divisor values as possible.
C
negdvrnrzero = minsgndvr
posdvrnrzero = maxsgndvr
maxdvr = minsgndvr
mindvr = maxsgndvr
call decc$gxprintf(%ref(
+ 'Selecting %u signed divisors from %Lu unique values...'
+ //char(10)//char(0)),
+ %val(NDIVISORS), %val(nsigned))
rtsn63 = 0
rtsneg = 0
rtspos = 0
if (asm('cmpule %a0, %a1, %v0', nsigned, NDIVISORS)) then
a_divisor = minsgndvr
do k = 1, NDIVISORS
do while (.not. int_invert(a_divisor,inverse(k),rtshift(k)))
a_divisor = a_divisor + 1
if (a_divisor .lt. 0) then
if (a_divisor .ge. -ubelow) a_divisor = ubelow + 1
else
if (a_divisor .gt. uabove) a_divisor = minsgndvr
endif
enddo
divisor(k) = a_divisor
if (a_divisor .lt. mindvr) mindvr = a_divisor
if (a_divisor .gt. maxdvr) maxdvr = a_divisor
if (a_divisor .lt. 0) then
if (a_divisor .gt. negdvrnrzero) negdvrnrzero = a_divisor
else
if (a_divisor .lt. posdvrnrzero) posdvrnrzero = a_divisor
endif
a_divisor = a_divisor + 1
if (a_divisor .lt. 0) then
if (a_divisor .ge. -ubelow) a_divisor = ubelow + 1
else
if (a_divisor .gt. uabove) a_divisor = minsgndvr
endif
if (rtshift(k) .ge. 0) then
rtspos = rtspos + 1
else if (rtshift(k) .eq. -63) then
rtsn63 = rtsn63 + 1
else
rtsneg = rtsneg + 1
endif
enddo
else
call seed_int_ran(bintim)
do k = 1, NDIVISORS
50000 a_divisor = int_ran(k,mindvrbits,maxdvrbits)
if (.not. int_invert(a_divisor,inverse(k),rtshift(k)))
+ goto 50000
i = 1
do while (i .lt. k)
if (divisor(i) .eq. a_divisor) then
50010 a_divisor = a_divisor + 1
if (a_divisor .lt. 0) then
if (a_divisor .ge. -ubelow) a_divisor = ubelow + 1
else
if (a_divisor .gt. uabove) a_divisor = minsgndvr
endif
if (.not. int_invert(a_divisor,inverse(k),rtshift(k)))
+ goto 50010
i = 1
else
i = i + 1
endif
enddo
divisor(k) = a_divisor
if (a_divisor .lt. mindvr) mindvr = a_divisor
if (a_divisor .gt. maxdvr) maxdvr = a_divisor
if (a_divisor .lt. 0) then
if (a_divisor .gt. negdvrnrzero) negdvrnrzero = a_divisor
else
if (a_divisor .lt. posdvrnrzero) posdvrnrzero = a_divisor
endif
if (rtshift(k) .ge. 0) then
rtspos = rtspos + 1
else if (rtshift(k) .eq. -63) then
rtsn63 = rtsn63 + 1
else
rtsneg = rtsneg + 1
endif
enddo
endif
call decc$gxprintf(%ref(
+ '# of divisors with right shift -63/negative/positive: %d/%d/%d'
+ //char(10)//char(0)),
+ %val(rtsn63), %val(rtsneg), %val(rtspos))
call decc$gxprintf(%ref(
+ 'Minimum signed divisor %016LX (%+Ld).'//char(10)//char(0)),
+ %val(mindvr), %val(mindvr))
call decc$gxprintf(%ref(
+ 'Negative dvr nearest zero %016LX (%+Ld).'//char(10)//char(0)),
+ %val(negdvrnrzero), %val(negdvrnrzero))
call decc$gxprintf(%ref(
+ 'Positive dvr nearest zero %016LX (%+Ld).'//char(10)//char(0)),
+ %val(posdvrnrzero), %val(posdvrnrzero))
call decc$gxprintf(%ref(
+ 'Maximum signed divisor %016LX (%+Ld).'//char(10)//char(0)),
+ %val(maxdvr), %val(maxdvr))
C
C Initialize the histogram arrays.
C
do i = -MAXHST, MAXHST
difhst(i) = 0
enddo
do i = 1, NDIVISORS
posdiv(i) = 0
negdiv(i) = 0
enddo
C
C Perform signed divides, both in the compiler-selected way,
C and in the Robert Gries way. Compare quotient and modulus values.
C
call decc$gxprintf(%ref(
+ 'Performing %d signed divides both ways...'//char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
nbadquo = 0
nbadmod = 0
sumdifquo = 0
sumdifmod = 0
a_dividend = kibset(0,63)
negdvdnrzero = a_dividend
maxdvd = a_dividend
a_dividend = a_dividend - 1
posdvdnrzero = a_dividend
mindvd = a_dividend
sumcmpquo = 0
sumcmpmod = 0
sumgriquo = 0
sumgrimod = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = int_ran(j,mindvdbits,maxdvdbits)
if (a_dividend .lt. mindvd) mindvd = a_dividend
if (a_dividend .gt. maxdvd) maxdvd = a_dividend
if (a_dividend .lt. 0) then
if (a_dividend .gt. negdvdnrzero) negdvdnrzero = a_dividend
else
if (a_dividend .lt. posdvdnrzero) posdvdnrzero = a_dividend
endif
do k = 1, NDIVISORS
a_divisor = divisor(k)
cmpquo = a_dividend / a_divisor
cmpmod = kmod(a_dividend,a_divisor)
griquo = int_div(a_dividend,a_divisor,inverse(k),rtshift(k))
grimod = int_mod(a_dividend,a_divisor,griquo)
sumcmpquo = sumcmpquo + cmpquo
sumcmpmod = sumcmpmod + cmpmod
sumgriquo = sumgriquo + griquo
sumgrimod = sumgrimod + grimod
difquo = griquo - cmpquo
difmod = grimod - cmpmod
sumdifquo = sumdifquo + difquo
sumdifmod = sumdifmod + difmod
if ((difquo .ne. 0) .or. (difmod .ne. 0)) then
if (difquo .ne. 0) nbadquo = nbadquo + 1
if (difmod .ne. 0) nbadmod = nbadmod + 1
call decc$gxprintf(%ref(char(10)//
+ 'Dividend #%4d = %016LX (%+Ld)'//char(10)//char(0)),
+ %val(j), %val(a_dividend), %val(a_dividend))
call decc$gxprintf(%ref(
+ 'Divisor #%4d = %016LX (%+Ld)'//char(10)//char(0)),
+ %val(k), %val(a_divisor), %val(a_divisor))
call decc$gxprintf(%ref(
+ 'Inverse = %016LX (%Lu)'//char(10)//char(0)),
+ %val(inverse(k)), %val(inverse(k)))
hi = asm('umulh %a0, %a1, %v0', a_dividend, inverse(k))
lo = asm('mulq %a0, %a1, %v0', a_dividend, inverse(k))
call decc$gxprintf(%ref(
+ 'UMULH and MULQ = %016LX %016LX'//char(10)//char(0)),
+ %val(hi), %val(lo))
call decc$gxprintf(%ref(
+ 'Right Shift = %+d'//char(10)//char(0)),
+ %val(rtshift(k)))
call decc$gxprintf(%ref(
+ 'Compiler Quo. = %016LX (%+Ld)'//char(10)//char(0)),
+ %val(cmpquo), %val(cmpquo))
call decc$gxprintf(%ref(
+ 'Compiler Mod. = %016LX (%+Ld)'//char(10)//char(0)),
+ %val(cmpmod), %val(cmpmod))
call decc$gxprintf(%ref(
+ 'Gries Quotient = %016LX (%+Ld)'//char(10)//char(0)),
+ %val(griquo), %val(griquo))
call decc$gxprintf(%ref(
+ 'Gries Modulus = %016LX (%+Ld)'//char(10)//char(0)),
+ %val(grimod), %val(grimod))
endif
if (difquo .le. -MAXHST) then
difhst(-MAXHST) = difhst(-MAXHST) + 1
else if (difquo .ge. MAXHST) then
difhst(MAXHST) = difhst(MAXHST) + 1
else
difhst(difquo) = difhst(difquo) + 1
endif
if (difquo .lt. 0) then
negdiv(k) = negdiv(k) + 1
else if (difquo .gt. 0) then
posdiv(k) = posdiv(k) + 1
endif
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d compiler signed divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpquo))
call decc$gxprintf(%ref(
+ '...for %d compiler signed modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpmod))
call decc$gxprintf(%ref(
+ '...for %d Gries signed divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
call decc$gxprintf(%ref(
+ '...for %d Gries signed modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgrimod))
call decc$gxprintf(%ref(
+ 'Minimum signed dividend %016LX (%+Ld).'//char(10)//char(0)),
+ %val(mindvd), %val(mindvd))
call decc$gxprintf(%ref(
+ 'Negative dvd nearest zero %016LX (%+Ld).'//char(10)//char(0)),
+ %val(negdvdnrzero), %val(negdvdnrzero))
call decc$gxprintf(%ref(
+ 'Positive dvd nearest zero %016LX (%+Ld).'//char(10)//char(0)),
+ %val(posdvdnrzero), %val(posdvdnrzero))
call decc$gxprintf(%ref(
+ 'Maximum signed dividend %016LX (%+Ld).'//char(10)//char(0)),
+ %val(maxdvd), %val(maxdvd))
if ((nbadquo .ne. 0) .or. (nbadmod .ne. 0)) then
call decc$gxprintf(%ref(
+ 'Sum of quotient differences = %+Ld; %d incorrect quotients.'
+ //char(10)//char(0)),
+ %val(sumdifquo), %val(nbadquo))
call decc$gxprintf(%ref(
+ 'Sum of modulus differences = %+Ld; %d incorrect modulii.'
+ //char(10)//char(0)),
+ %val(sumdifmod), %val(nbadmod))
call decc$gxprintf(%ref(
+ 'Histogram of quotient differences:'//char(10)//char(0)))
do i = -MAXHST, MAXHST
if (difhst(i) .ne. 0) then
call decc$gxprintf(%ref('%+03d: %d'//char(10)//char(0)),
+ %val(i), %val(difhst(i)))
endif
enddo
do i = 1, NDIVISORS
if ((negdiv(i) .ne. 0) .or. (posdiv(i) .ne. 0)) then
call decc$gxprintf(%ref(
+ 'divisor #%d: %d smaller, %d larger than correct quotient.'
+ //char(10)//char(0)),
+ %val(i), %val(negdiv(i)), %val(posdiv(i)))
endif
enddo
endif
C
C Do a timing run for signed compiler divides only.
C
call decc$gxprintf(%ref(
+ 'Performing %d compiler signed divides...'//char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumcmpquo = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = int_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
cmpquo = a_dividend / a_divisor
sumcmpquo = sumcmpquo + cmpquo
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d compiler signed divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpquo))
C
C Do a timing run for signed compiler quotient and modulus.
C
call decc$gxprintf(%ref(
+ 'Performing %d compiler signed quotient and modulus...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumcmpquo = 0
sumcmpmod = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = int_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
cmpquo = a_dividend / a_divisor
cmpmod = kmod(a_dividend,a_divisor)
sumcmpquo = sumcmpquo + cmpquo
sumcmpmod = sumcmpmod + cmpmod
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d compiler signed divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpquo))
call decc$gxprintf(%ref(
+ '...for %d compiler signed modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpmod))
C
C Do a timing run for signed Gries divides only.
C
call decc$gxprintf(%ref(
+ 'Performing %d Gries signed divides...'//char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumgriquo = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = int_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
griquo = int_div(a_dividend,a_divisor,inverse(k),rtshift(k))
sumgriquo = sumgriquo + griquo
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d Gries signed divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
C
C Do a timing run for signed Gries quotient and modulus.
C
call decc$gxprintf(%ref(
+ 'Performing %d Gries signed quotient and modulus...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumgriquo = 0
sumgrimod = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = int_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
griquo = int_div(a_dividend,a_divisor,inverse(k),rtshift(k))
grimod = int_mod(a_dividend,a_divisor,griquo)
sumgriquo = sumgriquo + griquo
sumgrimod = sumgrimod + grimod
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d Gries signed divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
call decc$gxprintf(%ref(
+ '...for %d Gries signed modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgrimod))
C
C Pre-select some unsigned divisors.
C Get their inverse and right shift values.
C Ensure as many unique divisor values as possible.
C
mindvr = uabove
maxdvr = 0
call decc$gxprintf(%ref(
+ 'Selecting %d unsigned divisors from %Lu unique values...'
+ //char(10)//char(0)),
+ %val(NDIVISORS), %val(nunique))
rtsn63 = 0
rtsneg = 0
rtspos = 0
if (asm('cmpule %a0, %a1, %v0', nunique, NDIVISORS)) then
a_divisor = ubelow + 1
do k = 1, NDIVISORS
do while (.not. uint_invert
+ (a_divisor,inverse(k),rtshift(k)))
a_divisor = a_divisor + 1
if (asm('cmpult %a0, %a1, %v0', uabove, a_divisor)) then
a_divisor = ubelow + 1
endif
enddo
divisor(k) = a_divisor
if (asm('cmpult %a0, %a1, %v0', a_divisor, mindvr)) then
mindvr = a_divisor
endif
if (asm('cmpult %a0, %a1, %v0', maxdvr, a_divisor)) then
maxdvr = a_divisor
endif
a_divisor = a_divisor + 1
if (asm('cmpult %a0, %a1, %v0', uabove, a_divisor)) then
a_divisor = ubelow + 1
endif
if (rtshift(k) .ge. 0) then
rtspos = rtspos + 1
else if (rtshift(k) .eq. -63) then
rtsn63 = rtsn63 + 1
else
rtsneg = rtsneg + 1
endif
enddo
else
call seed_int_ran(bintim)
do k = 1, NDIVISORS
60000 a_divisor = uint_ran(k,mindvrbits,maxdvrbits)
if (.not. uint_invert(a_divisor,inverse(k),rtshift(k)))
+ goto 60000
i = 1
do while (i .lt. k)
if (divisor(i) .eq. a_divisor) then
60010 a_divisor = a_divisor + 1
if (asm('cmpult %a0, %a1, %v0', uabove, a_divisor)) then
a_divisor = ubelow + 1
endif
if (.not. uint_invert(a_divisor,inverse(k),rtshift(k)))
+ goto 60010
i = 1
else
i = i + 1
endif
enddo
divisor(k) = a_divisor
if (asm('cmpult %a0, %a1, %v0', a_divisor, mindvr)) then
mindvr = a_divisor
endif
if (asm('cmpult %a0, %a1, %v0', maxdvr, a_divisor)) then
maxdvr = a_divisor
endif
if (rtshift(k) .ge. 0) then
rtspos = rtspos + 1
else if (rtshift(k) .eq. -63) then
rtsn63 = rtsn63 + 1
else
rtsneg = rtsneg + 1
endif
enddo
endif
call decc$gxprintf(%ref(
+ '# of divisors with right shift -63/negative/positive: %d/%d/%d'
+ //char(10)//char(0)),
+ %val(rtsn63), %val(rtsneg), %val(rtspos))
call decc$gxprintf(%ref(
+ 'Minimum unsigned divisor %016LX (%Lu).'//char(10)//char(0)),
+ %val(mindvr), %val(mindvr))
call decc$gxprintf(%ref(
+ 'Maximum unsigned divisor %016LX (%Lu).'//char(10)//char(0)),
+ %val(maxdvr), %val(maxdvr))
C
C Initialize the histogram arrays.
C
do i = -MAXHST, MAXHST
difhst(i) = 0
enddo
do i = 1, NDIVISORS
posdiv(i) = 0
negdiv(i) = 0
enddo
C
C Perform unsigned divides, both in the compiler-selected way,
C and in the Robert Gries way. Compare quotient and modulus values.
C
call decc$gxprintf(%ref(
+ 'Performing %d unsigned divides both ways...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
nbadquo = 0
nbadmod = 0
sumdifquo = 0
sumdifmod = 0
mindvd = max_uint_value(8, 0)
maxdvd = 0
sumcmpquo = 0
sumcmpmod = 0
sumgriquo = 0
sumgrimod = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = uint_ran(j,mindvdbits,maxdvdbits)
if (asm('cmpult %a0, %a1, %v0', a_dividend, mindvd)) then
mindvd = a_dividend
endif
if (asm('cmpult %a0, %a1, %v0', maxdvd, a_dividend)) then
maxdvd = a_dividend
endif
do k = 1, NDIVISORS
a_divisor = divisor(k)
cmpquo = ots$div_ul(%val(a_dividend),%val(a_divisor))
cmpmod = ots$rem_ul(%val(a_dividend),%val(a_divisor))
griquo = uint_div
+ (a_dividend,a_divisor,inverse(k),rtshift(k))
grimod = int_mod(a_dividend,a_divisor,griquo)
sumcmpquo = sumcmpquo + cmpquo
sumcmpmod = sumcmpmod + cmpmod
sumgriquo = sumgriquo + griquo
sumgrimod = sumgrimod + grimod
difquo = griquo - cmpquo
difmod = grimod - cmpmod
sumdifquo = sumdifquo + difquo
sumdifmod = sumdifmod + difmod
if ((difquo .ne. 0) .or. (difmod .ne. 0)) then
if (difquo .ne. 0) nbadquo = nbadquo + 1
if (difmod .ne. 0) nbadmod = nbadmod + 1
call decc$gxprintf(%ref(char(10)//
+ 'Dividend #%4d = %016LX (%Lu)'//char(10)//char(0)),
+ %val(j), %val(a_dividend), %val(a_dividend))
call decc$gxprintf(%ref(
+ 'Divisor #%4d = %016LX (%Lu)'//char(10)//char(0)),
+ %val(k), %val(a_divisor), %val(a_divisor))
call decc$gxprintf(%ref(
+ 'Inverse = %016LX (%Lu)'//char(10)//char(0)),
+ %val(inverse(k)), %val(inverse(k)))
hi = asm('umulh %a0, %a1, %v0', a_dividend, inverse(k))
lo = asm('mulq %a0, %a1, %v0', a_dividend, inverse(k))
call decc$gxprintf(%ref(
+ 'UMULH and MULQ = %016LX %016LX'//char(10)//char(0)),
+ %val(hi), %val(lo))
call decc$gxprintf(%ref(
+ 'Right Shift = %+d'//char(10)//char(0)),
+ %val(rtshift(k)))
call decc$gxprintf(%ref(
+ 'Compiler Quo. = %016LX (%Lu)'//char(10)//char(0)),
+ %val(cmpquo), %val(cmpquo))
call decc$gxprintf(%ref(
+ 'Compiler Mod. = %016LX (%Lu)'//char(10)//char(0)),
+ %val(cmpmod), %val(cmpmod))
call decc$gxprintf(%ref(
+ 'Gries Quotient = %016LX (%Lu)'//char(10)//char(0)),
+ %val(griquo), %val(griquo))
call decc$gxprintf(%ref(
+ 'Gries Modulus = %016LX (%Lu)'//char(10)//char(0)),
+ %val(grimod), %val(grimod))
endif
if (difquo .le. -MAXHST) then
difhst(-MAXHST) = difhst(-MAXHST) + 1
else if (difquo .ge. MAXHST) then
difhst(MAXHST) = difhst(MAXHST) + 1
else
difhst(difquo) = difhst(difquo) + 1
endif
if (difquo .lt. 0) then
negdiv(k) = negdiv(k) + 1
else if (difquo .gt. 0) then
posdiv(k) = posdiv(k) + 1
endif
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d compiler unsigned divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpquo))
call decc$gxprintf(%ref(
+ '...for %d compiler unsigned modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpmod))
call decc$gxprintf(%ref(
+ '...for %d Gries unsigned divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
call decc$gxprintf(%ref(
+ '...for %d Gries unsigned modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgrimod))
call decc$gxprintf(%ref(
+ 'Minimum unsigned dividend %016LX (%Lu).'//char(10)//char(0)),
+ %val(mindvd), %val(mindvd))
call decc$gxprintf(%ref(
+ 'Maximum unsigned dividend %016LX (%Lu).'//char(10)//char(0)),
+ %val(maxdvd), %val(maxdvd))
if ((nbadquo .ne. 0) .or. (nbadmod .ne. 0)) then
call decc$gxprintf(%ref(
+ 'Sum of quotient differences = %+Ld; %d incorrect quotients.'
+ //char(10)//char(0)),
+ %val(sumdifquo), %val(nbadquo))
call decc$gxprintf(%ref(
+ 'Sum of modulus differences = %+Ld; %d incorrect modulii.'
+ //char(10)//char(0)),
+ %val(sumdifmod), %val(nbadmod))
call decc$gxprintf(%ref(
+ 'Histogram of quotient differences:'//char(10)//char(0)))
do i = -MAXHST, MAXHST
if (difhst(i) .ne. 0) then
call decc$gxprintf(%ref('%+03d: %d'//char(10)//char(0)),
+ %val(i), %val(difhst(i)))
endif
enddo
do i = 1, NDIVISORS
if ((negdiv(i) .ne. 0) .or. (posdiv(i) .ne. 0)) then
call decc$gxprintf(%ref(
+ 'divisor #%d: %d smaller, %d larger than correct quotient.'
+ //char(10)//char(0)),
+ %val(i), %val(negdiv(i)), %val(posdiv(i)))
endif
enddo
endif
C
C Do a timing run for unsigned compiler divides only.
C
call decc$gxprintf(%ref(
+ 'Performing %d compiler unsigned divides...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumcmpquo = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = uint_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
cmpquo = ots$div_ul(%val(a_dividend),%val(a_divisor))
sumcmpquo = sumcmpquo + cmpquo
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d compiler unsigned divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpquo))
C
C Do a timing run for unsigned compiler quotient and modulus.
C
call decc$gxprintf(%ref(
+ 'Performing %d compiler unsigned quotient and modulus...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumcmpquo = 0
sumcmpmod = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = uint_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
cmpquo = ots$div_ul(%val(a_dividend),%val(a_divisor))
cmpmod = ots$rem_ul(%val(a_dividend),%val(a_divisor))
sumcmpquo = sumcmpquo + cmpquo
sumcmpmod = sumcmpmod + cmpmod
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d compiler unsigned divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpquo))
call decc$gxprintf(%ref(
+ '...for %d compiler unsigned modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumcmpmod))
C
C Do a timing run for unsigned Gries divides only.
C
call decc$gxprintf(%ref(
+ 'Performing %d Gries unsigned divides...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumgriquo = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = uint_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
griquo = uint_div
+ (a_dividend,a_divisor,inverse(k),rtshift(k))
sumgriquo = sumgriquo + griquo
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d Gries unsigned divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
C
C Do a timing run for unsigned Gries quotient and modulus.
C
call decc$gxprintf(%ref(
+ 'Performing %d Gries unsigned quotient and modulus...'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumgriquo = 0
sumgrimod = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = int_ran(j,mindvdbits,maxdvdbits)
do k = 1, NDIVISORS
a_divisor = divisor(k)
griquo = uint_div
+ (a_dividend,a_divisor,inverse(k),rtshift(k))
grimod = int_mod(a_dividend,a_divisor,griquo)
sumgriquo = sumgriquo + griquo
sumgrimod = sumgrimod + grimod
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d Gries unsigned divides: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
call decc$gxprintf(%ref(
+ '...for %d Gries unsigned modulii: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgrimod))
C
C Do a timing run for loop overhead only.
C
call decc$gxprintf(%ref(
+ 'Performing %d loop overhead operations...'//char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS))
sumgriquo = 0
call seed_int_ran(bintim)
call lib$init_timer(%val(0))
do j = 1, NDIVIDENDS
a_dividend = uint_ran(j,mindvdbits,maxdvdbits)
sumgriquo = sumgriquo + a_dividend
do k = 1, NDIVISORS
a_divisor = divisor(k)
sumgriquo = sumgriquo + a_divisor
enddo
enddo
call lib$show_timer(%val(0))
call decc$gxprintf(%ref(
+ '...for %d loop overhead operations: sum = %+Ld.'
+ //char(10)//char(0)),
+ %val(NDIVIDENDS*NDIVISORS), %val(sumgriquo))
99999 end
|