プログラム

angles.hs

#!/usr/bin/env cabal
{- cabal:
build-depends: base, vector
-}
module Main where

import Data.Foldable (foldl')
import Control.Monad (guard)

import qualified Data.Vector.Unboxed as V

type Table = [V.Vector Int]
type Marginal = V.Vector Int

partition :: Int -> Int -> [[Int]]
partition n len = partitionWithCap n (replicate len n)

partitionWithCap :: Int -> [Int] -> [[Int]]
partitionWithCap n [] = [ [] | n == 0]
partitionWithCap n [k] = [ [n] | n <= k ]
partitionWithCap n (k:ks) =
  do x <- [0 .. min n k]
     xs <- partitionWithCap (n-x) ks
     pure (x:xs)

searchOnMarginal :: Marginal -> [Table]
searchOnMarginal qs = go (zip [1..] (V.toList qs)) qs qs
    where
      len = V.length qs
      isEmpty = V.all (== 0)
      go [] bs cs = [ [] | isEmpty bs && isEmpty cs ]
      go ((i,a):as) bs cs =
        do row <- V.fromList <$> partitionWithCap a (V.toList bs)
           let bs' = V.zipWith (-) bs row
           cs' <- subCs i row cs
           rest <- go as bs' cs'
           pure (row : rest)
      
      subCs i row cs = [ cs' | V.all (>= 0) cs' ]
        where
          cs' = V.imap subCell cs
          subCell k c =
            let j = 2 * len - i - k
                r | 1 <= j && j <= len = row V.! (j-1)
                  | otherwise          = 0
            in c - r

isReduced :: Table -> Bool
isReduced = (== 1) . foldl' gcd 0 . map (V.foldl' gcd 0)

search :: Int -> [Table]
search size =
  do total <- [1..]
     qs <- partition total size
     table <- searchOnMarginal (V.fromList qs)
     guard (isReduced table)
     pure table

printTable :: Table -> IO ()
printTable table = do
    putStrLn header
    mapM_ printRow table
    putStrLn ""
  where
    sumTable = sum $ map V.sum table
    header = "# n=" ++ show sumTable
    printRow = putStrLn . unwords . map show . V.toList

main :: IO ()
main = mapM_ printTable (take 200 $ search 4)

angles.jl

# https://twitter.com/aoki_taichi/status/1469101798552391680

using JuMP
using GLPK

# number of possible angles
# θ = [1:n]/n * π
n = 4

model = Model(GLPK.Optimizer)

# p[i,j] is the probability the tuple of angles (i/n*π, j/n*π, (2*n-i-j)/n*π)
@variable(model, p[1:n, 1:n] >= 0)
# q[i] is the marginal probability of each angle (which should be same for the three angle)
@variable(model, q[1:n] >= 0)

@constraint(model, total, sum(q) == 1)
@constraint(model, [i in 1:n], sum(p[i,:]) == q[i])
@constraint(model, [j in 1:n], sum(p[:,j]) == q[j])
@constraint(model, [k in 1:n], sum(p[i,j] for i in 1:n for j in 1:n if 2 * n - i - j == k) == q[k])
for i in 1:n
    for j in 1:n
        if i + j < n || i + j >= 2*n
            @constraint(model, p[i,j] == 0)
        end
    end
end

#print(model)

optimize!(model)
@show value.(q)

print("p = ")
value.(p)