Short & Sweet Math Challenge #21: Powers that be

[Valentin Albillo]

 
Hi all:

About a week has elapsed since I posted S&SMC#21 and in that time some of you have contributed a number of interesting ideas and even actual code implementing them, as well as relevant discussion.

Thanks so much for your interest and contributions, now is the time to give my own original solution to the challenge, but first I'll present some theory which will help understand the algorithm I used.

The challenge asked for you to write a program which should find and output all the constants in the range >1 and <2 with the property that their increasing integer powers would be nearer and nearer to integer values, outputting both each constant and its minimal polynomial, limited to degree 8 or less and coefficients -1,0,+1, with the leading coefficient in particular being always +1. The Golden Ratio (Phi) was given as an example of such a constant.

J-F Garnier tried the natural approach of constructing all -1,0,+1 polinomials up to degree 8, finding a candidate root, then computing its increasing integer powers and examining the fractional parts to see if they were converging to either 0.0000... or 0.9999.....

While natural, this approach is doomed to fail because the exponentially growing powers quickly exceed the accuracy of any calculator and the fractional parts get inaccurate or even lost, thus the test eventually fails.

The proper way comes by the way of a little theory, which I'll succintly explain with just the minimum amount of detail to make it short. Given the integer coefficients a(n) of a monic polynomial (leading term 1), say:

        x^n + a(n-1)*x^(n-1) + ... + a(2)*x^2 + a(1)*x + a(0)

any symmetric function of the roots can be expressed in terms of the coefficients, where symmetric means that the value of the function is the same (invariant) for all permutations of the roots. For example, the function:

        Sum of the roots = x1 + x2 + x3+ ... + xN

is a symmetric function because its value is the same no matter how the roots are permuted, so it can be computed as an integer function of the coefficients, in this case:

        Sum of the roots = x1 + x2 + ... + xN = -a(n-1)

        (Example: Sum of the roots of [x^3-6*x^2+11*x-6 = 0] = -(-6) = +6. The roots are 1,2,3 and 1+2+3 = 6)

The same goes for the sum of the roots raised to any integer power, such as for example:

        Sum of the 10-th powers of the roots = x1^10 + x2^10 + ... + xN^10

which, being a symmetric funtion, can also be computed as an integer expression of the coefficients, the expression being more complicated but still delivering an integer result, thus we have:

        x1^10 + x2^10+ ... + xN^10 = integer value

no matter what the roots are and no matter if they are real, complex, or a mix. Now consider what happens when one of the roots is real with an absolute value > 1 while all other roots (real or complex) have an absolute value < 1, such as is the case for the Golden Ratio, namely:

        Min. Polynomial : x^2-x-1
        root x1 : +1.61803398875
        root x2 : -0.61803398875

        Sum of their Nth powers S(N) = x1^N + x2^N = 1.61803398875^N + (-0.61803398875)^N

If we the compute S(N) for N=10, 20, 30, .. on a 12-digit HP calc we get:

        1.61803398875^10 + (-0.61803398875)^10 = 122.991869381 + 0.0081306187558 = 123 exactly
        1.61803398875^20 + (-0.61803398875)^20 = 15126.9999339 + 0.0000661069613 = 15127 exactly
        1.61803398875^30 + (-0.61803398875)^30 = 1860498.00000 + 0.0000005374904 = 1860498 exactly

and you see what's happening: the powers of the root with absolute value > 1 grow larger and larger while the powers of the root with absolute value < 1 grow smaller and smaller. As their sum has to be an integer value and the smaller root's powers are contributing less and less to the sum, the powers of the larger root are forced to approach the integer value of the whole sum by themselves and thus be almost-integer.

This generalizes to all monic polynomials of any degree so that the criterium is that the monic polynomial of degree N must have one real root of absolute value > 1 and all its other N-1 roots, whether real or complex, must have absolute values <1. This way the sums of the powers of the smaller roots will tend to 0 and the powers of the one root with absolute value > 1 will tend to the integer sum and so will tend to being integer themselves.

