# Challenge

Given a left- or right-stochastic matrix where the limit as x approaches infinity of the matrix to the power of x approaches a matrix with all finite values, return the matrix to which the matrix converges. Basically, you want to keep multiplying the matrix by itself until the result no longer changes.

# Test Cases

[[7/10, 4/10], [3/10, 6/10]] -> [[4/7, 4/7], [3/7, 3/7]]
[[2/5, 4/5], [3/5, 1/5]] -> [[4/7, 4/7], [3/7, 3/7]]
[[1/2, 1/2], [1/2, 1/2]] -> [[1/2, 1/2], [1/2, 1/2]]
[[1/3, 2/3], [2/3, 1/3]] -> [[1/2, 1/2], [1/2, 1/2]]
[[1/10, 2/10, 3/10], [4/10, 5/10, 6/10], [5/10, 3/10, 1/10]] -> [[27/130, 27/130, 27/130], [66/130, 66/130, 66/130], [37/130, 37/130, 37/130]]
[[1/7, 2/7, 4/7], [2/7, 4/7, 1/7], [4/7, 1/7, 2/7]] -> [[1/3, 1/3, 1/3], [1/3, 1/3, 1/3], [1/3, 1/3, 1/3]]


# Rules

• Standard Loopholes Apply
• You may choose whether you want a right- or a left-stochastic matrix
• You may use any reasonable number type, such as floats, rationals, infinite precision decimals, etc
• This is , so the shortest submission in bytes for each language is declared the winner for its language. No answer will be accepted
• @FryAmTheEggman it seems like some earlier comments have been deleted, so this might be redundant, but reducible and periodic matrices are already excluded by "Given a left- or right-stochastic matrix where the limit as x approaches infinity of the matrix to the power of x approaches a matrix with all finite values", which I read as saying the input is guaranteed to converge to a unique solution. (i.e. the input matrix is guaranteed to be ergodic.) Mar 8 '18 at 1:49
• @Nathaniel That isn't quite true, as if the chain is reducible you can have a result (like for the identity matrix) that meets what you said but the chain it describes is not irreducible and therefore the input wouldn't be guaranteed to be ergodic (since it won't be positive recurrent). Guaranteeing ergodicity is what the OP wants, and I think they have that now, thanks to the additional constraint of all of the row values being identical. If you know a better way to constrain the input without needing to add an explanation of Markov chains, I'm sure HyperNeutrino would appreciate it! :) Mar 8 '18 at 2:06
• @FryAmTheEggman ah, you're right, sorry. I was thinking of power iteration rather than raising the matrix to a power. (So by "unique solution" I meant "one that's independent of the starting point of the iteration process", but that's not relevant here.) I agree that the 'all rows are identical' condition does the job. I suppose the OP could just say "the Markov chain is guaranteed to be ergodic," which would satisfy people like us who are likely to worry about it! Mar 8 '18 at 2:56
• Actually, if B is a solution to BA = B, then so is cB for any scalar constant c. So a non-zero solution cannot be strictly unique, unless you somehow fix the scale. (Requiring B to be stochastic would do it.) Also, obviously, whether its the rows or the columns of B that are equal will depend on whether A is left or right stochastic. Mar 8 '18 at 3:07
• For anyone else who never learned about matrices during Math-class in high school and how to multiply them: mathsisfun.com/algebra/matrix-multiplying.html. I had to look it up to understand what was being asked.. Maybe it's common knowledge elsewhere on Earth.. :S Mar 8 '18 at 8:49

# R,  44  43 bytes

function(m){X=m
while(any(X-(X=X%*%m)))0
X}


Try it online!

Just keeps multiplying until it finds a fixed matrix. Apparently X!=(X=X%*%m) does the comparison, then reassigns X, so that's neat.

Thanks to @Vlo for shaving off a byte; even though crossed out 44 is still regular 44.

• I wonder why function(m){ while(any(m!=(m=m%*%m)))0 m} doesn't work. Numerical inaccuracies preventing the termination condition from triggering? Mar 8 '18 at 9:19
• @CodesInChaos most likely it's a lack of precision. Switching to an arbitrary precision library doesn't help, either -- they either time out (Rmpfr) or fail (gmp) in the same manner, although I'm probably doing something wrong. Mar 8 '18 at 16:22
• @Giuseppe I guess the suggested approach is repeated squaring until no longer changes? (I can't read R) Mar 8 '18 at 16:35
• @user202729 yeah it is. R uses 64-bit floating point numbers and I know errors propogate pretty quickly. Mar 8 '18 at 16:41
• I think that algorithm is unstable. Jelly has the same problem too. (TODO prove this and find an alternative) Mar 8 '18 at 16:45

# Octave, 4542 35 bytes

@(A)([v,~]=eigs(A,1))/sum(v)*any(A)


Try it online!

Saved 3 bytes thanks to Giuseppe, and 7 more thanks to Luis Mendo!

This uses that the eigenvector corresponding to the eigenvalue 1 (also the maximum eigenvalue) is the column vector that is repeated for each value of the limiting matrix. We have to normalise the vector to have sum 1 for it to be stochastic, then we simply repeat it to fill out the matrix. I'm not too well versed in golfing Octave, but I haven't been able to find a functional way to do repeated multiplication, and a full program seems like it will always be longer than this.

We can use any(A) since from the restrictions we know that the matrix must describe an irreducible Markov chain, and so each state must be reachable from the other states. Therefore at least one value in each column must be non-zero.

• How does eigs always return the eigenvector corresponding to 1? My memory of markov chains is a bit fuzzy. Mar 7 '18 at 19:24
• 42 bytes Mar 7 '18 at 19:32
• @Giuseppe Because the matrix is stochastic, and a couple other things, its maximum eigenvalue is 1, and eigs returns starting from the largest eigenvalue. Also, thanks for the golf! Mar 7 '18 at 19:39
• Ah, right. Markov chains are on my next exam but since it's for actuaries, some of the details are entirely missing. Mar 7 '18 at 19:47

# Wolfram Language (Mathematica), 14 bytes

Output floats:

N@#//.x_:>x.#&


Try it online!

# Wolfram Language (Mathematica), 30 bytes

Output fractions:

Limit[#~MatrixPower~n,n->∞]&


Try it online!

æ×³$ÐL  This is a full program. Try it online! # APL (Dyalog), 18 6 bytes 12 bytes saved thanks to @H.PWiz +.×⍨⍣≡  Try it online! • +.×⍨⍣≡ for 6 bytes. That is, square until nothing changes Mar 7 '18 at 18:57 • @H.PWiz I believe it is. I have no idea why I hadn't done it in the first place. Thanks! Mar 7 '18 at 19:02 # k/q, 10 bytes k/q because the program is identical in both languages. The code is really just idiomatic k/q. {$[x]/[x]}


### Explanation

{..} is out lambda syntax, with x as an implicit parameter

$[x] has$ as the binary matrix multiplication operator, providing only one parameters creates a unary operator that left multiplies by the Markov Matrix

/[x] applies the left multiplication until convergence, starting with x itself.

# C (gcc), 207192190181 176 bytes + 2 flag bytes -lm

• Saved fifteen seventeen twenty-two bytes thanks to ceilingcat.
• Saved nine bytes; removing return A;.
float*f(A,l,k,j)float*A;{for(float B[l*l],S,M=0,m=1;fabs(m-M)>1e-7;wmemcpy(A,B,l*l))for(M=m,m=0,k=l*l;k--;)for(S=j=0;j<l;)m=fmax(m,fdim(A[k],B[k]=S+=A[k/l*l+j]*A[k%l+j++*l]));}


Try it online!

• @ceilingcat Counting the compiler flag bytes, this results in 192 again. Included your suggestion nonetheless. Jun 20 '18 at 21:08
• @ceilingcat Thank you. Dec 18 '18 at 21:24

# Python 3, 75 61 bytes

f=lambda n:n if allclose(n@n,n)else f(n@n)
from numpy import*


Try it online!

In test cases there are float inaccuracies, so the values may differ by a bit from the exact fractions.

PS. numpy.allclose() is used because numpy.array_equal() is longer and prone to float inaccuracies.

-14 bytes Thanks HyperNeutrino ;) Oh yes I forgot the @ operator ;P

• Use dot instead of matmul :D Mar 8 '18 at 4:01
• Actually, take numpy arrays as input and do x=n@n :P tio.run/… Mar 8 '18 at 4:03
• 59 bytes Mar 8 '18 at 4:04
• Added back the f= at front because it is recursively called ;) Mar 8 '18 at 4:07
• Oh yeah you're right :) good call! Mar 8 '18 at 4:09

