Julia

Needs JuMP and MiniZinc.jl. Takes ~0.5 seconds on my machine to solve all 3.

using JuMP
using MiniZinc

# Only need this for Windows.
ENV["JULIA_LIBMINIZINC_DIR"] = "C:\\Program Files\\MiniZinc"

function find_matrix(mat_sum)
    n = length(mat_sum[1, :])
    model = Model(() -> MiniZinc.Optimizer{Float64}("highs"))
    @variable(model, 1 ≤ x[1:n, 1:n] ≤ n^2, Int)
    @constraint(model, vec(x) in MOI.AllDifferent(n^2))
    for i in 1:n
        @constraint(model, sum(x[i, :]) == mat_sum[1, i])
        @constraint(model, sum(x[:, i]) == mat_sum[2, i])
    end
    optimize!(model)
    display(value.(x))
end


find_matrix([12 18 15; 9 18 18])
find_matrix([33 21 43 39; 24 36 41 35])
find_matrix([43 79 58 60 85; 76 66 44 70 69])
```

Java

(Simple) backtracking approach with random candidates:

import java.util.*;
import java.util.concurrent.atomic.AtomicLong;

public class Main {
    public static void main(String[] args) {
        long start = System.currentTimeMillis();
        System.out.println(Arrays.deepToString(calcSolution(new int[][]{{12, 18, 15}, {9, 18, 18}})));
        System.out.println(Arrays.deepToString(calcSolution(new int[][]{{33, 21, 43, 39}, {24, 36, 41, 35}})));
        System.out.println(Arrays.deepToString(calcSolution(new int[][]{{43, 79, 58, 60, 85}, {76, 66, 44, 70, 69}})));
        long end = System.currentTimeMillis();
        System.out.println("Time: " + (end - start) + "ms");

        int n = 20;
        long time1 = 0;
        for (int i = 0; i < n; i++) {
            start = System.currentTimeMillis();
            System.out.println(Arrays.deepToString(calcSolution(new int[][]{{43, 79, 58, 60, 85}, {76, 66, 44, 70, 69}})));
            end = System.currentTimeMillis();
            time1 += end - start;
        }
        time1 /= n;
        System.out.println("Time (average): " + time1 + "ms for " + n + " times");
    }

    public static int[][] calcSolution(int[][] grid) {
        final int n1 = grid[0].length;
        final int n2 = n1 * n1;
        final ArrayList<ArrayList<Integer>> randoms = new ArrayList<>();
        for (int i = 0; i < n2; i++) {
            ArrayList<Integer> list = new ArrayList<>();
            for (int j = 0; j < n2; j++) {
                list.add(j + 1);
            }
            Collections.shuffle(list);
            randoms.add(list);
        }

        int[][] solution = new int[n1][n1];
        while (calcSolution2(grid, n1, n2, solution, 0, new int[2][n1], new HashSet<>(), randoms, new AtomicLong()) == 2) {
            //System.out.println("new round");
            for (int i = 0; i < n2; i++) {
                Collections.shuffle(randoms.get(i));
            }
            solution = new int[n1][n1];
        }

        // (Only for testing/checking for me)
        test(grid, n1, solution);

        return solution;
    }

    private static int calcSolution2(int[][] grid, int n1, int n2, int[][] solution, int idx, int[][] sums, HashSet<Integer> seen, ArrayList<ArrayList<Integer>> randoms, AtomicLong callsCounter) {
        if (idx == n2) {
            return 1;
        }
        if (callsCounter.incrementAndGet() == 1_000_000) {
            return 2;
        }
        int r = idx / n1;
        int c = idx % n1;
        for (int i : randoms.get(idx)) {
            if (seen.contains(i)) {
                continue;
            }
            solution[r][c] = i;
            sums[0][r] += i;
            sums[1][c] += i;
            if (sums[0][r] > grid[0][r] || sums[1][c] > grid[1][c]) {
                sums[0][r] -= i;
                sums[1][c] -= i;
                continue;
            }
            if ((idx + 1) % n1 == 0 && sums[0][r] != grid[0][r]) {
                sums[0][r] -= i;
                sums[1][c] -= i;
                continue;
            }
            if (idx / n1 == n1 - 1 && sums[1][c] != grid[1][c]) {
                sums[0][r] -= i;
                sums[1][c] -= i;
                continue;
            }
            seen.add(i);
            int retVal = calcSolution2(grid, n1, n2, solution, idx + 1, sums, seen, randoms, callsCounter);
            if (retVal > 0) {
                return retVal;
            }
            seen.remove(i);
            sums[0][r] -= i;
            sums[1][c] -= i;
        }
        return 0;
    }

    // (Only for testing/checking for me)
    private static void test(int[][] grid, int n1, int[][] solution) {
        int[][] sums = new int[2][n1];
        for (int i = 0; i < n1; i++) {
            for (int j = 0; j < n1; j++) {
                sums[0][i] += solution[i][j];
                sums[1][i] += solution[j][i];
            }
        }
        for (int i = 0; i < n1; i++) {
            if (sums[0][i] != grid[0][i] || sums[1][i] != grid[1][i]) {
                System.out.println(false);
                return;
            }
        }
        System.out.println(true);
    }
}

Try it online!

Average time on my CPU:

< 200ms

Log:

• 2024-06-14: I use random candidates to fill now.
• 2024-06-17: Improved exit condition. Runs in under 200ms now!

Brachylog, 43 bytes

{hlL;Mṁ₍+ᵐ~h?∧M\+ᵐ~t?∧Mc≠{>0∧.^₂≥?∧}ᵛL∧M≜}ᵐ

Try it online!

Takes about 5 seconds for all 3 cases on TIO.

Explanation

Since Brachylog uses constraint logic arithmetic by default for all integer operations, we don’t even need to do anything special to be able to solve this kind of problem relatively efficiently. The performance of the solve will depend on the performance of the clp(fd) library of SWI-Prolog which is used by Brachylog.

Note that this code is not golfed.

{                                           }ᵐ  Map for each input list:
 hlL                                            L is the number of sums
   L;Mṁ₍                                        M is a square matrix of size L
     M  +ᵐ~h?                                   Sums of rows of M is the first element of the input
             ∧
               M\+ᵐ~t?                          Sums of columns of M is the last element of the input
                      ∧
                        Mc≠                     All elements of M are different
                           {          }ᵛL       For each element:
                            >0                  It is strictly greater than 0
                              ∧
                                .^₂≥?∧          It is greater or equal to L²
                                         ∧
                                          M≜    Find values for elements of M that match all these
                                                constraints, M is the output

Nekomata + -1, 21 bytes

ᵃ#*R↕v#Š:ᵐ∑ᵈv=¿:∑ᵈ§=¿

Attempt This Online!

I just wanted to see how fast my golfing language can go. The language uses a simple backtracking algorithm, and it's not optimized for speed.

The 3x3 and 4x4 grids can be solved in around 4 seconds on ATO. Better than I expected.

Do not try the 5x5 grid. It will eat up all your memory and will never finish.

ᵃ#*R↕v#Š:ᵐ∑ᵈv=¿:∑ᵈ§=¿
ᵃ#*                     Multiply the lengths of both inputs
   R                    Range from 1 to the product
    ↕                   Find a permutation of the range
     v#Š                Split into slices of the lengths of the inputs
        :ᵐ∑ᵈv=¿         Check if the sums of each slice are equal to the first input
               :∑ᵈ§=¿   Check if the (vectorized) sum of all slices is equal to the second input

-1 finds only the first solution.

If you want to run it on your machine, you need to install cabal and ghc, then run the following commands:

git clone https://github.com/AlephAlpha/Nekomata.git
cd Nekomata
cabal build
cabal run Nekomata -- -c 'ᵃ#*R↕v#Š:ᵐ∑ᵈv=¿:∑ᵈ§=¿' -1 -i '[12,18,15] [9,18,18]'
cabal run Nekomata -- -c 'ᵃ#*R↕v#Š:ᵐ∑ᵈv=¿:∑ᵈ§=¿' -1 -i '[33,21,43,39] [24,36,41,35]'

Picat

import cp.

puzzle(Rows, Cols, N) = A =>
	A = new_array(N,N),
	A :: 1..N*N,
	all_distinct([A[R,C]: R in 1..N, C in 1..N]),
	foreach(R in 1..N) sum(A[R]) #= Rows[R] end,
	foreach(C in 1..N) sum([A[R,C]: R in 1..N]) #= Cols[C] end,
	solve(A).

main =>
	println(puzzle([12, 18, 15], [9, 18, 18], 3)),
	println(puzzle([33, 21, 43, 39], [24, 36, 41, 35], 4)),
	println(puzzle([43, 79, 58, 60, 85], [76, 66, 44, 70, 69], 5)).

Output:

{{1,4,7},{3,6,9},{5,8,2}}
{{1,2,14,16},{3,6,8,4},{11,15,7,10},{9,13,12,5}}
{{1,2,3,12,25},{6,15,17,21,20},{22,14,4,10,8},{23,16,7,9,5},{24,19,13,18,11}}

Attempt This Online!

Solves all three puzzles in 0.1 second on ATO.

Python

To run the code, install the z3-solver package. It can usually solve all three cases in ~1s.

from z3 import *

def solve(rows, cols, n):
    grid = [Int(f'a{i}') for i in range(n * n)]

    s = Solver()

    s.add(Distinct(grid))
    for i in range(n * n):
        s.add(Or([grid[i] == j + 1 for j in range(n * n)]))
    
    for i in range(n):
        s.add(sum(grid[i*n:i*n+n]) == rows[i])
        s.add(sum(grid[i::n]) == cols[i])
    
    assert s.check() == sat
    model = s.model()

    sol = [model[g].as_long() for g in grid]
    return [sol[i*n:i*n+n] for i in range(n)]


print(solve([12, 18, 15], [9, 18, 18], 3))
print(solve([33, 21, 43, 39], [24, 36, 41, 35], 4))
print(solve([43, 79, 58, 60, 85], [76, 66, 44, 70, 69], 5))