S&SMC#18: My Original Solutions & Comments [LONG] Message #50 Posted by Valentin Albillo on 9 Apr 2007, 7:07 a.m., in response to message #1 by Valentin Albillo
Hi all,
Thanks to all lurkers and posters alike for your interest in this S&SMC#18,
there's been a healthy number of keen solutions and interesting comments to some
of the most difficult subchallenges and I'll now give my original solutions
plus assorted relevant comments:
My Original Solutions:
1) (HP-71B specific) Find a root X of:
Abs(Ln(x*x)-Ln(x2)) > 4.012007
As Jean-François Garnier early pointed out, this subchallenge seems at first
sight to have no solution at all, since you seem to be subtracting two mathematically equivalent quantities and requesting that the difference is to be a relatively
large value which can't possibly be brought about by minor rounding
errors.
The solution, which is HP-71B specific, lies in the fact that the
developing team for the internal ROMs and the Math ROM didn't have the
necessary allowances for ROM space and time to be able to fully implement
the Math IEEE proposals to consistently handle signed zero, most specially
when dealing with complex numbers and branch cuts. So there are times when
+0 and -0 aren't treated equivalently in a number of complex-valued functions
where they should be. This "x*x vs x^2" affair is but one such case.
Jean-François gave the simplest possible example for the case x*x versus x^2, where
an explicit zero component of a complex value had its sign explicitly changed.
This can be further masked out of user's sight by using instead some
function or operation that returns a -0 result from some nonzero argument,
such as, in degrees, Cos(270).
So, a complex number such as X = COS(270) + i will be a particular
solution for this challenge, and, in general:
X = COS( 270+360*k) + m*i,
X = COS(-270-360*k) + m*i,
with k integer >=0 and m#0 in both cases, will be a solution too. The
following short program demonstrates this by evaluating the expression LN(X*X)-LN(X^2)
for complex arguments formed using all integer values of the angle between 0 and 360 degrees:
10 DEF FNF(X)=LN(X*X)-LN(X^2)
20 OPTION ANGLE DEGREES @ STD @ COMPLEX X @ REAL Y,I
30 FOR I=0 TO 360 @ X=(COS(I),1) @ Y=ABS(FNF(X))>1
40 IF Y THEN DISP "A solution is: x = (COS(";STR$(I);" deg),1)"
50 NEXT I
>RUN
A solution is: x = (COS(270 deg),1)
and, as you can see, only COS(270) results in a difference greater
than 1 (or than 4.012007 for that matter, namely 2*Pi). So it's not
just a question of starting with a weird initial value, -0, and
not-so-unexpectedly getting weird results, the problem is that
you can start with such a nice integer value as 270 and, totally unaware
of this pitfall, get weird, outrageously out-of-range results where
all other integer (or even real) values in 0-360 would give you no
unexpected complications at all.
Even the neighboring values 269.999999999 and 270.000000001
will return the expected results, but there's a huge, abnormal discontinuity exactly
at 270 degrees in this case.
2) (any model) Write a program that takes as input three positive integers A,B,M
(A,B not being exact powers of 10), and outputs the difference
between the probabilities P(A,M) and P(B,M) that two random powers of A
and B (say Ai and Bj
for some random positive integers i,j) begin with M.
This subchallenge was also cleverly solved by Egan Ford, and indeed the solution is a program
which simply outputs 0, regardless of the inputs.
This is based on the proven fact that the probability that a power of any positive
integer A, whose logarithm to the base 10 is irrational (i.e.: any positive integer
other than a power of 10), begins with the digits that represent the number M
(in decimal notation) is
P(A,M) = LGT(1+1/M)
i.e.: independent of A. This is actually quite easy to prove, both formally and
informally, and due to its independence of A, the difference of both probabilities
is always 0, regardless of the values of A and B within the original constraints.
3) (any model) Let F(X) be a function defined thus for positive X:
F(1) = 1
F(X) = F(X-1) + 1/X
Write a program that, given any value N > 1 and N < 200 (for 10-digit
models) or N < 1000 (for 12-digit models), will find
with the maximum possible accuracy a value X such that F(X) = N. Use your program to find X such that F(X) = 4.012007
This definition is one of the many that can be formulated for the
Harmonic Number function,
H(N) = 1 + 1/2 +1/3 + ... + 1/N
The challenge resides in how to extend the function's domain
to arbitrary real values, as the given definition is only useful for
discrete positive integer values, which is perfectly Ok as a
theoretical definition is not necessarily an effective computation scheme.
This is quite similar to the well-known factorial function, N!,
usually defined as
F(0) = 1
F(X) = F(X-1) * X
which can only be used to compute the function for discrete positive
integer values, and its extension to arbitrary real (and even complex) values X, usually
known as Gamma function, requires a new definition which totally agrees
with the discrete case and further meets a number of continuity requirements, etc.
In our case, H(N), a simple but appropriate extension would be this:
/ 1 N
| 1 - (1 - X)
H(N) = | ------------ . dX
| X
/ 0
which coincides with the original definition for all positive integer
values of N, plus it also works for suitable positive real arguments.
An HP-71B program to find N for a given value of F(N) would then
simply be:
10 DESTROY X,N @ INPUT "F(N)=";N @ X=FNROOT(1,100,LN(FVAR)+.577-N)
20 DISP FNROOT(X-1,X+1,INTEGRAL(0,1,1E-12,(1-(1-IVAR)^FVAR)/IVAR)-N)
>RUN
F(N)=4.012007 [ENTER]
30.523595225 (25.37 sec. in Emu71)
which can be shortened to a single line, or directly evaluated from
the command line with no program necessary, but as written we first
compute a suitable first approximation by using a well-known asymptotic
approach to H(N) (which involves the EulerGamma constant 0.577...), and this initial approximation is then used to bracket
the two initial guesses for the FNROOT (Solve) built-in function when
applied to the integral extension in order to appreciably reduce
computation time. The single-line version works the same, but takes
longer unless you supply good initial guesses.
4) (any model) Write a program, command line or key sequence
to compute the exact symbolic value of (in radians mode):
/ Inf
|
I1 = | 1/((1+x2)*(1+x4.012007)) .dx
|
/ 0
/ Pi/2
|
I2 = | 1/(1+Tan(x)4.012007) .dx
|
/ 0
As some of you quickly found out, both integrals are one and the same upon a simple
change of variables. Furthermore, the value of
/ Pi/2
|
I = | 1/(1+Tan(x)S) .dx
|
/ 0
is independent of the particular value of S, i.e. it is the same whether
S = 4.012007 or S = 2007, or S = 0 for that matter, so we can take S = 0
and then the integral's value always will be:
/ Pi/2
|
I2 = | 1/(1+Tan(x)0) .dx
|
/ 0
/ Pi/2
|
= | 1/2 .dx = Pi/4.
|
/ 0
As for a program to compute this, it suffices to supply the integrand 1/2
to any integration program or built-in functionality, and either recognize
the computed value by comparing it with submultiples of Pi or else call QPI
or any other such program to do the recognizing, that is, assuming your HP
model doesn't feature some CAS because if it does, any CAS can symbolically
evaluate the definite integral of 1/2 and return Pi/4 at once. The details
are absolutely trivial.
The fact that the original general integral is independent of the value of
S is but a particular case of a much more general, extremely powerful, multiparameter
master formula whose value also happens to be independent
of one of its parameters. For instance, a particularized case for three
specific values of its parameters would get us this unusual exact symbolic computation:
/ Inf
| dx
I = | ----------------------------------------------- = Pi/(2*(1+Sqr(2))
| (x4 + (1+2*Sqr(2))*x2 + 1)*(x100 - x98 + ... + 1)
/ 0
5) (any 12-digit model) Given the succession defined thus:
X0 = 0
30*M2 - 89*M + 26
XN = FractionalPart( 24*XN-1 + ------------------------------------------ )
2^3*M4 - 26*M3 + 178*M2 - 206*M + 84
where M = 22*N
write a program to compute the first 9 terms (X1, X2, ..., etc.)
10 X=0 @ Z=0 @ FOR N=1 TO 9 @ M=2^2*N
20 X=FP(2^4*X+(30*M^2-89*M+2^6)/(2^3*M^4-2^6*M^3+178*M^2-206*M+84))
30 DISP INT(2^4*X); @ NEXT N @ DISP
>RUN
2 4 3 15 6 10 8 8 8
which happen to be the first 9 hexadecimal digits of the fractional part of Pi.
This HP-71B program, which of course uses 12-digit arithmetic, can't do
better and the very next terms are all wrong, starting with the 10th
which comes out as 2 instead of the correct 5. This is a problem of insufficient
computational precision, not of the formula itself. So the correct answer
to the first question is:
"The Nth term is the Nth hexadecimal digit of the fractional part of Pi"
Surprisingly, it is not known if this is true for the whole
infinite sequence, but it's been
proved that there can only exist a finite number of exceptions, and the
formula has been numerically checked correct for the first 1,000,000 terms.
The theoretical importance of this sequence lies in the following fact: if
it could be proved that the generated terms are uniformly distributed
in the range 0-16, then Pi would be proved normal to base 16 (and also to base 2).
Currently, Pi hasn't been proved to be normal to any base, despite the fact that
more than 1012 digits of Pi have been computed and they pass all statistical tests
for normality to all bases, and despite the proven fact that almost all real numbers
are normal to all bases.
6) (HP-15C specific) Assuming A and B are both NxN matrices and that their elements are integer
numbers less than 5,000,000,000, write a program that will exchange their
contents, for NxN up to and including 5x5.
The problem here is that, not being able to use any conditional instruction, it's quite
impossible to exchange the elements using a loop, and further there's not enough RAM to define an
auxiliary 5x5 matrix to help make the exchange. The correct solution is thus the
following 11-step program for the HP-15C:
LBL A
RCL MATRIX A
RCL MATRIX B
RESULT A
+
RCL MATRIX A
LAST X
RESULT B
-
RESULT A
+
which succeeds in exchanging the contents of A and B by arithmetic means, without
loss of precision because due to the original constraints on the elements, no intermediate
result ever exceeds the allowed range for exact integers. Paul Dale correctly saw the right arithmetic approach, with his quickly concocted 13-step solution.
7) A very interesting and commonsense-defying subchallenge which no one took (perhaps Calculus-themed challenges aren't that popular here), thus best left for a future Datafile article.
8) (any model) Write a program to find, compute, and print all square factorials N!
from N=1 to N=512.
Despite this being an "April's Fool" subchallenge, which frequently need some lateral thinking nevertheless, and despite my hints suggesting the need to think plainly about what a
"square" is (the kind of way Sesame Street's Grover would explain the concept to toddlers,
which probably wouldn't grasp irrationals, complex arguments, modular residues, and
existence proofs), noone stumbled upon the correct interpretation initially.
The alleged "squareness" of the sought-for factorials made actually no reference to
their *value* but to their *shape* when printed. For instance, 12! can be printed thus:
4 7 9
0 0 1
6 0 0
which, sure enough, it's pretty square on paper or screen. This is possible because
12! = 479001600 is a 9-digit integer, and thus, having a square number of digits, it
can be printed as a square array of digits, which is all perfect and good and quite
pleasant from an aesthetic point of view.
The concept itself is hardly new and I saw it
for the first time several decades ago in one of Martin Gardner's extremely popular Scientific
American columns, where he discussed not only "square" factorials, but "polygonal"
factorials in general, including a very nice-looking, pretty-printed "octogonal" factorial !
All we need then is some program which can detect which factorials from 1! to 512! do
have an square number of digits, compute them to full precision, and print them
as an square array of digits.
This short program (465 bytes) is my complete solution
for the HP-71B, which includes a multiprecision factorial-computing subprogram,
BIGFACT (which accepts two parameters, one by value which is the integer
argument, the other by reference, which is a string variable where the resulting
multiprecision factorial representation will be returned), and a simple "driver" program,
SQRFACT, which finds out which factorials do
have an square number of digits by using a simple Stirling approximation, then
calls BIGFACT to fully compute the ones found in order to subsequently print them (20 in all):
100 DESTROY ALL @ FOR N=1 TO 512
110 D=MAX(SQR(INT(N*LGT(N/EXP(1))+LGT(2*PI*N)/2+LGT(1+1/12/N))+1),1)
120 IF FP(D) THEN 140 ELSE DIM S$[D*D] @ CALL BIGFACT(N,S$) @ DISP N;LEN(S$)
130 FOR I=1 TO D @ DISP S$[1+(I-1)*D,D+(I-1)*D] @ NEXT I
140 NEXT N
150 !
160 SUB BIGFACT(N,S$) @ OPTION BASE 0 @ K=9 @ L=10^K @ STD
170 M=INT((N*LGT(N/EXP(1))+LGT(2*PI*N)/2)/K)+1 @ DIM F(M)
180 MAT F=ZER @ F(0)=1 @ FOR I=2 TO N @ MAT F=(I)*F @ FOR J=0 TO M-1
190 T=F(J) @ IF T>=L THEN F(J+1)=F(J+1)+T DIV L @ F(J)=MOD(T,L)
200 NEXT J @ NEXT I @ P=1 @ S$="" @ B$="" @ FOR I=M TO 0 STEP -1
210 IF F(I) AND P THEN P=0 ELSE IF P THEN 230 ELSE B$="0"
220 A$=STR$(F(I)) @ S$=S$&RPT$(B$,K-LEN(A$))&A$
230 NEXT I
>RUN
1 1
1
2 1
2
3 1
6
7 4
5 0
4 0
12 9
4 7 9
0 0 1
6 0 0
18 16
6 4 0 2
3 7 3 7
0 5 7 2
8 0 0 0
32 36
2 6 3 1 3 0
8 3 6 9 3 3
6 9 3 5 3 0
1 6 7 2 1 8
0 1 2 1 6 0
0 0 0 0 0 0
59 81
1 3 8 6 8 3 1 1 8
5 4 5 6 8 9 8 3 5
7 3 7 9 3 9 0 1 9
7 2 0 3 8 9 4 0 6
3 4 5 9 0 2 8 7 6
7 7 2 6 8 7 4 3 2
5 4 0 8 2 1 2 9 4
9 4 0 1 6 0 0 0 0
0 0 0 0 0 0 0 0 0
...
508 1156
5 1 1 9 9 ... 8 5 3
7 9 7 9 8 ... 6 0 1
1 1 7 4 6 ... 6 5 4
... ... ... ...
4 7 5 2 9 ... 8 4 9
8 0 9 0 4 ... 0 0 0
0 0 0 0 0 ... 0 0 0
0 0 0 0 0 ... 0 0 0
0 0 0 0 0 ... 0 0 0
That's all. Hope you enjoyed the ride, and thank you very much for your interest and outstandingly clever solutions, extensions, and comments.
Best regards from V.
|