Strange behavior related to value / reference

[Lambda]

I defined a function to raise a 2x2 matrix to nth power:

def matrix_power(m, n):
result = m[:]
print result is m

for i in range(n - 1):
result[0][0] = result[0][0] * m[0][0] + result[0][1] * m[1][0]
result[0][1] = result[0][0] * m[0][1] + result[0][1] * m[1][1]
result[1][0] = result[1][0] * m[0][0] + result[1][1] * m[1][0]
result[1][1] = result[1][0] * m[0][1] + result[1][1] * m[1][1]
print result == m
print "result:", result, "m:", m
return result

And the output:
False
True
result: [[1, 2], [2, 5]] m: [[1, 2], [2, 5]]
True
result: [[5, 20], [20, 425]] m: [[5, 20], [20, 425]]
True
result: [[425, 17000], [17000, 289180625]] m: [[425, 17000], [17000,
289180625]]

I find the matrix and m lists are always equal.
But these two lists are separate objects, right?
What's wrong with this code?

BTW, numpy has such a function, but it doesn't support really large
number.
Does numpy has some functions that support arbitrarily large number?

 

[Chris Rebert]

I defined a function to raise a 2x2 matrix to nth power:

def matrix_power(m, n):
 result = m[:]

Note that this only copies the *outer* list. It does NOT copy any of
the inner, nested lists, it just makes fresh *references* to them,
which is what's causing your problem; the inner lists get shared, thus
the 2D lists are almost the same list for practical purposes.

You want:
result = [row[:] for row in m]
Which explicitly copies the inner lists.

I find the matrix and m lists are always equal.
But these two lists are separate objects, right?

Right, but the *inner* list objects *are* the same in your code, and
lists with the same contents (in this case, the same inner lists)
compare as equal (==) but not identical (`is`). Remember that all
lists in Python are 1-dimensional; "multidimensional" lists are merely
by convention, hence you have to think a bit more carefully when
working with them.

Cheers,
Chris

 

[Terry Reedy]

I defined a function to raise a 2x2 matrix to nth power:

There is a much faster way to raise x to a count power n than the
definitional but naive method of multipling 1 by x n times. It is based
on the binary representation of n.

Example: x**29 = x**(16+8+4+1) = x**16 * x**8 * x**4 * x * 1
So square x 4 times and do 4 more multiplies (total 8) instead of 28!

General algorithm is something like (untested):

def expo(x, n): # assume n is a count
if n%2: res = x
else: res = 1 # or the identity for type(x)
base = x
n //= 2 # if n < 2, we are done
while n:
base *= base # or appropriate mul function
if n%2: res *= base # ditto
n //= 2
return res

Terry Jan Reedy

 

[Lambda]

I defined a function to raise a 2x2 matrix to nth power:

def matrix_power(m, n):
 result = m[:]

Note that this only copies the *outer* list. It does NOT copy any of
the inner, nested lists, it just makes fresh *references* to them,
which is what's causing your problem; the inner lists get shared, thus
the 2D lists are almost the same list for practical purposes.

You want:
result = [row[:] for row in m]
Which explicitly copies the inner lists.

I find the matrix and m lists are always equal.
But these two lists are separate objects, right?

Right, but the *inner* list objects *are* the same in your code, and
lists with the same contents (in this case, the same inner lists)
compare as equal (==) but not identical (`is`). Remember that all
lists in Python are 1-dimensional; "multidimensional" lists are merely
by convention, hence you have to think a bit more carefully when
working with them.

Cheers,
Chris
--http://blog.rebertia.com

Thank you!
Following is my final code:

def matrix_power(m, n):
"""
Raise 2x2 matrix m to nth power.
"""
if n == 0: return [[1, 0], [0, 1]]

x = matrix_power(m, n / 2)
square = matrix_mul(x, x)
if n % 2 == 0:
return square
else:
return matrix_mul(m, square)

def matrix_mul(a, b):
result = [row[:] for row in a]
result[0][0] = a[0][0] * b[0][0] + a[0][1] * b[1][0]
result[0][1] = a[0][0] * b[0][1] + a[0][1] * b[1][1]
result[1][0] = a[1][0] * b[0][0] + a[1][1] * b[1][0]
result[1][1] = a[1][0] * b[0][1] + a[1][1] * b[1][1]
return result

 

[Mark Dickinson]

Thank you!
Following is my final code:

Looks good, but are you sure about that word 'final'? ;-)

def matrix_power(m, n):
  """
  Raise 2x2 matrix m to nth power.
  """
  if n == 0: return [[1, 0], [0, 1]]

  x = matrix_power(m, n / 2)

I suggest using n // 2 instead of n / 2 here: n // 2 always
does integer divsion (aka 'floor division' in Python), while
the behaviour of n/2 depends on whether you're using Python
2.x or 3.1, and (in 2.x) on whether anyone's put a 'from
__future__ import division' at the top of the file.

  square = matrix_mul(x, x)
  if n % 2 == 0:
    return square
  else:
    return matrix_mul(m, square)

def matrix_mul(a, b):
  result = [row[:] for row in a]
  result[0][0] = a[0][0] * b[0][0] + a[0][1] * b[1][0]
  result[0][1] = a[0][0] * b[0][1] + a[0][1] * b[1][1]
  result[1][0] = a[1][0] * b[0][0] + a[1][1] * b[1][0]
  result[1][1] = a[1][0] * b[0][1] + a[1][1] * b[1][1]
  return result

