Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi all,
After nearly half a year has elapsed since my latest S&SMC, here's
a new one for you to try and conquer, namely my brandnew S&SMC#18:
April 1st, 2007 Spring's Special, similar in spirit to the one I posted
in 2006 but this time for April 1st, 2007, i.e.:
4.012007 in plain numeric "date format".
To allow for variety, this S&SMC is further subdivided into
8 subchallenges, all of them very short and sweet, and each one being
assigned 10 points so that if you manage to solve them all, you'll get
a perfect 100% score of 80 points. Some of them are intended for specific
models, others are general in nature so any HP calc model can be
used in principle to solve them though one or two might require suitably
advanced models due to RAM requirements if nothing else.
You're expected to get your favourite HP calc, your programming
ingenuity, experience and wits, and put everything to the task of solving them.
Notes:
 The modelspecific questions apply only to that model in particular. They
would probably be meaningless for any other models.
 If you use an HP41C, you can also use any of these:
Math ROM, Advantage ROM, Card Reader ROM, Printer ROM. No other ROMs or extra
routines allowed.
 If you use an HP71B, you can also use any of these: Math ROM, HPIL ROM,
STRNGLEX. No other ROMs or LEX files allowed.
The Challenge:
1) (HP71B specific) Find a root x of:
Abs(Ln(x*x)Ln(x^{2})) > 4.012007
where, of course, Abs is the "Absolute Value" function and
Ln is the "Natural Logarithm" function.
Score: Finding one root: 8 points
Finding additional roots: 2 points:
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 A^{ i} and B^{ j}
for some random positive integers i, j) begin with M.
For instance, if A=265, B=1387, M=123456 > RUN > P(A,M)  P(B,M)
where P(A,M) = probability that a random power of 265 begins with 123456
P(B,M) = probability that a random power of 1387 begins with 123456
Use your program to compute this difference for the case: A=314159, B=271828, M= 4012007
Score: Writing a working program: 5 points:
Getting the correct answer for the specific case: 5 points
3) (Any model): Let F(X) be a function defined thus for real positive X:
F(1) = 1
F(X) = F(X1)+1/X
for instance: F(2) = F(1)+1/2 = 1+1/2
F(3) = F(2)+1/3 = 1+1/2+1/3
F(4) = F(3)+1/4 = 1+1/2+1/3+1/4
...
F(100) = 1+1/2+1/3+1/4+1/5+...+1/98+1/99+1/100
... etc ...
 Write a program that, given any real value N > 1 and N < 200 (for 10digit
models) or N < 1000 (for 12digit models), it will find
the value of X (with full 1012 digit accuracy) such that F(X) = N.
 Use your program to find X such that F(X) = 4.012007
Score: Writing a working program: 5 points
Getting the correct answer for the specific case: 5 points
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+x^{2})*(1+x^{ 4.012007})) .dx

/ 0
/ Pi/2

I2 =  1/(1+Tan(x)^{4.012007}) .dx

/ 0
Score: 5 points for each correct symbolic value.
5) (Any 12digit model): Given the succession defined thus:
X_{0} = 0
30*M^{2}  89*M + 2^{6}
X_{N} = FractionalPart( 2^{4}*X_{N1} +  )
2^{3}*M^{4}  2^{6}*M^{3} + 178*M^{2}  206*M + 84
where M = 2^{2}*N
write a program to compute the first 9 terms (X_{1}, X_{2}, ...,
X_{9}) and print
not each term X_{k} but IntegerPart(2^{4}*X_{k}).
Once this is correctly done, you must answer these questions about the integer values thus obtained:
 Give a description as accurate and short as possible
of the printed values in this succession, i.e., what are they ?
 Compute and output the next 9 additional (integer) values and state
whether you deem them correct or not.
 State whether you think the description you gave in (1) applies to
your newly computed extra terms, and to the theoretically exact infinite
sequence.
Score: Writing a working program: 4 points
Correct answer to (1): 4 points
Correct answer to (2): 1 point
Correct answer to (3): 1 point
6) (HP15C specific): Assuming A and B are both NxN matrices, their elements being integer values less than 5,000,000,000, write a program that will exchange their
contents, for NxN up to and including 5x5.
The program must take no input but assume both matrices already exist.
None of the 21 conditional instructions available in the HP15C are allowed.
Score: Writing a working program: 10 points
7) (Any model): Calculus tells us that the first derivative of a function
such as x^{3} is D(x^{3}) = 3*x^{2}
and its second derivative is D_{2}(x^{3})
= D(D(x^{3})) = D(3*x^{2}) = 6*x,
its third derivative would be 6, and finally its fourth derivative would be
the first derivative of 6, but the first, second, etc, derivatives of 6,
being a constant, are zero.
Let's consider the "Derivative" operator, D_{n}(f(x)), such that, for instance:
D_{1}(x^{3}) = 3*x^{2} > at x_{0}= 2 , D_{1}(x^{3}) = 3*2^{2} = 12
D_{2}(x^{3}) = 6*x > at x_{0}= 5 , D_{2}(x^{3}) = 6*5 = 30
D_{3}(x^{3}) = 6 > at x_{0}=Pi , D_{3}(x^{3}) = 6
D_{4}(x^{3}) = 0 > at x_{0}=4 , D_{4}(x^{3}) = 0

Write a program that, given m, n, and x_{0},
it computes and outputs the
value of D_{m}(x^{n}) at x_{0}. In particular, use
your program to compute the value of
D_{4.012007}(x^{Pi}) at x_{0} = e
 Write another program that, given m, a constant c,
and x_{0}, it computes and outputs the
value of D_{m}(c) at x_{0}. In particular, use
your program to compute the value of
D_{e}(4.012007) at x_{0} = Pi
(Strictly speaking, there are some combinations of m, n, and x_{0}
than aren't allowed because they would result in out of range values or division by zero,
but your program needs not address them for this subchallenge)
Score: For writing working programs: 3 points each
For giving the correct answers to the specific cases: 2 points each
8) (Any model): Write a program to find, compute, and print all square factorials N!
from N=1 to N=512.
Score: For writing a working program: 6 points
For correctly printing all results: 4 points
That's all. As usual, I'll post my solutions and comments within a week or so.
This is the time to test what you're made of as far as HPcalcs programming goes, so go ahead
full steam, do your best, engage, and may the force be with you.
Best regards from V.
Posts: 412
Threads: 40
Joined: Mar 2006
Hi Valentin,
I was sure that your April 1st edition of your S&SMC was a joke, and I read it for the fun. Your first challenge just looked crazy:
Abs(Ln(x*x)Ln(x^2)) > 4.012007
As x^2 IS x*x by definition, it seems that there can be no solution. But I know you, Valentin, and you would not post nosense. So I thought that some special cases could give different results for x*x and x^2 due to floating point rounding effects. I quickly thought that only complex numbers could give such large differences.
After investigating a bit, what was my surprise to find THIS:
X=(1E11,11E11)
X*X > (9.9999999998E21,1.99999999998)
X^2 > (9.9999999998E21,174532925.196)
We can call it a bug of the HP71B Math ROM, no?
I never heard of that. Incredible!
JF
Edited: 1 Apr 2007, 11:11 a.m.
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hi, Valentin 
Quite an ambitious effort! I remember one like number 6, so it wasn't too difficult to develop a solution:
Quote:
6) (HP15C specific): Assuming A and B are both NxN matrices, their elements being integer values less than 5,000,000,000, write a program that will exchange their contents, for NxN up to and including 5x5.
The program must take no input but assume both matrices already exist. None of the 21 conditional instructions available in the HP15C are allowed.
Assume that the two matrices are identified as A and B. The user is responsible for ensuring compatibility of dimensions, which need not be square.
LBL A
MATRIX 1
LBL 0
RCL A
x<> B
STO A (user mode)
GTO 0
RTN
The matrix indices are automatically incremented after a STO or RCL instruction when "USER" mode is set, but otherwise are not. This program takes advantage of a thoughtful feature of the HP15C: Within a program, if a matrix index is incremented past the last element back to (1, 1), the next instruction is skipped. This allows escape from a loop without a conditional test.
The feature is documented on pages 176177 in the HP15C Owner's Handbook.
This program also uses the versatile openended "x< > " function, which can be used with any numbered memory register, with the I register, with (i) for indirection, and with AE for matrix elements. It was a significant improvement over the HP11C's "x< >I" and "x< >(i)" functions, and also made a needed keyboard position available as well. ("x< > " was also unavailable on the HP32S, but restored for the HP32SII.)
BTW, How do you figure "21 conditional instructions"? There are 12 instructions that compare stack x with zero or stack y, plus 10 flag tests, plus DSE and ISG. That's 24 operations having a "do if true" rule, although "F? 8" should be avoided as a generalpurpose instruction.
 KS
Edited: 1 Apr 2007, 2:58 p.m.
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hi, JF 
Quote:
X=(1E11,11E11)
X*X > (9.9999999998E21,1.99999999998)
X^2 > (9.9999999998E21,174532925.196)
We can call it a bug of the HP71B Math ROM, no?
Probably so. The HP71B Math ROM was, I suspect, HP's first effort at handheld complexnumber functionality for 12digit arguments.
The successor HP42S and HP28C with their 12digit arguments and builtin complexnumber support give matching correct answers.
The successor HP32S does not offer "COMPLXx^{2}", so no comparison is possible.
The predecessor HP15C with its 10digit arguments and complexnumber support give matching correct answers for 1E09 and 11E09.
 KS
