My Original Solutions and Comments - Part 3 (Re: S&SMC#20) [LONG and I mean it] Message #84 Posted by Valentin Albillo on 26 Apr 2008, 10:47 p.m., in response to message #78 by Valentin Albillo
Hi, all:
This is Part 3 of "S&SMC#20: My Original Solutions & Comments"
which at long last concludes with my solution for the very last (and most difficult) subchallenge #5.
Several of you bravely
attempted to conquer it with various degrees of success: both Eamonn and Egan Ford
posted actual results and finally Egan posted HP-71B code. However, their results do not
fully agree with mine in all cases for whatever reasons. That said, let's proceed:
The Challenge (concluded)
5 - Least Simple
(n + 1) * xn+1 = n*xnk + xn
with k = 2 and x1 = 3, continues with
x2 = 6, x3 = 16,
x4 = 76, x5 = 1216,
and x6 = 247456, all of which are nicely
integer. Regrettably, x7 = 8747993810.2857..., which is not. Write a program which, given k and x1, both integer
and > 1 , finds and outputs both the index of the first element of the sequence
which fails to be an integer, if any, as well as its nonzero fractional part."
The fun factor here is that the task would be absolutely trivial were it not for the
fact that the terms of the sequence grow unbelievably enormous extremely quickly for most
values of k and x1. Just for instance, for the "simplest"
case k=2 and x1=2,
the sequence merrily goes:
x2 = 3, x3 = 5, x4 = 10, x5 = 28, x6 = 154, x7 = 3520, x8 = 1551880,
x9 = 267593772160, x10 = 7160642690122633501504, ...
and so on till the last integer element at index 42, which has about 89,288,343,500 digits !!.
The next element, x43, is *not* integer and, using the asymptotic formula:
xn = ( n + 2 - 1/n + 4/n2 - 21/n3 + 138/n4 - 1091/n5 + ...) * C2n
where C = 1.04783144758..., approximately evaluates to:
x43 = 5.41 * 10 178,485,291,567
but a number more than 178 billion digits long isn't something that any current
computer can handle,
let alone handheld programmable calculators, not even HP's. Now, if this is the case for
the quite modest k=2, x1=2, and x43,
just try to imagine the humongous size of x601 (!!),
the first non-integer element
for k=3, x1=11. Even merely computing its
exact size is far from trivial, let alone the element itself !
Consider also the extremely disconcerting mathematical nature of what's happening here: a simple recurrence which unfailingly produces a sequence of integer elements for the first 600 cases suddenly and most unexpectedly results in a non-integer one for case 601, and the same goes for other values of the parameters.
What to do ? Well, I see three possible approaches. On one extreme you could try using
the full multiprecision route and proceed to compute each element till the first non-integer
one is found. This may succeed for a few values of k and x1
(k=2, x1=3 or k=8, x1=10,
for instance) but for the most part, as seen in the case of k=2, x1=2,
it's a totally hopeless dead end, now and in the foreseeable technological future.
The other extreme would be using as little multiprecision as possible or even none at all
and try to cope using modular arithmetic and solving congruences to go forward
seeking for the elusive non-integer element. This is mostly the route Eamonn and Egan
followed, but though they succeeded in the majority of cases at minimal computing cost,
they found some recalcitrant ones which, at the time of writing this, they haven't completely succeeded
in ironing out.
My original solution is an intermediate one. I don't use the full multiprecision approach,
as I quickly saw the task was impossible, but I'm using modular arithmetic to a sufficiently
high precision (up to 500 digits for the version I'm posting here) that the algorithm is
straightforward and no ambiguities
arise at any point, just a very simple loop till either the first
non-integer element is found, or the maximum precision is reached, or time becomes excessive, whichever comes first.
The only problem for me is that I wanted to implement my solution for the HP-71B which,
regrettably, lacks the multiprecision integer features of advanced RPL models. However, this
only means that besides the simple code to implement the main algorithm (6 lines),
I needed to
implement a multiprecision "library" which would allow calls from the main program to
perform the requested multiprecision operations on demand.
This isn't particularly easy, most specially as it must both run fast and be as small and compact
as possible since this is a simple S&S Math Challenge, not
a full-fledged article on how to implement multiprecision libraries in the HP-71B.
With that in mind, I created the following original solution for the HP-71B:
Main program (just 6 lines of HP-71B BASIC code implementing the base algorithm):
10 DESTROY ALL @ OPTION BASE 1 @ D=200 DIV 6+1 @ DIM F(D),X(D),Y(D),Z(D)
20 INPUT "K,X1=";K,X1 @ SETTIME 0 @ FOR N=10 TO 2000 STEP 10 @ CALL MFAC(N,F)
30 DISP N; @ CALL N2M(X1,Y) @ FOR I=2 TO N @ MAT Z=Y @ CALL MDIVK(Z,I-1)
40 CALL MPOW(Z,K,X) @ CALL MSUM(Y,X,Z) @ CALL MMOD(Z,F,Y) @ CALL MMODK(Y,I,M)
50 IF M THEN DISP @ DISP I;M/I;TIME @ END ALL
60 NEXT I @ NEXT N
As HP-71B's BASIC lacks any multiprecision operations, integer or otherwise, the main program
simply calls the relevant library functions (you'll find a reasonably detailed description
of the library in Apendix A below). The main program is executing the equivalent of the following
pseudo-code:
dimension multiprecision variables F, X, Y, and Z to hold a certain number of digits (200 above)
ask the user for K and X1
for values of N from 10 to 2000 (both included) and counting in steps of 10, do:
assign the factorial of N to F
assign the value of X1 to Y
for values of I from 2 to N, both included, do
assign the value of Y to Z
integer-divide Z by I-1
assign the value of Z raised to the kth power to X
assign the value of X + Y to Z
assign the value of Z Modulus F to Y
assign the value of Y Modulus I to M
If M is not zero, then we've found the first non-integer element, so
output index and fractional part and terminate
loop to process the next value of I
loop to process the next value of N
My original program is an implementation extensively relying in full multiprecision modular arithmetic, and
taking into account the ideas discussed in the book "Mainly Natural Numbers"
which you can freely download in PDF format from here. This
particular problem is featured in Chapter IV: "Some Sequences of Large Integers",
pp. 32-37., though the notation there differs slightly from mine; for instance my
case k=2,
x1=2, appears there as m=1, x1=2, etc.
You're advised to read it for the rigorous analysis and relevant details, as space
(and notational difficulty inherent to text-based HTML) does not permit to discuss the subject at lenght here.
Now this is a sample run for the case k=2, x1=41:
>RUN
k,x1=2,41
10 20 ... -> 17 .235294117647 1.53
where the running index is shown at steps of 10 while the computation is in progress, then the index
of the first non-integer element is output (17), as well as its fractional part
(0.235294117647)
and the time it took to find it (1.53 seconds under Emu71, a little over 6 min.
in a physical HP-71B).
Note:
The default precision listed above in line 10 is 200 digits, which is more than enough for all the cases
I asked except the two final, extra ones which require about 360 and 415 digits, respectively. If you get
an error (normally "Subscript") while checking some case, just increase the precision from 200 to some
suitable higher value less than 500 in line 10 above.
Running my program for the asked values, one gets the following results and times (for Emu71 @ 2.4 Ghz, single Pentium 4 CPU with
512 Mb RAM under Windows XP Home Edition 2002; for a physical HP-71B, multiply the given times by roughly 250):
k x1 index #digits fractional part time
-----------------------------------------------------
2 2 43 200 0.558139534884 19.10
2 3 7 200 0.285714285714 0.26
2 15 67 200 0.268656716418 77.60
2 41 17 200 0.235294117647 1.53
2 42 59 200 0.559322033898 44.40
2 50 34 200 0.117647058824 10.33
2 67 51 200 0.411764705882 41.56
2 71 17 200 0.235294117647 1.58
2 1958 17 200 0.235294117647 1.64
2 1992 17 0.235294117647 1.64
2 2008 61 200 0.983606557377 77.10
2 99 103 360 0.650485436893 593.41
3 2 89 415 0.460674157303 686.60
and further this is a short table of extra results for comparison purposes with the book and/or your
own programs (140 digits are more than enough so edit line 10 accordingly):
k x1 index #digits fractional part time
----------------------------------------------------
3 6 31 140 0.645161290323 16.40
4 7 13 140 0.692307692308 2.02
6 11 19 140 0.263157894737 5.19
8 3 7 140 0.285714285714 0.60
8 8 17 140 0.588235294118 6.01
Finally, this is an extensive table of results obtained with a different version of my
ad-hoc library,
which allows for much greater number of digits but takes at least one order of magnitude longer (the
posted version can only handle slightly less than 500 digits, being limited by the exponent range of the
built-in REAL precision, but it's smaller and much faster, and thus better suited for the
cases asked in this subchallenge).
For each assorted x1 and k, two numbers are given: the index of the
first non-integer element, and the approximate number of digits necessary to carry out the
computations,
which will appear as "?" if the program could not reach far enough before I gave up due to
timing considerations. Some of these results appear in the book I mentioned (and fully coincide with them),
but the rest I've never seen published before so this can be considered original research
publicly appearing for the first time in the World Wide Web and so independent corroboration would be most welcome.
x1 K=2 K=3 K=4 K=5 K=6 K=7 K=8 K=9 K=10
----------------------------------------------------------------------------------------------------------
2 43 103 89 395 97 598 214 2024 19 95 239 ? 37 331 79 1036 83 1222
3 7 6 89 403 17 54 43 255 83 733 191 2464 7 20 127 1901 31 323
4 17 27 89 403 23 84 139 1182 13 52 359 ? 23 167 158 2502 41 469
5 34 74 89 402 97 599 107 845 19 95 419 ? 37 331 79 1030 83 1224
6 17 27 31 97 149 1032 269 ? 13 47 127 1477 23 169 103 1449 71 1000
7 17 26 151 788 13 35 107 849 37 250 127 1475 37 332 103 1457 83 1226
8 51 128 79 343 13 35 214 2025 13 50 239 ? 17 107 163 ? 71 1000
9 17 25 89 402 83 491 139 1180 37 247 191 2464 23 167 103 1453 23 211
10 7 6 79 343 23 84 251 2460 347 ? 239 ? 7 19 163 ? 41 476
11 34 73 601 ? 13 34 107 848 19 88 461 ? 37 329 79 1035 31 323
12 17 27 197 1095 97 600 ? ? 37 248 191 2436 17 102 79 1034 41 478
13 17 26 151 788 23 85 ? ? 37 248 127 1480 37 327 158 2495 31 323
14 43 103 158 834 67 369 139 1184 37 245 191 2462 23 161 158 2503 41 473
15 67 186 197 1097 173 1244 139 1182 37 249 ? ? 37 332 127 1902 31 322
41 17 25 79 346 17 52 107 851 83 736 127 1479 23 168 79 1026 23 209
42 59 157 151 789 23 79 251 2462 37 244 191 2463 17 105 ? ? 31 324
50 34 74 193 1070 13 35 107 848 37 248 191 2463 23 169 79 1028 41 478
67 51 129 89 401 97 599 251 2463 37 250 191 2462 23 164 79 1031 83 1227
71 17 26 79 344 17 53 107 850 13 52 127 1477 46 446 ? ? 89 1342
99 103 324 31 97 13 34 139 1182 13 44 ? ? 37 332 151 2355 41 465
1958 17 27 158 834 13 34 251 2461 13 47 ? ? 149 2066 ? ? 41 474
1992 17 27 197 1098 17 53 107 850 109 1040 ? ? 23 168 103 1456 89 1342
2008 61 164 229 1321 167 1191 ? ? 13 52 191 2462 37 327 ? ? 31 322
That's all. This "Short & Sweet Math Challenges" series exclusive to
the MoHP ends its run now. The very
first S&SMC#1 was posted in May, 27th 2002 at 9:59 a.m.,
so it's been nearly six years sharing
what I hope has been interesting and rarely-seen math ideas and highly creative programming techniques for
an incredible variety of HP (and non-HP) calc (and non-calc) models.
Thanks
to all of you for your unfailing interest and for your extremely worthy contributions, it's been a pleasure
and a continued source of learning and awe to me.
Very Best Regards From V.
Appendix A: The multiprecision library
I'll briefly discuss here the small (just 32 lines of code) "Ad-hoc multiprecision library" I've
implemented specifically
for this particular subchallenge (the following listing is numbered as if to be included in
the same file as the main program without line numbers conflicting, but could reside in a separate
program file in RAM as well, if intended to be used by other programs).
Ad-hoc multiprecision library:
70 !
80 SUB MMUL(A(),B(),C()) @ MAT C=ZER @ CALL MTOP(A,U) @ CALL MTOP(B,V)
90 K=1000000 @ FOR I=1 TO V @ M=B(I) @ IF NOT M THEN 120 ELSE L=I
100 FOR J=1 TO U @ N=A(J) @ IF N THEN P=C(L)+M*N @ C(L)=RMD(P,K) @ C(L+1)=C(L+1)+P DIV K
110 L=L+1 @ NEXT J
120 NEXT I
130 SUB MPOW(A(),N,C()) @ IF NOT N THEN MAT C=ZER @ C(1)=1 @ END ELSE MAT C=A
140 IF N=1 THEN END ELSE DIM B(UBND(C,1)) @ U=LOG2(N) @ IF FP(U) THEN 160
150 FOR I=1 TO U @ MAT B=C @ CALL MMUL(B,B,C) @ NEXT I @ END
160 FOR I=2 TO N @ MAT B=C @ CALL MMUL(B,A,C) @ NEXT I
170 SUB MSUM(A(),B(),C()) @ MAT C=A+B @ CALL MNOR(C)
180 SUB MSBI(A(),B()) @ MAT A=A-B @ CALL MNOR(A)
190 SUB MMODK(A(),K,M) @ DIM B(UBND(A,1)) @ MAT B=A @ M=K @ CALL MDIVK(B,M)
200 SUB MDIVK(A(),K) @ FOR I=UBND(A,1) TO 2 STEP -1 @ P=A(I) @ IF P=0 THEN 220
210 A(I)=P DIV K @ A(I-1)=A(I-1)+MOD(P,K)*1000000
220 NEXT I @ P=A(1) @ A(1)=P DIV K @ K=MOD(P,K)
230 SUB MNOR(A()) @ FOR I=1 TO UBND(A,1) @ P=A(I)
240 IF P>999999 THEN A(I)=RMD(P,1000000) @ A(I+1)=A(I+1)+P DIV 1000000
250 IF P<0 THEN A(I)=P+1000000 @ A(I+1)=A(I+1)-1
260 NEXT I
270 SUB MTOP(A(),U) @ FOR U=UBND(A,1) TO 2 STEP -1 @ IF A(U) THEN END
280 NEXT U
290 SUB M2N(A(),N) @ N=0 @ FOR I=UBND(A,1) TO 1 STEP -1 @ N=N*1000000+A(I)
300 IF N>1E12 THEN N=SCALE10(N,6*I-6) @ END
310 NEXT I
320 SUB N2M(N,A()) @ MAT A=ZER @ I=0 @ M=N
330 IF M THEN I=I+1 @ A(I)=MOD(M,1000000) @ M=M DIV 1000000 @ GOTO 330
340 SUB MMOD(A(),B(),C()) @ U=UBND(A,1) @ DIM X(U),Z(U) @ CALL M2N(B,V) @ MAT C=A
350 CALL M2N(C,U) @ IF U<V THEN END
360 CALL N2M(MAX(INT(NEIGHBOR(NEIGHBOR(NEIGHBOR(NEIGHBOR(NEIGHBOR(U,0),0),0),0),0)/V),1),X)
370 CALL MMUL(X,B,Z) @ CALL MSBI(C,Z) @ GOTO 350
380 SUB MFAC(N,A()) @ CALL N2M(FACT(MIN(N,17)),A) @ IF N<18 THEN END
390 FOR I=18 TO N @ MAT A=(I)*A @ CALL MNOR(A) @ NEXT I
"Ad-hoc library" in this context means that it's tailored to do what the main
program needs to do
and absolutely nothing more, so it's mainly useful for the present purpose and not general-purpose by
any stretch of the imagination. Thus, it has the following disclaimer attached:
- the function set it's not complete or symmetric: only the indispensable
functionality for the
task at hand is included so many unneeded operations aren't implemented at all such as integer
division of two multiprecision arguments, for instance, while others are implemented only for a
multiprecision argument and a single precision one. Also, some of the operations return
the result in a variable different from the arguments while others are done in place. Note in
particular that no routines to input or output the multiprecision
values are provided.
- it's not optimized: all arguments are operated upon at maximum precision (say 200-digit
precision), even if the actual value stored in the multiprecision variable is only 10 digits long.
- it has no error checking at all: most errors will result
in a hard, untrapped error which will stop program execution. Also, debugging and checking hasn't been exhaustive so errors may have crept in.
- only positive integers are supported: any negative intermediate or final
results arising will result in a hard (untrapped) error.
Most of these shortcomings are the result of designing it to run fast, be small in size, and most importantly, take as
little of my scarce free time as possible to develop and debug. On the positive side, you'll find that:
- most operations are optimized for speed by
using assembly-language matrix operations, as well as one or
two neat tricks.
- as many as a full dozen routines are provided, all of them callable from your own programs (use
them at your own risk).
Available calls reference:
Note: In what follows,
- Multiprecision variables can hold multiprecision positive integer values up
to slightly less than 500 digits long (approximately), depending on the DIMension
of the underlying arrays.
Each of them should be DIMensioned (with OPTION BASE 1 in effect)
as a REAL precision
array with (N DIV 6)+1 elements where N is the maximum number
of digits, as each array element will hold up to 6 digits.
- Single-precision variables can hold positive integers up to 6 digits long.
- REAL variables can hold any positive integer value within the range
of built-in REAL precision.
SUB MMUL(A(),B(),C()) { Multiprecision MULtiplication }
SUB MPOW(A(),N,C()) { Multiprecision POWer }
SUB MSUM(A(),B(),C()) { Multiprecision SUM }
Takes two multiprecision variables A and B and returns their sum
in the multiprecision variable C. A and B are unaffected.
SUB MSBI(A(),B()) { Multiprecision SuBtraction In-place }
SUB MMODK(A(),K,M) { Multiprecision MODulus single-precision }
SUB MDIVK(A(),K) { Multiprecision DIV (and modulus) single-precision }
Takes a multiprecision variable A and a single-precision variable K and returns
the multiprecision value of the integer division operation A DIV K in the multiprecision variable A
(i.e., in place)
as well as the single-precision value of A modulus K in
the single-precision variable K
itself, which must be passed by reference. If you don't want or need the modulus
and want K to remain unaffected you can simply pass K by value instead.
SUB MNOR(A()) { Multiprecision NORmalization }
SUB MTOP(A(),U) { Multiprecision TOP }
Internal use only, though it can be useful for printing multiprecision variables, etc.
Returns in the single-precision variable U the index of the topmost
non-zero element of the multiprecision variable A. Variable U must be
passed by reference. Variable A is unaffected.
SUB M2N(A(),N) { Multiprecision to real-precision Numeric variable }
SUB N2M(N,A()) { Numeric Real-precision variable to Multiprecision }
SUB MMOD(A(),B(),C()) { Multiprecision MODulus }
SUB MFAC(N,A()) { Multiprecision FACtorial }
Sample use:
69! MOD 1368 = 161707119156747337147933149666645422128978599226831537632782504385470771962
This is how we would proceed (from the keyboard in this case):
1) DIMension several multiprecision variables capable of holding up to 100 digits (actually 102):
> DESTROY ALL @ OPTION BASE 1 @ P=100 DIV 6+1 @ DIM F(P),A(P),B(P)
2) Make the necessary calls to the multiprecision library in order to perform the computation:
> CALL MFAC(69,F) { compute the factorial of 69 }
@ CALL N2M(13,A) { assign the value 13 to multiprecision variable A }
@ CALL MPOW(A,68,B) { raise A (13) to the 68th power and assign the result to B }
@ CALL MMOD(F,B,A) { compute the modulus operation and assign the result to A }
3) Output the multiprecision result (watch your WIDTH setting):
> FOR I=P TO 1 STEP -1 @ DISP A(I); @ NEXT I @ DISP
0 0 0 0 161 707119 156747 337147 933149 666645 422128 978599 226831 537632 782504 385470 771962
which is the correct result as shown above, disregarding the non-significant 0's at the beginning.
|