Подсчитайте количество матриц Hankelable

12

Фон

Бинарная ганкелева матрица - это матрица с постоянными косыми диагоналями (положительными наклонными диагоналями), содержащая только 0s и 1s. Например, бинарная ганкелева матрица 5x5 выглядит следующим образом

a b c d e
b c d e f
c d e f g
d e f g h
e f g h i

где a, b, c, d, e, f, g, h, iлибо 0или 1.

Давайте определим матрицу M как ганкелеву, если есть перестановка порядка строк и столбцов M, так что M является ганкелевой матрицей. Это означает, что можно применить одну перестановку к порядку строк и, возможно, другую к столбцам.

Соревнование

Задача состоит в том, чтобы посчитать, сколько Hankelable n по nматрицам существует для всехn до максимально возможного значения.

Выход

Для каждого целого числа n от 1 и выше выведите число Hankelable n по nматрицам с записями, которые являются 0или 1.

Для n = 1,2,3,4,5ответов должны быть2,12,230,12076,1446672 . (Спасибо orlp за код для их создания.)

Лимит времени

Я запусту ваш код на моей машине и остановлю его через 1 минуту. Код, который выводит правильные ответы до наибольшего значения n побед. Ограничение по времени для всего от n = 1самой большой стоимостиn на которое вы даете ответ.

Победитель будет лучшим ответом к концу субботы 18 апреля.

Tie Breaker

В случае ничьей по максимуму nя определю, сколько времени потребуется, чтобы довести результаты до максимума и выиграть n+1самый быстрый. В случае, если они работают одновременно с точностью до секунды доn+1 , выигрывает первая заявка.

Языки и библиотеки

Вы можете использовать любой язык со свободно доступным компилятором / интерпретатором и т. Д. для Linux и любых библиотек, которые также свободно доступны для Linux.

Моя машина

Время будет работать на моей машине. Это стандартная установка Ubuntu на восьмиъядерный процессор AMD FX-8350 на материнской плате Asus M5A78L-M / USB3 (Socket AM3 +, 8 ГБ DDR3). Это также означает, что мне нужно иметь возможность запустить ваш код. Как следствие, используйте только легкодоступное бесплатное программное обеспечение и, пожалуйста, включите подробные инструкции по компиляции и запуску вашего кода.

Примечания

Я рекомендую не делать итерации по всем n по n матрицам и пытаться определить, обладает ли каждая из них свойством, которое я описываю. Во-первых, их слишком много, а во-вторых, кажется, нет быстрого способа сделать это обнаружение .

Ведущие записи пока

  • n = 8 Питер Тейлор. Джава
  • n = 5 по orlp. питон
Сообщество
источник
4
Мне очень нравится слово «Hankelable».
Алекс А.
3
Для n=6всего есть 260357434. Я думаю, что нехватка памяти является более важной проблемой, чем время процессора.
Питер Тейлор
Это потрясающий вопрос. Я был полностью обрезан ботаником.
Александр-Бретт

Ответы:

7

Java (n = 8)

import java.util.*;
import java.util.concurrent.*;

public class HankelCombinatorics {
    public static final int NUM_THREADS = 8;

    private static final int[] FACT = new int[13];
    static {
        FACT[0] = 1;
        for (int i = 1; i < FACT.length; i++) FACT[i] = i * FACT[i-1];
    }

    public static void main(String[] args) {
        long prevElapsed = 0, start = System.nanoTime();
        for (int i = 1; i < 12; i++) {
            long count = count(i), elapsed = System.nanoTime() - start;
            System.out.format("%d in %dms, total elapsed %dms\n", count, (elapsed - prevElapsed) / 1000000, elapsed / 1000000);
            prevElapsed = elapsed;
        }
    }

    @SuppressWarnings("unchecked")
    private static long count(int n) {
        int[][] perms = new int[FACT[n]][];
        genPermsInner(0, 0, new int[n], perms, 0);

        // We partition by canonical representation of the row sum multiset, discarding any with a density > 50%.
        Map<CanonicalMatrix, Map<CanonicalMatrix, Integer>> part = new HashMap<CanonicalMatrix, Map<CanonicalMatrix, Integer>>();
        for (int m = 0; m < 1 << (2*n-1); m++) {
            int density = 0;
            int[] key = new int[n];
            for (int i = 0; i < n; i++) {
                key[i] = Integer.bitCount((m >> i) & ((1 << n) - 1));
                density += key[i];
            }
            if (2 * density <= n * n) {
                CanonicalMatrix _key = new CanonicalMatrix(key);
                Map<CanonicalMatrix, Integer> map = part.get(_key);
                if (map == null) part.put(_key, map = new HashMap<CanonicalMatrix, Integer>());
                map.put(new CanonicalMatrix(m, perms[0]), m);
            }
        }

        List<Job> jobs = new ArrayList<Job>();
        ExecutorService pool = Executors.newFixedThreadPool(NUM_THREADS);

        for (Map.Entry<CanonicalMatrix, Map<CanonicalMatrix, Integer>> e : part.entrySet()) {
            Job job = new Job(n, perms, e.getKey().sum() << 1 == n * n ? 0 : 1, e.getValue());
            jobs.add(job);
            pool.execute(job);
        }

        pool.shutdown();
        try {
            pool.awaitTermination(1, TimeUnit.DAYS); // i.e. until it's finished - inaccurate results are useless
        }
        catch (InterruptedException ie) {
            throw new IllegalStateException(ie);
        }

        long total = 0;
        for (Job job : jobs) total += job.subtotal;
        return total;
    }