Posts: 536
Threads: 56
Joined: Jul 2005
using hplua on the 50g (counts as any calculator right?)
im thinking of using the digamma function instead of loop and also to get the realvalued harmonic function.
start with definitions:
harmonic(x) = sum(1,x,1/x),
euler = euler's constant,
digamma(x) = gamma'(x)/gamma(x) [prime = derivative]
and
harmonic(x) = euler + digamma(x+1)
since the digamma doesnt exist native, i try with a hacky numerical derivative like this,
function digamma(x)
local dx = 1e6
return (math.fact(x1+dx)  math.fact(x1dx))/2/dx/math.fact(x1)
end
which is probably not terribly accurate but could be improved.
then i get an answer for #3 using some functional programming and solve like this:
euler = 0.5772156649015328606
function harmonic(x)
return euler + digamma(x+1)
end
function inversefn(f, eps)
return function(x)
local ff = function(y) return f(y)  x end
return solve(ff, 1, 100, eps) end
end
transcript,
>hplua i solve.lua
Lua 5.1.1 Copyright (C) 19942006 Lua.org, PUCRio
HPLua version 0.2
dofile("valentin3.lua")
invdi = inversefn(harmonic, 1e10)
=invdi(4.012007)
> 30.523595225432980746
>
anywhere close?
Posts: 536
Threads: 56
Joined: Jul 2005
here's my program:
1
R/S
for there are no square factorials except 1.
proof:
Suppose n! is a perfect square, n != 1. let p be the largest prime dividing n!, then we must have p^2  n! and n! >= p^2 > 2p when p > 2.
But, there exist a prime q, 2p > q > p (bertran)
But, q  n! so p was not the largest prime. contradiction.
!
Posts: 1,619
Threads: 147
Joined: May 2006
#3 50G program:
<< 1 + Psi 1 Psi  >NUM >>
Enter
Psi(X+1)Psi(1)=4.012007
in NUM.SLV and press 'SOLVE' to get X=30.5235952258
#3 15C program:
1 LBL A
2 STO 0
3 0
4 ENTER
5 1
6 INTEGRAL 0
7 RTN
8 LBL 0
9 RCL 0
10 Y^X
11 CHS
12 1
13 +
14 X<>Y
15 CHS
16 1
17 +
18 /
19 RTN
Put 30.52359523 on stack, run A, and get: 4.012007.
Edited: 2 Apr 2007, 2:11 a.m.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, JeanFrançois:
Thanks for your interest and kind comments, and yes, your intuition was correct and as sharp as ever: complex numbers provide a correct solution for this subchallenge.
The key idea behind your solution is given in the May 1983 issue of the HP Journal, more specifically in the article "Scientific Pocket Calculator Extends Range of BuiltIn Functions" which in particular discusses the complex algorithms implemented for the HP15C (pp 2831). The concepts and compromises discussed there also apply somewhat to the HP71B implementation and probably to other HP calcs as well.
Now that you've found a correct solution, let's give you a new goal perfectly suited to your worthiness:
Your solution is indeed correct but on aesthetic grounds one could say your X is somewhat ungainly, what with its absolute value being very large (about 1E11), and its real and imaginary parts being 22 orders of magnitude apart. With so extreme a value, one would expect odd behaviour when using it in computations: it might exceed some limits here and there or incur in signigicantdigitcancellation problems.
That being so, see if you can find a solution which doesn't exceed 1 in absolute value and where such digitcancellations simply can't happen. If you succeed, that will tell you even more about the HP71B's complex number implementation. :)
Best regards from V.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Karl:
Karl posted:
"Quite an ambitious effort! I remember one like number 6, so it wasn't too difficult to develop a solution:"
Thanks for your interest and words of appreciation, but sorry, your alleged solution isn't, read below why.
The requirements for this subchallenge specifically state that "None of the 21 conditional instructions available in the HP15C are allowed." What's a "conditional instruction" ? I thought it was pretty obvious but perhaps it wasn't, so I'll make it clear here:
"A conditional instruction is any HP15C's programmable instruction that can perform a conditional jump over the next step depending on some condition checked by the instruction as part of its execution. Examples of such instructions are the TEST instructions, F?, etc."
Now, your alleged solution does include an STO A instruction in user mode. That's a conditional instruction as per the definition above, as it can skip the next step depending on some condition that it checks as part of its execution. So your alleged solution, most regrettably, isn't.
"BTW, How do you figure "21 conditional instructions"?"
As per the definition above. Also, you shouldn't count "10 flag tests" as ten different conditional instructions but merely one, F?.
Thanks again for your very didactic comments and worthy techniques, which certainly add value to the forum's community even if they didn't result in a valid solution in this case, and I hope you'll keep on contributing to this current S&SMC.
Best regards from V.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Hugh:
Hugh posted:
Best regards from V.
Posts: 412
Threads: 40
Joined: Mar 2006
Hello all,
I'm sorry, but there is no bug in the Math ROM, the behaviour I reported yesterday (on April 1st :) is only due to inherent limited accuracy of floating point calculations (and is not HP71B specific) that I tried to mask using a number chosen on purpose to show apparently very different results for x*x and x^2 (there are very close, actually).
And this number even doesn't satisfy Valentin's first problem (I didn't claim it does...):
X=(1e11,11e11)
ABS(LN(X*X)LN(X^2)) = 0 !
Sorry again...
A "solution" of the first problem is much simpler (and is specific to the HP71B):
X=(0,1)
ABS(LN(X*X)LN(X^2)) = 6.28 (about)
JF
Posts: 1,619
Threads: 147
Joined: May 2006
The above 15C code works for N > 200 too, however it could be off by 1e9.
Posts: 3,229
Threads: 42
Joined: Jul 2006
I'm pretty busy at the moment so won't be able to do much on this set of problems, However, how about something along these lines for the sixth:
RCL MATRIX A
RCL MATRIX B
RESULT B
+ // B = A+B
RCL MATRIX B
RCL MATRIX A
RESULT A
 // A = (A+B)A = B
RCL MATRIX B
RCL MATRIX A
RESULT B
 // B = (A+B)  B = A
I'm relying on the assumption that the matrix contents are integers and that a sum of two such won't overflow a register or lose any precision. This is justified given the problem statement.
 Pauli
Posts: 3,229
Threads: 42
Joined: Jul 2006
Quote:
Karl posted:
"BTW, How do you figure "21 conditional instructions"?"
As per the definition above. Also, you shouldn't count "10 flag tests" as ten different conditional instructions but merely one, F?.
Lets see:
x <= y?
x = 0 ?
TEST
F?
DSE
ISG
SOLVE
Integrate
matrix STO in USER mode
matrix RCL in USER mode
Even counting the TEST as the full 10, I'm two short.
 Pauli
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hi, Valentin 
Ah, there's the rub! I had several nagging doubts: Why you would essentially repeat an earlier challenge, and why integer elements of less than half the exactlydisplayable range were specified. I'd say that Paul Dale identified the reason, and that his approach is on track.
Of course, if a user wanted some HP15C program code to perform a matrix swap, why else would one use slow (and sometimes spacecostly) matrix arithmetic instead of "surgical elementwise replacement"?
To illustrate, here's an example from linear algebra. Suppose that a user were to calculate eigenvalues for small matrices on an HP15C (not that I'd advise it, due to computational slowness).
Using "r" to denote a right eigenvalue of a square matrix A with a corresponding right eigenvector x, we have the following familiar derivation:
Ax = rx
Ax  rx = 0
(A  rI)x = 0
Since the only eigenvector for a nonsingular (A  rI) would be the trivial solution x = 0, (A  rI) must be singular. Thus,
det (A  rI) = 0.
Calculation of (A  rI) from a known or assumed value of r could be performed by creating an identity matrix, multiplying it by r, then subtracting from A. Or, a simple loop program could subtract r from each value on the main diagonal of A. Which is more efficient in terms of computing speed and space requirements, both of which are rather limited on the HP15C?
I suspect that these are the reasons why no utility for generating identity matrices was provided on the HP15C. (Not to belabor the point, but the comparison of methods in your challenge offered a segue.)
Quote:
... Also, you shouldn't count "10 flag tests" as ten different conditional instructions but merely one, F?.
Point taken about conditional instructions, but I generally define instructions based on opcodes: If it has a unique code, it's a complete instruction. The basic command is a function. There are 700 programmable instructions on an HP15C, but far fewer functions.
 KS
Posts: 305
Threads: 17
Joined: Jun 2007
Quote:
"> 30.523595225432980746 anywhere close?"
Indeed. Your result agrees with my original oneline solution for the HP71B to 12 significant digits save for a few ulps.
The function to be used for this challenge was defined like this:
Quote: Let F(X) be a function defined thus for real positive X:
F(1) = 1
F(X) = F(X1)+1/X
for instance: F(2) = F(1)+1/2 = 1+1/2
F(3) = F(2)+1/3 = 1+1/2+1/3
F(4) = F(3)+1/4 = 1+1/2+1/3+1/4
...
F(100) = 1+1/2+1/3+1/4+1/5+...+1/98+1/99+1/100
... etc ...
The description "real positive X" includes nonintegers, although when I first looked at the definition, I couldn't see just how the recursion would terminate if X wasn't an integer. I suppose the process could be terminated just before the next term of the sequence becomes negative, but as I mention below, that stopping criterion doesn't always give a good result.
At any rate, if the value 30.523595225432980746 is substituted for X and the summation carried out, after adding 30 terms the sum is about 3.37647; after 31 terms, the sum is about 5.28634. With 32 terms, the sum begins to decrease, and its value is about 3.18724.
Using the function definition given in the original statement of the challenge, for X = 30.523595225432980746, the function is never equal to 4.012007 for any number of terms in the sum, so that for this part of the challenge:
Quote: 2.Use your program to find X such that F(X) = 4.012007
30.523595225432980746 is not a solution (if we stick to the given definition of the function).
However, if we let X = 29.98952530377276865 and add up 30 terms according to the given definition of the function, we will get a result very close to 4.012007. And, in fact, there are 30 values of X which will give 4.012007 with 30 terms added up. The values of X go something like this:
29.98952530377276865
28.425274519995517519
.
.
.
5.165990779772168
.
.
.
0.12209269005912
If we choose to add up 29 terms, there are 29 starting values of X that give 4.012007. If we choose to add up 28 terms, there are 28 values of X, etc.
Notice that for the small starting values of X and with 30 terms, there are quite a few terms in the sum that go negative.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Egan:
Egan posted:
01 LBL A
02 MATRIX 1
03 X<> 0
04 0
05 INTEGRAL 0
06 RTN
07 LBL 0
08 RCL 0
09 Y^X
10 RCL 1
11 1
12 Roll Up
13 
14 /
In your original program, even the use of register 0 can be avoided altogether if needed, so using no numbered registers at all.
Congratulations again and I look forward to further contributions from you.
Best regards from V.
Posts: 536
Threads: 56
Joined: Jul 2005
not sure how to get an exact symbolic value. i could "compute" pi/4 in symbolic mode, but i don't think that counts as evaluating the integral (although the problem doesn't say it requires the integrals to be evaluated).
progress made:
 I1 == I2, since in I1, let tan y = x to get I2 in y.
 replacing 4.012007 with k, and considering I as I(k), i can show I(2n) are all Pi/4, n >= 1. also I(2) == I(1).
 I(k) appears to be a constant k >= 1. haven't proved this. not even sure it's true. maybe k has to be rational?
