Policz liczbę matryc Hankelable

12

tło

Binarna macierz Hankela to macierz o stałych przekątnych (dodatnich ukośnych przekątnych) zawierająca tylko 0s i 1s. Np. Wygląda binarna macierz Hankela 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

gdzie a, b, c, d, e, f, g, h, isą albo 0albo 1.

Zdefiniujmy macierz M jako Hankelable, jeśli istnieje permutacja rzędu rzędów i kolumn M, tak że M jest macierzą Hankela. Oznacza to, że można zastosować jedną permutację do kolejności wierszy i ewentualnie inną do kolumn.

Wyzwanie

Wyzwanie polega na policzeniu, ile jest Hankelable n według nmacierzy dla wszystkich, naż do jak największej wartości.

Wynik

Dla każdej liczby całkowitej od 1 w górę, wypisz liczbę Hankelablen według nmacierzy z wpisami, które są 0lub 1.

Na n = 1,2,3,4,5odpowiedzi powinno być 2,12,230,12076,1446672. (Dzięki orlp za kod do ich wytworzenia.)

Limit czasu

Uruchomię twój kod na moim komputerze i zatrzymam go po 1 minucie. Kod, który wyświetla poprawne odpowiedzi do największej wartości n, wygrywa. Limit czasu dotyczy wszystkiego, od n = 1największej wartości, nna którą dajesz odpowiedź.

Zwycięzca będzie najlepszą odpowiedzią do końca soboty 18 kwietnia.

Tie Breaker

W przypadku remisu na maksymalny nczas określę, ile czasu zajmie oddanie wyników, n+1a najszybszy wygra. W przypadku, gdy biegną one w tym samym czasie do sekundy n+1, wygrywa pierwsze zgłoszenie.

Języki i biblioteki

Możesz używać dowolnego języka, który ma swobodnie dostępny kompilator / tłumacz / etc. dla systemu Linux i dowolnych bibliotek, które są również bezpłatnie dostępne dla systemu Linux.

Moja maszyna

Czasy będą uruchamiane na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350 na płycie głównej Asus M5A78L-M / USB3 (Socket AM3 +, 8 GB DDR3). Oznacza to również, że muszę być w stanie uruchomić Twój kod. W związku z tym korzystaj tylko z łatwo dostępnego bezpłatnego oprogramowania i dołącz pełne instrukcje, jak skompilować i uruchomić kod.

Notatki

Odradzam iterację wszystkich macierzy n przez n i próbę wykrycia, czy każda z nich ma właściwość, którą opisuję. Po pierwsze, jest ich zbyt wiele, a po drugie, wydaje się, że nie ma szybkiego sposobu na wykrycie tego .

Dotychczasowe wpisy

  • n = 8 autorstwa Petera Taylora. Jawa
  • n = 5 przez orlp. Pyton
Społeczność
źródło
4
Naprawdę podoba mi się słowo „Hankelable”.
Alex A.,
3
W n=6sumie jest 260357434. Myślę, że presja pamięci jest większym problemem niż czas procesora.
Peter Taylor
To niesamowite pytanie. Byłem bardzo ostry jak kujon.
Alexander-Brett

Odpowiedzi:

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;
        }
    }
}

Zapisz jako HankelCombinatorics.java, skompiluj jako javac HankelCombinatorics.java, uruchom jako java -Xmx2G HankelCombinatorics.

Z NUM_THREADS = 4moją czterordzeniową maszyną trwa to 20420819767436od n=850 do 55 sekund, z dużą zmiennością między seriami; Oczekuję, że powinno to łatwo zarządzać tym samym na twoim komputerze z ośmiordzeniowym rdzeniem, ale zajmie to godzinę lub dłużej n=9.

Jak to działa

Biorąc pod uwagę n, istnieją macierze 2^(2n-1)binarne nx nHankela. Wiersze można permutować na n!różne sposoby, a kolumny na n!różne sposoby. Wszystko, co musimy zrobić, to uniknąć podwójnego liczenia ...