    private static int genPermsInner(int idx, int usedMask, int[] a, int[][] perms, int off) {
        if (idx == a.length) perms[off++] = a.clone();
        else for (int i = 0; i < a.length; i++) {
            int m = 1 << (a[idx] = i);
            if ((usedMask & m) == 0) off = genPermsInner(idx+1, usedMask | m, a, perms, off);
        }
        return off;
    }

    static class Job implements Runnable {
        private volatile long subtotal = 0;
        private final int n;
        private final int[][] perms;
        private final int shift;
        private final Map<CanonicalMatrix, Integer> unseen;

        public Job(int n, int[][] perms, int shift, Map<CanonicalMatrix, Integer> unseen) {
            this.n = n;
            this.perms = perms;
            this.shift = shift;
            this.unseen = unseen;
        }

        public void run() {
            long result = 0;
            int[][] perms = this.perms;
            Map<CanonicalMatrix, Integer> unseen = this.unseen;
            while (!unseen.isEmpty()) {
                int m = unseen.values().iterator().next();
                Set<CanonicalMatrix> equiv = new HashSet<CanonicalMatrix>();
                for (int[] perm : perms) {
                    CanonicalMatrix canonical = new CanonicalMatrix(m, perm);
                    if (equiv.add(canonical)) {
                        result += canonical.weight() << shift;
                        unseen.remove(canonical);
                    }
                }
            }

            subtotal = result;
        }
    }

    static class CanonicalMatrix {
        private final int[] a;
        private final int hash;

        public CanonicalMatrix(int m, int[] r) {
            this(permuteRows(m, r));
        }

        public CanonicalMatrix(int[] a) {
            this.a = a;
            Arrays.sort(a);

            int h = 0;
            for (int i : a) h = h * 37 + i;
            hash = h;
        }

        private static int[] permuteRows(int m, int[] perm) {
            int[] cols = new int[perm.length];
            for (int i = 0; i < perm.length; i++) {
                for (int j = 0; j < cols.length; j++) cols[j] |= ((m >> (perm[i] + j)) & 1L) << i;
            }
            return cols;
        }

        public int sum() {
            int sum = 0;
            for (int i : a) sum += i;
            return sum;
        }

        public int weight() {
            int prev = -1, count = 0, weight = FACT[a.length];
            for (int col : a) {
                if (col == prev) weight /= ++count;
                else {
                    prev = col;
                    count = 1;
                }
            }
            return weight;
        }

        @Override public boolean equals(Object obj) {
            // Deliberately unsuitable for general-purpose use, but helps catch bugs faster.
            CanonicalMatrix that = (CanonicalMatrix)obj;
            for (int i = 0; i < a.length; i++) {
                if (a[i] != that.a[i]) return false;
            }
            return true;
        }

        @Override public int hashCode() {
            return hash;
        }
    }
}

Сохранить как HankelCombinatorics.java, скомпилировать как javac HankelCombinatorics.java, запустить как java -Xmx2G HankelCombinatorics.

На NUM_THREADS = 4моей четырехъядерной машине он длится 20420819767436от n=850 до 55 секунд, с достаточной вариативностью между запусками; Я ожидаю, что он легко справится с тем же на вашей восьмиъядерной машине, но потребуется час или больше, чтобы получитьn=9 .

Как это устроено

Учитывая n, существуют 2^(2n-1)бинарные nх nганкелевы матрицы. Строки можно переставлять n!разными способами, а столбцы - n!разными способами. Все, что нам нужно сделать, это избежать двойного счета ...

Если вы вычисляете сумму каждой строки, то ни перестановка строк, ни перестановка столбцов не изменяют мультимножество сумм. Например

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

имеет многозначную сумму строк {3, 3, 2, 2, 2} , как и все матрицы Ханкеля, полученные из него. Это означает, что мы можем сгруппировать матрицы Ханкеля по этим мультимножествам суммы строк, а затем обрабатывать каждую группу независимо, используя несколько процессорных ядер.