idea:
prove I(4.012007) == I(1) or I(2) and find some dandy way to sym integrate the latter to Pi/4.
anybody?
Posts: 1,619
Threads: 147
Joined: May 2006
50g:
I1:
Put on the stack:
4: 0
3: Inf
2: 1/((1+X^2)*(1+X^4.012007))
1: X
Press Integral (Right Shift TAN), then press EVAL, results:
.785398163394
To get symbolic press:
CONVERT (Left Shift 6), F4 (REWRITE), NXT, F6 (>QPI), results:
(1/4)*PI
I2:
Put on the stack:
4: 0
3: PI/2
2: 1/(1+Tan(x)^4.012007)
1: X
Press Integral (Right Shift TAN), then press EVAL, results:
.785398163394
To get symbolic press:
CONVERT (Left Shift 6), F4 (REWRITE), NXT, F6 (>QPI), results:
(1/4)*PI
Posts: 1,619
Threads: 147
Joined: May 2006
I forgot to add a program for I1:
%%HP: T(3)A(R)F(.);
\<< 0. \oo 1. X 2. ^ + 1. X 4.012007 ^ + * INV X \.S EVAL \>Q\pi
\>>
I2 is equally as trival:
%%HP: T(3)A(R)F(.);
\<< 0. \pi 2. / X TAN 4.012007 ^ 1. + INV X \.S EVAL \>Q\pi
\>>
Edited: 3 Apr 2007, 8:39 p.m.
Posts: 1,619
Threads: 147
Joined: May 2006
Posts: 536
Threads: 56
Joined: Jul 2005
hi egan,
that's what i did first, but i wouldn't call that a sure answer because you have no proof that your answer is exactly pi/4. my interpretation of symbolic in this context is one in which approximations have not been used.
for example, maybe the answer is 22/83 + 11/34 + 37/188 which if you calculate on the 50g and press >QPI you also get pi/4.
my motivation for doubt is that if you replace 4.012007 with 4, straightforward integration is possible to give pi/4 and how can you prove the perturbation of .012007 does not, in fact, change the answer very very slightly away from pi/4.
what i've been trying to do is show replacing 4.012007 with 4 (or something else) does not alter the result and then go with that. peforming >QPI on the simplified integrand is then more credible.
any ideas?
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Hugh !
I do appreciate your mathematically orthodox reasoning but taking into account the particular nature of these "Spring Specials" I suggest some lateral thinking would be in order to deal with #8.
Best regards from V.
Posts: 182
Threads: 17
Joined: Oct 2005
Perhaps I'm thinking too lateral now, but your suggestion convinces me that there are exactly 512 square factorials that would do for an answer, each being the square of its square root. You aren't asking for integers, so ......
It kind of depends on the definition of "square", doesn't it?
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Bram:
"It kind of depends on the definition of "square", doesn't it? "
Yes, it does, and I'd further suggest that you go to the simplest definition (actually more like a mental image) you can think of for "square", actually what a young child would do; think of "Sesame Street"'s character Grover trying to explain the concept to a toddler audience :)
BTW, I've been missing your contributions, you frequently posted very interesting solutions and ideas in past S&SMC's.
Best regards from V.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Egan:
Congratulations, absolutely correct if a bit terse.
Would you care to elaborate the basis and genesis of your solution for the benefit of the gentle forum readers interested in this thread ?
Thanks for your interest and contributions and
Best regards from V.
Posts: 536
Threads: 56
Joined: Jul 2005
unless you mean all factorials of squares by the phrase "square factorials".
and indeed it does depend what you mean by a "square". i think it's up to the problem poster to define the problem properly if this is the case.
so there's bram's suggestion that square can apply to reals. in which case it might as well apply to rationals, complex and so on.
a square could also refer to being a quadratic residue. so n! is a square if there's an a, with a^2 = n! (mod m). for some m?
who knows?
Posts: 1,619
Threads: 147
Joined: May 2006
IANS, Benford's law. Simply put, for natural (random) events the probability for any starting n is log((n+1)/n)).
Seeing is believing. I generated as many 2^x and 3^x numbers as I could in a reasonable amount of time. (I used Perl, but with any 12 digit model you should be able to generate 1023 2^x numbers and 646 3^x numbers). After generating the 2 tables of numbers I computed the ratio of starting 1s, 2s, 3s, etc... within each set. Incredibly it matched log((n+1/n)), i.e. the number of 1s/1023 (for the 2^x case) and the number of 1s/646 (for the 3^x case) equaled log((1+1/1)). The same was true for the 2s, 3s, etc... Benford's law is not limited to only the starting digit, but any sequence of starting digits. After removing the decimals from my tables (the larger numbers were formatted n.nnn....ennn) I checked various random starting numbers, such as 12, 123, 55, etc... All ratios in both tables matched log((n+1)/n)), i.e., log(13/12), log(124/123), log(56/55)...
I concluded that if A<>B and i<>j, then P(A,M) = P(B,M), and therefore P(A,M)  P(B,M) = 0 for any M.
Edited: 4 Apr 2007, 10:43 a.m.
Posts: 1,619
Threads: 147
Joined: May 2006
You are probably right >Qpi should not be used. I also noticed that if I change 4.012007 to 4 or 9 for that matter that I still get 4/PI (in approx mode).
Edited: 4 Apr 2007, 10:58 a.m.
Posts: 182
Threads: 17
Joined: Oct 2005
Hi Valentin,
you remarked:
"BTW, I've been missing your contributions, you frequently posted very interesting solutions and ideas in past S&SMC's",
which is very kind to say. Be assured that I do read your S&SMC's and I do think about them, but all too often a properly worked out answer would require the time I must spend on other things. Yet I'm sure that some contributions will emerge in the future.
For now I seem to kind of have unlearnt the ability to think easily, so it's difficult to get these things square ;)
(I keep thinking about it, though)
Posts: 230
Threads: 11
Joined: Jan 1970
Hello, just to see if I understand your point.
First bizarre thing :
Any number starts with either 1, 2, 3, 4, 5, 6, 7, 8 or 9.
So since the probability of all these cases must be 1 (one), I get with that Benford's law :
1 = log(2)/1 + log(3)/2 + log(4)/3 + log(5)/4 + log(6)/5 + log(7)/6 + log(8)/7 + log(9)/8 + log(10)/9
That equality is very surprising.
Second bizarre thing :
 Is P(1,1) = 1 ?
 Is P(A,M) = P(B,M) for any A, B and M (true of course if A=B) ?
 If both are true, does it mean that for any B : P(B,1) = 1 ?
That last theorem is very powerful and contradicts badly my guts feelings.
As Valentin can testify, I can make mistakes...
Please put my mind at rest !
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, GE:
GE posted:
"As Valentin can testify, I can make mistakes..."
Yes, I testify. Here you are, a few more. Let's see:
"I get with that Benford's law :
1 = log(2)/1 + log(3)/2 + log(4)/3 + log(5)/4 + log(6)/5 + log(7)/6 + log(8)/7 + log(9)/8 + log(10)/9
That equality is very surprising."
First mistake: it is not:
1 = log(2)/1 + log(3)/2 + log(4)/3 + [...] + log(10)/9
but:
1 = log(2/1) + log(3/2) + log(4/3) + [...] + log(10/9)
which taking into account that, informally speaking, log(a)+log(b) equals log(a*b), then the expression is equivalent to:
2 3 4 5 6 7 8 9 10
1 = log( ********) = log(10)
1 2 3 4 5 6 7 8 9
which is hardly surprising because your "log" should be base10 logarithms, i.e. "lgt" in HP71B BASIC. "Log", unqualified, is normally understood to be natural logarithms, also notated as "ln".
"Second bizarre thing :
 Is P(1,1) = 1 ?
 Is P(A,M) = P(B,M) for any A, B and M (true of course if A=B) ?
 If both are true, does it mean that for any B : P(B,1) = 1 ?
That last theorem is very powerful and contradicts badly my guts feelings."
The problem's description specifies that neither A nor B can be a power of 10, and as it happens, 1 is a power of 10, namely 10^{0}, so Benford's law does not apply in that case. Any nonnegative integer power of 10 (1 included) begins exclusively by 1, 10, 100, ..., obviously.
I hope your mind is at rest now. Thanks for your interest and
Best regards from V.
Posts: 230
Threads: 11
Joined: Jan 1970
Thanks !
It was not obvious from E.Ford's message that he was using base 10 logarithms. I hope he did.
Also, it is quite strange to take into account 1 as a leading number both in its "1" form and "10" form. This accounts for my misunderstanding.
Last, it was not stated that Benford's law doesn't apply for powers of the base.
Again, let me say that your challenges are a real pleasure to bang one's head at, and this one is **especially** mindboggling. I can't even say which one is the most novel, intriguing piece of math (well #1 is a bit too machinespecific to my taste...).
Best regards
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, GE:
Thank you very much for your continued appreciation of my humble work.
Despite the long 6month pause since my latest S&SMC#17, this "Spring Special" (aka "April's Fool") S&SMC#18 hasn't been as successful as I hoped it would be.
Perhaps people are too tired of these mostly mathematical ramblings or maybe people are out enjoying the Easter holidays or whatever, but there's been less inputs than, in my humble opinion, these subchallenges deserved.
In particular, subchallenge #5 is a real marvel, both from the mathematical and the strictly computational points of view, with its amazing result and the hard fact that our calculators may 'deceive' us without we being even aware. It's also very important from a theoretical point of view, and if some facts about it could be proved, it would make worldwide headlines in the mathematical research world.
Yet, despite the additional fact that it's very easy to program, noone seems interested in it and it has received exactly zero inputs. A real pity. :(
Enough 'rants' :) Thanks and
Best regards from V.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Valentin,
Here is an attempt to number 5:
A simple 32SII program:
LBL A *
0 84
STO i +
STO X RCL M
LBL B x^2
1 30
STO+ i *
4 89
RCL* i RCL* M
STO M 
ENTER 64
ENTER +
ENTER x<>y
8 /
* 16
64 RCL* X
 +
* FP
178 STO X
+ LASTx
* IP
206 R/S
 GTO B
A: CK=57E9 006.0
B: CK=E781 063.0
XEQ A > 3
R/S > 2
R/S > 4
R/S > 3
R/S > 15
R/S > 6
R/S > 10
R/S > 8
R/S > 8
1. These are the first hexadecimal digits of pi:
3.243F6A88
2. The next 9 digits given by the program are:
8, 2, 3, 9, 6, 15, 5, 7 and 11
Only the first digit is this sequence is correct, due to rounding errors.
3. This would require further analysis. All I can say, the formulas gives at least 15 correct digits on sixteendigit machines.
Best regards,
Gerson.
Edited: 6 Apr 2007, 3:14 p.m.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Valentin,
Quote:
Yet, despite the additional fact that it's very easy to program, noone seems interested in it and it has received exactly zero inputs. A real pity. :(
I was tired and sleepy this afternoon, but since you have called our attention I gave it a try. No indepth analysis, just a quick program and the outputs but at least a start. I am looking forward to more solutions and, of course, to your final comments and solutions.
Best regards,
Gerson.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Gerson:
Your solution is correct, congratulations ! :)
I'll post my original solutions and comments within 48 hours, so you've got plenty of time to add whatever you'd like to.
By the way, I've noticed that you last edited your post exactly at 3:14 p.m. I assume this was fully intentional, right ? ;)
Thanks for your continued interest and kind appreciation and
Best regards from V.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi all,
It would be a real pity to left #8 "unattended" just because of some misunderstanding about the meaning of "square factorials".
I've given some advices and subtle hints about the proper meaning but either they weren't understood as well or simply there's no interest in #8.
Assuming the former, I'll give one last hint before I post my original solutions within 48 hours from now: "square" doesn't necessarily have to apply to the factorial's value. It may also apply to some other physical attribute, such as its shape when output.
Have you never heard of "triangular numbers', 'pentagonal numbers', 'polygonal numbers' in general ? Why do you think that square numbers were called precisely 'squares' back in ancient times ? Just try to dispose 16 marbles in your table to form some most regular pattern and you'll clearly see why ancient mathematicians would say that 16 is a square number ! Same for these factorials I'm asking for.
'Nuff said. Thanks for your interest and
Best regards from V.
Posts: 1,619
Threads: 147
Joined: May 2006
Quote:
Assuming the former, I'll give one last hint before I post my original solutions within 48 hours from now: "square" doesn't necessarily have to apply to the factorial's value. It may also apply to some other physical attribute, such as its shape when output.
I kind of figured that would be the case after the sesame street reference. Like the others I was confused about the term "square factorial". My first assumption matched Hugh's original solution. After that was shot down I turned to http://www.research.att.com/~njas/sequences/A100777 (Squarefactorial numbers). And wrote the following Perl prototype:
#!/usr/bin/perl
use Math::BigInt;
my $a = Math::BigInt>new("1");
foreach(1..512) {
print $_;
print ":\t";
$a = $a * l($_);
print $a;
print "\n";
}
sub l
{
my $n=shift;
for(my $i=int(sqrt($n));$i>0;$i) {
my $s=$i**2;
if($n % $s == 0) {
return($s);
}
}
}
However, I never got around to coding it in my 50G.
Your last two statements made me think about the shape as if each digit was a standalone unit (e.g. marbles).
The following 15C program computes the number of digits for N! for the range 1 to 512 then checks to see if the number of digits is a perfect square, if so print N.S, where N is 1 to 512 and S is the number of digits.
Output:
1.1
2.1
3.1
7.4
12.9
18.16
32.36
59.81
81.121
105.169
132.225
228.441
265.529
284.576
304.625
367.784
389.841
435.961
483.1089
508.1156
Code:
1 LBL A
2 1
3 .
4 5
5 1
6 2
7 0
8 1
9 STO 2
10 LBL 0
11 RCL 2
12 INT
13 GSB 1
14 STO 1
15 SQRT
16 FRAC
17 X=0
18 GSB 2
19 0
20 f ISG 2
21 GTO 0
22 RTN
23 LBL 1
24 STO 0
25 1
26 EXP
27 /
28 LOG
29 RCL* 0
30 RCL 0
31 RCL+ 0
32 PI
33 *
34 SQRT
35 LOG
36 +
37 1
38 +
39 INT
40 X=0?
41 1
42 RTN
43 LBL 2
44 RCL 2
45 INT
46 RCL 1
47 RCL 1
48 LOG
49 INT
50 1
51 +
52 STO I
53 10x
54 /
55 +
56 f FIX I
57 R/S
58 RTN
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hi, Valentin!
Quote:
By the way, I've noticed that you last edited your post exactly at 3:14 p.m. I assume this was fully intentional, right ? ;)
You're right! After I posted I noticed I had skipped one line in the 32SII program listing. When I finish to edit my posting I saw it was 3:11 p.m. so I decided to wait three more minutes before saving the edited post :)
Best regards,
Gerson.
P.S.: I should have used the extra three minutes to correct "the formulas gives" to "the formula gives" :)
Edited: 7 Apr 2007, 11:45 a.m.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Egan:
Congratulations, your interpretation and solution are both fully correct, there are 20 "square factorials" in the 1512 range. Two comments:
 It would be extranice if you would also print the full "square factorials" in a square arrangement to futher make the point and for visual impact. My original solution for the HP71B, which is quite short, actually does it.
 Your HP15C program admits a number of easy improvements. Just for instance, in the very first steps you set a constant for the counter, 1.51201. The final "01" is not necessary at all, if absent the "01" increment is assumed by default, saving 2 steps. Also, it would be better to use a decrementing counter, going from 512 down to 1 so the constant would be simply 512, additionally saving two extra steps.
