delacorte_number.py

distance_cache = {}
gcd_cache = {}
def DelacorteInit():
  # TODO: This half of the initialization could be done faster.
  for row1 in range(28):
    for col1 in range(28):
      for row2 in range(28):
        row_distance = (row1 - row2) ** 2
        for col2 in range(28):
          distance_cache[row1, col1, row2, col2] = (
            row_distance + (col1 - col2) ** 2)

  for a in range(1, 27**2+1):
    for b in range(a+1, 27**2+1):
      gcd_cache[a, b] = fractions.gcd(a, b)
      gcd_cache[b, a] = fractions.gcd(a, b)

def DelacorteScore(values):
  """Calculates the Delacorte number of a grid.

  Uses a cache to save distance and gcd values from previous runs.

  Args:
    values: a list of lists of integers. Each contained list is a row, and
    all of the rows must be the same size.

  Returns:
    int, as defined at <http://www.azspcs.net/Contest/DelacorteNumbers>
  """
  if not distance_cache or not gcd_cache:
    DelacorteInit()

  width, height = len(values[0]), len(values)
  num = 0
  for row1 in range(height):
    for col1 in range(width):
      # First, the rest of the numbers in this row.
      for col2 in range(col1 + 1, width):
        num += distance_cache[row1, col1, row1, col2] * gcd_cache[values[row1][col1], values[row1][col2]]
      # Then the rest of the rows below this one.
      for row2 in range(row1 + 1, height):
        for col2 in range(width):
          num += distance_cache[row1, col1, row2, col2] * gcd_cache[values[row1][col1], values[row2][col2]]
  return num

print DelacorteScore( [[4,1,3],[6,5,2]] )  # 53

Look up functools.lru_cache https://docs.python.org/3/library/functools.html#functools.lru_cache

If you are not using cython to speed up your code, you can use this and not have to write your own caching code.

Also, if you make your matrix linear instead of two-dimensional and convert an ordinal position to row-column, your code to calculate the delacorte square number can be much simpler:

def delacorte_score(matrix, coordinates):
    return np.sum(
        distance(coordinates[i], coordinates[j]) * gcd(p, q)
        for i, p in enumerate(matrix)
        for j, q in enumerate(matrix[i + 1:], start=i + 1)
    )

or:

def delacorte_score(matrix, coordinates):
    return np.sum(
        distance(coordi, coordj) * gcd(p, q)
        for coordi, p in zip(coordinates, matrix)
        for coordj, q in zip(coordinates, matrix)
    ) // 2

depending on whether generating the slices is expensive enough to warrant iterate n x n over the matrix and divide by two.

For one thing, the distance from X to X is 0 so including the diagonal has no effect on the computation.

Also, you need to import fractions.