Short & Sweet Math Challenge #18: April 1st, 2007 Spring Special



#62

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 brand-new 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 sub-challenges, 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 model-specific questions apply only to that model in particular. They
      would probably be meaningless for any other models.

    • If you use an HP-41C, 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 HP-71B, you can also use any of these: Math ROM, HP-IL ROM,
      STRNGLEX. No other ROMs or LEX files allowed.

    The Challenge:


    1) (HP-71B specific) Find a root x of:

            Abs(Ln(x*x)-Ln(x2)) > 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(X-1)+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 ...

    1. Write a program that, given any real value N > 1 and N < 200 (for 10-digit
      models) or N < 1000 (for 12-digit models), it will find
      the value of X (with full 10-12 digit accuracy) such that F(X) = N.

    2. 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+x2)*(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 12-digit model): Given the succession defined thus:

          X0 = 0
    30*M2 - 89*M + 26
    XN = FractionalPart( 24*XN-1 + ------------------------------------ )
    23*M4 - 26*M3 + 178*M2 - 206*M + 84

    where M = 22*N

    write a program to compute the first 9 terms (X1, X2, ...,
    X9) and print
    not each term Xk but IntegerPart(24*Xk).
    Once this is correctly done, you must answer these questions about the integer values thus obtained:

    1. Give a description as accurate and short as possible
      of the printed values in this succession, i.e., what are they ?

    2. Compute and output the next 9 additional (integer) values and state
      whether you deem them correct or not.

    3. 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) (HP-15C 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 HP-15C are allowed.

      Score: Writing a working program: 10 points


    7) (Any model): Calculus tells us that the first derivative of a function
    such as x3 is D(x3) = 3*x2
    and its second derivative is D2(x3)
    = D(D(x3)) = D(3*x2) = 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, Dn(f(x)), such that, for instance:

             D1(x3) = 3*x2   ->  at x0= 2 , D1(x3) = 3*22 = 12
    D2(x3) = 6*x -> at x0= 5 , D2(x3) = 6*5 = 30
    D3(x3) = 6 -> at x0=Pi , D3(x3) = 6
    D4(x3) = 0 -> at x0=-4 , D4(x3) = 0
    1. Write a program that, given m, n, and x0,
      it computes and outputs the
      value of Dm(xn) at x0. In particular, use
      your program to compute the value of
                        D4.012007(xPi) at x0 = e
    2. Write another program that, given m, a constant c,
      and x0, it computes and outputs the
      value of Dm(c) at x0. In particular, use
      your program to compute the value of
                        De(4.012007) at x0 = Pi

      (Strictly speaking, there are some combinations of m, n, and x0
      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 sub-challenge)

      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 HP-calcs programming goes, so go ahead
    full steam, do your best, engage, and may the force be with you.

Best regards from V.

#63

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=(1E-11,1-1E11)
X*X --> (-9.9999999998E21,-1.99999999998)
X^2 --> (-9.9999999998E21,-174532925.196)

We can call it a bug of the HP-71B Math ROM, no?

I never heard of that. Incredible!

J-F

Edited: 1 Apr 2007, 11:11 a.m.


#64

Hi, J-F --

Quote:
X=(1E-11,1-1E11)
X*X --> (-9.9999999998E21,-1.99999999998)
X^2 --> (-9.9999999998E21,-174532925.196)

We can call it a bug of the HP-71B Math ROM, no?


Probably so. The HP-71B Math ROM was, I suspect, HP's first effort at handheld complex-number functionality for 12-digit arguments.

The successor HP-42S and HP-28C with their 12-digit arguments and built-in complex-number support give matching correct answers.

The successor HP-32S does not offer "COMPLXx2", so no comparison is possible.

The predecessor HP-15C with its 10-digit arguments and complex-number support give matching correct answers for 1E-09 and 1-1E09.

-- KS

#65

Hi, Jean-Franç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 Built-In Functions" which in particular discusses the complex algorithms implemented for the HP-15C (pp 28-31). The concepts and compromises discussed there also apply somewhat to the HP-71B 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 signigicant-digit-cancellation problems.

      That being so, see if you can find a solution which doesn't exceed 1 in absolute value and where such digit-cancellations simply can't happen. If you succeed, that will tell you even more about the HP-71B's complex number implementation. :-)

Best regards from V.
#66

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 HP-71B 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=(1e-11,1-1e11)
ABS(LN(X*X)-LN(X^2)) = 0 !

Sorry again...

A "solution" of the first problem is much simpler (and is specific to the HP-71B):

X=-(0,1)
ABS(LN(X*X)-LN(X^2)) = 6.28 (about)

J-F

#67

Hi, Valentin --

