S&SMC#19: My Original Solutions & Comments [LONG] Message #72 Posted by Valentin Albillo on 28 June 2007, 12:38 p.m., in response to message #1 by Valentin Albillo
Hi all,
Thanks to all of you who were interested in this S&SMC#19, most specially to the ones
who worked really hard at the individual Grades and posted solutions. I'll give now my original
solutions for the HP71B (very easy to convert to other similarly capable models), plus additional comments.
Grade F:
Write a program to compute the N^{th} term for the following sequences:
Sequence 1: u_{1} = 0, u_{2} = 0, u_{n+2} = u_{n+1} + u_{n}/n
Sequence 2: u_{1} = 0, u_{2} = 0, u_{n+2} = u_{n+1}/n + u_{n}
where the number of terms N (which can be distinct por each sequence) is either
hardcoded or given, your choice. After computing the N^{th} term, u_{N},
your program must then compute and output the following values:
Sequence 1: L = N*1 / u_{N}^{1}
Sequence 2: L = N*2 / u_{N}^{2}
Now use your program to help you try and predict, for each sequence, the corresponding
theoretical exact limit of L when the number of terms N tends to infinity.
My original solution for the HP71B is the following trivial 3line (135 byte) program
which computes and outputs the corresponding values of L for N = 13 and N = 13000, respectively.
1 DESTROY ALL @ STD @ K1=13 @ K2=13000
2 U=0 @ V=1 @ FOR N=1 TO K1 @ W=U/N+V @ U=V @ V=W @ NEXT N @ DISP N*1/U^1
3 U=0 @ V=1 @ FOR N=1 TO K2 @ W=V/N+U @ U=V @ V=W @ NEXT N @ DISP N*2/U^2
>RUN
2.71828182845
3.14195515673
As you may see, the sequences (very quickly) converge to
e = 2.71828182846 and (very slowly) to Pi = 3.14159265359,
respectively.
The Surprise ! factor comes from the fact that both
sequences are very, very similar (as can be seen both from their definitions and
from the program listing above, where lines 2 and 3 are almost a slightly
particularized copypaste version of one another), yet one converges to the
uberimportant constant Pi, the other to the
nolessimportant constant e, suggesting a certain
symmetry, a relationship, a "common origin" for both constants.
Some of you went to also compute 13,000 terms for the first sequence. This does no
good, as 13 terms are all that is needed to achieve full 12digit precision. After
that, rounding errors do accumulate and by the time the 13,000th term is reached,
most of the precision is gone.
Grade D:
Write a program which, given N, will autoscale,
label, and plot the ISin(x, N) function, for arguments X in radians belonging to the
fixed interval [0, 2*Pi] at increments of Pi/40, say, where
ISin(x, N) = Sin(Sin(Sin( ... Sin(x) ... )))
i.e., the iterated Sine of X, where N, a positive integer, is the number of iterations
Looking at the plot itself and the value of Maxim, try and predict what happens
for increasing values of N and in the limit when N tends to infinite.
Given this challenge's nature, I was absolutely expecting that HP graphics models would
be put to good use here, much more so than nongraphics ones, so I was utterly
surprised that not a single RPL solution was posted at all. Go figure ! ...
My original solution for the (nongraphics) HP71B is the following simple 8line (268byte) program:
10 DESTROY ALL @ OPTION ANGLE RADIANS @ D=3 @ FIX D @ INPUT "N=";N
20 K=FNS(PI/2,N) @ DISP "Maxim =";K @ FOR X=0 TO 2*PI STEP PI/20
30 Y=FNR(FNS(X,N)/K,D) @ U=45 @ V=FNR(Y*25+45,0) @ DISP X;Y;
40 IF V>U THEN DISP TAB(U);":";TAB(V);"*" ELSE DISP TAB(V);"*";
50 IF V<U THEN DISP TAB(U);":" ELSE IF V=U THEN DISP
60 NEXT X
70 DEF FNR(X,N)=INT(X*10^N+.5)/10^N
80 DEF FNS(X,N) @ FOR N=N TO 1 STEP 1 @ X=SIN(X) @ NEXT N @ FNS=X
where FNS(X, N) is a userdefined function that can be used from the program itself or
directly from the command line to compute the Iterated Sine of X with N iterations, for instance
>FNS(1,1)
.841470984808 ( = SIN(1) )
>FNS(PI/3,5)
.594535945122 ( = SIN(SIN(SIN(SIN(SIN(PI/3))))) )
and FNR(X, N) is a convenience function that returns X rounded to N digits
(where N can be negative to achieve rounding to the left of the decimal
point), for instance:
>FNR(PI,3)
3.142
>FNR(100*PI,1)
310
Running the program for increasing values of N we get plots like the following for N = 10:
Maxim = 0.481
0.000 0.000 *
0.157 0.314 : *
0.314 0.565 : *
0.471 0.737 : *
0.628 0.845 : *
0.785 0.911 : *
0.942 0.951 : *
1.100 0.976 : *
1.257 0.990 : *
1.414 0.998 : *
1.571 1.000 : *
1.728 0.998 : *
1.885 0.990 : *
2.042 0.976 : *
2.199 0.951 : *
2.356 0.911 : *
2.513 0.845 : *
2.670 0.737 : *
2.827 0.565 : *
2.985 0.314 : *
3.142 0.000 *
3.299 0.314 * :
3.456 0.565 * :
3.613 0.737 * :
3.770 0.845 * :
3.927 0.911 * :
4.084 0.951 * :
4.241 0.976 * :
4.398 0.990 * :
4.555 0.998 * :
4.712 1.000 * :
4.869 0.998 * :
5.027 0.990 * :
5.184 0.976 * :
5.341 0.951 * :
5.498 0.911 * :
5.655 0.845 * :
5.812 0.737 * :
5.969 0.565 * :
6.126 0.314 * :
and for N = 1000:
Maxim = 0.055
0.000 0.000 *
0.157 0.946 : *
0.314 0.987 : *
0.471 0.995 : *
0.628 0.998 : *
0.785 0.999 : *
0.942 0.999 : *
1.100 1.000 : *
1.257 1.000 : *
1.414 1.000 : *
1.571 1.000 : *
1.728 1.000 : *
1.885 1.000 : *
2.042 1.000 : *
2.199 0.999 : *
2.356 0.999 : *
2.513 0.998 : *
2.670 0.995 : *
2.827 0.987 : *
2.985 0.946 : *
3.142 0.000 *
3.299 0.946 * :
3.456 0.987 * :
3.613 0.995 * :
3.770 0.998 * :
3.927 0.999 * :
4.084 0.999 * :
4.241 1.000 * :
4.398 1.000 * :
4.555 1.000 * :
4.712 1.000 * :
4.869 1.000 * :
5.027 1.000 * :
5.184 1.000 * :
5.341 0.999 * :
5.498 0.999 * :
5.655 0.998 * :
5.812 0.995 * :
5.969 0.987 * :
6.126 0.946 * :
and we see that the limiting waveform of nested SIN is identically 0, but if we keep
scaling the function up,
we get in the limit a square wave. If you've got a copy of
Mathematica
you can see it very clearly by looking at the "animation" created by this oneliner:
Do[ Plot[ Nest[Sin, x, n], {x, 2Pi, 2Pi}, PlotRange > {1, 1}, PlotLabel > n ], {n, 1, 50} ]
where it's clear that what's happening is that the whole plot slowly collapses towards zero but
different parts collapse at different rates, making the plot get more and more square. This is not
difficult to prove by using a wellknown theorem from real analysis that states that every bounded
monotonic sequence approaches a finite limit.
The convergence to 0 is dreadfully slow, because for small arguments, we have with
high accuracy:
Sin(x) => x  x^{3}/6
which means that when we replace x by Sin(x), the only thing that
gets us any closer to 0 is the term x^{3}/6, which, if x
itself is small, it is a much smaller decrement indeed.
Grade C:
Find and output all values of N < 10000 such that both N is prime and
P(N) = S(N)
where P(N) gives the number of primes <= N and S(N) gives the sum
of the factorials of the digits of N.
Your program must be optimized primarily for speed, then for program size.
The following is my original solution for the barebones HP71B, a plainvanilla 6line (220 byte) program
which doesn't require any ROMs or LEX addons and trades RAM usage for speed to deliver
the goods rather quickly:
10 DESTROY ALL @ OPTION BASE 0 @ M=5000 @ INTEGER P(M) @ P(1)=1
20 FOR N=2 TO (SQR(2*M)+1)/2 @ IF P(N) THEN 40
30 FOR K=2*N*(N1)+1 TO M STEP 2*N1 @ P(K)=1 @ NEXT K
40 NEXT N @ T=1 @ FOR N=2 TO M @ IF P(N) THEN 60 ELSE K=2*N1 @ S=0 @ T=T+1
50 S=S+FACT(MOD(K,10)) @ K=K DIV 10 @ IF K THEN 50 ELSE IF S=T THEN DISP 2*N1
60 NEXT N @ DISP "OK"
>RUN
6521
OK
Lines 1040 perform a fast sieving procedure to flag all odd numbers < 10000 as
either prime or composite, recording the results in P, an integer array of adequate
dimension.
This allows lines 4060 to quickly scan the array, skipping all
nonprime entries, and computing the Pi(N) function effortlessly by simply
updating the prime count each time an element is flagged as a prime. This
running value of Pi(N) is compared to the sum of the factorials of the
digits of the current N (which is computed in a tight loop at line 50) and any
solution found is then immediately output, resuming the search afterwards.
As it happens, 6521 is the only solution < 10,000, but there's just another one
(and no more) in the range from 10,000 to infinity which some of you correctly
and cleverly located, namely 5224903.
If composite numbers were also allowed as solutions, there would be 23 solutions in
all, starting at 6500 and ending at 11071599. Among these,
I find it interesting that 5854409, 5854419, and 5854429
happen to be both consecutive and in arithmetic progression.
This problem can be generalized to numerical bases other than 10, for example
6719_{10} is a solution in base 7 and
5378437_{10} is a solution in base 15.
Grade B:
Regrettably, this very interesting challenge, which has a truly surprising result,
didn't seem to catch the eye of the vast majority of you for whatever reason and
was left largely untouched.
PeterP made a brave attempt at dealing with it, but his final version still has
a number of fatal flaws (such as the fact that it doesn't properly separate the
roots from the poles and mixes them at times while computing the sum, etc), and
so the numerical result isn't accurate. As is my policy in these cases, when a challenge
gets essentially no acceptable responses I simply forfeit it as well, lest its appeal would be wasted.
The "surprise" factor of this challenge is the extremely unexpected result you get.
Having an infinite sum of the roots of a transcendental trigonometric equation with
arguments in radians, one would expect to get some weird noname irrational, possibly transcendental number related to Pi in some way, say PI^{3}/310 or some such value.
But nothing could be further from the truth because what you actually get absolutely defies intuition, and paradoxically enough, it is fairly easy to concoct a purely formal proof "a la Euler"
to prove that this unexpected result is mathematically exact and correct. The value
you get for the B+ variant isn't that big a surprise, but is also most
unexpected in its very unique way.
By the way, one of the main difficulties of this
challenge lies in the fact that you must compute a lot of roots (possibly thousands), very quickly,
and very accurately, in the sense that all roots should be accounted for with no
missing entries, in their proper order, and that no poles are mistaken for roots.
This requires a very precise separation between poles and roots and, as it happens,
they tend to become extremely close for larger values, which enormously increases
the difficulty of separating them. Besides, as you know, HP builtin solvers have
a very hard time telling apart a pole from a root !
Grade A:
The limit of the sum of the reciprocals of the natural numbers 1 to N does diverge when N tends to infinite:
N>oo

\ 1
S = >  > oo
/ N

N = 1
However, if we simply leave out of the sum just those terms whose denominators
include the digit 9 (such as 199, 34343900132, ...), the infinite sum then
nicely converges to a finite value.
You must write a program to compute as accurately as possible
this finite value of the infinite sum. You must optimize your program for both
accuracy and speed.
This has been the grade of choice for many of you, with excellent results. This series
converges far too slowly for direct summation, as has been extensively commented in the thread,
so better algorithms are definitely required.
As some of you pointed out, there's a PDF document titled "SUMMING A CURIOUS, SLOWLY CONVERGENT
SERIES" by Schmelzer and Baillie, freely downloadable from various web locations, which discusses the problem at length and gives in sufficient detail a couple
of algorithms: the simpler Fischer algorithm, which can give the sum for the original Grade A
problem to high accuracy in very fast times (but which can't be used for the A+ variant as
it does not generalize to other digits or sequences of digits), and a more convoluted but
extremely general algorithm by the authors, which can be used for the generalized sum. See
the paper for the details.
I'll give now my two original A solutions, both implementing the Fischer algorithm but
one of them optimized mainly for size, the other for speed.
This is my shorter original solution, a 6line (265byte) HP71B program:
1 DESTROY ALL @ OPTION BASE 0 @ M=11 @ D=10 @ DIM B(M) @ B(0)=D @ FOR N=2 TO M
2 S=0 @ FOR K=N TO 2 STEP 1 @ S=S+FNC(N,K)*(D^(NK+1)D^K+1)*B(NK) @ NEXT K
3 B(N1)=(D*(11^ND^N)S)/N/(D^N9) @ NEXT N @ H=D*LN(D)
4 FOR K=2 TO M1 @ H=HB(K1)*FNZ(K)/D^K @ NEXT K @ DISP H
5 DEF FNC(N,K)=FACT(N)/FACT(K)/FACT(NK)
6 DEF FNZ(X)=INTEGRAL(0,100,0,IVAR^(X1)/(EXP(IVAR)1))/GAMMA(X)
>RUN
22.9206766192 ( 6.83 seconds, Emu71 @ 2.4 Ghz, about 20 min in a real HP71B )
and this is my longer (but 60+ times faster !) original solution, a 10line (443byte)
solution for the HP71B:
1 DESTROY ALL @ OPTION BASE 0 @ M=11 @ D=10 @ DIM B(M) @ B(0)=D @ FOR N=2 TO M
2 S=0 @ FOR K=N TO 2 STEP 1 @ S=S+FNC(N,K)*(D^(NK+1)D^K+1)*B(NK) @ NEXT K
3 B(N1)=(D*(11^ND^N)S)/N/(D^N9) @ NEXT N @ H=D*LN(D)
4 FOR K=2 TO M1 @ H=HB(K1)*FNZ(K)/D^K @ NEXT K @ DISP H
5 DEF FNC(N,K)=FACT(N)/FACT(K)/FACT(NK)
6 DEF FNZ(Z) @ IF Z=2 THEN FNZ=PI*PI/6 ELSE IF Z=3 THEN FNZ=1.20205690316
7 IF Z=4 THEN FNZ=PI^4/90 ELSE IF Z=5 THEN FNZ=1.03692775514
8 IF Z=6 THEN FNZ=PI^6/945 ELSE IF Z=7 THEN FNZ=1.00834927738
9 IF Z<8 THEN END ELSE S=1 @ T=0 @ N=2
10 S=S+N^(Z) @ N=N+1 @ IF S<>T THEN T=S @ GOTO 10 ELSE FNZ=S
>RUN
22.9206766192 ( 0.1 seconds, Emu71 @ 2.4 Ghz, about 20 seconds in a real HP71B)
Both results are identical and correct to 12 digits, but the longer version runs 60+ times
faster. The only difference between them lies in the definition of the Riemann Zeta function:
 The shorter version implements it as a userdefined function, FNZ(X), which can
compute the Zeta function for real arguments, not only integers, using integration.
As such, it does require the Math ROM, and it's somewhat slow. For instance:
>FNZ(10), FNZ(PI)
1.00099457513 1.17624173838
 On the other hand, the longer version implements the same userdefined function but
only for integer arguments, using theoretical constant values for arguments from 2 to 7
and the power series for integer arguments greater than 7, which converges quickly
enough for such arguments. This results in much faster computation times when compared
to the integral version and further does not require the Math ROM, a bare bones HP71B
will do nicely.
As for the A+ variant, Egan Ford has worked very hard and very successfully to
implement the quite complicated SchmelzerBaillie's general algorithm, with excellent
results, my most sincere congratulations for his outstanding contribution and
thanks a lot for his interest in this challenge.
That's all. I hope you've enjoyed this S&SMC (and even learned a math thingy and/or programming trick or two), I certainly did enjoy seeing your
truly excellent posts, as clever and insightful as always, you're the best !
And last but not least, a fond welcome to PeterP, which finally got to try and solve some of the
grades like a pro. See, Peter ? It wasn't that difficult, man ! :)
Best regards from V.
