michael@0: Division
michael@0: 
michael@0: This describes the division algorithm used by the MPI library.
michael@0: 
michael@0: Input:    a, b; a > b
michael@0: Compute:  Q, R; a = Qb + R
michael@0: 
michael@0: The input numbers are normalized so that the high-order digit of b is
michael@0: at least half the radix.  This guarantees that we have a reasonable
michael@0: way to guess at the digits of the quotient (this method was taken from
michael@0: Knuth, vol. 2, with adaptations).
michael@0: 
michael@0: To normalize, test the high-order digit of b.  If it is less than half
michael@0: the radix, multiply both a and b by d, where:
michael@0: 
michael@0:              radix - 1
michael@0: 	d = -----------
michael@0:               bmax + 1
michael@0: 
michael@0: ...where bmax is the high-order digit of b.  Otherwise, set d = 1.
michael@0: 
michael@0: Given normalize values for a and b, let the notation a[n] denote the
michael@0: nth digit of a.  Let #a be the number of significant figures of a (not
michael@0: including any leading zeroes).
michael@0: 
michael@0: 	Let R = 0
michael@0: 	Let p = #a - 1
michael@0: 
michael@0: 	while(p >= 0)
michael@0: 	  do
michael@0: 	    R = (R * radix) + a[p]
michael@0: 	    p = p - 1
michael@0: 	  while(R < b and p >= 0)
michael@0: 
michael@0: 	  if(R < b)
michael@0: 	    break
michael@0: 
michael@0: 	  q = (R[#R - 1] * radix) + R[#R - 2]
michael@0: 	  q = q / b[#b - 1]
michael@0: 
michael@0: 	  T = b * q
michael@0: 
michael@0: 	  while(T > L)
michael@0: 	    q = q - 1
michael@0: 	    T = T - b
michael@0: 	  endwhile
michael@0: 
michael@0: 	  L = L - T
michael@0: 
michael@0: 	  Q = (Q * radix) + q
michael@0: 
michael@0: 	endwhile
michael@0: 
michael@0: At this point, Q is the quotient, and R is the normalized remainder.
michael@0: To denormalize R, compute:
michael@0: 
michael@0: 	R = (R / d)
michael@0: 
michael@0: At this point, you are finished.
michael@0: 
michael@0: ------------------------------------------------------------------
michael@0:  This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  # License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  # file, You can obtain one at http://mozilla.org/MPL/2.0/.