Thanks for your contribution and
Best regards from V.
Posts: 1,619
Threads: 147
Joined: May 2006
Quote:
It would be extranice if you would also print the full "square factorials" in a square arrangement to futher make the point and for visual impact. My original solution for the HP71B, which is quite short, actually does it.
I figured as much. I'll post a 50G solution (no point in two 71B solutions) once I figure out how to display small fonts. The output will look like this (sans the extra space I added so that it looked more square here):
1
2
6
5 0
4 0
4 7 9
0 0 1
6 0 0
6 4 0 2
3 7 3 7
0 5 7 2
8 0 0 0
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
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
5 7 9 7 1 2 6 0 2 0 7
4 7 3 6 7 9 8 5 8 7 9
7 3 4 2 3 1 5 7 8 1 0
9 1 0 5 4 1 2 3 5 7 2
4 4 7 3 1 6 2 5 9 5 8
7 4 5 8 6 5 0 4 9 7 1
6 3 9 0 1 7 9 6 9 3 8
9 2 0 5 6 2 5 6 1 8 4
5 3 4 2 4 9 7 4 5 9 4
0 4 8 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
1 0 8 1 3 9 6 7 5 8 2 4 0
2 9 0 9 0 0 5 0 4 1 0 1 3
0 5 8 0 0 3 2 9 6 4 9 7 2
0 6 4 6 1 0 7 7 7 4 9 0 2
5 7 9 1 4 4 1 7 6 6 3 6 5
7 3 2 2 6 5 3 1 9 0 9 9 0
5 1 5 3 3 2 6 9 8 4 5 3 6
5 2 6 8 0 8 2 4 0 3 3 9 7
7 6 3 9 8 9 3 4 8 7 2 0 2
9 6 5 7 9 9 3 8 7 2 9 0 7
8 1 3 4 3 6 8 1 6 0 9 7 2
8 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
1 1 1 8 2 4 8 6 5 1 1 9 6 0 0
4 3 0 7 4 4 9 9 6 3 0 7 6 0 7
6 1 6 9 0 2 9 9 7 5 6 2 4 7 5
5 7 1 8 4 2 6 3 3 8 3 8 4 1 2
1 6 7 5 6 8 3 6 1 1 6 9 6 7 2
8 2 0 1 1 8 4 5 4 0 4 5 7 3 0
2 6 0 6 8 8 5 1 0 0 8 7 9 9 0
9 2 7 1 9 6 1 0 4 9 6 2 6 8 5
4 6 2 5 9 5 8 3 7 3 6 0 3 3 6
0 9 4 2 6 7 2 0 5 1 3 4 9 4 8
2 5 0 3 8 9 0 3 2 4 6 1 9 2 4
9 0 9 7 6 6 6 0 7 7 1 5 9 2 4
0 8 6 4 8 9 2 9 7 7 1 5 2 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 0 0 0 0 0
1 4 7 3 0 5 6 2 5 6 8 2 2 9 9 7 0 0 0 5 3
4 1 1 3 2 1 4 1 4 3 3 0 3 8 6 2 2 3 9 2 1
6 2 7 7 0 0 6 9 3 6 2 7 5 9 2 6 3 2 7 0 4
0 8 8 5 6 9 2 8 8 2 5 9 7 5 6 7 8 2 6 0 1
1 6 1 0 4 0 4 5 2 7 1 5 3 2 1 0 0 6 4 3 2
7 7 0 1 2 1 8 8 1 3 7 5 2 1 0 2 0 0 0 0 4
7 5 6 1 6 0 3 9 5 4 5 1 9 7 6 3 5 1 7 3 1
3 4 5 5 9 1 6 7 2 1 9 1 1 0 1 2 1 5 4 8 0
1 7 7 5 1 6 4 0 9 5 5 3 6 3 5 3 1 7 0 5 8
1 0 5 2 7 8 6 2 4 6 3 6 8 7 4 0 4 2 6 6 4
7 4 8 9 9 7 0 8 7 1 5 4 7 5 5 2 5 7 8 0 5
5 0 2 2 3 8 3 4 5 1 3 2 8 3 7 6 3 3 6 2 8
9 2 5 0 6 4 4 0 3 5 7 9 6 6 1 7 3 7 0 0 0
4 1 7 6 5 5 7 3 6 0 0 4 3 0 7 1 8 9 2 6 1
7 6 5 8 2 1 7 7 6 2 3 1 7 4 3 2 2 8 8 2 0
6 4 3 7 6 6 4 1 7 9 3 4 0 7 2 4 4 7 7 2 1
9 0 4 7 5 4 6 9 3 3 4 7 7 8 3 9 4 7 7 9 2
9 7 2 7 3 9 0 8 4 9 0 7 4 9 4 8 8 0 5 1 2
4 8 9 5 5 3 9 2 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 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 0
4 8 1 9 1 2 6 5 4 5 1 8 9 1 1 5 2 4 9 8 5 7 9
9 6 7 3 0 9 6 0 9 8 6 2 5 1 5 3 3 7 7 0 3 6 6
3 3 1 8 5 7 3 6 8 9 5 7 6 1 5 6 6 1 2 2 0 7 9
8 6 2 8 5 8 5 0 7 3 7 8 3 9 1 4 4 8 3 6 4 6 7
2 5 8 6 7 7 3 9 6 8 7 6 2 2 5 0 6 3 5 7 3 2 4
0 1 6 0 5 3 1 4 5 5 6 2 7 7 0 0 9 4 2 5 1 3 7
3 1 3 8 3 4 2 2 4 9 8 2 3 7 9 6 0 3 7 6 6 5 4
3 7 6 0 1 7 9 8 0 8 1 1 3 1 8 5 5 8 9 5 3 2 2
0 3 0 4 8 4 1 7 7 7 4 1 4 6 8 3 0 0 0 0 7 0 3
2 6 6 3 1 5 0 8 1 9 7 4 4 1 8 4 9 5 5 5 5 9 6
5 1 7 6 7 5 4 2 4 9 6 4 7 5 3 8 9 3 5 5 1 1 8
2 8 5 5 1 3 9 7 4 7 3 6 0 7 2 9 2 2 3 8 9 4 5
6 0 0 1 3 2 5 5 6 3 3 1 6 2 0 7 8 9 2 3 9 8 1
0 1 9 1 6 5 8 2 8 7 9 7 4 3 1 0 8 1 0 7 9 3 9
1 4 8 1 3 2 0 7 7 7 9 5 8 3 8 1 4 5 7 6 0 0 3
6 7 4 7 7 4 1 5 4 0 4 5 4 0 7 3 7 8 0 3 6 9 0
2 9 3 5 3 4 9 0 6 1 0 4 3 9 6 7 3 3 2 1 2 3 9
0 5 6 2 1 4 4 6 5 8 4 7 8 5 9 0 0 9 8 8 4 9 4
8 4 9 2 5 9 5 9 9 2 2 0 9 2 0 1 4 2 0 7 9 5 6
4 1 3 1 2 8 4 3 6 6 3 2 5 2 1 0 8 2 3 8 8 0 8
6 7 8 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 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 0 0 0 0 0 0 0 0 0 0 0
1 0 6 8 1 9 9 2 0 5 9 8 2 4 2 1 2 4 2 1 0 4 8 9
7 7 0 5 6 0 5 0 9 2 3 9 3 4 0 1 2 5 0 6 8 8 7 3
1 2 9 9 4 0 9 1 3 4 7 7 7 7 1 2 6 4 7 9 3 2 2 5
4 9 2 6 2 7 0 8 5 9 5 9 5 2 2 3 2 5 2 3 2 3 8 8
4 4 5 6 9 4 1 0 1 6 1 5 4 8 8 1 5 6 2 7 4 8 6 8
6 3 2 1 4 1 4 9 8 3 2 5 3 7 5 4 0 7 6 4 8 2 3 5
2 3 8 8 4 3 1 2 5 9 0 8 7 4 1 5 2 8 3 4 8 3 5 2
4 4 1 6 4 1 0 9 7 4 5 2 3 3 8 5 7 9 6 4 6 9 6 9
8 5 7 2 2 9 4 7 4 6 1 8 5 4 5 1 6 0 2 9 7 7 9 4
3 8 1 1 7 2 0 7 6 4 0 8 7 8 1 8 2 7 8 2 4 0 7 3
1 8 1 4 5 1 1 9 1 2 1 4 8 9 0 2 4 4 0 0 2 6 6 3
9 9 4 5 1 4 3 0 0 9 9 7 0 4 2 8 2 0 6 1 7 4 3 2
9 0 6 8 0 2 7 8 3 1 8 9 3 2 3 1 2 0 5 0 3 2 5 9
6 0 5 0 7 0 2 3 0 3 8 5 2 7 2 0 7 7 4 0 3 5 7 0
4 2 3 3 9 6 8 3 4 6 9 5 4 4 3 9 2 9 4 5 4 0 3 9
2 4 1 4 8 7 2 9 3 9 1 9 9 9 4 3 7 8 0 7 7 5 2 9
1 9 4 0 4 1 0 1 3 8 1 7 8 9 5 5 0 2 7 0 2 4 1 4
1 3 3 8 4 0 7 4 5 5 6 2 2 8 9 5 0 7 5 9 1 4 5 9
9 7 9 1 6 5 9 9 6 8 6 7 7 4 4 6 4 8 3 3 5 1 9 3
8 0 9 0 3 5 1 0 1 4 6 6 1 1 7 0 0 6 0 8 0 7 6 1
7 2 7 5 2 9 4 4 9 1 5 6 7 3 5 8 1 4 2 1 6 9 7 4
3 3 6 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 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 5 6 2 6 7 0 0 5 6 6 2 3 1 5 3 4 5 1 4 2 3 6 7 7
9 7 3 4 1 7 0 3 5 3 5 4 3 2 8 9 5 0 2 7 9 8 2 4 1
8 0 1 5 6 4 8 7 4 2 3 4 0 8 7 3 0 5 4 6 6 1 6 1 4
4 0 6 1 2 8 5 8 1 3 0 9 6 1 6 4 2 7 5 8 5 6 4 2 4
5 9 8 1 0 4 9 8 5 6 3 0 0 4 0 7 4 0 0 0 0 1 9 2 7
7 7 1 5 5 9 5 3 2 4 7 4 3 2 0 7 7 8 7 8 9 6 5 8 9
2 4 3 2 9 1 0 1 2 0 3 4 4 1 5 6 7 3 0 2 7 2 0 1 7
0 9 8 3 9 5 6 5 8 1 0 0 9 0 2 4 5 5 7 3 7 0 3 9 9
1 8 3 9 4 6 6 9 7 8 0 1 1 9 9 8 6 3 0 0 9 4 0 7 6
5 2 5 6 8 7 7 5 8 7 7 0 4 5 9 6 3 6 2 1 2 6 0 8 0
7 3 8 1 2 7 4 6 9 2 1 2 0 7 3 4 0 4 2 0 8 0 2 1 6
2 0 3 8 4 6 0 6 0 6 4 2 6 5 2 4 2 4 6 6 5 5 7 0 5
6 4 9 5 2 8 9 0 2 6 6 6 3 9 2 8 5 4 7 7 0 8 3 1 6
8 6 4 4 6 6 7 5 0 1 5 4 3 3 4 1 6 6 0 8 4 9 3 7 1
7 6 9 5 2 6 3 4 7 3 8 2 0 9 7 8 2 1 6 2 1 8 9 3 1
6 0 5 1 2 6 8 7 8 1 1 3 4 8 4 5 3 9 2 6 9 6 6 9 9
3 6 7 5 0 8 4 0 5 1 1 1 3 4 8 1 9 5 8 9 6 9 3 2 3
2 7 4 2 5 3 4 0 0 1 1 7 1 2 3 7 0 9 5 9 9 2 3 0 7
2 0 9 6 2 7 2 2 8 4 2 8 9 2 6 7 0 2 6 7 8 1 4 5 5
0 5 6 4 9 4 5 2 2 1 4 5 6 5 6 0 2 6 6 1 0 4 1 0 8
3 3 0 7 6 5 4 9 8 8 8 7 9 7 9 1 3 3 4 9 1 0 0 2 0
0 0 2 7 1 9 6 6 4 9 7 1 8 7 0 4 2 5 3 5 9 3 1 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 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3 3 7 2 0 3 6 7 7 1 9 5 8 4 0 0 0 3 9 0 5 0 8 6 8 7 4 3
2 2 7 9 2 9 1 8 5 4 1 3 2 4 3 0 1 3 6 9 0 2 9 4 7 8 6 4
4 5 8 6 0 7 3 1 6 3 5 1 8 4 3 7 8 6 9 8 1 1 2 7 5 6 0 8
4 3 7 8 8 7 1 0 0 4 9 3 2 1 4 0 5 3 8 9 4 6 2 1 7 6 0 9
5 9 5 8 6 8 3 5 0 3 6 3 6 2 9 5 2 9 4 3 3 3 0 2 4 3 3 8
7 4 4 8 8 7 0 8 4 9 4 8 4 9 4 4 6 3 0 9 8 7 9 4 8 8 4 4
3 8 6 8 8 2 8 2 3 0 7 9 6 6 0 5 3 1 1 7 8 5 3 5 5 7 7 6
6 6 0 0 5 5 0 6 9 8 5 8 3 5 0 3 3 2 9 2 2 1 6 5 3 9 1 7
0 5 6 7 2 3 4 9 5 6 9 7 4 7 4 0 3 7 3 5 4 6 1 4 2 7 8 5
7 7 0 3 1 7 4 5 7 9 8 6 6 5 5 5 6 1 6 1 6 0 4 3 6 1 2 7
0 0 2 5 9 2 2 6 5 6 2 1 1 8 4 4 1 2 8 6 3 0 0 2 1 4 3 9
7 1 2 5 8 7 9 0 4 5 7 4 8 8 4 1 8 1 1 2 0 0 0 6 4 7 7 4
2 0 2 1 8 8 6 7 9 3 4 0 9 3 2 3 1 8 3 6 8 7 5 4 2 5 7 2
3 0 3 7 9 9 4 3 6 8 1 0 1 8 8 7 0 6 1 2 0 3 7 5 7 4 4 0
1 1 8 6 6 4 9 0 6 4 6 1 1 5 1 0 7 2 1 1 6 0 9 5 6 1 4 1
6 2 4 0 1 1 1 9 2 3 4 9 2 3 0 0 5 2 3 4 9 8 9 3 5 6 9 1
5 7 6 4 1 4 2 5 9 6 2 4 4 0 4 6 1 9 6 2 6 6 5 7 5 1 1 6
2 0 7 6 6 9 7 1 3 9 3 5 4 4 3 7 4 7 6 6 9 0 0 8 3 2 7 8
2 8 4 5 7 1 6 8 2 1 2 7 7 4 1 9 4 1 2 2 6 6 1 0 1 4 9 1
7 0 7 0 1 4 6 7 0 4 3 0 5 8 1 9 3 7 2 1 0 4 4 8 1 5 2 0
0 3 5 1 1 9 2 1 9 4 1 6 1 8 2 6 3 0 2 8 0 8 3 0 5 7 9 5
9 8 3 5 4 1 9 0 8 0 3 2 8 8 0 9 5 2 3 5 1 3 3 9 3 9 9 6
9 5 0 8 3 8 4 7 2 3 1 3 6 2 0 4 3 1 8 6 7 6 8 6 7 1 4 1
6 0 6 0 4 7 4 5 2 3 0 9 0 9 0 8 0 0 9 0 8 9 7 9 1 5 6 2
4 5 4 7 4 8 0 7 1 3 0 6 3 6 7 9 8 6 5 6 5 1 2 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 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 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 0 0 0 0 0 0 0 0
1 7 5 3 8 7 7 3 7 8 8 1 5 9 5 0 4 1 6 3 8 6 2 0 2 4 7 4 0
2 6 1 6 6 5 3 2 3 8 8 4 1 3 6 6 0 1 9 7 5 0 9 0 6 4 2 7 2
0 8 1 5 1 8 0 2 1 1 4 7 2 7 1 5 4 2 6 0 1 4 4 3 6 5 6 3 0
2 9 6 9 5 7 8 4 0 1 0 2 7 4 6 0 4 4 8 8 7 7 9 8 0 0 7 9 4
7 3 8 6 2 0 8 1 0 6 3 5 5 5 3 7 7 3 1 7 8 4 3 3 1 9 1 5 3
0 5 4 6 8 2 3 6 6 5 8 7 1 5 6 1 0 5 5 4 6 2 0 6 3 8 6 1 4
6 2 9 6 7 7 2 5 2 0 3 6 1 0 5 4 5 9 6 0 1 3 1 1 2 2 2 8 3
9 5 4 4 9 1 9 5 2 5 1 5 8 7 2 3 8 0 3 1 3 5 0 6 7 6 6 5 6
2 5 8 5 4 8 9 6 2 5 1 7 1 0 6 8 8 3 9 4 0 7 7 2 9 2 7 3 4
7 9 8 1 9 2 1 7 2 5 8 7 6 7 1 1 4 0 5 4 8 5 4 3 6 7 6 8 4
8 7 8 5 1 7 8 3 7 6 3 3 4 6 5 0 0 8 7 6 0 6 5 0 3 9 3 0 9
1 2 1 8 5 4 5 7 1 1 4 8 1 0 9 5 0 1 7 7 4 7 7 7 2 0 2 4 2
6 6 7 9 5 4 9 9 8 8 0 5 2 9 3 6 5 7 4 8 0 4 6 3 6 1 5 8 7
6 9 1 3 4 5 8 7 9 2 9 0 1 8 5 8 4 5 5 2 2 4 8 4 0 1 8 7 8
4 8 4 6 3 4 8 9 5 6 2 1 8 6 1 9 8 1 2 8 4 2 6 4 3 1 5 0 6
5 6 2 7 8 2 5 5 1 7 9 6 2 6 9 0 4 8 4 9 5 6 2 2 7 3 5 0 8
6 3 5 7 7 7 2 8 0 4 4 7 2 7 6 0 1 4 9 1 1 6 5 2 9 4 9 0 6
2 7 1 1 0 2 7 4 2 3 4 3 5 7 6 4 5 9 3 9 9 5 1 2 7 2 3 2 2
6 3 6 5 5 2 4 3 9 3 9 3 0 6 1 8 6 4 8 2 1 6 5 6 0 2 6 4 3
2 9 8 8 7 5 9 1 7 8 3 2 6 6 0 2 4 0 7 2 5 1 3 4 1 1 4 2 2
6 3 4 8 6 0 6 7 7 6 2 2 1 1 6 0 5 6 4 1 7 8 0 4 9 8 3 9 3
9 5 9 2 9 2 3 5 9 5 2 0 3 8 9 8 5 9 1 9 6 0 2 0 1 1 2 7 0
8 4 0 3 4 5 5 9 2 0 6 0 9 7 9 6 0 3 3 6 0 8 4 0 0 7 4 5 3
4 3 4 5 3 6 9 1 9 6 4 6 6 6 8 9 1 1 6 0 1 4 2 9 6 7 5 9 7
8 3 4 7 1 2 5 7 0 9 4 9 8 9 2 6 3 3 2 0 1 7 9 9 6 8 4 7 3
1 7 7 2 1 5 2 1 6 3 9 5 6 7 8 1 2 1 9 8 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 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 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3 4 9 2 5 6 8 3 0 2 5 8 6 6 6 0 8 0 6 8 4 2 3 8 6 8 5 6 4 4 4
5 0 3 6 6 2 5 8 3 8 3 1 1 0 3 1 5 1 9 6 3 2 7 0 4 2 6 6 7 4 6
2 1 6 3 5 6 6 0 8 0 5 2 9 9 5 6 0 1 8 3 8 5 7 8 4 8 4 7 2 6 0
6 2 3 5 7 0 9 1 7 6 3 2 4 9 2 8 1 3 5 7 8 5 5 6 6 2 7 2 1 0 5
1 3 5 8 3 7 3 6 6 5 9 4 2 0 9 8 3 3 5 8 5 5 5 8 8 1 4 3 6 3 0
7 0 3 5 8 1 8 2 9 7 7 6 1 2 8 5 7 4 5 5 9 1 7 1 5 7 2 8 6 2 9
7 6 0 1 9 3 0 4 7 9 3 8 7 4 8 5 6 0 9 0 0 1 4 7 6 5 1 2 3 7 4
3 5 1 0 8 3 0 8 8 3 6 5 3 7 7 7 8 4 1 6 9 9 1 5 3 5 2 9 7 5 3
3 5 0 8 4 3 1 9 6 3 5 3 7 0 9 0 3 9 7 5 6 3 2 7 3 9 9 8 3 7 4
1 2 4 2 6 4 9 4 7 6 8 4 3 8 5 3 8 6 8 7 0 1 9 6 1 2 6 1 2 5 3
1 1 6 7 7 5 2 6 4 2 7 4 3 4 8 9 8 1 9 8 7 9 5 7 1 1 9 0 3 9 5
0 7 3 6 0 3 0 0 2 2 8 3 5 7 0 3 5 9 5 5 1 7 3 8 8 5 9 9 6 1 5
7 1 1 5 5 2 6 0 9 8 7 2 9 8 2 8 4 0 4 7 2 1 4 6 2 9 8 5 3 0 6
0 5 0 4 0 1 7 9 9 8 6 2 5 7 3 6 4 2 0 2 1 4 4 6 6 9 1 2 3 8 5
2 5 9 9 9 4 3 5 4 6 0 3 5 0 0 4 2 1 0 6 7 5 1 2 4 0 9 0 7 3 3
4 1 7 1 8 0 4 2 3 8 7 3 5 4 6 4 8 1 0 9 1 0 9 2 3 8 1 0 9 1 3
9 4 5 0 5 0 1 6 3 8 0 2 6 4 2 3 6 7 6 2 1 1 3 3 3 4 0 6 9 1 9
9 8 2 6 7 4 7 7 0 9 8 4 3 5 7 0 6 8 9 5 3 2 5 9 7 4 8 9 1 1 9
8 3 7 0 2 9 9 7 1 2 2 6 1 1 9 5 5 5 1 0 1 2 0 3 1 4 7 4 4 3 1
2 4 4 2 1 8 3 6 2 1 7 2 0 7 2 5 8 8 2 1 3 3 7 1 1 6 8 3 5 4 5
5 4 3 5 0 8 1 8 6 0 1 8 9 1 0 6 5 2 9 8 7 1 8 9 9 7 3 9 8 9 5
5 1 6 3 8 2 0 7 8 7 1 3 3 9 9 2 1 4 1 6 8 4 6 6 8 1 1 0 4 6 6
9 6 6 8 1 6 2 0 1 6 5 5 9 3 6 2 8 4 7 1 3 6 7 7 3 7 9 1 4 3 4
5 6 1 8 5 5 4 9 1 9 1 6 9 5 0 4 9 2 9 5 4 8 9 7 9 9 1 3 1 2 4
2 9 4 0 8 5 4 5 9 4 6 6 6 5 0 3 6 2 5 1 4 3 0 7 5 7 6 1 8 4 3
6 7 6 6 9 4 2 7 7 9 1 7 4 4 0 3 4 0 0 6 0 4 2 1 4 6 1 9 6 9 9
4 5 1 6 8 8 7 7 5 4 5 0 5 7 3 2 7 7 2 1 1 1 9 5 2 7 5 2 4 5 9
6 3 5 9 8 6 0 9 8 2 0 6 8 0 1 9 2 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 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 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 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
2 1 0 5 5 6 0 9 7 7 9 2 8 6 4 0 6 8 6 0 7 7 0 2 0 4 0 9 5 9 8 4 2
9 2 5 6 2 0 3 4 0 9 0 4 2 4 3 4 7 1 3 1 7 0 6 2 9 7 1 8 7 0 6 4 3
7 2 2 6 3 8 1 6 7 2 1 5 2 2 0 1 8 6 7 7 4 8 2 5 7 4 8 1 3 0 6 2 7
9 2 6 7 5 1 1 5 5 2 4 3 3 5 1 1 2 5 9 9 9 9 3 3 2 3 6 0 9 9 0 7 6
4 6 4 0 5 7 8 7 8 7 9 8 9 4 5 1 9 0 7 0 5 1 9 7 1 9 1 8 2 3 2 7 7
2 0 3 1 6 6 9 0 1 4 8 8 5 0 8 3 2 1 8 1 7 9 6 1 7 1 4 3 4 5 0 3 3
0 2 0 9 1 5 7 2 3 0 8 7 0 3 4 0 1 5 7 8 5 2 2 2 3 6 7 7 6 4 5 3 3
2 1 0 3 6 5 8 1 6 7 8 6 1 0 0 5 0 8 5 7 6 1 5 9 5 6 6 5 1 5 8 3 2
8 0 1 2 9 9 7 4 0 2 0 5 5 3 0 3 7 0 2 4 8 0 4 2 0 0 9 9 6 7 3 1 6
0 0 6 4 7 4 5 2 0 1 0 4 4 9 2 3 4 3 6 4 5 4 4 1 1 8 6 5 6 4 8 2 5
0 3 1 8 3 9 1 8 6 0 9 6 0 8 4 6 9 6 2 9 5 9 0 8 3 1 4 7 2 1 7 0 8
9 7 3 6 9 6 0 0 3 3 6 3 8 1 8 7 9 0 8 4 4 4 3 2 2 5 0 9 4 9 7 5 4
0 0 7 7 2 1 7 4 9 4 8 6 8 6 3 7 3 7 3 3 0 7 8 5 9 0 7 1 9 4 1 7 6
4 8 9 3 9 3 5 3 3 2 7 3 3 2 6 5 3 7 8 9 1 2 7 2 1 7 7 3 9 8 0 5 9
4 1 6 0 7 7 1 2 3 6 6 7 0 4 6 5 2 8 7 1 9 8 8 8 8 3 9 1 8 9 7 5 1
1 4 6 8 7 5 1 6 9 6 1 3 1 1 4 8 1 4 4 3 5 9 4 6 0 3 6 7 3 1 2 3 8
0 2 4 8 6 6 0 2 4 2 9 7 2 6 1 2 6 4 4 7 0 6 2 4 5 5 6 1 8 7 8 0 5
6 1 8 6 5 9 5 3 1 3 3 9 0 5 6 6 0 7 3 5 9 8 2 3 6 5 4 9 8 6 6 3 4
6 4 5 4 8 5 4 4 2 6 1 3 4 4 4 4 9 4 0 0 6 5 4 5 3 2 5 3 2 8 9 9 9
0 2 5 1 4 0 9 5 0 6 4 7 4 0 0 6 5 3 0 4 2 2 5 9 4 0 6 9 6 0 5 7 4
8 2 8 1 4 7 7 3 7 1 9 7 0 0 9 8 6 2 4 8 7 9 0 0 1 7 3 3 4 4 2 8 8
0 0 5 4 1 0 3 0 9 0 8 2 3 6 5 9 1 6 2 2 2 7 8 2 9 2 3 0 4 4 6 7 2
3 7 8 6 4 7 2 0 3 3 6 6 3 4 1 4 1 4 3 3 1 1 7 9 4 7 5 3 6 5 1 6 0
9 8 9 1 5 0 8 8 7 5 4 7 9 4 3 9 4 8 7 6 8 9 7 2 5 4 6 2 0 2 9 2 8
7 8 4 0 7 1 3 0 5 8 2 3 3 7 8 1 4 8 7 8 7 3 1 5 9 7 4 0 0 2 1 4 8
9 4 0 6 1 0 1 7 5 3 5 1 1 9 8 5 9 8 1 7 3 4 5 2 5 1 5 0 5 1 9 1 9
2 0 6 8 2 9 2 3 4 5 1 7 6 6 4 1 1 8 6 4 7 4 9 9 5 8 0 7 1 8 4 1 1
1 0 3 2 6 0 1 4 6 3 0 8 1 2 3 8 5 1 1 5 3 3 5 1 1 5 6 8 2 8 6 7 0
3 9 0 3 7 2 8 2 0 9 5 6 2 6 5 0 0 9 6 1 6 9 1 3 1 8 1 8 9 0 7 7 9
0 6 5 7 2 3 9 7 5 1 0 6 5 6 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 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 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 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 0 0 0 0 0 0 0 0 0 0
5 1 1 9 9 0 6 9 2 7 7 5 5 8 7 9 2 6 6 0 0 3 6 1 5 2 5 8 1 9 1 8 5 3
7 9 7 9 8 4 3 6 0 6 7 7 2 9 8 4 7 0 1 3 3 9 5 8 9 0 6 7 1 4 4 6 0 1
1 1 7 4 6 3 3 9 6 4 3 9 8 5 8 3 9 1 1 2 2 3 3 1 6 5 7 7 2 9 5 6 5 4
8 4 9 6 1 6 6 2 5 4 9 3 5 5 1 6 7 9 5 1 4 5 6 5 0 7 9 5 2 2 5 8 8 6
7 7 6 0 8 0 1 2 6 4 2 3 4 8 9 0 4 5 6 6 2 1 4 7 4 5 3 1 2 6 3 4 9 8
2 5 7 9 0 0 3 6 4 3 7 1 5 8 6 4 3 2 6 6 4 8 2 0 0 2 8 8 1 1 3 5 0 5
6 9 4 9 1 6 9 2 4 2 4 3 9 2 9 1 2 1 6 3 9 7 9 9 5 1 2 3 3 2 0 6 8 0
2 0 5 3 8 8 1 4 9 8 2 9 5 3 6 7 2 0 6 9 7 5 4 6 5 8 9 3 3 8 1 0 5 1
2 0 0 2 0 0 0 5 6 7 4 7 0 5 1 4 5 2 8 6 4 1 4 0 9 9 7 8 9 7 8 9 5 6
6 3 1 6 6 4 6 0 8 4 5 2 2 5 3 9 2 2 2 1 8 2 1 3 9 3 2 2 0 9 1 2 6 0
8 8 9 7 1 1 7 1 0 2 1 7 5 0 0 9 3 4 5 9 8 6 5 9 5 4 6 4 8 7 9 2 9 4
5 9 2 1 4 7 3 5 0 0 7 2 0 0 7 6 9 1 0 5 6 6 7 7 3 5 5 4 0 7 4 2 8 9
5 4 8 6 5 5 6 5 9 9 7 7 2 2 6 2 0 0 5 4 0 1 6 0 3 3 5 0 5 8 1 3 1 8
3 6 5 3 8 4 2 3 5 5 1 0 7 1 4 0 7 1 4 9 1 0 9 8 8 3 5 8 1 2 7 3 6 5
8 8 9 2 2 7 9 5 5 1 1 4 5 6 4 6 1 4 2 1 2 5 4 7 7 3 8 0 4 9 0 7 8 5
3 0 7 3 3 8 4 4 8 4 8 8 8 7 8 4 0 9 0 7 5 0 3 0 9 6 2 8 7 5 9 1 2 5
0 9 5 2 1 9 9 9 5 2 5 2 9 2 5 9 8 3 5 9 8 8 0 8 4 6 4 2 3 9 5 2 3 9
3 1 2 0 4 1 1 1 8 1 8 2 8 0 9 7 9 2 1 3 5 4 4 7 7 7 6 4 4 7 5 1 5 3
8 4 3 5 2 0 8 7 7 4 6 0 3 0 8 8 4 7 7 1 1 6 0 3 2 2 2 3 6 5 1 1 6 4
4 3 9 4 1 9 2 2 0 0 0 2 0 7 3 5 6 7 3 2 5 1 8 0 1 5 1 9 5 8 3 5 3 5
4 7 2 8 8 9 7 6 0 4 9 0 5 2 6 9 2 8 9 0 1 5 3 0 7 7 9 7 6 1 8 9 8 4
4 6 4 6 5 4 0 4 2 9 3 4 9 1 2 7 8 8 2 7 3 3 4 7 9 8 2 5 6 1 6 9 5 5
5 3 1 2 1 6 1 0 7 0 5 0 2 7 1 4 0 1 2 5 9 4 5 9 8 7 5 2 4 9 5 0 8 1
6 9 4 4 0 0 1 3 3 2 7 3 9 5 3 1 6 8 8 7 0 0 0 8 3 3 9 1 1 7 6 4 4 8
3 2 8 4 9 8 7 6 1 9 0 7 5 0 8 8 3 4 3 7 9 7 7 8 6 4 7 3 7 1 9 4 5 1
5 7 9 1 8 0 4 6 2 5 2 2 2 6 9 6 9 5 4 6 6 1 6 8 1 1 4 3 4 0 3 5 4 6
1 8 1 5 7 9 2 9 6 8 2 7 3 1 9 8 2 5 4 5 6 2 5 6 1 3 7 0 5 0 4 9 8 3
4 2 3 8 5 4 4 5 5 7 7 0 2 6 9 4 5 3 6 3 8 5 2 9 2 1 4 5 3 4 6 0 8 0
3 3 6 0 7 1 4 2 4 2 8 9 1 6 0 1 1 1 7 2 0 8 4 9 0 1 8 9 0 3 2 4 9 0
4 7 5 2 9 1 2 8 4 2 2 8 8 6 4 6 7 7 6 4 2 6 7 8 7 7 8 6 1 5 6 8 4 9
8 0 9 0 4 2 9 6 4 4 8 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 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 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 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Quote:
Your HP15C program admits a number of easy improvements. Just for instance, in the very first steps you set a constant for the counter, 1.51201. The final "01" is not necessary at all, if absent the "01" increment is assumed by default, saving 2 steps. Also, it would be better to use a decrementing counter, going from 512 down to 1 so the constant would be simply 512, additionally saving two extra steps.
Yeah, I know, but I wanted to count up.
Posts: 182
Threads: 17
Joined: Oct 2005
Ah! Right. So thát's what is meant with square numbers. I was not even near this solution.
It does make some sense, I must say.
And immediately a new question arises. The 20 answers can be connected to these numbers:
1 1 1 2 3 4 6 9 11 13 15 21 23 24 25 28 29 31 33 34
Quite irregular. Are there any known properties and behaviour of this series?
Posts: 1,619
Threads: 147
Joined: May 2006
Quote:
I'll post a 50G solution (no point in two 71B solutions) once I figure out how to display small fonts.
Ok, here it is:
%%HP: T(3)A(R)F(.);
\<< 1 \> N
\<< "" 1 512
FOR I I N * 'N' STO N \>STR DUP SIZE \v/ DUP FP
IF 0 >
THEN DROP DROP
ELSE \> S M
\<< 0 M 1 
FOR J S J M * 1 + J 1 + M * SUB 10 CHR + +
NEXT
\>> 10 CHR +
END
NEXT
\>> 1 \>GROB
\>>
Output (as seen on 50G):
Posts: 536
Threads: 56
Joined: Jul 2005
ha! well done Egan.
but valentin, you lateral dudes are so base 10, why stop there? a quick hack gives me,
function lngamma(z)
return (z.5)*math.ln(z)z+math.ln(math.pi*2)/2+1/12/z
end
function isSquare(x)
r = math.floor(math.sqrt(x))
return (r*r == x)
end
function findSquares(b, list)
c = 0
s = math.ln(b)
for i = 1, 512 do
y = math.ceil(lngamma(i+1)/s)
if (isSquare(y)) then
print(i, y)
list[i]=true
c = c + 1
end
end
print("total", c, "squares base", b)
end
function findAllSquares()
list = {}
for i = 2,10 do
findSquares(i, list)
end
c = 0
for k in pairs(list) do
io.write(k, ",")
c = c + 1
end
print("\ntotal squares in any base", c)
end
with transcript,
hplua i valentin8.lua
Lua 5.1.1 Copyright (C) 19942006 Lua.org, PUCRio
HPLua version 0.3
> findAllSquares()
1 1
8 16
17 49
37 144
162 961
216 1369
298 2025
309 2116
391 2809
416 3025
total 10 squares base 2
1 1
2 1
11 16
19 36
34 81
107 361
116 400
166 625
188 729
359 1600
405 1849
421 1936
437 2025
470 2209
487 2304
total 15 squares base 3
1 1
2 1
5 4
17 25
28 49
34 64
98 256
118 324
151 441
163 484
175 529
400 1444
436 1600
473 1764
total 14 squares base 4
1 1
2 1
14 16
19 25
31 49
38 64
45 81
53 100
70 144
133 324
145 361
257 729
289 841
358 1089
395 1225
414 1296
494 1600
total 17 squares base 5
1 1
2 1
6 4
10 9
15 16
41 64
49 81
97 196
120 256
145 324
158 361
186 441
298 784
316 841
353 961
372 1024
432 1225
453 1296
496 1444
total 19 squares base 6
1 1
2 1
3 1
6 4
11 9
16 16
22 25
52 81
61 100
71 121
92 169
128 256
155 324
169 361
199 441
319 784
358 900
378 961
441 1156
463 1225
total 20 squares base 7
1 1
2 1
3 1
6 4
11 9
23 25
30 36
46 64
55 81
97 169
109 196
122 225
135 256
149 289
194 400
210 441
total 16 squares base 8
1 1
2 1
3 1
7 4
17 16
24 25
31 36
39 49
57 81
67 100
78 121
114 196
171 324
237 484
255 529
312 676
353 784
396 900
441 1024
total 19 squares base 9
1 1
2 1
3 1
7 4
12 9
18 16
32 36
59 81
81 121
105 169
132 225
228 441
265 529
284 576
304 625
367 784
389 841
435 961
483 1089
508 1156
total 20 squares base 10
1,2,3,5,6,7,8,10,11,12,14,15,16,122,378,508,128,257,132,133,135,
391,265,395,396,17,145,19,22,23,405,28,30,31,32,414,34,416,163,37,
38,166,421,41,169,298,45,46,175,49,304,432,52,53,435,309,437,57,312,
186,441,188,316,319,194,70,71,453,400,389,367,78,284,55,81,463,210,
105,59,18,255,436,216,171,199,473,67,39,24,149,97,98,353,470,228,
483,155,358,359,487,107,118,109,237,61,92,494,114,496,116,158,372,
151,120,289,162
total squares in any base 116
Posts: 1,755
Threads: 112
Joined: Jan 2005
Awesome !
Best regards from V.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Hugh!
Hugh posted:
"[...] but valentin, you lateral dudes are so base 10, why stop there?"
Yes, why ? :)
Thanks for your interesting extension to bases other than 10, Hugh, it's nice to see people getting so involved with a challenge that they'll want to try and explore new possibilities beyond the original goals, as you've just superbly done. Way to go !
Thanks again and
Best regards from V.
Posts: 536
Threads: 56
Joined: Jul 2005
i have a symbolic solution for #4.
because the calculator can't solve the problem itself symbolically, my approach is to prove the problem is equivalent one that the calculator *can* solve symbolically. consider the following lemma,
now, by putting lambda = 4.012007, the lemma tells us we only have
to calculate the definite integral of 1/(1+x^2) between 0 and 1
so now the calculator part (using the 50g):
enter the integral (rpn mode):
ALPHA X 2 ^ 1 + 1/x
integrate symbolically,
CALC INTVX
substitute the limits and simplify,
0 1 CALC PREVAL
leaves the symbolic answer pi/4 on the screen.
for the second part of problem #4, the substituion tan y = x
transforms the integral into the first problem. therefore the answer is the same again.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Quote:
1 1 1 2 3 4 6 9 11 13 15 21 23 24 25 28 29 31 33 34
Quite irregular. Are there any known properties and behaviour of this series?
I'd say this sequence is just what it is: the sequence of the side lengths of square factorials. Anyway, if this can be of help, more terms can be obtained with the following HP50G program:
%%HP: T(3)A(R)F(.);
\<< 1 \> m n
\<< { } { }
DO 'FLOOR(LOG(2*\pi*n)/2+n*LOG(n/e)+LOG(1+1/(12*n)+1/(288*n^2)139/(51840*
n^3)571/(2488320*n^4)+163879/(209018880*n^5)))' \>NUM 1 + DUP \v/ FP 0 ==
IF
THEN \v/ + SWAP n + SWAP
ELSE DROP
END 1 'n' STO+ DUP SIZE
UNTIL m ==
END
\>>
\>>
Output for m=50:
{ 1 2 3 7 12 18 32 59 81 105 132 228 265 284 304 367 389 435 483 508 697 726 944 1011 1045 1080 1115 1187 1454 1494 1617
1659 1788 1921 2012 2105 2200 2248 2395 2445 2861 2915 3192 3480 3539 3902 3964 4476 4542 4675 }
{ 1. 1. 1. 2. 3. 4. 6. 9. 11. 13. 15. 21. 23. 24. 25. 28. 29. 31. 33. 34. 41. 42. 49. 51. 52. 53. 54. 56. 63. 64. 67. 68.
71. 74. 76. 78. 80. 81. 84. 85. 93. 94. 99. 104. 105. 111. 112. 120. 121. 123. }
That's the same formula used in this OEIS sequence:
http://www.research.att.com/~njas/sequences/A006488
The formula needn't be that accurate for this application. 'FLOOR(LOG(2*\pi*n)/2+n*LOG(n/e)+LOG(1+1/(12*n)))' would be accurate enough and faster, except the first number in the second sequence would be 0 instead of 1, unless properly handled. Actually, nothing new here, this appears to be equivalent to Egan Ford's HP15C program.
Regards,
Gerson.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi all,
Abs(Ln(x*x)Ln(x^{2})) > 4.012007
As JeanFranç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 HP71B 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 complexvalued functions
where they should be. This "x*x vs x^2" affair is but one such case.
JeanFranç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(270360*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
notsounexpectedly 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 outofrange results where
all other integer (or even real) values in 0360 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 A^{i} and B^{j}
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(X1) + 1/X
Write a program that, given any value N > 1 and N < 200 (for 10digit
models) or N < 1000 (for 12digit 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 wellknown factorial function, N!,
usually defined as
F(0) = 1
F(X) = F(X1) * 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 HP71B 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)+.577N)
20 DISP FNROOT(X1,X+1,INTEGRAL(0,1,1E12,(1(1IVAR)^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 wellknown 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) builtin function when
applied to the integral extension in order to appreciably reduce
computation time. The singleline 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

I_{1} =  1/((1+x^{2})*(1+x^{4.012007})) .dx

/ 0
/ Pi/2

I_{2} =  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

I_{2} =  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 builtin 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))
 (x^{4} + (1+2*Sqr(2))*x^{2} + 1)*(x^{100}  x^{98} + ... + 1)
