Extended precision library

[G.E.]

Hello,
In the 11th issue of the french journal "La Gazette des Pocketicaires", I published an extended precision package for the HP Prime (starts on page 22).
Would you mind give me your opinion on this ?
There are still some subtle issues on precision for logarithms.

The journal is seeking contributions at the moment, they don't need to be in French.
You're welcome if you want to drop a couple of pages or more !
Thank you
G.E.

[G.E.]

Sorry here is a link...
http://silicium.org/site/index.php/telec...keticaires

[Dan]

Nice, what language did you write them in, what precision and can it be used on other platforms (e.g. microcontrollers)?

[G.E.]

Hello,
It is simple PPL code.
Precision is chosen by the user but more than a hundred digits will be slow.
Transcendental functions have a precision problem (half the digits only will be correct usually), I hope to fix it, don't hesitate to post suggestions.
Ideas welcome !
G.E.

[rprosperi]

Thanks for sharing info on this newsletter. In addition to your interesting looking article, there are additional articles on the 16C, 75C, Sharp PC-1211, and others. While translation tools will help somewhat, I really wish I could read French.

[StephenG1CMZ]

Precision

Dan, as I understand it the digits are stored in a list, with 6 digits per item in the list.
The HP Prime PPL allows up to 10 000 items in a list.
So that suggests an upper limit of 60 000 digits that can be stored in that format - less about a dozen digits because a couple of items in the list are used for exponent etc.

[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.
 

[G.E.]

Hello,
Thank you rprosperi for your kind words.

Valentin, first of all I really appreciate all your work (SSMC, "Long Live...", lots...), always interesting and enlightening to read all your ideas.
On your points I have followed these almost exact steps, and have these remarks :
1. Fixed format is nice but there is always a case when you "need" more digits. Also due to precision loss on higher functions, it is simpler to raise the number of digits rather than fix the damn algorithm.
9. Yes, absolutely right, the division algorithm is an horrible kludge ! I tried to find an efficient way, and will try the Newton method you show.
11. The precision of the constant PI must be adjustable, the right way I think is as 4*atan(1).
12. I don't like series very much for calculations, they tend to be slow and imprecise, being a CORDIC fan :-)
Can't seem to find the SSMC#4 at the moment, will see.

I posted the programs in the Prime Library for interested people to try.
Regards,
G.E.

[pier4r]

(05-25-2019 10:30 AM)G.E. Wrote:  Can't seem to find the SSMC#4 at the moment, will see.

https://www.hpmuseum.org/cgi-sys/cgiwrap...read=19014

and tribute to Eric outstanding work. https://archived.hpcalc.org/museumforum/...order=desc

edit. I forgot. Thanks for sharing and nice that the french community is so active on calculators! I wish it would be so also elsewhere!

[Valentin Albillo]

.
Hi again, G.E.;

(05-25-2019 10:30 AM)G.E. Wrote:  Valentin, first of all I really appreciate all your work (SSMC, "Long Live...", lots...), always interesting and enlightening to read all your ideas.

Thank you very much for your interest, much appreciated.

Quote:1. Fixed format is nice but there is always a case when you "need" more digits. Also due to precision loss on higher functions, it is simpler to raise the number of digits rather than fix the damn algorithm.

Of course, you're right. As I said, this was written 17 years ago and was intended to help HP-41C owners to implement sometihng in order to have a try at my challenge.

I'd use a different approach now, and even did at the time. For instance, I created a Windows app which implemented a fully programmable RPN calculator with lots of extras, such as the XYZTL stack being visible at all times, 999 registers, 999 program steps, 999 labels, full alphanumerics, full mathematical functions, constants, the works, and of interest for this discussion, multiprecision up to 256 decimal digits, which was all I cared for at the time.

It was nice to be able to instantly compute Gamma of a 256-decimals argument to 256 correct digits, finding roots of equations to 50, 100, 200 or 256 digits (the precision was selectable), all in perfect, compatible RPN with many extensions (full stack and register manipulations, for instance).

I wrote the user interface in VB, for simplicity, which called a DLL I wrote in Turbo Pascal 5.0, implementing all the multiprecision functionalities at blinding speed thanks to optimized, compiled code. My DLL had full arithmetic, powers, trigs, exponential, logarithm, Gamma, Factorials, Lambert W and common multiprecision constants defined up to 256 digits (Pi, e, Phi, EulerGamma, etc) for instant recall.

I never released it, it was written before I had decent Internet access, must be there somewhere.

Quote:9. Yes, absolutely right, the division algorithm is an horrible kludge ! I tried to find an efficient way, and will try the Newton method you show.

Please do, it's what I used for my TP DLL library and it did well for me.

Quote:11. The precision of the constant PI must be adjustable, the right way I think is as 4*atan(1).

But there's no use in having to recompute it each time or even once at the beginning. Simply have it inside the initialization part of your library, defined as an explicit constant (say 1024 digits) which you'll simply trim to the specified precision (say 200 digits) without having to call the arctan function. If the result is fixed and known, what's the use in computing it ? You already know what result the Atan will give ! Don't use time and memory to compute a known constant, simply define and use the constant.

Quote:Can't seem to find the SSMC#4 at the moment, will see.

I see that pier4r (thanks, Pier!) already pointed you to two versions of it. You'll find the Guidelines there and several examples which you can use for testing purposes.

Have a nice weekend.
V.
.

[G.E.]

Hello,
I see someone uploaded my library to hpcalc.org.
This was done very well and faithfully, so thank you !
G.E.

[pier4r]

Likely that was Eric, that is doing a great job behind the scenes.