Matrix for the HP-48S/SX/G/GX

Overview

 1°) Magic Squares & Cubes

  a)  Magic Squares
  b)  Magic Square ?
  c)  Panmagic Squares
  d)  Panmagic Square ?
  e)  Bimagic Squares
  f)  Magic Multiplication
  g) Magic & Panmagic Cubes
  h) Magic Cube ?
  i) Primality ?

 2°)  Knight's Tour
 3°)  Hartley Transform

   a)  Discrete Hartley Transform ( Dim1 )
   b)  Discrete Hartley Transform ( Dim2 )
   c)  Discrete Hartley Transform ( Dim3 )
   d)  Hartley Transform with Filon Integration
   e)  Continous Data - Integration

 4°)  Eigenvalues & Eigenvectors - Jacobi's Method
 5°)  Linear Systems

   a)  Trangularization ( Householder Method )
   b)  Overdetermined Systems
   c)  Underdetermined Systems

 6°)  Linear Programming - Simplex Method
 7°)  Creating an Array
 8°)  Random Matrix
 9°)  Trace of a Square Matrix
10°) Rank of a Matrix
11°) Pseudo-Inverse ( Greville's Method )
12°) Power of a Square Matrix
13°) Matrix Square-Root
14°) Matrix N-th Root
15°) Exponential of a Square Matrix
16°) Logarithm of a Square Matrix
17°)  GETC
18°)  eps

1°)  Magic Squares & Cubes
 

     a) Magic Squares
 

-The elements of a normal magic square are  1 , 2 , ............. , n²  and it has the following property:
-All its row sums, column sums and both diagonal sums are equal to the magic constant   M = n.(n² + 1)/2

 
  • With MAG, the element a_i,j of a magic square of odd order n  is calculated by  a_i,j = 1 + ( i + 2 j - 2 ) mod n + n [ ( i + j - 1 + Int(n/2) ) mod n ]

    via the following method:

-Let  x' = 0  if  1 <= x <= n/2
 and  x' = 1  if  n/2 < x <= n

-Let  x" = n + 1 - x

  Let   t = n/2

    u = ( i' - j' ) mod t                  where           x' = min { x , n+1-x }
    v = 2 i" + j"                              and            x" = 0  if  0 < x <= n/2  &  x" = 1  if  n/2 < x <= n

  Then,  d(u,v) is defined by the table ( in the last 2 rows, x > 0 so that u > 2 ):
 
 

	d(u,0)	d(u,1)	d(u,2)	d(u,3)
d(0,v)	3	1	2	0
d(1,v)	1	3	0	2
d(2,v)	0	1	3	2
d(2x+1,v)	0	1	2	3
d(2x+2,v)	3	2	1	0

 
   Finally,  a_i,j =  t² d(u,v)  +  t  [  ( i' - j' + (t-1)/2 )  mod t ]  +  [  ( i' + j' - (t+1)/2 )  mod t ] + 1
 

       8   1   6
       3   5   7
       4   9   2

-Magic Constant = 15

  •  The element a_i,j of a panmagic square of order n where n is neither a mutiple of 2 nor a multiple of 3 is calculated by

    [ ( i - 1 ) + ( n/2 ) ( j - 1 ) ] mod n = x       if   x >= n/2   x' = 3n/2 - x - 1  &  if  x < n/2 ,  x' = x
    [ ( n/2 ) ( i - 1 ) + ( i - 1 ) ] mod n = y       if   y >= n/2   y' = 3n/2 - y - 1  &  if  y < n/2 ,  y' = y

   5   PMAG  gives:

   01   14   22   10   18
   07   20   03   11   24
   13   21   09   17   05
   19   02   15   23   06
   25   08   16   04   12

-Magic Constant = 65

   •   n = 4

   4   PMAG  returns:

   01   14   04   15
   08   11   05   10
   13   02   16   03
   12   07   09   06