# Java 8, 356 339 bytes

import java.math.*;m->{BigDecimal t[][],q;RoundingMode x=null;for(int l=m.length,f=1,i,k;f>0;m=t.clone()){for(t=new BigDecimal[l][l],i=l*l;i-->0;)for(f=k=0;k<l;t[i/l][i%l]=(q!=null?q:q.ZERO).add(m[i/l][k].multiply(m[k++][i%l])))q=t[i/l][i%l];for(;++i<l*l;)f=t[i/l][i%l].setScale(9,x.UP).equals(m[i/l][i%l].setScale(9,x.UP))?f:1;}return m;}


-17 bytes thanks to @ceilingcat.

Definitely not the right language.. Damn floating point precision..

Explanation:

Try it online.

import java.math.*;     // Required import for BigDecimal and RoundingMode
m->{                    // Method with BigDecimal-matrix as both parameter and return-type
BigDecimal t[][],q;   //  Temp BigDecimal-matrix
RoundingMode x=null;  //  Static RoundingMode value to reduce bytes
for(int l=m.length,   //  The size of the input-matrix
f=1,          //  Flag-integer, starting at 1
i,k;          //  Index-integers
f>0;              //  Loop as long as the flag-integer is still 1
m=t.clone()){     //    After every iteration: replace matrix m with t
for(t=new BigDecimal[l][l],
//   Reset matrix t
i=l*l;i-->0;)   //   Inner loop over the rows and columns
for(f=k=0;        //    Set the flag-integer to 0
k<l           //    Inner loop to multiply
;             //      After every iteration:
//       Sum the value at the current cell in matrix t with:
m[i/l][k]  //        the same row, but column k of matrix m,
.multiply(m[k++][i%l])))
//        multiplied with the same column, but row k of matrix m
q=t[i/l][i%l];  //     Set temp q to the value of the current cell of t
for(;++i<l*l;)      //   Loop over the rows and columns again
f=t[i/l][i%l].setScale(9,x.UP).equals(m[i/l][i%l].setScale(9,x.UP))?
//    If any value in matrices t and m are the same:
f              //     Leave the flag-integer unchanged
:               //    Else (they aren't the same):
1;}            //     Set the flag-integer to 1 again
return m;}            //  Return the modified input-matrix m as our result

• Why is the main function so verbose? Mar 8 '18 at 16:31
• @user202729 Because float/double don't have the proper floating point precision, java.math.BigDecimal should be used instead. And instead of simply +-*/, BigDecimals use .add(...), .subtract(...), .multiply(...), .divide(...). So something as simply as 7/10 becomes new BigDecimal(7).divide(new BigDecimal(10)). In addition, the ,scale,RoundingMode in the divide are required for values with 'infinite' decimal values (like 1/3 being 0.333...). The main method can of course be golfed, but I didn't bother when I did a search&replace to convert the floats to BigDecimals. Mar 8 '18 at 17:45
• @ceilingcat Thanks! Jan 4 '20 at 13:07