The algorithm is now quite clear:

        - iterate through all monic polynomials of coefficients 1,0,+1 and degrees up to 8 (per the challenge)
        - for each polynomial find all its roots and check:
                - that there's only one root of absolute value >1
                - that all other roots have absolute values < 1
        - if such a root is found, check if its absolute value is < 2 (per the challenge requirement)
                - if it is, record both the root and its minimal polynomial
        - after all polynomials have been examined, sort and ouput the recorded roots/polynomials.

and this is my 11-line generic implementation for the HP-71B (39 statements, 587 bytes), readily adaptable to any RPL and RPN HP models with a polynomial root-finding capability:

>CAT

SSMC21 BASIC 587 bytes

>LIST

         1  DESTROY ALL @ OPTION BASE 1 @ INPUT L @ DIM P(L+1),F(L),T$[60] @ COMPLEX R(L)
         2  INPUT "Sort by C,P ? ","C";X$ @ K=16*(X$="P") @ N=0 @ P(1)=1 @ H=0 @ SETTIME 0
         3  H=H+1 @ F(H)=-1 @ REPEAT @ P(H+1)=F(H) @ IF H<L THEN 3 ELSE MAT R=PROOT(P)
         4  X=0 @ FOR I=1 TO L @ IF ABS(R(I))<1 THEN 5 ELSE IF X THEN 6 ELSE X=REPT(R(I))
         5  NEXT I @ IF X>1 AND X<2 THEN N=N+1 @ DIM S$(N)[60] @ Z=L+2 @ GOSUB 9
         6  F(H)=F(H)+1 @ UNTIL F(H)>(H#1 AND H#L) @ H=H-1 @ IF H THEN 6
         7  FOR I=1 TO N @ FOR J=I+1 TO N @ IF S$(I)[K]>S$(J)[K] THEN VARSWAP S$(I),S$(J)
         8  NEXT J @ NEXT I @ FOR I=1 TO N @ DISP S$(I) @ NEXT I @ DISP N;TIME @ END
         9  Z=Z-1 @ IF NOT P(Z) THEN 9 ELSE FIX 11 @ X$=" "&STR$(X)&" " @ STD @ T$=""
        10  FOR I=1 TO L+1 @ IF P(I) THEN J=2-(P(I)>0) @ T$=T$&"+-"[J,J]&"x^"&STR$(Z-I)
        11  NEXT I @ S$(N)=REPLACE$(REPLACE$(X$&T$[2],"^1",""),"x^0","1") @ RETURN

(Lines 1-2 do the initialization, lines 3-6 do the search, lines 7-8 do the sorting and output, lines 9-11 format the output.)

It uses the Math ROM (PROOT) and the JPCROM for structured statements (REPEAT, UNTIL) and assorted utility functions (VARSWAP, REPLACE$). The JPCROM functions are mere conveniences and can be dispensed with but PROOT is mandatory or else the program would be much longer and slower.

My program, as written, is generic for monic -1,0,+1 polynomials of any degree, not just up to 8, and can display its output sorted by either constant or (for degree < 10) by minimal polynomial, as well as providing a count of the constants found and timing.

For example, for polynomials up to degree 4, sorting by constant:

>RUN

? 4
Sort by C,P ? C

1.32471795724 x^3-x-1
1.38027756910 x^4-x^3-1
1.46557123188 x^3-x^2-1
1.61803398875 x^2-x-1
1.83928675521 x^3-x^2-x-1
1.92756197548 x^4-x^3-x^2-x-1

6 .76

i.e., it found 6 constant in 0.76 seconds (Emu71 on an old 2.4 GHz single-core).
The same run sorted by polynomial (same count and timing):

>RUN

? 4
Sort by C,P ? P

1.61803398875 x^2-x-1
1.32471795724 x^3-x-1
1.46557123188 x^3-x^2-1
1.83928675521 x^3-x^2-x-1
1.38027756910 x^4-x^3-1
1.92756197548 x^4-x^3-x^2-x-1


For polynomials up to degree 6, sorted by constant:

>RUN

? 6
Sort by C,P ? C

1.32471795724 x^3-x-1
1.38027756910 x^4-x^3-1
1.44326879127 x^5-x^4-x^3+x^2-1
1.46557123188 x^3-x^2-1
1.50159480354 x^6-x^5-x^4+x^2-1
1.53415774491 x^5-x^3-x^2-x-1
1.57014731220 x^5-x^4-x^2-1
1.61803398875 x^2-x-1
1.66040772247 x^6-x^5-x^3-x^2-1
1.70490277604 x^5-x^4-x^3-1
1.74370016590 x^6-x^5-x^4-x-1
1.77847961614 x^5-x^4-x^3-x^2+1
1.81240361927 x^5-x^4-x^3-x-1
1.83928675521 x^3-x^2-x-1
1.88851884548 x^5-x^4-x^3-x^2-1
1.91118343669 x^6-x^5-x^4-x^3-x-1
1.92756197548 x^4-x^3-x^2-x-1
1.96594823665 x^5-x^4-x^3-x^2-x-1
1.98358284342 x^6-x^5-x^4-x^3-x^2-x-1

i.e.: 19 constants found in 15.25 seconds.

Finally, for polynomials up to degree 8 sorted by constant, as per the challenge requirements:

>RUN

? 8
Sort by C,P ? C

1.32471795724 x^3-x-1
1.38027756910 x^4-x^3-1
1.44326879127 x^5-x^4-x^3+x^2-1
1.46557123188 x^3-x^2-1
1.50159480354 x^6-x^5-x^4+x^2-1
1.53415774491 x^5-x^3-x^2-x-1
1.54521564973 x^7-x^6-x^5+x^2-1
1.57014731220 x^5-x^4-x^2-1
1.57367896839 x^8-x^7-x^6+x^2-1
1.59000537390 x^7-x^5-x^4-x^3-x^2-x-1
1.60134733379 x^7-x^6-x^4-x^2-1
1.61803398875 x^2-x-1
1.65363458677 x^7-x^6-x^5-1
1.66040772247 x^6-x^5-x^3-x^2-1
1.67752174784 x^7-x^6-x^5-x^2+1
1.70490277604 x^5-x^4-x^3-1
1.71428532914 x^6-x^5-x^4-x^2+1
1.72778030821 x^8-x^7-x^5-x^4-x^3-x^2-1
1.74370016590 x^6-x^5-x^4-x-1
1.74745742449 x^7-x^6-x^5-x^4+x^3-1
1.77452005864 x^7-x^6-x^5-x^3-1
1.77847961614 x^5-x^4-x^3-x^2+1
1.80750202302 x^6-x^5-x^4-x^3+1
1.81240361927 x^5-x^4-x^3-x-1
1.83524591514 x^8-x^7-x^6-x^4-x^3-x-1
1.83928675521 x^3-x^2-x-1
1.84771602509 x^8-x^7-x^6-x^5-1
1.85454747658 x^7-x^6-x^5-x^4-1
1.86221396917 x^8-x^7-x^6-x^5-x-1
1.88004410997 x^7-x^6-x^5-x^4-x-1
1.88670294847 x^8-x^7-x^6-x^5-x^2-x-1
1.88851884548 x^5-x^4-x^3-x^2-1
1.91118343669 x^6-x^5-x^4-x^3-x-1
1.92212800436 x^7-x^6-x^5-x^4-x^2-x-1
1.92756197548 x^4-x^3-x^2-x-1
1.96594823665 x^5-x^4-x^3-x^2-x-1
1.97061782036 x^8-x^7-x^6-x^5-x^4-x^3-1
1.97504243425 x^7-x^6-x^5-x^4-x^3-x^2-1
1.98358284342 x^6-x^5-x^4-x^3-x^2-x-1
1.99196419661 x^7-x^6-x^5-x^4-x^3-x^2-x-1
1.99603117974 x^8-x^7-x^6-x^5-x^4-x^3-x^2-x-1


i.e.: it found 41 constants in 239.96 seconds, just shy of 4 min. For fun, a run for polynomials up to degree 10 finds 85 constants in less than an hour.

That's all. Thanks again for your contributions, I hoped you enjoyed the challenge, see you all in S&SMC#22.

Best regards.
v.