The provided implementation in the post is flawed and fails to port the Pari/GP code. In this statement

Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)

the actual intentions of rem(1, x**r-1) and (1%n) are to create symbolic representations of
% (x**r - 1) and % n for each matrix element, which is going to take effect during the computation of matrix power.

However, under this implementation, rem(1, x**r-1) and (1%n) are both evaluated right away:

r = 3
n = 11
print(rem(1, x**r-1), (1 % n))

Output:

1 1

Therefore, these are not no longer symbolic representations and would be eliminated immediately after being multiplied to the matrix.

All in all I do not think the overall algorithm could be implemented easily using the sympy library. It would be more convenient to implement the matrix power-modulo operation yourself with the help of numpy to work with polynominal computations (see doc).

The provided implementation in the post is flawed and fails to port the Pari/GP code. In this statement

Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)

the actual intentions of rem(1, x**r-1) and (1%n) are to create symbolic representations of
% (x**r - 1) and % n for each matrix element, which is going to take effect during the computation of matrix power.

However, under this implementation, rem(1, x**r-1) and (1%n) are both evaluated right away:

r = 3
n = 11
print(rem(1, x**r-1), (1 % n))

Output:

1 1

Therefore, these are not no longer symbolic representations and would be eliminated immediately after being multiplied to the matrix.

All in all I do not think the overall algorithm could be implemented easily using the sympy library. It would be more convenient to implement the matrix power-modulo operation yourself with the help of numpy to work with polynominal computations (see doc). By the way, I am not sure whether that helps but numpy actually have a chebyshev module.

The provided implementation in the post is flawed and fails to port the Pari/GP code. In this statement

Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)

The actual intentions of rem(1, x**r-1) and (1%n) are to create symbolic representations of % (x**r - 1) and % n for each matrix element, which is going to take effect during the computation of matrix power.

However, under this implementation, rem(1, x**r-1) and (1%n) are both evaluated right away:

r = 3
n = 11
print(rem(1, x**r-1), (1 % n))

Output:

1 1

Therefore, these are not no longer symbolic representations and would be eliminated immediately after being multiplied to the matrix.

All in all I do not think the overall algorithm could be implemented easily using the sympy library. It would be more convenient to implement the matrix power-modulo operation yourself with the help of numpy to work with polynominal computations (see doc).

The provided implementation in the post is flawed and fails to port the Pari/GP code. In this statement

Matrix([[2*x, -1], [1, 0]])*rem(1,x**r-1)*(1%n)

the actual intentions of rem(1, x**r-1) and (1%n) are to create symbolic representations of 
% (x**r - 1) and % n for each matrix element, which is going to take effect during the computation of matrix power.

However, under this implementation, rem(1, x**r-1) and (1%n) are both evaluated right away:

r = 3
n = 11
print(rem(1, x**r-1), (1 % n))

Output:

1 1

Therefore, these are not no longer symbolic representations and would be eliminated immediately after being multiplied to the matrix.

All in all I do not think the overall algorithm could be implemented easily using the sympy library. It would be more convenient to implement the matrix power-modulo operation yourself with the help of numpy to work with polynominal computations (see doc).