Split a 2d array creating array from row to row values with unitary displacement

The condition for the subarrays to exactly fit is (M+1)*(M-1) divides (N+1)*(N-1), the number of subarrays you can put is the quotient of these numbers. Note that these numbers are equal to M*M-1 and N*N-1. In this form the rule also applies to non square matrices.

Examples

M      N      M*M-1  N*N-1  Y
-----------------------------
3      5      8      24     3
3      7      8      48     6
5      7      24     48     2
4      11     15     120    8

Implementation: Please note that this returns overlapping views into the original array. If you want to modify them you may want to make a copy. Also note that this implementation fits as many subsquares as fit, any leftover elements in the larger matrix are dropped.

Update: I've added two functions which calculate given N all possible M and vice versa.

Output:

# Testing predictions ...
# ok

# M = 105
# solutions: [105, 1273, 1377, 4135, 4239, 5407, 5511, 5513]
# this pattern repeats at offsets 5512, 11024, 16536, ...

# N = 1000001
# solutions: [2, 3, 4, 5, 7, 9, 11, 31, 49, 1000001]

# example N, M = (5, 3)
# [[[ 0  1  2]
#   [ 3  4  5]
#   [ 6  7  8]]

#  [[ 8  9 10]
#   [11 12 13]
#   [14 15 16]]

#  [[16 17 18]
#   [19 20 21]
#   [22 23 24]]]

Code:

import numpy as np
import sympy
import itertools as it
import functools as ft
import operator as op

def get_subsquares(SqNN, M0, M1=None):
    M1 = M0 if M1 is None else M1
    N0, N1 = SqNN.shape
    K = (N0*N1-1) // (M0*M1-1)
    SqNN = SqNN.ravel()
    s, = SqNN.strides
    return np.lib.stride_tricks.as_strided(SqNN, (K, M0, M1),
                                           (s*(M0*M1-1), s*M1, s))


def get_M_for_N(N):
    """Given N return all possible M
    """
    assert N >= 2
    f = 1 + (N & 1)
    factors = sympy.factorint((N+1)//f)
    factors.update(sympy.factorint((N-1)//f))
    if f == 2:
        factors[2] += 2
    factors = [ft.reduce(op.mul, fs) for fs in it.product(
        *([a**k for k in range(n+1)] for a, n in factors.items()))]
    return [fs + 1 for fs in sorted(set(factors) & set(fs - 2 for fs in factors)) if (N*N-1) % (fs * (fs+2)) == 0]

def get_N_for_M(M):
    """Given M return all possible N in the form period, smallest

    smallest is a list of all solutions between M and M*M if M is even
    and between M and (M*M+1) / 2 if M is odd, all other solutions can be
    obtained by adding multiples of period
    """
    assert M >= 2
    f = 1 + (M & 1)
    factors = sympy.factorint((M+1)//f)
    factors.update(sympy.factorint((M-1)//f))
    factors = [k**v for k, v in factors.items()]
    rep = (M+1)*(M-1) // f
    f0 = [ft.reduce(op.mul, fs) for fs in it.product(*zip(it.repeat(1), factors))]
    f1 = [rep // (f*a) for a in f0]
    inv = [f if b==1 else f*b + 2 if a==1 else 2 * sympy.mod_inverse(a, b)
           for a, b in zip(f1, f0)]
    if f==1:
        inv[1:-1] = [a%b for a, b in zip(inv[1:-1], f0[1:-1])]
    return rep, sorted(a*b - 1 for a, b in zip(f1, inv))

def test_predict(N):
    def brute_force(a, b):
        return [i for i in range(a, b) if (i*i-1) % (a*a-1) == 0]
    for x in range(2, N+1):
        period, pred = get_N_for_M(x)
        assert brute_force(x, period*4+2) \
            == [a + b for b in range(0, 4*period, period) for a in pred]
    def brute_force(b):
        return [i for i in range(2, b+1) if (b*b-1) % (i*i-1) == 0]
    for x in range(2, N+1):
        assert brute_force(x) == get_M_for_N(x)
    print('ok')

# test
print("Testing predictions ...")
test_predict(200)
print()
# examples
M = 105
period, pred = get_N_for_M(M)
print(f"M = {M}")
print(f"solutions: {pred}")
print(f"this pattern repeats at offsets {period}, {2*period}, {3*period}, ...")
print()
N = 1000001
pred = get_M_for_N(N)
print(f"N = {N}")
print(f"solutions: {pred}")
print()
N, M = 5, 3
SqNN = np.arange(N*N).reshape(N, N)
print(f"example N, M = ({N}, {M})")
print(get_subsquares(SqNN, M))