Note:

  [[ 01   14   04   15 ]
   [ 08   11   05   10 ]
   [ 13   02   16   03 ]
   [ 12   07   09   06 ]     PMG?    ->     34

 BMAG  constructs a bimagic square ( a magic square that remains magic after squaring all its elements ) of order p² provided  p is a prime  p > 2

  If p > 3 , the bimagic square is also panmagic.

-The method uses 4 key-numbers:  r = ( a b c d )_p  r' = ( a' b' c' d' )_p  r" = ( a" b" c" d" )_p  r''' = ( a''' b''' c''' d''' )_p  witten in base p
-The element a_i,j  is given by

   a_i,j =  p³ [  { a [ ( i - 1 ) mod p ]  + a' int [ ( i - 1 ) / p ] + a" [ ( j - 1 ) mod p ] + a''' int [ ( j - 1 ) / p ] } mod p ]
        +  p² [  { b [ ( i - 1 ) mod p ]  + b' int [ ( i - 1 ) / p ] + b" [ ( j - 1 ) mod p ] + b''' int [ ( j - 1 ) / p ] } mod p ]
        +  p  [  { c [ ( i - 1 ) mod p ]  + c' int [ ( i - 1 ) / p ] + c" [ ( j - 1 ) mod p ] + c''' int [ ( j - 1 ) / p ] } mod p ]
        +     [  { d [ ( i - 1 ) mod p ]   + d' int [ ( i - 1 ) / p ] + d" [ ( j - 1 ) mod p ] + d''' int [ ( j - 1 ) / p ] } mod p ] + 1
 

-This formula will give a bimagic square if the following determinant is different from 0  ( mod p )

      a   b   c   d
      a'  b'  c'  d'
     a"  b"  c"  d"
     a''' b''' c''' d'''

 and if   a b' - a' b # 0 ,  a c' - a' c # 0 , a d' - a' d # 0 , b c' - b' c # 0 , ...... , c d' - c' d # 0    ( mod p )
 and similar relations with the keys  r" & r'''

-The square will be panmagic if the keys

       r + r" &  r' + r'''  and
       r - r"  &  r' - r'''  satisfy the same similar relations too.

-In BMAG,  the 4 keys are:

  ( 2 1 0 1 )   ( 0 1 1 2 )  ( 1 0 2 1 )  ( 2 1 2 0 )    if p = 3
  ( 0 1 1 1 )   ( 1 0 1 2 )  ( 3 3 2 1 )  ( 2 1 0 1 )    if p > 3 , p # 11
  ( 0 1 1 1 )   ( 1 0 1 2 )  ( 3 3 2 0 )  ( 2 1 0 1 )    if p = 11

  where  p is a prime,  p > 2

Example:     p = 3    We'll have a bimagic square of order 9

   3   BMAG   returns:

   01   35   60   70   14   39   49   74   27
   65   18   40   53   78   19   05   30   61
   48   79   23   09   31   56   69   10   44               Magic Constant = 369
   15   37   71   75   25   50   36   58   02
   76   20   54   28   62   06   16   41   66               Magic Constant of the squared elements = 20049
   32   57   07   11   45   67   80   24   46
   26   51   73   59   03   34   38   72   13
   63   04   29   42   64   17   21   52   77
   43   68   12   22   47   81   55   08   33

 MxM   takes 2 magic squares and calculates a "product" of these matrices that is also a magic square.
 

  Formula:     A is an mxm matrix   B is an nxn matrix  then  C = A*B = [ C_i,j ]  is defined as

   C_i,j = n² A + b_i,j 1_mxm          where       1_mxm  is the mxm matrix whose all elements equal 1

   01   14   22   10   18
   07   20   03   11   24
   13   21   09   17   05   
   19   02   15   23   06
   25   08   16   04   12

-The "product" C = AxB is the following 20x20 matrix:
 
 

26	351	101	376	39	364	114	389	47	372	122	397	35	360	110	385	43	368	118	393
201	276	126	251	214	289	139	264	222	297	147	272	210	285	135	260	218	293	143	268
326	51	401	76	339	64	414	89	347	72	422	97	335	60	410	85	343	68	418	93
301	176	226	151	314	189	239	164	322	197	247	172	310	185	235	160	318	193	243	168
32	357	107	382	45	370	120	395	28	353	103	378	36	361	111	386	49	374	124	399
207	282	132	257	220	295	145	270	203	278	128	253	211	286	136	261	224	299	149	274
332	57	407	82	345	70	420	95	328	53	403	78	336	61	411	86	349	74	424	99
307	182	232	157	320	195	245	170	303	178	228	153	311	186	236	161	324	199	249	174
38	363	113	388	46	371	121	396	34	359	109	384	42	367	117	392	30	355	105	380
213	288	138	263	221	296	146	271	209	284	134	259	217	292	142	267	205	280	130	255
338	63	413	88	346	71	421	96	334	59	409	84	342	67	417	92	330	55	405	80
313	188	238	163	321	196	246	171	309	184	234	159	317	192	242	167	305	180	230	155
44	369	119	394	27	352	102	377	40	365	115	390	48	373	123	398	31	356	106	381
219	294	144	269	202	277	127	252	215	290	140	265	223	298	148	273	206	281	131	256
344	69	419	94	327	52	402	77	340	65	415	90	348	73	423	98	331	56	406	81
319	194	244	169	302	177	227	152	315	190	240	165	323	198	248	173	306	181	231	156
50	375	125	400	33	358	108	383	41	366	116	391	29	354	104	379	37	362	112	387
225	300	150	275	208	283	133	258	216	291	141	266	204	279	129	254	212	287	137	262
350	75	425	100	333	58	408	83	341	66	416	91	329	54	404	79	337	62	412	87
325	200	250	175	308	183	233	158	316	191	241	166	304	179	229	154	312	187	237	162

 
-Here, the magic constant = 4510
-But if you subtract 25 to all these elements, you get a normal magic square, with all the integers between 1 and 400, magic constant = 4010

-Since A & B are panmagic squares,  C is a panmagic square too !

-An Andrews magic cube is a cubic array where the sums of the elements equal the magic constant n(n³+1)/2  in the 3 directions and the 4 space diagonals.
-The cube is perfect pandiagonal if all the orthogonal sections are panmagic squares.
-In a panmagic cube , all the orthogonal or diagonal sections - broken or unbroken - are panmagic squares.
 

 "MGCB" builds:

    Andrews magic cubes if n = 3 or 5 and - more generally - if n is odd.
    Perfect pandiagonal magic cube if n = 7
    Panmagic cubes if n is a prime > 7 

Formulae:

    a_i,j,k =  n²  { [ a ( i - 1 ) + a' ( j - 1 ) + a" ( k - 1 ) ] mod n }
              +  n  { [ b ( i - 1 ) + b' ( j - 1 ) + b" ( k - 1 ) ] mod n }
               +  [  c ( i - 1 ) + c' ( j - 1 ) + c" ( k - 1 ) ] mod n + 1

    with

   ( a b c )   =  ( 1 -1  1 )
  ( a' b' c' )  =  ( -1  1  1 )      to create Andrews Magic Cubes of odd orders
 ( a" b" c" ) =  ( 1  1  -1 )

   ( a b c )   =  ( 6  5  4 )
  ( a' b' c' )  =  ( 5  4  5 )       if  n = 7
 ( a" b" c" ) =  ( 4  6  6 )

-Let  x' = 0  if  1 <= x <= n/2
 and  x' = 1  if  n/2 < x <= n

-Let  x" = n + 1 - x

    01  17  24
    23  03  16       15  19  08
    18  22  02       07  14  21      26  06  10
                           20  09  13      12  25  05
                                                 04  11  27

-The magic constant = 42 =  1+17+24 = 1+23+18 = .....  = 4+11+27 = 10+5+27

  = 1+14+27 = 18+14+10 = 24+14+4 =2+14+10  ( the 4 space diagonals )

  where   { i  j }  is the starting square                                M = m x n Matrix = a solution
   and    { m  n } is the dimension of the chessboard           N = number of visited squares

-The program fails if  N <  n m  ( in this case, M contains one or several zero(s) )

Example1:          m = 3 ,  n = 4  ,  starting square = { 1  1 }

   { 1  1 }   ENTER
   { 3  4 }   KNIGHT  returns  N = 12         ( in 41 seconds with an HP-48S/SX , in 27 seconds with an HP-48G/GX )

                   [ [  1   4   7  10 ]
   and   M =   [  8  11  2   5  ]
                     [  3   6   9  12 ] ]

Example2:          m = n = 8  - standard chessboard -  starting square = { 1  1 }

   { 1  1 }   ENTER
   { 3  4 }   KNIGHT  returns  N = 64         ( in 4mn57s with an HP-48S/SX , in 3mn17s with an HP-48G/GX )

                    [ [  1    26   15   24   29   34   13   32 ]
                      [ 16   23   28   37   14   31   64   35 ]
                      [ 27    2    25   30   61   36   33   12 ]
   and   M  =   [ 22   17   40   49   38   55   60   63 ]
                      [  3    48   21   56   41   62   11   54 ]
                      [ 18   45   42   39   50   53    8    59 ]
                      [ 43    4    47   20   57    6    51   10 ]
                      [ 46   19   44    5    52    9    58    7  ] ]

-If m = n = 8 , all the starting squares produce a complete tour, except  { i  j } = { 6  4 }
 

               [ [  1  3  4  10  ]
 Let  V =   [  4  5  7  14  ]       DHT2   gives 
                 [  2  9  6  11  ] ]

                       [ [    76             -28         -28          8         ]
     DHT (A) =   [  -9.2679    5.9282    5.4641   -3.1962  ]        ( rounded to 4 decimals )
                         [ -12.7321  -7.9282   -1.4641    7.1962  ] ]

Notes:

-In this case, the DHT of a  mxnxp  hypermatrix [ a_i,j,k ]  ( 1 <= i <= m ; 1 <= j <= n ; 1 <= k <= p ) is a  mxnxp hypermatrix  [ b_i,j,k ]

  where  b_i,j,k =  SUM_{q=1,2,...,m  ; r=1,2,....,n ; s=1,2,....,p}   a_q,r,s. cas 360°[(q-1)(i-1)/m+(r-1)(j-1)/n+(s-1)(k-1)/p]         ( cas = sin + cos )

    Let L =  { [ [   1  25  40   5    2  ]                  
                       [   4  36  18  32  37 ]                 
                       [   9   8   39  20  33 ]
                       [ 16  23  21  10  31 ] ]               
 
                             [ [  31  10  21  23  16 ]
                               [  33  20  39   8    9  ]
                               [  37  32  18  36   4  ]
                               [   2    5   40  25   1  ] ] 

                                       [ [   0  16  23  21  10 ]
                                         [   1  25  40   5    2  ]
                                         [   4  36  18  32  37 ]
                                         [   9   8   39  20  33 ] ] }        DHT3     returns a list of three  4x5-matrices

-For instance, the (2;1)-element of the 3rd matrix became:  b_2,1,3 = 67.0070415521

Notes:

-Actually, Hartley has defined the real transform  H(µ) =  (2.pi) ^-1/2  §_-infinity^+infinity  f(t) ( cos µ.t + sin µ.t ) dt      ( § = integral )
-This transform is strictly reciprocal.

     We take  f(t) =  e ^-t/2  for  t >= 0  and  f(t) = 0  if  t < 0

- H(µ) may be expressed analytically and  H(µ) = 2(1+2µ)/(1+4µ²) / sqrt ( 2.Pi )

-We can choose a & b  so that f(x) is enough small to be neglected in the integral

-If we use f(0) , f(1) , .............. , f(16) to represent the function,

                     0     ENTER
                    16    ENTER
[ 1  0.60653...  0.36787....    ........................   0.000335... ]
                     2      HFIL                                                             ->      0.2366

 we obtain the following results ( rounded to 4 decimals )

 
 

x	-1	-0.5	-0.25	0	0.2071	1	2
HFIL	-0.1584	0.0013	0.3197	0.7979	0.9636	0.4797	0.2366
exact	-0.1596	0	0.3192	0.7979	0.9631	0.4787	0.2347

-The graph of the Hartley transform looks like this:

                                                              H(x)
                                                                 |
                                                                 |  *                               H(x) is maximum for x = 0.2071
                                                                 |*   *
                                                                 |          *
                                                                 |                       *
                                                                 |                                                          *
  ---------------------------------------* -|----------------------------------------------- x
            *                                                   |
                                     *                 *       |
                                                                 |
 

Example:    The same one:    f(t) =  e ^-t/2  for  t >= 0  and  f(t) = 0  if  t < 0

                    4  FIX

                    0             ENTER
                   16            ENTER
              'EXP(-X/2)'  ENTER
                   2               HRTL     ->      0.2346

 
 we obtain the following results ( rounded to 4 decimals )

 
 

x	-1	-0.5	-0.25	0	0.2071	1	2
HFIL	-0.1597	0.0003	0.3190	0.7976	0.9635	0.4789	0.2346
exact	-0.1596	0	0.3192	0.7979	0.9631	0.4787	0.2347

-The precision is better in 5  FIX  ,  6 FIX  ....   but the execution will increase a lot.
-You can use  a = - E499 & b = + E499,  the change of variable  t = Tan u  will change them into ( about )  - Pi/2 & + Pi/2

 JCB uses a variant of the Jacobi's algorithm:

*** The eigenvalues of a matrix A are the diagonal-elements of an upper triangular matrix T equal to the infinite product   ...U_k^-1.U_k-1^-1....U₂^-1.U₁^-1.A.U₁.U₂....U_k-1.U_k....
 where the U_k are unitary matrices. The eigenvectors are the columns of  U = U₁.U₂....U_k-1.U_k.... if A is Hermitian ( i-e if A equals its transconjugate )
 Actually, T is diagonal if A is Hermitian.

-"JCB" finds the greatest element a_{i j} below the main diagonal
-Then,  U₁  is determined so that it places a zero in position ( i , j )  in U₁^-1.A.U₁

  U₁  has the same elements as the Identity matrix, except that  u_{i i} = u_{j j} = x  and  u_{i j} = y + i.z ,  u_{j i} = -y + i.z

     with  y + i.z  = C.x  ,  x = ( 1 + | C |² ) ^-1/2  and   C = 2.a_{i j} / [ a_{i i} - a_{i j} + ( ( a_{i i} - a_{i j} )² + 4 a_{i j} a_{j i} )^1/2 ]

-The process is repeated until the greatest sub-diagonal element is smaller than eps

-The successive greatest | a_{i j} |  ( i > j ) are displayed when the routine is running.

*** If the matrix is non-Hermitian, and if all the eigenvalues are distinct,  we define an upper triangular matrix W = [ w_{i j} ]  by

   w_{i j} = 0    if  i > j
   w_{i i} = 1
   w_{i j} = Sum_{k = i+1, i+2 , .... , j}      t_{i k} w_{k j} / ( t_{j j} - t_{i i} )      if  i < j

            1  2  4
    A =  2  7  3     JCB   returns    the 3-vector in level 1    [ -0.730676198691    4.91074121335   12.8199349853  ]
            4  3  9

-So, the eigenvalues are:

                 k₁ = -0.730676198691
                 k₂ =  4.91074121335
                 k₃ = 12.8199349853

and the 3 corresponding eigenvectors are the 3 columns of the matrix in level 2:

Example2:

             1  2  4
     A =  4  3  5    JCB   gives the following results;
             7  4  7

         k₁ =  -2.09209759346
         k₂ =  0.185167648589
         k₃ = 12.90692994 48

and the 3 corresponding eigenvectors:

          V₁ ( 0.800454174265 ; -0.0416510801675 ; -0.597945066394 )
          V₂ ( -0.0958180743763 ; 0.907396310041 ; -0.434179238312 ) 
          V₃ ( 0.360058723032 ; 0.548018519593 ;  0.797789237995 ) 

-The eigenvectors are not unit vectors except the first one if the A is not hermitian ( or not symmetric if all the coefficients are real ).
-Divide V₂ & V₃ by  || V₂ || & || V₃ || respectively if you want unit eigenvectors.

Notes:
 
-If A is non-symmetric ( or not hermitian for complex matrices ), this program only works if all the eigenvalues are distinct:
-The eigenvalues are correctly computed, but not all the eigenvectors.
-The eigenvalues of a symmetric real matrix are always real.
-The eigenvalues of a hermitian matrix are always real too.

-The Householder method uses orthogonal symmetries to gradually triangularize a matrix.

-If ( a₁ , ...... , a_n ) is the first column of the matrix A,  the other columns  ( x₁ , ...... , x_n )  are replaced by  ( y₁ , ...... , y_n )  defined as

   ( y₁ , ...... , y_n ) = ( x₁ , ...... , x_n ) - C [ a₁ + sgn(a₁) ( a₁² + ....... + a_n² )^1/2 , a₂ , ...... , a_n ]

 where   C = [ SUM_i=1,...,n a_i x_i + x₁ sgn(a₁) ( SUM_i=1,...,n a_i² )^1/2 ] / [ SUM_i=1,...,n a_i² + a₁ sgn(a₁) ( SUM_i=1,...,n a_i² )^1/2 ]

-The first column becomes  ( A₁ , 0 , .......... , 0 )   with   A₁ = -sgn(a₁) ( SUM_i=1,...,n a_i² )^1/2
-Then, the same process is repeated with the 2nd column, without taking into account the 1st column and the 1st raw, and so on...

-An overdetermined system  A.X = B  has often no solution at all
-But the following programs calculate the least-squares solution:
-It minimizes || A.X - B ||

     5x + y + z  =   8          
     4x - y + 2z =  13
       x + 2y - z =  -5
     7x - 4y + 5z =  32
     2x + 5y - 9z = -20

-The purpose is to find  n  non-negative real numbers:   x₁ , ..... , x_n   satisfying:

                          b_i >= a_i;1x₁ + ....... + a_i;n x_n            i = 1 , ....... , m                              ( m inequations )           the b_i  , b_i'  may be positive or not !
                          b_i'  =  a_i';1x₁ + ....... +a_i';n x_n            i' = 1 , ...... , m'                              ( m' equations )

and such that      F  =  c₁x₁ + ....... +  c_nx_n       is maximized.    ( to find the minimum of F, seek the maximum of  -F )
 

-Here, we first treat each equation by pivoting on the largest coefficient to develop a complete basic variable set

-Then, the dual simplex pivoting rule is used to transform all the negative b_i into non-negative numbers ( if possible )

       •  pivot = a_i,j < 0  with the largest  c_j / a_i,j

-Then, the primal simplex pivoting rule is applied to transform all the positive c_i into non-positive numbers

       •  pivot = a_i,j > 0  with the smallest  b_i / a_i,j

Find  4 non-negative numbers  x, y, z and t satisfying:

             56 >=  x + 2y + 3z +4t
             41 >= 5x +  y - 3z - t
             61 >= 4x + 7y +2z - 3t  
               7 <=   x + y + z + t
             17 <= 2x - y + z + 3t
             25 = 2x + 3y - z + 2t 

such that F = 2x + 4y +3z + 5t is maximized.

1-Re-write the LP as follows:  the >= inequations first ( if any ), then the = equations ( if any )

             56 >=  x + 2y + 3z +4t
             41 >= 5x +  y - 3z - t
             61 >= 4x + 7y +2z - 3t
             -7 >=  -x  -  y  -  z  -  t                          the "<=" inequations must be multiplied by (-1) to become  ">=" inequations
           -17 >= -2x + y  -  z - 3t
             25 =  2x + 3y - z + 2t
 

2-Then, place the following matrix in the stack

-We also get the 5 slack variables:  0  52  0  14  0

Notes:

Formula:

    B₀   =  ( A A^T )  /   Trace ( A A^T )
   B_k+1 = 2 B_k - B_k²

  A⁺ is the unique mxn matrix ( m rows, n columns ) such that:  A A⁺ A = A ,  A⁺ A A⁺  = A⁺ ,  A A⁺ and  A⁺ A  symmetric

-The method of Greville gradually builds the pseudo-inverse row by row:

-Let  a_k = the k-th column of A
        A_k = the matrix built with the k first columns of A
       A⁺_k  = the pseudoinverse at step k

-We start with  A⁺₁ = ^ta₁ / ( ^ta₁ a₁ )   if   a₁ # 0
                  or  A⁺₁ = ^ta₁                   if   a₁ = 0     and then,  for k = 2 , 3 , ..... , m

     Let    d_k = A⁺_k-1  a_k
              c_k = a_k - A_k-1 d_k
              b_k = ^tc_k / ( ^tc_k c_k )   if   c_k # 0      or      b_k = ^td_k A⁺_k-1 / ( 1 +  ^td_k d_k )   if  c_k = 0

   >>>   A⁺_k  = [ A⁺_k-1 - d_k b_k ]           In other words,  A⁺_k  is obtained by placing  the row   b_k  under the matrix  A⁺_k-1 - d_k b_k
                          [      b_k      ]

                 1  1  4  2          
      A =     0  1  2  3       GREV       gives
                 3  2  6  7           

                      -21/112   -85/112   43/112
      A⁺ =           7/112      23/112   -9/112             ( with some roundoff errors of course ).
                       49/112      1/112   -15/112
                      -35/112    29/112    13/112

-In this example,  A A⁺ = Id₃

          1  4  9    
  A =  3  5  7     ENTER    7     POWM    gives:
          2  1  8    

-Given a non-singular nxn matrix A  ( det A # 0 ), "SQRTM" uses the iterations:

    M₀ = A      M_k+1 = (1/2) [ I + ( M_k + M_k^-1 ) / 2 ]      where  I  is the Identity matrix
    X₀ = A      X_k+1 = X_k ( I + M_k^-1 ) / 2

-The principal n-th root of a non-singular matrix A ( det A # 0 )  may be computed by the algorithm:

    M₀ = A      M_k+1 = M_k [ ( 2.I + ( n - 2 ) M_k ) ( I + ( n - 1) M_k ) ^-1 ] ^p      where  I  is the Identity matrix
    X₀ =  I       X_k+1  = X_k  ( 2.I + ( n - 2 ) M_k ) ^-1 ( I + ( n - 1) M_k )

    M_k  tends to  I        as k tends to infinity
    X_k   tends to A ^1/p  as k tends to infinity

15°) Exponential of a Square Matrix

-The exponential of a nxn matrix A is defined as

   exp(A) = I + A + A²/2! + A³/3! + ..... + A^k/k! + ....                ( I = the Unit matrix = the Identity matrix:  1 in the diagonal, 0 elsewhere )

16°)  Logarithm of a Square Matrix

-If  || A - I || < 1 , the logarithm of a nxn matrix A is defined by the series expansion:

               1.2   0.1   0.3   
      A =   0.1   0.8   0.1       LNM     returns:
               0.1   0.2   0.9   

                   0.16708339584        0.0695779233559   0.287707998737  
    Ln(A) =  0.0977830050234   -0.240971673255     0.103424021357
                   0.0865009723563    0.23505312438      -0.131906635565

-In this example,   || A - I || = 0.5099...  < 1

Notes:

References: