e2lstr - Online Python3 Interpreter & Debugging Tool

from math import gcd
import random

def _coprime(a, n):
    return gcd(a, n) == 1

def _find_linear_coeffs(N, max_tries=5000, rng_seed=None):
    """
    Find (a1,b1,c1,a2,b2,c2) so that:
      - a1,b1,a2,b2 are coprime to N
      - det = a1*b2 - a2*b1 is nonzero mod N  (orthogonality)
      - a1±b1 and a2±b2 are coprime to N      (diagonals are permutations)
    c1=c2=0 is fine; we keep them 0.
    """
    if rng_seed is not None:
        random.seed(rng_seed)

    # search space: units modulo N (coprime to N)
    units = [x for x in range(1, N) if _coprime(x, N)]
    for _ in range(max_tries):
        a1, b1 = random.choice(units), random.choice(units)
        a2, b2 = random.choice(units), random.choice(units)
        if (a1 * b2 - a2 * b1) % N == 0:
            continue
        # diagonals: need slopes (a±b) to be units too
        if not (_coprime((a1 + b1) % N or N, N) and _coprime((a1 - b1) % N or N, N)):
            continue
        if not (_coprime((a2 + b2) % N or N, N) and _coprime((a2 - b2) % N or N, N)):
            continue
        return a1, b1, 0, a2, b2, 0

    raise RuntimeError(f"Could not find suitable coefficients for N={N} (try increasing max_tries).")

def build_magic_square(N, base=None, rng_seed=None):
    """
    Build an N×N 'digit-balanced' magic square for odd N.
    Values are B*L1 + L2 where L1,L2 are orthogonal Latin squares
    with diagonal-permutation property.
    Returns the matrix and the base B.
    """
    if N <= 0 or N % 2 == 0:
        raise ValueError("N must be a positive odd integer.")
    B = max(10, N) if base is None else base
    a1, b1, c1, a2, b2, c2 = _find_linear_coeffs(N, rng_seed=rng_seed)

    A = [[0]*N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            l1 = (a1*i + b1*j + c1) % N
            l2 = (a2*i + b2*j + c2) % N
            A[i][j] = B*l1 + l2
    return A, B, (a1, b1, a2, b2)

def verify_properties(A, B):
    """
    Verify that all rows, columns, main & anti-diagonals have the same sum.
    Also checks 'inverse' (swap tens/ones -> swap L1 and L2) keeps the same sum.
    """
    N = len(A)
    # row/col sums
    row_sums = [sum(r) for r in A]
    col_sums = [sum(A[i][j] for i in range(N)) for j in range(N)]
    target = row_sums[0]
    ok_rows = all(s == target for s in row_sums)
    ok_cols = all(s == target for s in col_sums)
    main = sum(A[i][i] for i in range(N))
    anti = sum(A[i][N-1-i] for i in range(N))
    ok_diags = (main == target) and (anti == target)

    # inverse: swap digits (tens<->ones), i.e., A' = (A % B)*B + (A // B)
    inv = [[(x % B)*B + (x // B) for x in row] for row in A]
    inv_row = [sum(r) for r in inv]
    ok_inverse = all(s == inv_row[0] == target for s in inv_row) and \
                 all(sum(inv[i][j] for i in range(N)) == target for j in range(N)) and \
                 sum(inv[i][i] for i in range(N)) == target and \
                 sum(inv[i][N-1-i] for i in range(N)) == target

    return {
        "target_sum": target,
        "rows_ok": ok_rows,
        "cols_ok": ok_cols,
        "diagonals_ok": ok_diags,
        "inverse_ok": ok_inverse
    }

if __name__ == "__main__":
    # demo
    N = 5                      # <- change this to any odd N (3,5,7,9,...)
    A, B, coeffs = build_magic_square(N, rng_seed=42)
    props = verify_properties(A, B)

    print(f"N={N}, base B={B}, coefficients (a1,b1,a2,b2)={coeffs}")
    print("Magic sum:", props["target_sum"])
    print("Rows/Cols/Diagonals OK? ", props["rows_ok"], props["cols_ok"], props["diagonals_ok"])
    print("Inverse (swap digits) also OK? ", props["inverse_ok"])
    print("\nSquare:")
    for row in A:
        print(" ".join(f"{x:>3d}" for x in row))