DivideError in affine_transform_mpo

[MFrankenbach]

I did some unit tests involving affine transform and found the following. In this affine transform example with rational entries:

using ITensors
using Quantics

function main()
    A = [ 1      0      0;
         -1//7   6//7   6//7;
        -10//21  6//7  -1//7]
    b =[0,27//7,61//14]

    R = 6
    D = 3
    y_sites = [ Index(2, "d=$(d),r=$(r)") for d in 1:D for r in 1:R ]
    x_sites = prime.(y_sites)
    y_sites_mat = reshape(y_sites, R, D)
    x_sites_mat = reshape(x_sites, R, D)
    display(y_sites_mat)
    display(x_sites_mat)
    # CRASHES (DivisionError)
    trafo_mpo = Quantics.affine_transform_mpo(
        y_sites_mat,
        x_sites_mat,
        A,
        b,
        Quantics.OpenBoundaryConditions()
    )
end

The code crashes because affine_transform_core returns a value_out whose first size is 0. The error message and stacktrace are:

ERROR: DivideError: integer division error
Stacktrace:
  [1] div
    @ ./int.jl:295 [inlined]
  [2] divrem
    @ ./div.jl:203 [inlined]
  [3] divrem
    @ ./div.jl:181 [inlined]
  [4] _reshape_uncolon
    @ ./reshapedarray.jl:138 [inlined]
  [5] reshape
    @ ./reshapedarray.jl:129 [inlined]
  [6] reshape
    @ ./reshapedarray.jl:128 [inlined]
  [7] affine_transform_tensors(R::Int64, A::StaticArraysCore.SMatrix{3, 3, Int64, 9}, b::StaticArraysCore.SVector{3, Int64}, s::Int64; boundary::Quantics.OpenBoundaryConditions)
    @ Quantics ~/Tensor4Fields/Quantics.jl/src/affine.jl:169
  [8] affine_transform_tensors(R::Int64, A::Matrix{Rational{Int64}}, b::Vector{Rational{Int64}}, boundary::Quantics.OpenBoundaryConditions)
    @ Quantics ~/Tensor4Fields/Quantics.jl/src/affine.jl:99
  [9] affine_transform_mpo(y::Matrix{Index{Int64}}, x::Matrix{Index{Int64}}, A::Matrix{Rational{Int64}}, b::Vector{Rational{Int64}}, boundary::Quantics.OpenBoundaryConditions)
    @ Quantics ~/Tensor4Fields/Quantics.jl/src/affine.jl:63
 [10] main()
    @ Main ~/Tensor4Fields/Quantics.jl/affinetransform_bug_MWE.jl:19
 [11] top-level scope
    @ REPL[13]:1

Printing the sizes of a few objects in affine_transform_tensors gives:

size(cap_matrix) = (1, 0)
size(tensors[1]) = (0, 0, 8, 8)
size.(tensors) = [(0, 0, 8, 8), (0, 0, 8, 8), (0, 0, 8, 8), (0, 0, 8, 8), (0, 0, 8, 8), (0, 1, 8, 8)]
size.(tensors_ext) = [(1, 0, 1, 1), (0, 0, 1, 1)]

demonstrating the zero size problem.

Branch: main, latest commit
Julia Version: 1.11.7

[shinaoka]

Analysis

The root cause is that this particular affine transformation has no valid integer mappings at all.

Why no valid mappings exist

The common denominator for the rational entries is s = lcm(1, 7, 21, 14) = 42 (even). After scaling:

A_int = 42 * A = [42  0   0;
                  -6  36  36;
                 -20  36  -6]

b_int = 42 * b = [0, 162, 183]

For the transformation y = (A_int * x + b_int) / s to yield integers, we need A_int * x + b_int ≡ 0 (mod 42).

Looking at the 3rd component: A_int[3,:] = [-20, 36, -6] (all even) and b_int[3] = 183 (odd). So A_int[3,:]*x + 183 is always odd, and can never be divisible by 42 (even). The transformation maps no integer grid points to integer grid points.

What happens in the code

When s is even, affine_transform_core requires z = A*x + b + carry to be all even. Since b_int[3] is odd and all entries of A_int are even, z[3] is always odd → no valid transitions are found → the carry set becomes empty (length 0) → downstream reshape(tensor, 0, :) causes 0 ÷ 0 → DivideError.

Fix

PR #47 handles empty carry sets by returning zero tensors (representing a zero MPO) with a warning, instead of crashing. This correctly represents a transformation with no valid mappings.

[MFrankenbach]

Oh, of course. Sorry, I forgot about that requirement while coding. Thanks for pointing this out.