|Short & Sweet Math Challenges for HP fans #4 [LONG!]|
Message #1 Posted by Ex-PPC member on 16 June 2002, 9:35 p.m.
Welcome to Short & Sweet Math Challenges for HP fans #4.
Not so "short" this time, but rather interesting, I hope ...
Last week we were searching for integer numbers having some interesting property.
This time, we depart from the integer realm and dive deep into the real numbers realm. As you know, many amazing things lurk in the depths, be
it marine depths, interstellar space depths, whatever. Same thing for mathematics,
once you stop scratching the surface and really go down, you may find lots of
unexpected wonders waiting to be found. This week I'm challenging you to
leave the comfortable built-in precission of your HP calculator, be it 10, 12,
or 20 digits, and really go for multiprecission calculations, I'm talking 50
decimal digits ...
Sounds frightening ? Too difficult ? It shouldn't be. Back in the good old times
of PPC, people were overcoming the built-in accuracy limitations of their HP-65s,
HP-67s, and HP-41s, and writing wonderfully clever programs to perform computations
to high-precission. Perhaps some of you will remember programs to compute the
constant Pi to 1000+ digits and more on a 41C. Outside of PPC, you could find
programs to compute multi-precission square roots and other functions in the
HP User's Library. If people could do that with their 65s, 67s, and 41Cs back in
the early eighties, there's no reason why you shouldn't try now, specially using
the newer HP models, which boast much faster CPUs and much larger RAM.
So, I'm proposing you a challenge, were I will ask you to perform several computations
to a prescribed high precission, then explain the amazing results you'll get. That will be
the second half of the challenge. The first is you'll need to develop a multiprecission
program, library or package to be able to perform said compuations on your HP.
You may beg, steal or borrow such a program, library or package, or you can write one
yourself. For those of you wanting to try, I'll give some guidelines at the bottom of
this posting. Be assured that this challenge is perfectly feasible, and even easy !
But you'll need a high performance HP calc. I think the barest minimum would be an
HP-41CV, and you'll be much advised to use instead an HP-42S, HP-71B, or HP-48/49
for this task. FOr ease of presentation and explanation, all guidelines, examples and
snippets of code will be given for an HP-71B's native BASIC language. They can be
easily adapted/converted to RPN or RPL.
The Challenge (#4):
- If you compute the square root of 308600 and 308699 to 41 decimal places, you'll get:
sqrt(308600) = 555.51777649324598368631537686142297732063495
sqrt(308699) = 555.60687540742330162311819242204090515640356
which look as expected, long sequences of random-looking decimal digits, with no
periodicity or pattern at all.
But now compute the square root of 308642 to that same precission, 41 decimal places.
Explain the amazing result. Can you find other numbers which behave like this ?
- Compute 987654321/123456789 to 25 decimal places. Explain the result. In particular,
is it just sheer coincidence that 729 (=9^3) which appears isolated at the 8th decimal
- Compute E^(PI*Sqrt(163)) to 12 decimal places (i.e: to 12 digits after the decimal point),
where E = 2.718..., PI = 3.14159..., and Sqrt is the square root function. Explain the
amazing result. Can you find any other values (instead of 163) behaving like this ?
- This doesn't really require high precission, but it's funny: compute a real root for
the following equation between 6 an 8 [Note: Log(x) is base-10 logarithm, not Ln(x)].
Give the root accurate to 12 decimal places:
Log(x) - x - Sqrt(x-2) + Sqrt(19*x) = PI
Optional: Guidelines for writing a multiprecission program/library/package for your HP calculator:
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:
- Chose a fixed-format representation for your multiprecission numbers, for instance
i.e: 50 digits for the integer part, 50 digits for the decimal part, plus sign.
- 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.
- 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 multiprecission.
- Write a procedure to output a number in internal representation to a user-readable format.
Similarly, write a procedure to convert a multiprecission value to the 10- or 12-digit
format used natively by your calculator.
- Write a procedure to change the sign of a multiprecission value
- Write a procedure to add two multiprecission values and give the result as another
multiprecission 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.
- Write a procedure to subtract two multiprecission values. You just need to call the
procedure for changing the sign of the second one, then call the add procedure.
- Write a procedure to multiply two multiprecission 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.
- Write a procedure to divide one multiprecission 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 x0 to be a good approximation to 1/B then a better
approximation is x1, computed as
x1 = x0*(2-B*x0)
which means we can compute reciprocals without division, just using multiplication.
For instance, if x0=0.3 is an approximation to the reciprocal of B=3.7, then
x1 = 0.3*(2-3.7*0.3) = 0.267
x2 = 0.267*(2-3.7*0.267) = 0.2702307
x3 = 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 multiprecission 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 precission, then A/B reduces to A*(1/B), which you
can compute with just one extra multiprecission multiplication.
- Write a procedure to compute the square root of a multiprecission value (A). Just use
Newton's method once again:
x1 = (x0 + A/x0) / 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.
- Include in your program the constant PI accurate to 50 decimal places. You can use this value:
PI = 3.14159265358979323846264338327950288419716939937511
- Finally, write a procedure to compute E^X where X is a multiprecission value.
One possible way is to use its Taylor's Series Expansion, like this:
E^X = 1 + X + X^2/2! + X^3/3! + X^4/4! + ... + X^n/n! + ...
where 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))^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, etc.