Существует также эксплуатируемая симметрия: матрицы с большим числом нулей, чем единицы, находятся в биекции с матрицами с большим числом нулей.

Двойной учет имеет место , когда Ганкель матрица M_1с перестановкой строк r_1и перестановкой столбцов c_1соответствует матрице ганкелева M_2с перестановкой строк r_2и перестановкой столбцов c_2(до двух , но не все три M_1 = M_2, r_1 = r_2, c_1 = c_2). Строки и столбцы перестановки являются независимыми, так что если мы применим строки перестановку r_1к M_1и строке перестановку r_2в M_2столбцах как мультимножества должны быть равны. Поэтому для каждой группы я вычисляю все мультимножества столбцов, полученные путем применения перестановки строк к матрице в группе. Простой способ получить каноническое представление мультимножеств - это отсортировать столбцы, что также полезно на следующем шаге.

Получив различные мультимножества столбцов, нам нужно выяснить, сколько n!перестановок каждого из них уникально. На этом этапе двойной подсчет может иметь место только в том случае, если в заданном мультимножестве столбцов есть повторяющиеся столбцы: нам нужно подсчитать количество вхождений каждого отдельного столбца в мультимножестве, а затем вычислить соответствующий коэффициент многочлена. Поскольку столбцы отсортированы, подсчет легко выполнить.

Наконец мы добавляем их все.

Асимптотическая сложность не является тривиальной для вычисления с полной точностью, потому что мы должны сделать некоторые предположения о множествах. Мы оцениваем порядок 2^(2n-2) n!столбцов мультимножеств, принимая n^2 ln nвремя для каждого (включая сортировку); если группировка занимает не больше, чем ln nфактор, у нас есть сложность времени Theta(4^n n! n^2 ln n). Но поскольку экспоненциальные факторы полностью доминируют над полиномиальными, это так Theta(4^n n!) = Theta((4n/e)^n).

Питер Тейлор
источник
Это очень впечатляет. Не могли бы вы рассказать что-нибудь об используемом алгоритме?
3

Python2 / 3

Довольно наивный подход, на медленном языке:

import itertools

def permute_rows(m):
    for perm in itertools.permutations(m):
        yield perm

def permute_columns(m):
    T = zip(*m)
    for perm in itertools.permutations(T):
        yield zip(*perm)

N = 1
while True:
    base_template = ["abcdefghijklmnopqrstuvwxyz"[i:i+N] for i in range(N)]

    templates = set()
    for c in permute_rows(base_template):
        for m in permute_columns(c):
            templates.add("".join("".join(row) for row in m))

    def possibs(free, templates):
        if free == 2*N - 1:
            return set(int(t, 2) for t in templates)

        s = set()
        for b in "01":
            new_templates = set(t.replace("abcdefghijklmnopqrstuvwxyz"[free], b) for t in templates)
            s |= possibs(free + 1, new_templates)

        return s

    print(len(possibs(0, templates)))
    N += 1

Запустите, набрав python script.py.

orlp
источник
У вас есть язык, указанный как Python 2/3, но для того, чтобы он работал в Python 2, вам не нужен from __future__ import print_function(или что-то в этом роде)?
Алекс А.
2
@AlexA. Обычно да, но не в этом случае. Учитывайте поведение Python2 при вводе return(1). Теперь замените returnна print:)
orlp
Здорово! Я узнаю что-то новое каждый день. :)
Алекс А.
2

Haskell

import Data.List
import Data.Hashable
import Control.Parallel.Strategies
import Control.Parallel
import qualified Data.HashSet as S

main = mapM putStrLn $ map (show.countHankellable) [1..]

a§b=[a!!i|i<-b]

hashNub :: (Hashable a, Eq a) => [a] -> [a]
hashNub l = go S.empty l
    where
      go _ []     = []
      go s (x:xs) = if x `S.member` s then go s xs
                                    else x : go (S.insert x s) xs

pmap = parMap rseq

makeMatrix :: Int->[Bool]->[[Bool]]
makeMatrix n vars = [vars§[i..i+n-1]|i<-[0..n-1]]

countHankellable :: Int -> Int
countHankellable n = let
    s = permutations [0..n-1]
    conjugates m = concat[permutations[r§q|r<-m]|q<-s]
    variableSets = sequence [[True,False]|x<-[0..2*(n-1)]]
 in
    length.hashNub.concat.pmap (conjugates.makeMatrix n ) $ variableSets

Нигде не так быстро, как у Питера - это довольно впечатляющая установка! Теперь с гораздо большим количеством кода, скопированного из Интернета. Использование:

$ ghc -threaded hankell.hs
$ ./hankell
александр-Brett
источник
Ответ на Haskell всегда приветствуется. Спасибо.
@Lembik - как твои дела на твоей машине?
Александр-Бретт