# Compute the perimeter density matrix

## Introduction

The perimeter density matrix is an infinite binary matrix M defined as follows. Consider a (1-based) index (x, y), and denote by M[x, y] the rectangular sub-matrix spanned by the corner (1, 1) and (x, y). Suppose that all values of M[x, y] except Mx, y, the value at index (x, y), have already been determined. Then the value Mx, y is whichever of 0 or 1 that puts the average value of M[x, y] closer to 1 / (x + y). In case of a tie, choose Mx, y = 1.

This is the sub-matrix M[20, 20] with zeros replaced by dots for clarity:

1 . . . . . . . . . . . . . . . . . . .
. . . . . 1 . . . . . . . . . . . . . .
. . 1 . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . 1 . . . . . . . . . . . . . . .
. 1 . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . 1 . .
. . . . . . . . . . . . . . 1 . . . . .
. . . . . . . . . . . . 1 . . . . . . .
. . . . . . . . . . 1 . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 1 . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . 1 . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . 1 . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .


For example, we have M1, 1 = 1 at the upper left corner, since 1 / (1 + 1) = ½, and the average of the 1 × 1 sub-matrix M[1, 1] is either 0 or 1; that's a tie, so we choose 1.

Consider then the position (3, 4). We have 1 / (3 + 4) = 1/7, and the average of the sub-matrix M[3, 4] is 1/6 if we choose 0, and 3/12 if we choose 1. The former is closer to 1/7, so we choose M3, 4 = 0.

Here is the sub-matrix M[800, 800] as an image, showing some of its intricate structure.

Given a positive integer N < 1000, output the N × N sub-matrix M[N, N], in any reasonable format. The lowest byte count wins.

# R, 158 154 141 bytes

Edit: Because the only 1 in the upper 2x2 submatrix is the top left M[1,1] we can start the search for 1s when {x,y}>1 so no need for the if statement.

M=matrix(0,n<-scan(),n);M[1]=1;for(i in 2:n)for(j in 2:n){y=x=M[1:i,1:j];x[i,j]=0;y[i,j]=1;d=1/(i+j);M[i,j]=abs(d-mean(x))>=abs(d-mean(y))};M


The solution is highly inefficient as the matrix is duplicated twice for each iteration. n=1000 took just under two and a half hours to run and produces a matrix of 7.6 Mb.

Ungolfed and explained

M=matrix(0,n<-scan(),n);                        # Read input from stdin and initialize matrix with 0s
M[1]=1;                                         # Set top left element to 1
for(i in 2:n){                                  # For each row
for(j in 2:n){                              # For each column
y=x=M[1:i,1:j];                         # Generate two copies of M with i rows and j columns
x[i,j]=0;                               # Set bottom right element to 0
y[i,j]=1;                               # Set bottom right element to 1
d=1/(i+j);                              # Calculate inverse of sum of indices
M[i,j]=abs(d-mean(x))>=abs(d-mean(y))   # Returns FALSE if mean(x) is closer to d and TRUE if mean(y) is
}
};
M                                               # Print to stdout