Jeśli obliczysz sumę każdego wiersza, to ani permutacja wierszy, ani permutacja kolumn nie zmieni multiset sum. Na przykład

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

ma {3, 3, 2, 2, 2}multiset sumy wierszy , podobnie jak wszystkie macierze Hankelable z niego pochodzące. Oznacza to, że możemy pogrupować macierze Hankela według tych multisets sumy wierszy, a następnie obsługiwać każdą grupę niezależnie, wykorzystując wiele rdzeni procesora.

Istnieje również użyteczna symetria: macierze z większą liczbą zer niż jedynki są wstrząsane z macierzami z większą liczbą zer niż zera.

Podwójnego zliczania występuje podczas Hankel matrycy M_1z rzędu permutacji r_1i kolumnę permutacji c_1odpowiada Hankel matrycę M_2z rzędu permutacji r_2i permutacji kolumnowej c_2(w liczbie do dwóch, lecz nie wszystkich trzech M_1 = M_2, r_1 = r_2, c_1 = c_2). Wiersz i kolumna permutacje są niezależne, więc jeśli zastosujemy wiersza permutacji r_1do M_1i wiersza permutacji r_2do M_2, kolumny jak multisets muszą być równe. Tak więc dla każdej grupy obliczam wszystkie wielokrotności kolumn uzyskane przez zastosowanie permutacji wiersza do macierzy w grupie. Prostym sposobem na uzyskanie kanonicznej reprezentacji multisets jest sortowanie kolumn, co jest również przydatne w następnym kroku.

Po uzyskaniu odrębnych multisets kolumnowych musimy dowiedzieć się, ile n!kombinacji każdego z nich jest unikalnych. W tym momencie podwójne liczenie może wystąpić tylko wtedy, gdy dany multiset kolumny ma zduplikowane kolumny: musimy policzyć liczbę wystąpień każdej odrębnej kolumny w multisecie, a następnie obliczyć odpowiedni współczynnik wielomianowy. Ponieważ kolumny są posortowane, łatwo jest policzyć.

Na koniec dodajemy je wszystkie.

Złożoność asymptotyczna nie jest łatwa do obliczenia z pełną precyzją, ponieważ musimy poczynić pewne założenia dotyczące zbiorów. Oceniamy według kolejności 2^(2n-2) n!multisets kolumn, zabierając n^2 ln nczas dla każdego (łącznie z sortowaniem); jeśli grupowanie nie zajmuje więcej niż ln nczynnik, mamy złożoność czasu Theta(4^n n! n^2 ln n). Ale ponieważ czynniki wykładnicze całkowicie dominują w przypadku wielomianów, tak właśnie jest Theta(4^n n!) = Theta((4n/e)^n).

Peter Taylor
źródło
To bardzo imponujące. Czy możesz powiedzieć coś o użytym algorytmie?
3

Python2 / 3

Dość naiwne podejście w wolnym języku:

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

Uruchom, wpisując python script.py.

orlp
źródło
Masz język wymieniony jako Python 2/3, ale żeby działał w Python 2, nie potrzebujesz from __future__ import print_function(czy coś takiego)?
Alex A.
2
@AlexA. Zwykle tak, ale nie w tym przypadku. Podczas pisania rozważ zachowanie Python2 return(1). Teraz zastąpić returnz print:)
orlp
Chłodny! Codziennie uczę się czegoś nowego. :)
Alex A.,
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

Nigdzie nie tak szybko jak u Petera - to imponująca konfiguracja, którą on tam ma! Teraz z większą ilością kodu skopiowanego z Internetu. Stosowanie:

$ ghc -threaded hankell.hs
$ ./hankell
Alexander-Brett
źródło
Odpowiedź Haskella jest zawsze mile widziana. Dziękuję Ci.
@Lembik - jak mi idzie na twoim komputerze?
Alexander-Brett