Extended precision library

[Valentin Albillo]

 
Hi, G.E.:
          

(05-23-2019 08:34 AM)G.E. Wrote:  [...]I published an extended precision package for the HP Prime [...] Would you mind give me your opinion on this ?

I've done a cursory read of the "basic arithmetic" part, up to the Square Root code, but not having a Prime I can't really comment on the actual running performance.

That said, at first sight the division code seems to me somewhat inefficient, there are better ways to implement multiprecision division, IMHO, (see point 9 below) and as many additional functions fully depend on it (Square Root, transcendental functions, etc.) it's essential for its code to be highly efficient and in particular as fast as possible, lest it'll be the weakest link on the chain, the bottleneck so to say.

Many years ago (16 June 2002) I posted my "Short & Sweet Math Challenge #4" to the old MoHP forum, which dealt with some multiprecision tasks and to ease the burden on the challenge-takers I included there what I titled as "12 Guidelines for writing a multiprecision program/library/package for your HP calculator". Quoting from said challenge, with some expanded bits (keep in mind that this was written for the 10-digit HP-41C, scale it to use the extra digits the Prime allows instead):

"
There are many ways to do it, but you can follow these simple guidelines. They aren't optimal, but will get the job done easily:

• 1. Chose a fixed-format representation for your multiprecision numbers, for instance
  
            (sign) 50 digits . 50 digits
  
  i.e: 50 digits for the integer part, 50 digits for the decimal part, plus sign.
   
• 2. Store those digits internally in blocks of, say, 5 digits, in an array or consecutive storage registers. Such a number will occupy 20 elements, the position of the decimal point is implicitly assumed to be between blocks 10th and 11th.
   
• 3. Write a procedure to allow the user to input such a number, converting the user's input to the internal representation. Along the same line, write a procedure to convert a normal, 10- or 12-digit number as used in your HP calculator, to the internal representation for multiprecision.
   
• 4. Write a procedure to output a number in internal representation to a user-readable format. Similarly, write a procedure to convert a multiprecision value to the 10- or 12-digit format used natively by your calculator.
   
• 5. Write a procedure to change the sign of a multiprecision value
   
• 6. Write a procedure to add two multiprecision values and give the result as another multiprecision value. You just need to:
   
  • stablish a loop which will sum all corresponding blocks (taking into account the numbers' signs). As you are using 5-digit blocks, they will never overflow, as adding up two 5-digits numbers will only result in a 6-digit number at most.
     
  • stablish a second loop which will normalize all blocks to be 5-digits again, propagating carries from one block to the previous if its size was 6 digits.
   
• 7. Write a procedure to subtract two multiprecision values. You just need to call the procedure for changing the sign of the second one, then call the add procedure.
   
• 8. Write a procedure to multiply two multiprecision values. You just need to perform the multiplication as you would by hand, but using blocks of 5 digits at a time, and taking proper care of the carries, and of the final result. No intermediate overflows are possible, because multiplying two 5-digit numbers together will give a 10-digit result at most. Assuming your elements can hold 10 to 12 digits, there's no problem.
   
• 9. Write a procedure to divide one multiprecision number (A) by another (B). Here you can reduce the problem to multiplication, by using Newton's method to compute the reciprocal of the second number (B): assuming x₀ to be a good approximation to 1/B then a better approximation is x₁, computed as
  
             x₁ = x₀ * (2 - B * x₀)
  
  which means we can compute reciprocals without division, just using multiplication. For instance, if x₀ = 0.3 is an approximation to the reciprocal of B = 3.7, then
  
            X₁ = 0.3 * (2 - 3.7 * 0.3) = 0.267
            x₂ = 0.267 * (2 - 3.7 * 0.267) = 0.2702307
            x₃ = 0.2702307 * (2 - 3.7 * 0.2702307) = 0.27027026
  
  are better approximations. The iterations converge quadratically, i.e: you get double the number of exact digits after each. As your initial approximation, you can use 1/B, which you can get to 10 or 12 digits by using the 1/X function of your HP calc ! Write a procedure to convert that value to the internal representation of multiprecision numbers, then compute 1/B using the iterative method above. As your initial 1/B is already accurate to 10 or 12 digits, just 2 or 3 iterations will give you 1/B accurate to in excess of 50 digits.
  
  Once you've got 1/B to 50-digit precision, then A/B reduces to A*(1/B), which you can compute with just one additional multiprecision multiplication.
   
• 10. Write a procedure to compute the square root of a multiprecision value (A). Just use Newton's method once again:
  
            x₁ = (x₀ + A/x₀) / 2
  
  As your initial approximation, just use the value of SQRT(A) given by the square root built-in function. This will give you a value accurate to 10 digits, and then just 2 or 3 iterations of the procedure above will give you a square root accurate to in excess of 50 digits. Notice also that there's an alternative Newton procedure which will instead converge to the reciprocal of Sqrt(A) without using divisions, and once you've obtained 1/Sqrt(A) then you get Sqrt(A) proper by doing a single additional multiplication (by A itself, of course).
   
• 11. Include in your program the constant Pi accurate to 50 decimal places. You can use this value:
  
            Pi = 3.14159265358979323846264338327950288419716939937511
   
• 12. Finally, write a procedure to compute e^X where X is a multiprecision value. One possible way is to use its Taylor's Series Expansion, like this:
  
            e^X = 1 + X + X²/2! + X³/3! + X⁴/4! + ... + Xⁿ/n! + ...
  
  where each term involves just one additional multiplication by X and division by a small integer constant, and you should stop when the next term to add would be smaller than 1E-50. Of course, n! is the factorial function = 1*2*3*4*...*n. There are a number of tricks you could try to accelerate convergence. For instance, you could use the identity:
  
            e^X = (e^X/2)²
  
  i.e: you can first divide X by 2, then compute e^X/2 and square the result. As the argument X is smaller, convergence of the Taylor Series will be faster. Similarly, you could divide X by 8, then square the result three times in a row, or make use of the fact that e^X = e^Int(X) * e^Frac(X) where the first multiplicand is computed using repeated multiplications and the second one using the Taylor Series Expansion, which will converge very quickly because Abs(Frac(X)) is less than 1.
  
  Other transcendental functions (logs, trigonometrics, Gamma, Lambert W) can be implemented similarly.

"

Perhaps some of this could be of help to you or some other people reading this thread. Also, all the multiprecision operations and results found in S&SMC#4 might be useful to you for testing your routines.

Have a nice weekend.
V.