07-31-2017, 07:16 AM

You should avoid converting from BCD to binary and back. This guarantees that you'll not get correct results for things you should. Arithmetic directly in decimal should be used. It will save a lot of multiplications and divisions in the conversion.

As for the initial estimate, you don't need to handle such a wide range. You only need to handle a mantissa between 0.1 and 1 and normalisation gives this automatically. For an even exponent, the exponent of the result is half that of the original and this doesn't alter the mantissa of the result. For an odd exponent, it is more troublesome.

This is the method described in Properly Rounded Variable Precision Square Root by T. Hull and A. Abrham from ACM Transactions on Mathematical Software, Vol 11 #3, September 1985:

Apply Newton's method: \( approximation = 0.5 \left( approximation + \frac{f}{approximation} \right) \) until convergence. The final result is \( approximation \times 10^{\frac{e}{2}} \). I.e. the approximation with the exponent set to e / 2. Correct rounding requires an extra step: you need to square the result and the result ±1 ULP to see which value is indeed most accurate.

There are others approaches if this isn't suitable.

Pauli

As for the initial estimate, you don't need to handle such a wide range. You only need to handle a mantissa between 0.1 and 1 and normalisation gives this automatically. For an even exponent, the exponent of the result is half that of the original and this doesn't alter the mantissa of the result. For an odd exponent, it is more troublesome.

This is the method described in Properly Rounded Variable Precision Square Root by T. Hull and A. Abrham from ACM Transactions on Mathematical Software, Vol 11 #3, September 1985:

Code:

`f = mantissa of x (which must be normalised 0.1 <= x < 1)`

e = exponent of x

if e is even then

approximation = 0.259 + 0.819 * f

else

f = f / 10 (i.e. a one digit shift, in BCD four bits)

e = e + 1

approximation = 0.0819 + 2.59 * f

endif

Apply Newton's method: \( approximation = 0.5 \left( approximation + \frac{f}{approximation} \right) \) until convergence. The final result is \( approximation \times 10^{\frac{e}{2}} \). I.e. the approximation with the exponent set to e / 2. Correct rounding requires an extra step: you need to square the result and the result ±1 ULP to see which value is indeed most accurate.

There are others approaches if this isn't suitable.

Pauli