Two solutions.
One is not to use numpy, since numpy numbers are limited. So just to it with pure python's int (that are big int)
# That's your matrix. But I had to add all the comas myself. That is painful.
# Please, next time, post a minimal reproducible example. You get less help
# when you make the life of people who might help you more difficult
B=[[ 96, 205, 325, 460, 599], [ 13, 109, 207, 307, 416], [ 11, 24, 122, 224, 333], [ 15, 32, 53, 155, 268], [ 4, 17, 37, 68, 191]]
Bpow17 = [[x**17 for x in line] for line in B]
Bpow17mod703 = [[x%703 for x in line] for line in Bpow17]
# Or, in one line
Res=[[(x**17)%703 for x in line] for line in B]
Or, if you insist on using numpy, then you must ensure that you never overflow the size of an int.
(If using floats, overflow is even easier, even it it takes another form. It is not that you exactly overflow, since maximum value is way bigger that int type of same size. But it is accuracy that goes down when values go up. Up to a point where accuracy cannot even ensure integer values. And the number of exact integer float64 is able to encode is lower than the number of exact integer you can encode with int64. So don't use float just because theoretical maximum value is bigger, when you need exact int values. That is not, in your case the tradeoff "bigger values in exchange for smaller accuracy". Is is "bigger values in exchange of it doesn't work".)
But even with ints, **17 will overflow.
So the trick is never compute such values as 96**17 % 703 directly.
The trick with modulo arithmetics, is that a*b≡c[703] iff a'≡a[703] and b'≡b[703] and a'*c'≡c[703].
For example 1051*2344%703 is 232.
But 1051%703 is 348 and 2344%703 is 235. So 348*235%703 is also 232.
So, same result. But operation is done with smaller values. So less risk of overflow (in the direct case, it uses 2463544 intermediary result, if computing %703 before, bigger intermediary result is 81780)
So, in practice, with numpy, you could compute this
B=np.array([[ 96, 205, 325, 460, 599], [ 13, 109, 207, 307, 416], [ 11, 24, 122, 224, 333], [ 15, 32, 53, 155, 268], [ 4, 17, 37, 68, 191]], dtype=np.uint32)
# Note: I could skip the dtype since by default, numpy would create a `int64` array.
# But I prefer explicit. Plus, I use `int32` to further demontrate how this avoid
# yet even easier overflow
B2=B**2%703
B4=B2**2%703
B8=B4**2%703
B16=B8**2%703
B17=B16*B%703
Result is the expected one.
You can generalize that to other values
def powmod(M, p, m):
if p==1: return M%m
X=powmod(M, p//2, m)
if p%2:
return X*X*M%m
else:
return X*X%m