/ 0
5) (any 12digit model) Given the succession defined thus:
X_{0} = 0
30*M^{2}  89*M + 2^{6}
X_{N} = FractionalPart( 2^{4}*X_{N1} +  )
2^3*M^{4}  2^{6}*M^{3} + 178*M^{2}  206*M + 84
where M = 2^{2}*N
write a program to compute the first 9 terms (X_{1}, X_{2}, ..., etc.)
10 X=0 @ Z=0 @ FOR N=1 TO 9 @ M=2^2*N
20 X=FP(2^4*X+(30*M^289*M+2^6)/(2^3*M^42^6*M^3+178*M^2206*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 HP71B program, which of course uses 12digit arithmetic, can't do
better and the very next terms are all wrong, starting with the 10^{th}
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 N^{th} term is the N^{th} 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 016, 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 10^{12} 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) (HP15C 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 11step program for the HP15C:
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 13step solution.
7) A very interesting and commonsensedefying subchallenge which no one took (perhaps Calculusthemed 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 soughtfor 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 9digit 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 nicelooking, prettyprinted "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 HP71B, which includes a multiprecision factorialcomputing 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+(I1)*D,D+(I1)*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 M1
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$,KLEN(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.
Posts: 536
Threads: 56
Joined: Jul 2005
hello valentin,
can you reveal the "master" formula?
Edited: 9 Apr 2007, 8:20 a.m.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Hugh:
Hugh posted:
Best regards from V.
Posts: 412
Threads: 40
Joined: Mar 2006
I think the unexpected result for LN(X*X)LN(X^2) is not directly linked to the signed 0 implemented in the HP71B.
It is perfectly consistent that:
LN((1,0)) = (1,PI)
and
LN((1,0)) = (1,PI)
even if the two numbers (1,0) and (1,0) are considered as equal:
((1,0)=(1,0)) is 1 (i.e. true)
BTW, this is an example of a violation of the rule that states that if X=Y then f(X)=f(Y). Another obvious violation is that 1/0 is different from 1/(0). Maybe this is why HP abandonned the signed 0 in the next Saturnbased machines.
We could argue that the branch cut for LN complex function should turn PI into +PI. But the HP71B numeric constant PI is an approximation and the cut is not exactly PI.
The root cause is, in my opinion, that
(0,1)*(0,1)=(1,0)
but
(0,1)^2 = (1,0)
instead of (1,0) due to rounding effects. If it was (1,0), no abnormal result would occur.
As a final point, this problem was specific to the HP71B in the handhed calculator domain, but if we include classic desktop HP calculators, we can find other examples, that have nothing to do with the signed 0:
For example, the 9816/26, and other HP9000 serie 200/300 calculators using the HP Basic (I used version 6.2 running on a Viper PC card):
10 COMPLEX X
20 X=CMPLX(0,1)
30 DISP X*X,X^2,LOG(X*X),LOG(X^2)
40 END
results are:
1 0 1 1.22E16 0 3.141592653589 0 3.141592653589
Anyway, congrat to Valentin for proposing this interesting and unexpected challenge!
JF
Edited: 9 Apr 2007, 9:34 a.m.
Posts: 68
Threads: 1
Joined: Jul 2005
Hi Valentin,
Thanks again for the entertaining challenge. I read through your solution and you point out that no one attempted #7. I searched online to see how noninteger derivatives are defined and found an explaination at the Mathworld site that led to the following two solutions that I hope are not too late for submission
For the first part of the challenge, here is a program that runs on the HP32S/32SII/33S:
LBL A
INPUT X
INPUT M
INPUT N
X!
RCL N
RCL M
X!
/
RCL X
RCL N
RCL M
Y^X
*
RTN
For M=4.012007, N=pi and X=e, this returns the value 0.414943481568.
For the second part of the challenge, I believe that the following program should suffice:
LBL B
INPUT C
INPUT M
INPUT X
RCL M
+/
Y^X
RCL C
*
RCL M
+/
X!
/
RTN
With M=e, C=4.012007 and X=pi, this program returns the result 0.068980430225.
Best Regards,
Eamonn.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Eamonn:
Eamonn posted:
"Thanks again for the entertaining challenge. I read through your solution and you point out that no one attempted #7. I searched online to see how noninteger
derivatives are defined and found an explaination at the Mathworld site that led to the following two solutions that I hope are not too late for submission"
Thanks a lot for your interest and kind words and no, it's never too late for submissions, the more (and for the most models) the better, we all can always learn something new or find some worthwhile routine, technique, or idea.
As for your programs and results, they're indeed fully correct and perfectly agree with my own original solutions and numeric values to full 12digit precision, congratulations.
Just for the record, this is my original solution for the HP71B, a pair of userdefined functions, callable from either programs or directly from the keyboard alone or as part of other expressions:
10 DEF FND(M,N,X)=GAMMA(N+1)/GAMMA(NM+1)*X^(NM)
20 DEF FND2(M,C,X)=C*X^(M)/GAMMA(1M)
The results for part 1, part 2, respectively:
>FND(4.012007,PI,EXP(1))
.414943481568
>FND2(EXP(1),4.012007,PI)
.068980430225
which fully agree with yours. The thrill of this particular subchallenge was twofold:
 On the one hand, the weird, unexpected concept of noninteger order derivatives, where most people, even the ones who feel comfortable with Calculus, have never seen nor heard about anything other than integerorder ones: first derivative, second derivative, ...
The mere fact that a noninteger order derivative can be defined, how to go about computing it, how to physically interpret it (if at all) and what its reallife uses might be (if any) are mindboggling questions to most mathematicallyaware people.
 On the other hand, the paradoxical concept that a constant function (which always has zero derivative for all positive integer orders), can have nonzero, nonlinear derivatives for positive noninteger orders, as can be seen in Part 2 of this subchallenge, where a noninteger derivative of the constant 4.012007 is evaluated for a given X value. This also seems to defy common mathematical sense, not only on computational grounds but on the very meaning of it all.
As you can see, this is a very deep subject that can't be easily dealt with here, so I'll leave it for a future Datafile article, most likely.
For now, it might be comforting to know that noninteger order derivatives can indeed be soundly defined, they can be efficiently computed, they do have physical interpretation, and most definitely, they do have important reallife applications.
The topic is not new (actually it's as old as Calculus itself), but it is now, with the advent of computers, when it finally gets to be studied and applied in full.
Thanks for your excellent solutions and
Best regards from V.
Edited: 10 Apr 2007, 10:50 a.m.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, JeanFrançois:
JeanFrançois posted:
"The root cause is, in my opinion, that (0,1)*(0,1)=(1,0) but (0,1)^2 = (1,0) instead of (1,0) due to rounding effects."
It has nothing to do with rounding, JF. Your statement seems to imply that multiplication of such small constants has no rounding errors whatsoever, while raising a complex number to the 2^{nd} power does entail some complicated operations which are prone to suffer from some rounding errors, however small, and this is what causes the incorrectly signed zero.
If that's what you mean, while plausible, let me assure you that this is not so, rounding has nothing to do with the differently signed zero, the same problem of incorrectly signed zeros does appear in operations not involving raising a complex number to a power but just a simple multiplication of small constants, say, where no rounding errors are possible.
"Anyway, congrat to Valentin for proposing this interesting and unexpected challenge!"
Thank you very much, I always appreciate your excellent inputs to my challenges and minichallenges and always will be indebted to you for Emu71, which is a most awesome HP71B emulator without which most of my articles and challenges would have never seen the light.
Best regards from V.
Posts: 536
Threads: 56
Joined: Jul 2005
hi valentin,
thanks for posting the pdf paper with the master formula, it is interesting. i am hoping you (or someone) can help me understand something about lemma 2.1. which says, if f(1/x) = x^2*f(x) then integral f(x)/(1+x^b) from 0 to infinity is independent of b. a detailed proof is not given only the phrase; the proof follows easily from differentiation wrt b and put x>1/x.
This is what i originally tried with your problem #4. the the end i managed to prove it another way because i got stuck with trying to show dI/db == 0 which is what i presume is done in the paper.
from lemma 2.1, if i try the "follows easily" bit, i get this:
but how is this necessarily zero, or have i made a mistake?
any help?
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Hugh:
Hugh posted:
"i am hoping you (or someone) can help me understand something about lemma 2.1 [...] a detailed proof is not given only the phrase; the proof follows easily from differentiation wrt b and put x>1/x."
Well, I assume you're experienced enough in all things mathematical to know that
"... the proof follows easily ..."
actually means
"... only to PhD's who specialize in that field, or to instructors who have taught the course 100 times ..."
i.e. the resulting integral after differentiation under the integral sign indeed is identically zero, but the above actual meaning fully applies.
Fortunately with a little ingenuity and experience in formal manipulation of mathematical entities it's perfectly possible to get away with the awkward differentiation process and concoct instead a much nicer demonstration of your Lemma 2.1, in the sense that all steps are crystalclear and I'll get you directly to a thistimefullyobvious result. Bear with me:
 Assuming that f(x) is a function that satisfies the functional equation f(1/x) = x^2 . f(x), I'll try and demonstrate that the value of the definite integral:
/Inf
 f(x)
 .dx
 1 + x^b
/0
does not depend on b.
 To that effect, I first split the integral into two parts:
/Inf /1 /Inf
 f(x)  f(x)  f(x)
 .dx =  .dx +  .dx
 1 + x^b  1 + x^b  1 + x^b
/0 /0 /1
 Now I left the first part alone and perform the change
of variable:
x = 1/t (so dx = 1/t^2)
on the second part, so that (after cosmetically replacing on the fly the dummy variable t for x again in order to keep the same variable in both parts) the sum now becomes:
/Inf /1 /1
 f(x)  f(x)  f(1/x) 1
 .dx =  .dx +  . .dx
 1 + x^b  1 + x^b  1 + (1/x)^b x^2
/0 /0 /0
 Now, still focusing our attention on the second part, I'll multiply both numerator and denominator by x^b so that I'll have:
/Inf /1 /1
 f(x)  f(x)  x^b.f(1/x) 1
 .dx =  .dx +  . .dx
 1 + x^b  1 + x^b  x^b + 1 x^2
/0 /0 /0
 Now it's time to make use of the initial assumption, namely that f(1/x) = f(x) . x^2 . Substituting this equivalence in the second part and arranging its denominator in the proper order, I finally get:
/Inf /1 /1
 f(x)  f(x)  x^b. f(x) x^2
 .dx =  .dx +  . .dx
 1 + x^b  1 + x^b  1 + x^b x^2
/0 /0 /0
/1 /1
 f(x)  x^b. f(x)
=  .dx +  .dx
 1 + x^b  1 + x^b
/0 /0
/1 /1
 f(x) + x^b.f(x)  f(x).(1 + x^b)
=  .dx =  .dx
 1 + x^b  1 + x^b
/0 /0
/1

=  f(x).dx

/0
which does not depend on b, Q.E.D.
Nice, uh ? :)
Thanks for your interest in these mathematical ramblings o'mine and
Best regards from V.
Posts: 536
Threads: 56
Joined: Jul 2005
hi valentin,
thanks very much for your detailed explanation. breaking up the integral is exactly what i did in my original post for the solution to #4, but i didn't have the general f(x) bit, only your original integrand. differentiating did not lead me anywhere helpful.
the paper, then, is a bit misleading because, the claimed, differentiation does not lead to a simple answer. i wondered if the authors knew of some "standard" result that allowed the differentiated version to be dismissed as zero quickly.
thanks again and thanks for your interesting challenges.
Posts: 230
Threads: 11
Joined: Jan 1970
Thanks to all (and Valentin in the first place), this challenge was very hard (and thus very interesting), IMHO this explains the lower than usual number of replies, not at all any lack of interest.
I know I have been reading this forum morning and evening for news on the SSMC. Very much enjoyed !
There are still lots of thing not fully explained...
Posts: 68
Threads: 1
Joined: Jul 2005
Quote: The thrill of this particular subchallenge was twofold:
* On the one hand, the weird, unexpected concept of noninteger order derivatives, where most people, even the ones who feel comfortable with Calculus, have never seen nor heard about anything other than integerorder ones: first derivative, second derivative, ...
The mere fact that a noninteger order derivative can be defined, how to go about computing it, how to physically interpret it (if at all) and what its reallife uses might be (if any) are mindboggling questions to most mathematicallyaware people.
* On the other hand, the paradoxical concept that a constant function (which always has zero derivative for all positive integer orders), can have nonzero, nonlinear derivatives for positive noninteger orders, as can be seen in Part 2 of this subchallenge, where a noninteger derivative of the constant 4.012007 is evaluated for a given X value. This also seems to defy common mathematical sense, not only on computational grounds but on the very meaning of it all.
I completely agree. I had not come across such definitions before, so it was a very interesting to see that it is possible to define noninteger order derivatives. And the fact that the noninteger order derivative of a constant is not necessarily zero was also counterintuitive, but it does follow from the definition of the noninteger order derivitaves.
Thanks again Valentin.
Eamonn.