Quite an ambitious effort! I remember one like number 6, so it wasn't too difficult to develop a solution:

Quote:
6) (HP-15C 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 HP-15C 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 HP-15C: 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 176-177 in the HP-15C Owner's Handbook.

This program also uses the versatile open-ended "x< > " function, which can be used with any numbered memory register, with the I register, with (i) for indirection, and with A-E for matrix elements. It was a significant improvement over the HP-11C's "x< >I" and "x< >(i)" functions, and also made a needed keyboard position available as well. ("x< > " was also unavailable on the HP-32S, but restored for the HP-32SII.)

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 general-purpose instruction.

-- KS


Edited: 1 Apr 2007, 2:58 p.m.


#68

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 HP-15C 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 HP-15C'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.

#69

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

#70

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 exactly-displayable 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 HP-15C program code to perform a matrix swap, why else would one use slow (and sometimes space-costly) matrix arithmetic instead of "surgical element-wise replacement"?

To illustrate, here's an example from linear algebra. Suppose that a user were to calculate eigenvalues for small matrices on an HP-15C (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 HP-15C?

I suspect that these are the reasons why no utility for generating identity matrices was provided on the HP-15C. (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 op-codes: If it has a unique code, it's a complete instruction. The basic command is a function. There are 700 programmable instructions on an HP-15C, but far fewer functions.

-- KS


#71

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

#72

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 real-valued 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 = 1e-6
return (math.fact(x-1+dx) - math.fact(x-1-dx))/2/dx/math.fact(x-1)
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) 1994-2006 Lua.org, PUC-Rio
HPLua version 0.2

dofile("valentin3.lua")
invdi = inversefn(harmonic, 1e-10)
=invdi(4.012007)
> 30.523595225432980746
>

anywhere close?


#73

Hi, Hugh:

Hugh posted:

    "using hplua on the 50g (counts as any calculator right?)"

      Certainly, it does. The 50g is an HP calc model and anything that runs on it is valid and legal for this S&SMC. Matter of fact, it's most welcome, as it provides diversity and further shows the forum's community which tools are available and what can be done with them in a practical, not contrived situation.

    "> 30.523595225432980746 anywhere close?"

      Indeed. Your result agrees with my original one-line solution for the HP-71B to 12 significant digits save for a few ulps.

    Thanks for your interest, congratulations for your correct and truly novel solution and I hope you'll consider trying your hand at some other subchallenges as well (or all!) :-)
Best regards from V.

#74

Quote:

"> 30.523595225432980746 anywhere close?"

    Indeed. Your result agrees with my original one-line solution for the HP-71B 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(X-1)+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 non-integers, 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.

#75

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.

!


#76

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.


#77

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?


#78

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.

#79

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?

#80

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)

#81

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.

#82

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 (Square-factorial 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

#83

Hi, Egan:

    Congratulations, your interpretation and solution are both fully correct, there are 20 "square factorials" in the 1-512 range. Two comments:

    • It would be extra-nice 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 HP-71B, which is quite short, actually does it.

    • Your HP-15C 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.


#84

Quote:
It would be extra-nice 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 HP-71B, 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 HP-15C 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.


#85

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?


#86

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) 1994-2006 Lua.org, PUC-Rio
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


#87

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.
#88

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 HP-50G 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 HP-15C program.

Regards,

Gerson.

#89

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):


#90

Awesome !

Best regards from V.

#91

#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.


#92

The above 15C code works for N > 200 too, however it could be off by 1e-9.

#93

Hi, Egan:

Egan posted:

    "Put 30.52359523 on stack, run A, and get: 4.012007."

      Yes, very nice and short solution, congratulations ! It agrees with mine to 12 digits save a few ulps, and your HP-15C program does indeed return 4.012007 in FIX 6, taking not too long.

      By the way, you can easily make your 19-step HP-15C program a full five steps shorter (and slightly faster as well) by taking advantage of the HP-15C's extensive instruction set, as follows:

           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.
#94

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?


#95

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

#96

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.

#97

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?


#98

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.


#99

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.

#2

15C:

LBL A
0

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.

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.


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 !


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 base-10 logarithms, i.e. "lgt" in HP-71B 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 100, so Benford's law does not apply in that case. Any non-negative 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.

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** mind-boggling. I can't even say which one is the most novel, intriguing piece of math (well #1 is a bit too machine-specific to my taste...).



Best regards


Hi, GE:

    Thank you very much for your continued appreciation of my humble work.

    Despite the long 6-month 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.

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 in-depth 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.

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 sixteen-digit machines.

Best regards,

Gerson.

Edited: 6 Apr 2007, 3:14 p.m.


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.


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.

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.)


      The following routine for the HP-71B implements the above recurrence:

          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.


hello valentin,

can you reveal the "master" formula?

Edited: 9 Apr 2007, 8:20 a.m.


Hi, Hugh:

Hugh posted:

    "Can you reveal the "master" formula ?"

    Thanks for your interest and
Best regards from V.

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?


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 crystal-clear and I'll get you directly to a this-time-fully-obvious result. Bear with me:

      1. 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.

      2. 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
      3. 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
      4. 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
      5. 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.

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.

I think the unexpected result for LN(X*X)-LN(X^2) is not directly linked to the signed 0 implemented in the HP-71B.
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 Saturn-based machines.

We could argue that the branch cut for LN complex function should turn -PI into +PI. But the HP-71B 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 HP-71B 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.22E-16 0 3.141592653589 0 -3.141592653589

Anyway, congrat to Valentin for proposing this interesting and unexpected challenge!

J-F


Edited: 9 Apr 2007, 9:34 a.m.


Hi, Jean-François:

Jean-Franç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, J-F. Your statement seems to imply that multiplication of such small constants has no rounding errors whatsoever, while raising a complex number to the 2nd 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 mini-challenges and always will be indebted to you for Emu71, which is a most awesome HP-71B emulator without which most of my articles and challenges would have never seen the light.

Best regards from V.

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 non-integer 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 HP-32S/32S-II/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.


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 non-integer
    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 12-digit precision, congratulations.

      Just for the record, this is my original solution for the HP-71B, a pair of user-defined 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(N-M+1)*X^(N-M)
      20 DEF FND2(M,C,X)=C*X^(-M)/GAMMA(1-M)

      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 two-fold:

      • On the one hand, the weird, unexpected concept of non-integer order derivatives, where most people, even the ones who feel comfortable with Calculus, have never seen nor heard about anything other than integer-order ones: first derivative, second derivative, ...

        The mere fact that a non-integer order derivative can be defined, how to go about computing it, how to physically interpret it (if at all) and what its real-life uses might be (if any) are mind-boggling questions to most mathematically-aware people.

      • On the other hand, the paradoxical concept that a constant function (which always has zero derivative for all positive integer orders), can have non-zero, non-linear derivatives for positive non-integer orders, as can be seen in Part 2 of this subchallenge, where a non-integer 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 non-integer order derivatives can indeed be soundly defined, they can be efficiently computed, they do have physical interpretation, and most definitely, they do have important real-life 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.


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...

Quote:
The thrill of this particular subchallenge was two-fold:

* On the one hand, the weird, unexpected concept of non-integer order derivatives, where most people, even the ones who feel comfortable with Calculus, have never seen nor heard about anything other than integer-order ones: first derivative, second derivative, ...

The mere fact that a non-integer order derivative can be defined, how to go about computing it, how to physically interpret it (if at all) and what its real-life uses might be (if any) are mind-boggling questions to most mathematically-aware people.

* On the other hand, the paradoxical concept that a constant function (which always has zero derivative for all positive integer orders), can have non-zero, non-linear derivatives for positive non-integer orders, as can be seen in Part 2 of this subchallenge, where a non-integer 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 non-integer order derivatives. And the fact that the non-integer order derivative of a constant is not necessarily zero was also counter-intuitive, but it does follow from the definition of the non-integer order derivitaves.

Thanks again Valentin.

Eamonn.


Possibly Related Threads...
Thread Author Replies Views Last Post
  Need help understanding math.... cyrille de Brébisson 9 1,210 12-13-2013, 02:23 AM
Last Post: Didier Lachieze
  HP Prime - Short "learning" modules CR Haeger 1 450 11-27-2013, 02:13 PM
Last Post: Jonathan Cameron
  I have written a short introduction to the HP Prime Michael Carey 7 985 11-18-2013, 08:04 PM
Last Post: Michael Carey
  HP-65 short circuit Ignacio Sánchez 2 519 10-22-2013, 08:27 AM
Last Post: Ignacio Sánchez Reig
  OT: a math competition site Pier Aiello 0 324 09-16-2013, 06:03 AM
Last Post: Pier Aiello
  Simple Math Question Namir 2 463 08-09-2013, 06:13 PM
Last Post: Eddie W. Shore
  Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand 3 512 06-29-2013, 03:39 PM
Last Post: Gjermund Skailand
  Cool math clock Bruce Bergman 28 2,292 04-10-2013, 03:13 AM
Last Post: Siegfried (Austria)
  NCTM, Denver April 17-20 Tim Wessman 1 298 03-16-2013, 06:41 AM
Last Post: bluesun08
  FRAM71 for HP-71B, short update #3 Hans Brueggemann 15 1,247 01-20-2013, 10:22 AM
Last Post: Jerry Raia

Forum Jump: