[VA] SRC#003- New Year 2019 Special

[Valentin Albillo]

      
Hi all:

First of all, many thanks to all of you who contributed your various RPN/RPL solutions and valuable comments. As Thomas Okken suspected and Thomas Klemm explained, the reason this procedure works and converges to the cubic root of 2019 has all to do with the eigenvalues of the matrix M.

My original code for the HP-71B (easiest to understand), which exactly follows the steps given in my OP (i.e.: computing the powers of M instead of repeatedly multiplying by a vector) is the following 5-line, 144-byte snippet of code, which also keeps count of the number of iterations needed:

      1   DESTROY ALL @ OPTION BASE 1 @ DIM M(3,3),B(3,3)
      2   DATA PI,2019,2019,1,PI,2019,1,1,PI
      3   READ M @ MAT B=M @ R=0 @ I=O
      4   REPEAT @ MAT B=B*M @ L=R @ R=B(1,3)/B(2,3)
      5   I=I+1 @ UNTIL L=R @ DISP I;R;R^3

      >RUN

            183      12.6389823194      2019

so after 183 iterations the limit is found to be 12.6389823194, which is 2019^(1/3), the cubic root of 2019. Now for a few comments:

The procedure can be generalized in many ways. For instance:

1) My example used Pi in the main diagonal just for aesthetics but actually the procedure will converge for other positive values K in the main diagonal, resulting always in the same limit but greatly affecting the number if iterations needed for convergence. For instance:

      K  Iterations
      1      198
      2      200
      Pi     183
      10     130
      20      83
      40      51
      80      39
      120     34
      160     33
      200     37

as you may see in the table above, the value of K which results in the lowest number of iterations needed seems to be around 120-140. In fact, the theoretically optimum value for K in the main diagonal which results in the minimum number of iterations to converge for a given number N (2019 in my OP) is:

      K = (N^(1/3)+N) / (N^(1/3)+1)

which for N=2019 would be

      (2019^(1/3)+2019)/(2019^(1/3)+1) = 148.958253244

thus placing, say, 148 or 149 in the main diagonal instead of Pi will result in the lowest number of iterations needed, about 33-34 instead of the 183 needed when K=Pi.


2) The procedure will converge for numbers N other than 2019, be they integer, real (or even complex !), and the limit will be N^(1/3), the cubic root of N. For instance:

     - using  K = 1, N = 5:

            1      5      5
      M =   1      1      5
            1      1      1

will converge in 25 iterations to 1.70997594668 = 5^(1/3), the cube root of 5

      - using K = 2, N = 2:

            2      2      2
      M =   1      2      2
            1      1      2

will converge in 14 iterations to 1.25992104989 = 2^(1/3), the cube root of 2.

      - using K = 1, N = 1 + 2 i: (remember to define M, B, L and R as COMPLEX)

            1    (1,2)  (1,2)
      M =   1      1    (1,2)
            1      1      1

will converge in 23 iterations to 1.21961650797 + 0.471711267789 i = (1 + 2 i)^(1/3), the cube root of 1 + 2 i.


3) Though my OP specified the ratio B(1,3)/B(2,3), the procedure will converge using many other different ratios. For instance;

      B(1,1)/B(2,1), B(2,1)/B(3,1), B(1,2)/B(2,2), B(2,2)/B(3,2), ..., B(3,3)/B(3,2), B(3,2)/B(3,1)

will converge to N^(1/3), the cubic root of N, while the ratios;

      B(1,1)/B(3,1), ... , B(3,3)/B(3,1)

will converge to N^(2/3), the square of the cubic root of N (or the cubic root of the square of N, your choice)


4) My OP used a 3x3 matrix and the limit was the cube root of N, but using DxD matrices with the same pattern will make the various ratios converge to N^(1/D), N^(2/D), N^(3/D), ..., i.e., the Dth root of N and its powers.


5) Besides computing cubic roots, the procedure can also be made to converge to the root of a given cubic equation of a particular form. For instance:

      Find a root of:    x^3 - 1.2*x - 2.1 = 0

We'll first create the following 3x3 initial matrix M:

            K      P      Q
      M =   1      K      0
            0      1      K

where we'll use K=1, P=1.2, Q=2.1, so the initial matrix will be:

            1      1.2    2.1
      M =   1      1      0
            0      1      1

which converges in 24 iterations to 1.58816816249, which is the real root of the given cubic equation. This can be checked like this:

      >FNROOT(1,2,FVAR^3-1.2*FVAR-2.1)

            1.5881681625

6) As was the case for Dth roots, it's possible to create an initial DxD matrix to have the limit converge to a root of a given polynomial of degree D. The details are somewhat lengthy and will not be discussed here.


As a final observation, notice the fact that for integer N and K this procedure will produce rational approximations (integer numerator and denominator) that will converge arbitrarily close to the floating point real value of the Dth root using just integer additions and multiplications and nothing else: no subtractions, no divisions except the optional one at the end if you wish to convert the fraction to a floating point value. Same for polynomial equations with integer coefficients.

Thanks again to all of you for your interest and valuable contributions, much appreciated, keep them coming !

Regards.
V.