It's slightly more natural to build the list up as you go,
instead of creating a list of the right size and assigning
to its entries. E.g.,

def matrix_mul(a, b):
row0 = [a[0][0]*b[0][0] + a[0][1]*b[1][0],
a[0][0]*b[0][1] + a[0][1]*b[1][1]]
row1 = <similar>
return [row0, row1]

Better still, use for loops to loop over the elements and
accumulate the sums, and then your code won't need rewriting
when you want to use it to take the power of 3-by-3 matrices.
And then check out the builtin 'sum' function, and consider
replacing some or all of the for loops with list
comprehensions (depending on personal taste)...

Not-so-final'ly yours,

Mark

 

[Peter Otten]

I defined a function to raise a 2x2 matrix to nth power:

BTW, numpy has such a function, but it doesn't support really large
number.
Does numpy has some functions that support arbitrarily large number?

You can tell it to use Python instead of C integers:

import numpy
m = numpy.matrix([[1,2],[3,4]], dtype=object)
m**100

matrix([[2477769229755647135316619942790719764614291718779321308132883358976984999,
3611168042210422417525189130504685706095268999850398812464239094782237250],
[5416752063315633626287783695757028559142903499775598218696358642173355875,
7894521293071280761604403638547748323757195218554919526829242001150340874]],
dtype=object)

Peter

 

[Dave Angel]

Thank you!
Following is my final code:

Looks good, but are you sure about that word 'final'? ;-)

def matrix_power(m, n):
"""
Raise 2x2 matrix m to nth power.
"""
if n =0: return [[1, 0], [0, 1]]

x =atrix_power(m, n / 2)

I suggest using n // 2 instead of n / 2 here: n // 2 always
does integer divsion (aka 'floor division' in Python), while
the behaviour of n/2 depends on whether you're using Python
2.x or 3.1, and (in 2.x) on whether anyone's put a 'from
__future__ import division' at the top of the file.

square =atrix_mul(x, x)
if n % 2 =0:
return square
else:
return matrix_mul(m, square)

def matrix_mul(a, b):
result =row[:] for row in a]
result[0][0] =[0][0] * b[0][0] + a[0][1] * b[1][0]
result[0][1] =[0][0] * b[0][1] + a[0][1] * b[1][1]
result[1][0] =[1][0] * b[0][0] + a[1][1] * b[1][0]
result[1][1] =[1][0] * b[0][1] + a[1][1] * b[1][1]
return result

It's slightly more natural to build the list up as you go,
instead of creating a list of the right size and assigning
to its entries. E.g.,

def matrix_mul(a, b):
row0 =a[0][0]*b[0][0] + a[0][1]*b[1][0],
a[0][0]*b[0][1] + a[0][1]*b[1][1]]
row1 =similar>
return [row0, row1]

Better still, use for loops to loop over the elements and
accumulate the sums, and then your code won't need rewriting
when you want to use it to take the power of 3-by-3 matrices.
And then check out the builtin 'sum' function, and consider
replacing some or all of the for loops with list
comprehensions (depending on personal taste)...

Not-so-final'ly yours,

exp = 600
mat = [ [3, 0], [0, 3]]
print mat, exp, matrix_power(mat, exp)
mat = [ [3.0, 0], [0, 3.0]]
print mat, exp, matrix_power(mat, exp); print power(mat, exp)

[[3, 0], [0, 3]] 600
[[18739277038847939886754019920358123424308469030992781557966909983211910963157763678726120154469030856807730587971859910379069087693119051085139566217370635083384943613868029545256897117998608156843699465093293765833141309526696357142600866935689483770877815014461194837692223879905132001L,
0L], [0L,
18739277038847939886754019920358123424308469030992781557966909983211910963157763678726120154469030856807730587971859910379069087693119051085139566217370635083384943613868029545256897117998608156843699465093293765833141309526696357142600866935689483770877815014461194837692223879905132001L]]
[[3.0, 0], [0, 3.0]] 600 [[1.8739277038847955e+286, 0.0], [0.0,
1.8739277038847955e+286]]
[[1.8739277038847965e+286, 0.0], [0.0, 1.8739277038847965e+286]]

Mark

That square/multiply technique brought back fond memories. I used it
when I did the microcode for floating point arithmetic in 1975 or so.
Exponentiation used logs in the general case, but for integers below 100
used this technique. (At the time, the hardware instruction set was
capable of multiplying only single digit numbers, the rest was loops in
my microcode)

Probably I should also mention that, if the values had been floats
instead of integers, the square/multiply technique generally improves
accuracy, not just speed. While repeatedly multiplying int/long values
will eventually get the same answer, the quantization errors if they're
floats do add up faster. If we use the matrix

mat = [ [3, 0], [0, 3]] n = 600
The first bunch of digits of the upper left corner are:
1873927703884793988675401...
If we give mat floating values, and run it again:
1.8739277038847955e+286
and if we use simple loop, multiplying 600 times, we get:
1.8739277038847965e+286

DaveA

 

[Aahz]

def matrix_power(m, n):
result = m[:]
print result is m

Use copy.deepcopy()
--
Aahz ([email protected]) <*> http://www.pythoncraft.com/

"You could make Eskimos emigrate to the Sahara by vigorously arguing --
at hundreds of screens' length -- for the wonder, beauty, and utility of
snow." --PNH to rb in r.a.sf.f