Output for n=20

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,]     1    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[2,]     0    0    0    0    0    1    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[3,]     0    0    1    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[4,]     0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[5,]     0    0    0    0    1    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[6,]     0    1    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[7,]     0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[8,]     0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     1     0     0
[9,]     0    0    0    0    0    0    0    0    0     0     0     0     0     0     1     0     0     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0     0     0     0     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     1     0     0     0     0     0     0     0     0     0
[12,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[13,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0     0     0     0     0     0     0
[14,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[15,]    0    0    0    0    0    0    0    0    1     0     0     0     0     0     0     0     0     0     0     0
[16,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[17,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[18,]    0    0    0    0    0    0    0    1    0     0     0     0     0     0     0     0     0     0     0     0
[19,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0
[20,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0


# BQN, 36 bytesSBCS

+⁼{𝔽𝔽˘>(⊣-⌊∘-)<˘⌊0.5+(×÷+)⌜˜1+↕𝕩}


Run online!

Two observations:

• Minimizing distance of mean to 1/(x+y) is the same as minimizing distance of sum to x×y/(x+y).
• At least up to 1000×1000, which is maximum size required in the challenge, there is no row or column with two ones.

The function performs three major steps:

⌊0.5+(×÷+)⌜˜1+↕𝕩 Generate matrix round(x×y/(x+y)).

(⊣-⌊∘-)<˘ Calculate matrix sum(M[x, y]).

+⁼+⁼˘ Convert that to the perimeter density matrix by performing inverse cumulative sum (+⁼) on rows and columns.

This runs reasonably fast, generating the full 1000x1000 in ~60ms locally, but some less golfed version were a bit faster.

# Python, 129119117 115 bytes

• Replaced the loop with a list comprehension so everything can go on one line and can be lambda
• Replaced .append() with :=z+[]
• -2 bytes thanks to @ceilingcat for pointing out that (1+a)*(1+c) = (-a-1)*(-b-1)
• -1 bytes thanks to @ceilingcat for removing the 0 from .5
• -1 byte thanks to @thonnu for a more efficinet way to do and z at the end.
lambda x,z=[]:[z:=z+[~(i//x)*~(i%x)/(i//x+i%x+2)>=sum(z[k]&(k%x<=i%x)for k in range(i))+.5]for i in range(x*x)][-1]


Attempt This Online!

Rust is no longer better than python. Port of my Rust answer. Outputs a flattened matrix

# Python, 118 111 bytes

def f(x,z=[]):
for i in range(x*x):z+=[~(i//x)*~(i%x)/(i//x+i%x+2)>=sum(z[k]&(k%x<=i%x)for k in range(i))+0.5]


Attempt This Online!

Thanks to @xnor and @ceilingcat. Outputs my modifying the input in place.

• 117 bytes I think I/O defaults also allow taking in z=[] as input and modifying it without returning.
– xnor
May 10, 2023 at 17:47
• 115 bytes (includes @ceilingcat's suggestion) May 14, 2023 at 8:16
• z[k]>(k%x>i%x) saves a byte, switching z from list to tuple saves another (lambda x,*z:/z:=z+(...,))
– ovs
May 14, 2023 at 11:27

# Scala 3, 565476460458447 440 bytes

Golfed version. Attempt this online!

Saved 125 bytes (565->440) thanks to the comment of @ceilingcat

@main def m={val n=scala.io.StdIn.readInt();val M=Array.ofDim[Int](n,n);M(0)(0)=1;for(i<-1 until n;j<-1 until n){val z=M.slice(0,i+1);val x=z.map(_.slice(0,j+1));val y=z.map(_.slice(0,j+1));x(i)(j)=0;y(i)(j)=1;val d=1.0/(i+j+2);M(i)(j)=if(a(x,d)<a(y,d))0 else 1};M.foreach(o=>{o.foreach(l=>print(if(l<1)"."else 1));println("")});def a(m:Array[Array[Int]],d:Double)=math.abs(d*m.length*m.headOption.map(_.length).getOrElse(0)-m.flatten.sum)}


Ungolfed version. Attempt this online!

import scala.io.StdIn.readInt

object Main extends App {
val M = Array.ofDim[Int](n, n)
M(0)(0) = 1

for (i <- 1 until n) {
for (j <- 1 until n) {
val x = M.slice(0, i + 1).map(_.slice(0, j + 1))
val y = M.slice(0, i + 1).map(_.slice(0, j + 1))
x(i)(j) = 0
y(i)(j) = 1
val d = 1.0 / (i + j + 2)

M(i)(j) = if (math.abs(d - mean(x)) >= math.abs(d - mean(y))) 1 else 0
}
}

for (row <- M) {
for (ele <- row) {
if (ele == 0) print(".") else print(1)
}
println("")
}

def mean(matrix: Array[Array[Int]]): Double = {
val rowCount = matrix.length
matrix.flatten.sum.toDouble / (rowCount * colCount)
}
}


• Can you do import io.StdIn._ instead?
– user
May 11, 2023 at 3:14

# Python 2, 189 Bytes

There are no crazy tricks in here, it is just calculating as described in the introduction. It isn't particularly quick but I don't need to create any new matrices to do this.

n=input()
k=[n*[0]for x in range(n)]
for i in range(1,-~n):
for j in range(1,-~n):p=1.*i*j;f=sum(sum(k[l][:j])for l in range(i));d=1./(i+j);k[i-1][j-1]=0**(abs(f/p-d)<abs(-~f/p-d))
print k


Explanation:

n=input()                                     # obtain size of matrix
k=[n*[0]for x in range(n)]                    # create the n x n 0-filled matrix
for i in range(1,-~n):                        # for every row:
for j in range(1,-~n):                      # and every column:
p=1.*i*j                                  # the number of elements 'converted' to float
f=sum(sum(k[l][:j])for l in range(i))     # calculate the current sum of the submatrix
d=1./(i+j)                                # calculate the goal average
k[i-1][j-1]=0**(abs(f/p-d)<abs(-~f/p-d))  # decide whether cell should be 0 or 1
print k                                       # print the final matrix


For those curious, here are some timings:

 20 x  20 took 3 ms.
50 x  50 took 47 ms.
100 x 100 took 506 ms.
250 x 250 took 15033 ms.
999 x 999 took 3382162 ms.


"Pretty" output for n = 20:

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0


# Rust, 152 bytes

|z|{let mut b=vec![];for i in 0..z*z{b.push(((1+i/z)*(1+i%z))as f32/(i/z+i%z+2)as f32>=(0..i).map(|d|(d%z<=i%z&&b[d])as u8 as f32).sum::<f32>()+0.5)};b}


Attempt This Online!

Lots of wasted space due to needing as f32 everywhere :(

## Racket 294 bytes

(define(g x y)(if(= 1 x y)1(let*((s(for*/sum((i(range 1(add1 x)))(j(range 1(add1 y)))#:unless(and(= i x)(= j y)))
(g i j)))(a(/ s(* x y)))(b(/(add1 s)(* x y)))(c(/ 1(+ x y))))(if(<(abs(- a c))(abs(- b c)))0 1))))


Ungolfed:

(define(f a b)
(define (g x y)
(if (= 1 x y) 1
(let* ((s (for*/sum ((i (range 1 (add1 x)))
#:unless (and (= i x) (= j y)))
(g i j)))
(a (/ s (* x y)))
(b (/ (add1 s) (* x y)))
(c (/ 1 (+ x y))))
(if (< (abs(- a c))
(abs(- b c)))
0 1))))
(for ((i (range 1 (add1 a))))
(for ((j (range 1 (add1 b))))
(print (g i j)))
(displayln ""))
)


Testing:

(f 8 8)


Output:

10000000
00000100
00100000
00000000
00001000
01000000
00000000
00000000


# Perl, 151 + 1 = 152 bytes

Run with the -n flag. The code will only work correctly every other iteration within the same instance of the program. To get it to work correctly every time, add 5 bytes by prepending my%m; to the code.

for$b(1..$_){for$c(1..$_){$f=0;for$d(1..$b){$f+=$m{"$d,$_"}/($b*$c)for 1..$c}$g=1/($b+$c);print($m{"$b,$c"}=abs$f-$g>=abs$f+1/($b*$c)-$g?1:_).$"}say""}''  Readable: for$b(1..$_){ for$c(1..$_){$f=0;
for$d(1..$b){
$f+=$m{"$d,$_"}/($b*$c)for 1..$c }$g=1/($b+$c);
print($m{"$b,$c"}=abs$f-$g>=abs$f+1/($b*$c)-$g?1:_).$"
}
say""
}


Output for input of 100:

1___________________________________________________________________________________________________
_____1______________________________________________________________________________________________
__1_________________________________________________________________________________________________
___________________________1________________________________________________________________________
____1_______________________________________________________________________________________________
_1__________________________________________________________________________________________________
_________________________1__________________________________________________________________________
_________________1__________________________________________________________________________________
______________1_____________________________________________________________________________________
____________1_______________________________________________________________________________________
__________1_________________________________________________________________________________________
____________________________________________________________________________________________________
_________1__________________________________________________________________________________________
____________________________________________________________________________________________________
________1___________________________________________________________________________________________
______________________________________________________________________________________1_____________
_________________________________________________________________1__________________________________
_______1____________________________________________________________________________________________
_____________________________________________________________1______________________________________
____________________________________________________1_______________________________________________
______________________________________________1_____________________________________________________
__________________________________________1_________________________________________________________
_______________________________________1____________________________________________________________
____________________________________1_______________________________________________________________
__________________________________1_________________________________________________________________
______1_____________________________________________________________________________________________
____________________________________________________________________________________________________
___1________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
________________________________1___________________________________________________________________
____________________________________________________________________________________________________
________________________1___________________________________________________________________________
____________________________________________________________________________________________________
_______________________1____________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
______________________1_____________________________________________________________________________
____________________________________________________________________________________________________
___________________________________________________________________________________________________1
_____________________1______________________________________________________________________________
____________________________________________________________________________________________________
_____________________________________________________________________________________1______________
__________________________________________________________________________________1_________________
____________________1_______________________________________________________________________________
____________________________________________________________________________________________________
________________________________________________________________________________1___________________
______________________________________________________________________________1_____________________
___________________________________________________________________________1________________________
_________________________________________________________________________1__________________________
___________________1________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
_______________________________________________________________________1____________________________
______________________________________________________________________1_____________________________
____________________________________________________________1_______________________________________
___________________________________________________________1________________________________________
__________________________________________________________1_________________________________________
_________________________________________________________1__________________________________________
__________________1_________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
________________1___________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
________________________________________________________1___________________________________________
_______________________________________________________1____________________________________________
____________________________________________________________________________________________________
___________________________________________________1________________________________________________
____________________________________________________________________________________________________
__________________________________________________1_________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
_________________________________________________1__________________________________________________
____________________________________________________________________________________________________
________________________________________________1___________________________________________________
____________________________________________________________________________________________________
_____________________________________________1______________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________1_______________________________________________________
_______________1____________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
____________________________________________________________________________________________________
_________________________________________1__________________________________________________________


# Wolfram Language (Mathematica), 180 bytes

(a=SparseArray[{1,1}->1,{#,#},-1];({x,y}={##};If[(t=a[[##]])<0,s=0;Do[If[{i,j}!={##},s+=#0[i,j]],{i,#1},{j,#2}];l=s(q=1/x/y)-1/(x+y);r=l+q;a[[##]]=If[Abs@r>Abs@l,0,1],t])&[#,#];a)&


Try it online!

Recursive approach: start from array filled with undefined -1 items and 1 for first item.
For n 20-30 as faster, as Python.
Full version self-explained:

Clear[m];
n = 15;
arr = SparseArray[{1, 1} -> 1, {n, n}, -1];
m[x_, y_] :=
If[(current = arr[[x, y]]) < 0,

total = 0; Do[total += m[i, j], {i, x}, {j, If[i == x, y - 1, y]}];
tmp1 = total/x/y - 1/(x + y); tmp2 = tmp1 + 1/x/y;
arr[[x, y]] = If[Abs@tmp2 > Abs@tmp1, 0, 1]

, current
];
m[n, n];
StringJoin[Append[#, "\n"] & /@ (Normal@arr /. {1 -> "*", 0 -> " "})]