Wyzwanie polega na napisaniu najszybszego możliwego kodu do obliczenia stałej macierzy .

Stała n-by- nmatrix A= ( ai,j) jest zdefiniowana jako

Tutaj S_nreprezentuje zestaw wszystkich permutacji [1, n].

Jako przykład (z wiki):

W tym pytaniu macierze są kwadratowe i mają tylko wartości -1i 1na nich.

Przykłady

Wkład:

[[ 1 -1 -1  1]
 [-1 -1 -1  1]
 [-1  1 -1  1]
 [ 1 -1 -1  1]]

Stały:

-4

Wkład:

[[-1 -1 -1 -1]
 [-1  1 -1 -1]
 [ 1 -1 -1 -1]
 [ 1 -1  1 -1]]

Stały:

0

Wkład:

[[ 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 -1  1 -1  1  1  1  1]]

Stały:

192

Wkład:

[[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, 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, 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, -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, -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, -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, -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, -1]]

Stały:

1021509632

Zadanie

Powinieneś napisać kod, który podany nprzez nmacierz, wyprowadza swój stały.

Ponieważ będę musiał przetestować kod, pomocne byłoby podanie prostego sposobu na przekazanie macierzy jako danych wejściowych do kodu, na przykład poprzez czytanie ze standardowego w.

Ostrzegamy, że permanent może być duży (skrajna macierz wszystkich 1s).

Wyniki i krawaty

Przetestuję twój kod na losowych matrycach + -1 o coraz większym rozmiarze i zatrzymam się, gdy twój kod po raz pierwszy zajmie więcej niż 1 minutę na moim komputerze. Macierze oceny będą spójne dla wszystkich wniosków, aby zapewnić uczciwość.

Jeśli dwie osoby otrzymają ten sam wynik, zwycięzcą jest ta, która jest najszybsza dla tej wartości n. Jeśli są one w odległości 1 sekundy od siebie, to jest to pierwszy.

Języki i biblioteki

Możesz użyć dowolnego dostępnego języka i bibliotek, które ci się podobają, ale nie ma wcześniejszej funkcji do obliczenia stałej. Tam, gdzie jest to wykonalne, dobrze byłoby móc uruchomić kod, więc proszę podać pełne wyjaśnienie, jak uruchomić / skompilować kod w systemie Linux, jeśli to w ogóle możliwe.

Referencyjne implementacje

Jest już pytanie o kodegolfa z dużą ilością kodu w różnych językach do obliczania stałej dla małych matryc. Zarówno Mathematica, jak i Maple mają również stałe implementacje, jeśli masz do nich dostęp.

Moja maszyna Czasy będą działać na moim 64-bitowym komputerze. Jest to standardowa instalacja ubuntu z 8 GB pamięci RAM, ośmiordzeniowym procesorem AMD FX-8350 i Radeon HD 4250. Oznacza to również, że muszę mieć możliwość uruchomienia kodu.

Informacje niskiego poziomu o mojej maszynie

cat /proc/cpuinfo/|grep flags daje

flagi: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl fqqm_sfs_sf_sfs_sfs f16c lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs xop skinit wdt lwp fma4 tce nodeid_msr tbm topoext perfctr_core perfctr_nb cpb hw_psts vmcall vmmcall bmw

Zadam ściśle powiązane, wielojęzyczne pytanie, które nie cierpi z powodu dużego problemu Int, więc miłośnicy Scali , Nima , Julii , Rusty i Basha mogą również pochwalić się swoimi językami.

Tablica liderów

• n = 33 (45 sekund. 64 sekundy dla n = 34). Ton Hospel w C ++ z g ++ 5.4.0.
• n = 32 (32 sekundy). Dennis w C z gcc 5.4.0 przy użyciu flag gcc Ton Hospel.
• n = 31 (54 sekund). Christian Sievers w Haskell
• n = 31 (60 sekund). primo w rpython
• n = 30 (26 sekund). ezrast in Rust
• n = 28 (49 sekund). xnor z Python + pypy 5.4.1
• n = 22 (25 sekund). Shebang z Python + pypy 5.4.1

Uwaga . W praktyce terminy Dennisa i Tona Hospela są bardzo różne z tajemniczych powodów. Na przykład wydają się być szybsze po załadowaniu przeglądarki internetowej! Podane czasy są najszybsze we wszystkich testach, które przeprowadziłem.

gcc C ++ n ≈ 36 (57 sekund w moim systemie)

Używa formuły Glynn z szarym kodem do aktualizacji, jeśli wszystkie sumy kolumn są parzyste, w przeciwnym razie używa metody Rysera. Gwintowane i wektoryzowane. Zoptymalizowany dla AVX, więc nie oczekuj dużo na starszych procesorach. Nie przejmuj n>=35się matrycą zawierającą tylko +1, nawet jeśli twój system jest wystarczająco szybki, ponieważ podpisany 128-bitowy akumulator się przepełni. W przypadku macierzy losowych prawdopodobnie nie trafisz do przepełnienia. Dla n>=37wewnętrznych mnożników zacznie się przepełniać dla całej 1/-1macierzy. Dlatego używaj tego programu tylko do n<=36.

Po prostu podaj elementy macierzy na STDIN oddzielone dowolnym rodzajem białych znaków

permanent
1 2
3 4
^D

permanent.cpp:

/*
  Compile using something like:
    g++ -Wall -O3 -march=native -fstrict-aliasing -std=c++11 -pthread -s permanent.cpp -o permanent
*/

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cstdint>
#include <climits>
#include <array>
#include <vector>
#include <thread>
#include <future>
#include <ctgmath>
#include <immintrin.h>

using namespace std;

bool const DEBUG = false;
int const CACHE = 64;

using Index  = int_fast32_t;
Index glynn;
// Number of elements in our vectors
Index const POW   = 3;
Index const ELEMS = 1 << POW;
// Over how many floats we distribute each row
Index const WIDTH = 9;
// Number of bits in the fraction part of a floating point number
int const FLOAT_MANTISSA = 23;
// Type to use for the first add/multiply phase
using Sum  = float;
using SumN = __restrict__ Sum __attribute__((vector_size(ELEMS*sizeof(Sum))));
// Type to convert to between the first and second phase
using ProdN = __restrict__ int32_t __attribute__((vector_size(ELEMS*sizeof(int32_t))));
// Type to use for the third and last multiply phase.
// Also used for the final accumulator
using Value = __int128;
using UValue = unsigned __int128;

// Wrap Value so C++ doesn't really see it and we can put it in vectors etc.
// Needed since C++ doesn't fully support __int128
struct Number {
    Number& operator+=(Number const& right) {
        value += right.value;
        return *this;
    }
    // Output the value
    void print(ostream& os, bool dbl = false) const;
    friend ostream& operator<<(ostream& os, Number const& number) {
        number.print(os);
        return os;
    }

    Value value;
};

using ms = chrono::milliseconds;

auto nr_threads = thread::hardware_concurrency();
vector<Sum> input;

// Allocate cache aligned datastructures
template<typename T>
T* alloc(size_t n) {
    T* mem = static_cast<T*>(aligned_alloc(CACHE, sizeof(T) * n));
    if (mem == nullptr) throw(bad_alloc());
    return mem;
}

// Work assigned to thread k of nr_threads threads
Number permanent_part(Index n, Index k, SumN** more) {
    uint64_t loops = (UINT64_C(1) << n) / nr_threads;
    if (glynn) loops /= 2;
    Index l = loops < ELEMS ? loops : ELEMS;
    loops /= l;
    auto from = loops * k;
    auto to   = loops * (k+1);

    if (DEBUG) cout << "From=" << from << "\n";
    uint64_t old_gray = from ^ from/2;
    uint64_t bit = 1;
    bool bits = (to-from) & 1;

    Index nn = (n+WIDTH-1)/WIDTH;
    Index ww = nn * WIDTH;
    auto column = alloc<SumN>(ww);
    for (Index i=0; i<n; ++i)
        for (Index j=0; j<ELEMS; ++j) column[i][j] = 0;
    for (Index i=n; i<ww; ++i)
        for (Index j=0; j<ELEMS; ++j) column[i][j] = 1;
    Index b;
    if (glynn) {
        b = n > POW+1 ? n - POW - 1: 0;
        auto c = n-1-b;
        for (Index k=0; k<l; k++) {
            Index gray = k ^ k/2;
            for (Index j=0; j< c; ++j)
                if (gray & 1 << j)
                    for (Index i=0; i<n; ++i)
                        column[i][k] -= input[(b+j)*n+i];
                else
                    for (Index i=0; i<n; ++i)
                        column[i][k] += input[(b+j)*n+i];
        }
        for (Index i=0; i<n; ++i)
            for (Index k=0; k<l; k++)
                column[i][k] += input[n*(n-1)+i];

        for (Index k=1; k<l; k+=2)
            column[0][k] = -column[0][k];

        for (Index i=0; i<b; ++i, bit <<= 1) {
            if (old_gray & bit) {
                bits = bits ^ 1;
                for (Index j=0; j<ww; ++j)
                    column[j] -= more[i][j];
            } else {
                for (Index j=0; j<ww; ++j)
                    column[j] += more[i][j];
            }
        }

        for (Index i=0; i<n; ++i)
            for (Index k=0; k<l; k++)
                column[i][k] /= 2;
    } else {
        b = n > POW ? n - POW : 0;
        auto c = n-b;
        for (Index k=0; k<l; k++) {
            Index gray = k ^ k/2;
            for (Index j=0; j<c; ++j)
                if (gray & 1 << j)
                    for (Index i=0; i<n; ++i)
                        column[i][k] -= input[(b+j)*n+i];
        }

        for (Index k=1; k<l; k+=2)
            column[0][k] = -column[0][k];

        for (Index i=0; i<b; ++i, bit <<= 1) {
            if (old_gray & bit) {
                bits = bits ^ 1;
                for (Index j=0; j<ww; ++j)
                    column[j] -= more[i][j];
            }
        }
    }

    if (DEBUG) {
        for (Index i=0; i<ww; ++i) {
            cout << "Column[" << i << "]=";
            for (Index j=0; j<ELEMS; ++j) cout << " " << column[i][j];
            cout << "\n";
        }
    }

    --more;
    old_gray = (from ^ from/2) | UINT64_C(1) << b;
    Value total = 0;
    SumN accu[WIDTH];
    for (auto p=from; p<to; ++p) {
        uint64_t new_gray = p ^ p/2;
        uint64_t bit = old_gray ^ new_gray;
        Index i = __builtin_ffsl(bit);
        auto diff = more[i];
        auto c = column;
        if (new_gray > old_gray) {
            // Phase 1 add/multiply.
            // Uses floats until just before loss of precision
            for (Index i=0; i<WIDTH; ++i) accu[i] = *c++ -= *diff++;

            for (Index j=1; j < nn; ++j)
                for (Index i=0; i<WIDTH; ++i) accu[i] *= *c++ -= *diff++;
        } else {
            // Phase 1 add/multiply.
            // Uses floats until just before loss of precision
            for (Index i=0; i<WIDTH; ++i) accu[i] = *c++ += *diff++;

            for (Index j=1; j < nn; ++j)
                for (Index i=0; i<WIDTH; ++i) accu[i] *= *c++ += *diff++;
        }

        if (DEBUG) {
            cout << "p=" << p << "\n";
            for (Index i=0; i<ww; ++i) {
                cout << "Column[" << i << "]=";
                for (Index j=0; j<ELEMS; ++j) cout << " " << column[i][j];
                cout << "\n";
            }
        }

        // Convert floats to int32_t
        ProdN prod32[WIDTH] __attribute__((aligned (32)));
        for (Index i=0; i<WIDTH; ++i)
            // Unfortunately gcc doesn't recognize the static_cast<int32_t>
            // as a vector pattern, so force it with an intrinsic
#ifdef __AVX__
            //prod32[i] = static_cast<ProdN>(accu[i]);
            reinterpret_cast<__m256i&>(prod32[i]) = _mm256_cvttps_epi32(accu[i]);
#else   // __AVX__
            for (Index j=0; j<ELEMS; ++j)
                prod32[i][j] = static_cast<int32_t>(accu[i][j]);
#endif  // __AVX__

        // Phase 2 multiply. Uses int64_t until just before overflow
        int64_t prod64[3][ELEMS];
        for (Index i=0; i<3; ++i) {
            for (Index j=0; j<ELEMS; ++j)
                prod64[i][j] = static_cast<int64_t>(prod32[i][j]) * prod32[i+3][j] * prod32[i+6][j];
        }
        // Phase 3 multiply. Collect into __int128. For large matrices this will
        // actually overflow but that's ok as long as all 128 low bits are
        // correct. Terms will cancel and the final sum can fit into 128 bits
        // (This will start to fail at n=35 for the all 1 matrix)
        // Strictly speaking this needs the -fwrapv gcc option
        for (Index j=0; j<ELEMS; ++j) {
            auto value = static_cast<Value>(prod64[0][j]) * prod64[1][j] * prod64[2][j];
            if (DEBUG) cout << "value[" << j << "]=" << static_cast<double>(value) << "\n";
            total += value;
        }
        total = -total;

        old_gray = new_gray;
    }

    return bits ? Number{-total} : Number{total};
}

// Prepare datastructures, Assign work to threads
Number permanent(Index n) {
    Index nn = (n+WIDTH-1)/WIDTH;
    Index ww = nn*WIDTH;

    Index rows  = n > (POW+glynn) ? n-POW-glynn : 0;
    auto data = alloc<SumN>(ww*(rows+1));
    auto pointers = alloc<SumN *>(rows+1);
    auto more = &pointers[0];
    for (Index i=0; i<rows; ++i)
        more[i] = &data[ww*i];
    more[rows] = &data[ww*rows];
    for (Index j=0; j<ww; ++j)
        for (Index i=0; i<ELEMS; ++i)
            more[rows][j][i] = 0;

    Index loops = n >= POW+glynn ? ELEMS : 1 << (n-glynn);
    auto a = &input[0];
    for (Index r=0; r<rows; ++r) {
        for (Index j=0; j<n; ++j) {
            for (Index i=0; i<loops; ++i)
                more[r][j][i] = j == 0 && i %2 ? -*a : *a;
            for (Index i=loops; i<ELEMS; ++i)
                more[r][j][i] = 0;
            ++a;
        }
        for (Index j=n; j<ww; ++j)
            for (Index i=0; i<ELEMS; ++i)
                more[r][j][i] = 0;
    }

    if (DEBUG)
        for (Index r=0; r<=rows; ++r)
            for (Index j=0; j<ww; ++j) {
                cout << "more[" << r << "][" << j << "]=";
                for (Index i=0; i<ELEMS; ++i)
                    cout << " " << more[r][j][i];
                cout << "\n";
            }

    // Send work to threads...
    vector<future<Number>> results;
    for (auto i=1U; i < nr_threads; ++i)
        results.emplace_back(async(DEBUG ? launch::deferred: launch::async, permanent_part, n, i, more));
    // And collect results
    auto r = permanent_part(n, 0, more);
    for (auto& result: results)
        r += result.get();

    free(data);
    free(pointers);

    // For glynn we should double the result, but we will only do this during
    // the final print. This allows n=34 for an all 1 matrix to work
    // if (glynn) r *= 2;
    return r;
}

// Print 128 bit number
void Number::print(ostream& os, bool dbl) const {
    const UValue BILLION = 1000000000;

    UValue val;
    if (value < 0) {
        os << "-";
        val = -value;
    } else
        val = value;
    if (dbl) val *= 2;

    uint32_t output[5];
    for (int i=0; i<5; ++i) {
        output[i] = val % BILLION;
        val /= BILLION;
    }
    bool print = false;
    for (int i=4; i>=0; --i) {
        if (print) {
            os << setfill('0') << setw(9) << output[i];
        } else if (output[i] || i == 0) {
            print = true;
            os << output[i];
        }
    }
}

// Read matrix, check for sanity
void my_main() {
    Sum a;
    while (cin >> a)
        input.push_back(a);

    size_t n = sqrt(input.size());
    if (input.size() != n*n)
        throw(logic_error("Read " + to_string(input.size()) +
                          " elements which does not make a square matrix"));

    vector<double> columns_pos(n, 0);
    vector<double> columns_neg(n, 0);
    Sum *p = &input[0];
    for (size_t i=0; i<n; ++i)
        for (size_t j=0; j<n; ++j, ++p) {
            if (*p >= 0) columns_pos[j] += *p;
            else         columns_neg[j] -= *p;
        }
    std::array<double,WIDTH> prod;
    prod.fill(1);

    int32_t odd = 0;
    for (size_t j=0; j<n; ++j) {
        prod[j%WIDTH] *= max(columns_pos[j], columns_neg[j]);
        auto sum = static_cast<int32_t>(columns_pos[j] - columns_neg[j]);
        odd |= sum;
    }
    glynn = (odd & 1) ^ 1;
    for (Index i=0; i<WIDTH; ++i)
        // A float has an implicit 1. in front of the fraction so it can
        // represent 1 bit more than the mantissa size. And 1 << (mantissa+1)
        // itself is in fact representable
        if (prod[i] && log2(prod[i]) > FLOAT_MANTISSA+1)
            throw(range_error("Values in matrix are too large. A subproduct reaches " + to_string(prod[i]) + " which doesn't fit in a float without loss of precision"));

    for (Index i=0; i<3; ++i) {
        auto prod3 = prod[i] * prod[i+3] * prod[i+6];
        if (log2(prod3) >= CHAR_BIT*sizeof(int64_t)-1)
            throw(range_error("Values in matrix are too large. A subproduct reaches " + to_string(prod3) + " which doesn't fit in an int64"));
    }

    nr_threads = pow(2, ceil(log2(static_cast<float>(nr_threads))));
    uint64_t loops = UINT64_C(1) << n;
    if (glynn) loops /= 2;
    if (nr_threads * ELEMS > loops)
        nr_threads = max(loops / ELEMS, UINT64_C(1));
    // if (DEBUG) nr_threads = 1;

    cout << n << " x " << n << " matrix, method " << (glynn ? "Glynn" : "Ryser") << ", " << nr_threads << " threads" << endl;

    // Go for the actual calculation
    auto start = chrono::steady_clock::now();
    auto perm = permanent(n);
    auto end = chrono::steady_clock::now();
    auto elapsed = chrono::duration_cast<ms>(end-start).count();

    cout << "Permanent=";
    perm.print(cout, glynn);
    cout << " (" << elapsed / 1000. << " s)" << endl;
}

// Wrapper to print any exceptions
int main() {
    try {
        my_main();
    } catch(exception& e) {
        cerr << "Error: " << e.what() << endl;
        exit(EXIT_FAILURE);
    }
    exit(EXIT_SUCCESS);
}

C99, n ≈ 33 (35 sekund)

#include <stdint.h>
#include <stdio.h>

#define CHUNK_SIZE 12
#define NUM_THREADS 8

#define popcnt __builtin_popcountll
#define BILLION (1000 * 1000 * 1000)
#define UPDATE_ROW_PPROD() \
    update_row_pprod(row_pprod, row, rows, row_sums, mask, mask_popcnt)

typedef __int128 int128_t;

static inline int64_t update_row_pprod
(
    int64_t* row_pprod, int64_t row, int64_t* rows,
    int64_t* row_sums, int64_t mask, int64_t mask_popcnt
)
{
    int64_t temp = 2 * popcnt(rows[row] & mask) - mask_popcnt;

    row_pprod[0] *= temp;
    temp -= 1;
    row_pprod[1] *= temp;
    temp -= row_sums[row];
    row_pprod[2] *= temp;
    temp += 1;
    row_pprod[3] *= temp;

    return row + 1;
}

int main(int argc, char* argv[])
{
    int64_t size = argc - 1, rows[argc - 1];
    int64_t row_sums[argc - 1];
    int128_t permanent = 0, sign = size & 1 ? -1 : 1;

    if (argc == 2)
    {
        printf("%d\n", argv[1][0] == '-' ? -1 : 1);
        return 0;
    }

    for (int64_t row = 0; row < size; row++)
    {
        char positive = argv[row + 1][0] == '+' ? '-' : '+';

        sign *= ',' - positive;
        rows[row] = row_sums[row] = 0;

        for (char* p = &argv[row + 1][1]; *p; p++)
        {
            rows[row] <<= 1;
            rows[row] |= *p == positive;
            row_sums[row] += *p == positive;
        }

        row_sums[row] = 2 * row_sums[row] - size;
    }

    #pragma omp parallel for reduction(+:permanent) num_threads(NUM_THREADS)
    for (int64_t mask = 1; mask < 1LL << (size - 1); mask += 2)
    {
        int64_t mask_popcnt = popcnt(mask);
        int64_t row = 0;
        int128_t row_prod = 1 - 2 * (mask_popcnt & 1);
        int128_t row_prod_high = -row_prod;
        int128_t row_prod_inv = row_prod;
        int128_t row_prod_inv_high = -row_prod;

        for (int64_t chunk = 0; chunk < size / CHUNK_SIZE; chunk++)
        {
            int64_t row_pprod[4] = {1, 1, 1, 1};

            for (int64_t i = 0; i < CHUNK_SIZE; i++)
                row = UPDATE_ROW_PPROD();

            row_prod *= row_pprod[0], row_prod_high *= row_pprod[1];
            row_prod_inv *= row_pprod[3], row_prod_inv_high *= row_pprod[2];
        }

        int64_t row_pprod[4] = {1, 1, 1, 1};

        while (row < size)
            row = UPDATE_ROW_PPROD();

        row_prod *= row_pprod[0], row_prod_high *= row_pprod[1];
        row_prod_inv *= row_pprod[3], row_prod_inv_high *= row_pprod[2];
        permanent += row_prod + row_prod_high + row_prod_inv + row_prod_inv_high;
    }

    permanent *= sign;

    if (permanent < 0)
        printf("-"), permanent *= -1;

    int32_t output[5], print = 0;

    output[0] = permanent % BILLION, permanent /= BILLION;
    output[1] = permanent % BILLION, permanent /= BILLION;
    output[2] = permanent % BILLION, permanent /= BILLION;
    output[3] = permanent % BILLION, permanent /= BILLION;
    output[4] = permanent % BILLION;

    if (output[4])
        printf("%u", output[4]), print = 1;
    if (print)
        printf("%09u", output[3]);
    else if (output[3])
        printf("%u", output[3]), print = 1;
    if (print)
        printf("%09u", output[2]);
    else if (output[2])
        printf("%u", output[2]), print = 1;
    if (print)
        printf("%09u", output[1]);
    else if (output[1])
        printf("%u", output[1]), print = 1;
    if (print)
        printf("%09u\n", output[0]);
    else
        printf("%u\n", output[0]);
}

Wejście jest obecnie nieco uciążliwe; jest pobierany z wierszami jako argumentami wiersza poleceń, gdzie każdy wpis jest reprezentowany przez jego znak, tj. + oznacza 1, a - oznacza -1 .

Testowe uruchomienie

$ gcc -Wall -std=c99 -march=native -Ofast -fopenmp -fwrapv -o permanent permanent.c
$ ./permanent +--+ ---+ -+-+ +--+
-4
$ ./permanent ---- -+-- +--- +-+-
0
$ ./permanent +-+----- --++-++- +----+++ ---+-+++ +--++++- -+-+-++- +-+-+-+- --+-++++
192
$ ./permanent +-+--+++----++++-++- +-+++++-+--+++--+++- --+++----+-+++---+-- ---++-++++++------+- -+++-+++---+-+-+++++ +-++--+-++++-++-+--- +--+---+-++++---+++- +--+-++-+++-+-+++-++ +-----+++-----++-++- --+-+-++-+-++++++-++ -------+----++++---- ++---++--+-++-++++++ -++-----++++-----+-+ ++---+-+----+-++-+-+ +++++---+++-+-+++-++ +--+----+--++-+----- -+++-++--+++--++--++ ++--++-++-+++-++-+-+ +++---+--++---+----+ -+++-------++-++-+--
1021509632
$ time ./permanent +++++++++++++++++++++++++++++++{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,}     # 31
8222838654177922817725562880000000

real    0m8.365s
user    1m6.504s
sys     0m0.000s
$ time ./permanent ++++++++++++++++++++++++++++++++{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,}   # 32
263130836933693530167218012160000000

real    0m17.013s
user    2m15.226s
sys     0m0.001s
$ time ./permanent +++++++++++++++++++++++++++++++++{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,} # 33
8683317618811886495518194401280000000

real    0m34.592s
user    4m35.354s
sys     0m0.001s

Python 2, n ≈ 28

from operator import mul

def fast_glynn_perm(M):
    row_comb = [sum(c) for c in zip(*M)]
    n=len(M)

    total = 0
    old_grey = 0 
    sign = +1

    binary_power_dict = {2**i:i for i in range(n)}
    num_loops = 2**(n-1)

    for bin_index in xrange(1, num_loops + 1):  
        total += sign * reduce(mul, row_comb)

        new_grey = bin_index^(bin_index/2)
        grey_diff = old_grey ^ new_grey
        grey_diff_index = binary_power_dict[grey_diff]

        new_vector = M[grey_diff_index]
        direction = 2 * cmp(old_grey,new_grey)      

        for i in range(n):
            row_comb[i] += new_vector[i] * direction

        sign = -sign
        old_grey = new_grey

    return total/num_loops

Używa formuły Glynn z szarym kodem do aktualizacji. n=23Za chwilę podbiega do mojej maszyny. Z pewnością można to zrobić lepiej, stosując szybszy język i lepsze struktury danych. Nie korzysta to z tego, że macierz ma wartość ± 1.

Implementacja formuły Ryser jest bardzo podobna, sumując wszystkie wektory 0/1 współczynników zamiast ± 1 wektorów. To zajmuje około dwa razy więcej niż formuła Glynna, ponieważ dodaje ponad 2 ^ n takich wektorów, podczas gdy Glynn dzieli połowę tego, że używa symetrii tylko na te, które zaczynają się od +1.

from operator import mul

def fast_ryser_perm(M):
    n=len(M)
    row_comb = [0] * n

    total = 0
    old_grey = 0 
    sign = +1

    binary_power_dict = {2**i:i for i in range(n)}
    num_loops = 2**n

    for bin_index in range(1, num_loops) + [0]: 
        total += sign * reduce(mul, row_comb)

        new_grey = bin_index^(bin_index/2)
        grey_diff = old_grey ^ new_grey
        grey_diff_index = binary_power_dict[grey_diff]

        new_vector = M[grey_diff_index]
        direction = cmp(old_grey, new_grey)

        for i in range(n):
            row_comb[i] += new_vector[i] * direction

        sign = -sign
        old_grey = new_grey

    return total * (-1)**n

Haskell, n = 31 (54s)

Z wieloma nieocenionymi wkładami @Angs: używaj Vector, używaj produktów zwarciowych, spójrz na nieparzyste n.

import Control.Parallel.Strategies
import qualified Data.Vector.Unboxed as V
import Data.Int

type Row = V.Vector Int8

x :: Row -> [Row] -> Integer -> Int -> Integer
x p (v:vs) m c = let c' = c - 1
                     r = if c>0 then parTuple2 rseq rseq else r0
                     (a,b) = ( x p                  vs m    c' ,
                               x (V.zipWith(-) p v) vs (-m) c' )
                             `using` r
                 in a+b
x p _      m _ = prod m p

prod :: Integer -> Row -> Integer
prod m p = if 0 `V.elem` p then 0 
                           else V.foldl' (\a b->a*fromIntegral b) m p

p, pt :: [Row] -> Integer
p (v:vs) = x (foldl (V.zipWith (+)) v vs) (map (V.map (2*)) vs) 1 11
           `div` 2^(length vs)
p [] = 1 -- handle 0x0 matrices too  :-)

pt (v:vs) | even (length vs) = p ((V.map (2*) v) : vs ) `div` 2
pt mat                       = p mat

main = getContents >>= print . pt . map V.fromList . read

Moje pierwsze próby paralelizmu w Haskell. W historii zmian można zobaczyć wiele kroków optymalizacji. O dziwo, były to głównie bardzo małe zmiany. Kod oparty jest na formule w sekcji „Formuła Balasubramanian-Bax / Franklin-Glynn” w artykule na temat obliczania stałego .

poblicza stałą. Nazywa się to, za pomocą ptktórego przekształca macierz w sposób, który jest zawsze poprawny, ale szczególnie przydatny w przypadku macierzy, które tu otrzymujemy.

Kompiluj z ghc -O2 -threaded -fllvm -feager-blackholing -o <name> <name>.hs. Aby uruchomić z parallelisation, dać to runtime parametry jak poniżej: ./<name> +RTS -N. Dane wejściowe pochodzą ze standardowego wejścia z zagnieżdżonymi listami oddzielonymi przecinkami w nawiasach, tak [[1,2],[3,4]]jak w poprzednim przykładzie (znaki nowej linii są dozwolone wszędzie).

Rdza + extprim

Ten prosty Ryser z implementacją Graya zajmuje około 65 90 sekund, aby uruchomić n = 31 na moim laptopie. Wyobrażam sobie, że twoja maszyna dotrze tam znacznie poniżej 60 lat. Używam extprim 1.1.1 dla i128.

Nigdy nie korzystałem z Rust i nie mam pojęcia, co robię. Żadnych opcji kompilatora innych niż cokolwiek cargo build --release. Komentarze / sugestie / optymalizacje są mile widziane.

Inwokacja jest identyczna z programem Dennisa.

use std::env;
use std::thread;
use std::sync::Arc;
use std::sync::mpsc;

extern crate extprim;
use extprim::i128::i128;

static THREADS : i64 = 8; // keep this a power of 2.

fn main() {
  // Read command line args for the matrix, specified like
  // "++- --- -+-" for [[1, 1, -1], [-1, -1, -1], [-1, 1, -1]].
  let mut args = env::args();
  args.next();

  let mat : Arc<Vec<Vec<i64>>> = Arc::new(args.map( |ss|
    ss.trim().bytes().map( |cc| if cc == b'+' {1} else {-1}).collect()
  ).collect());

  // Figure how many iterations each thread has to do.
  let size = 2i64.pow(mat.len() as u32);
  let slice_size = size / THREADS; // Assumes divisibility.

  let mut accumulator : i128;
  if slice_size >= 4 { // permanent() requires 4 divides slice_size
    let (tx, rx) = mpsc::channel();

    // Launch threads.
    for ii in 0..THREADS {
      let mat = mat.clone();
      let tx = tx.clone();
      thread::spawn(move ||
        tx.send(permanent(&mat, ii * slice_size, (ii+1) * slice_size))
      );
    }

    // Accumulate results.
    accumulator = extprim::i128::ZERO;
    for _ in 0..THREADS {
      accumulator += rx.recv().unwrap();
    }
  }
  else { // Small matrix, don't bother threading.
    accumulator = permanent(&mat, 0, size);
  }
  println!("{}", accumulator);
}

fn permanent(mat: &Vec<Vec<i64>>, start: i64, end: i64) -> i128 {
  let size = mat.len();
  let sentinel = std::i64::MAX / size as i64;

  let mut bits : Vec<bool> = Vec::with_capacity(size);
  let mut sums : Vec<i64> = Vec::with_capacity(size);

  // Initialize gray code bits.
  let gray_number = start ^ (start / 2);

  for row in 0..size {
    bits.push((gray_number >> row) % 2 == 1);
    sums.push(0);
  }

  // Initialize column sums
  for row in 0..size {
    if bits[row] {
      for column in 0..size {
        sums[column] += mat[row][column];
      }
    }
  }

  // Do first two iterations with initial sums
  let mut total = product(&sums, sentinel);
  for column in 0..size {
    sums[column] += mat[0][column];
  }
  bits[0] = true;

  total -= product(&sums, sentinel);

  // Do rest of iterations updating gray code bits incrementally
  let mut gray_bit : usize;
  let mut idx = start + 2;
  while idx < end {
    gray_bit = idx.trailing_zeros() as usize;

    if bits[gray_bit] {
      for column in 0..size {
        sums[column] -= mat[gray_bit][column];
      }
      bits[gray_bit] = false;
    }
    else {
      for column in 0..size {
        sums[column] += mat[gray_bit][column];
      }
      bits[gray_bit] = true;
    }

    total += product(&sums, sentinel);

    if bits[0] {
      for column in 0..size {
        sums[column] -= mat[0][column];
      }
      bits[0] = false;
    }
    else {
      for column in 0..size {
        sums[column] += mat[0][column];
      }
      bits[0] = true;
    }

    total -= product(&sums, sentinel);
    idx += 2;
  }
  return if size % 2 == 0 {total} else {-total};
}

#[inline]
fn product(sums : &Vec<i64>, sentinel : i64) -> i128 {
  let mut ret : Option<i128> = None;
  let mut tally = sums[0];
  for ii in 1..sums.len() {
    if tally.abs() >= sentinel {
      ret = Some(ret.map_or(i128::new(tally), |n| n * i128::new(tally)));
      tally = sums[ii];
    }
    else {
      tally *= sums[ii];
    }
  }
  if ret.is_none() {
    return i128::new(tally);
  }
  return ret.unwrap() * i128::new(tally);
}

Mathematica, n ≈ 20

p[m_] := Last[Fold[Take[ListConvolve[##, {1, -1}, 0], 2^Length[m]]&,
  Table[If[IntegerQ[Log2[k]], m[[j, Log2[k] + 1]], 0], {j, n}, {k, 0, 2^Length[m] - 1}]]]

Za pomocą Timingpolecenia macierz 20x20 wymaga około 48 sekund w moim systemie. Nie jest to tak wydajne jak inne, ponieważ opiera się na fakcie, że stałe można znaleźć jako współczynnik iloczynu wielomianów z każdego rzędu macierzy. Efektywne mnożenie wielomianowe odbywa się poprzez tworzenie list współczynników i wykonywanie splotów za pomocą ListConvolve. Wymaga to około O (2 ⁿ n ² ) założenia, że ​​splot jest wykonywany przy użyciu szybkiej transformaty Fouriera lub podobnej, która wymaga czasu O ( n log n ).

Python 2, n = 22 [Odwołanie]

Jest to implementacja „referencyjna”, którą wczoraj podzieliłem z Lembikiem, brakuje jej n=23wykonania o kilka sekund na jego komputerze, na mojej maszynie robi to w około 52 sekund. Aby osiągnąć te prędkości, musisz uruchomić to przez PyPy.

Pierwsza funkcja oblicza stałą podobnie do sposobu wyznaczania wyznacznika, przechodząc przez każdą podwielokrotkę, aż pozostanie 2x2, do którego można zastosować podstawową regułę. Jest niesamowicie wolny .

Druga funkcja to ta implementująca funkcję Rysera (drugie równanie wymienione w Wikipedii). Zestaw Sjest zasadniczo zestawem liczb {1,...,n}(zmiennym s_listw kodzie).

from random import *
from time import time
from itertools import*

def perm(a): # naive method, recurses over submatrices, slow 
    if len(a) == 1:
        return a[0][0]
    elif len(a) == 2:
        return a[0][0]*a[1][1]+a[1][0]*a[0][1]
    else:
        tsum = 0
        for i in xrange(len(a)):
            transposed = [zip(*a)[j] for j in xrange(len(a)) if j != i]
            tsum += a[0][i] * perm(zip(*transposed)[1:])
        return tsum

def perm_ryser(a): # Ryser's formula, using matrix entries
    maxn = len(a)
    n_list = range(1,maxn+1)
    s_list = chain.from_iterable(combinations(n_list,i) for i in range(maxn+1))
    total = 0
    for st in s_list:
        stotal = (-1)**len(st)
        for i in xrange(maxn):
            stotal *= sum(a[i][j-1] for j in st)
        total += stotal
    return total*((-1)**maxn)


def genmatrix(d):
    mat = []
    for x in xrange(d):
        row = []
        for y in xrange(d):
            row.append([-1,1][randrange(0,2)])
        mat.append(row)
    return mat

def main():
    for i in xrange(1,24):
        k = genmatrix(i)
        print 'Matrix: (%dx%d)'%(i,i)
        print '\n'.join('['+', '.join(`j`.rjust(2) for j in a)+']' for a in k)
        print 'Permanent:',
        t = time()
        p = perm_ryser(k)
        print p,'(took',time()-t,'seconds)'

if __name__ == '__main__':
    main()

RPython 5.4.1, n ≈ 32 (37 sekund)

from rpython.rlib.rtime import time
from rpython.rlib.rarithmetic import r_int, r_uint
from rpython.rlib.rrandom import Random
from rpython.rlib.rposix import pipe, close, read, write, fork, waitpid
from rpython.rlib.rbigint import rbigint

from math import log, ceil
from struct import pack

bitsize = len(pack('l', 1)) * 8 - 1

bitcounts = bytearray([0])
for i in range(16):
  b = bytearray([j+1 for j in bitcounts])
  bitcounts += b


def bitcount(n):
  bits = 0
  while n:
    bits += bitcounts[n & 65535]
    n >>= 16
  return bits


def main(argv):
  if len(argv) < 2:
    write(2, 'Usage: %s NUM_THREADS [N]'%argv[0])
    return 1
  threads = int(argv[1])

  if len(argv) > 2:
    n = int(argv[2])
    rnd = Random(r_uint(time()*1000))
    m = []
    for i in range(n):
      row = []
      for j in range(n):
        row.append(1 - r_int(rnd.genrand32() & 2))
      m.append(row)
  else:
    m = []
    strm = ""
    while True:
      buf = read(0, 4096)
      if len(buf) == 0:
        break
      strm += buf
    rows = strm.split("\n")
    for row in rows:
      r = []
      for val in row.split(' '):
        r.append(int(val))
      m.append(r)
    n = len(m)

  a = []
  for row in m:
    val = 0
    for v in row:
      val = (val << 1) | -(v >> 1)
    a.append(val)

  batches = int(ceil(n * log(n) / (bitsize * log(2))))

  pids = []
  handles = []
  total = rbigint.fromint(0)
  for i in range(threads):
    r, w = pipe()
    pid = fork()
    if pid:
      close(w)
      pids.append(pid)
      handles.append(r)
    else:
      close(r)
      total = run(n, a, i, threads, batches)
      write(w, total.str())
      close(w)
      return 0

  for pid in pids:
    waitpid(pid, 0)

  for handle in handles:
    strval = read(handle, 256)
    total = total.add(rbigint.fromdecimalstr(strval))
    close(handle)

  print total.rshift(n-1).str()

  return 0


def run(n, a, mynum, threads, batches):
  start = (1 << n-1) * mynum / threads
  end = (1 << n-1) * (mynum+1) / threads

  dtotal = rbigint.fromint(0)
  for delta in range(start, end):
    pdelta = rbigint.fromint(1 - ((bitcount(delta) & 1) << 1))
    for i in range(batches):
      pbatch = 1
      for j in range(i, n, batches):
        pbatch *= n - (bitcount(delta ^ a[j]) << 1)
      pdelta = pdelta.int_mul(pbatch)
    dtotal = dtotal.add(pdelta)

  return dtotal


def target(*args):
  return main

Aby skompilować, pobierz najnowsze źródło PyPy i wykonaj następujące czynności:

pypy /path/to/pypy-src/rpython/bin/rpython matrix-permanent.py

Wynikowy plik wykonywalny zostanie nazwany matrix-permanent-club podobny w bieżącym katalogu roboczym.

Począwszy od PyPy 5.0, prymitywy wątków RPython są znacznie mniej prymitywne niż kiedyś. Nowo odrodzone wątki wymagają GIL, który jest mniej lub bardziej bezużyteczny w obliczeniach równoległych. Użyłem forkzamiast tego, więc może nie działać zgodnie z oczekiwaniami w systemie Windows, chociaż nie przetestowałem, aby się nie skompilować ( unresolved external symbol _fork).

Plik wykonywalny akceptuje maksymalnie dwa parametry wiersza poleceń. Pierwszy to liczba wątków, drugi opcjonalny parametr to n. Jeśli zostanie podana, zostanie wygenerowana losowa macierz, w przeciwnym razie zostanie odczytana ze standardowego wejścia. Każdy wiersz musi być oddzielony znakiem nowej linii (bez końcowego znaku nowej linii), a każda przestrzeń wartości oddzielona. Trzeci przykładowy przykład będzie podany jako:

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 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 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 -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 -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 -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 -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 -1

Przykładowe użycie

$ time ./matrix-permanent-c 8 30
8395059644858368

real    0m8.582s
user    1m8.656s
sys     0m0.000s

metoda

Użyłem formuły Balasubramanian-Bax / Franklin-Glynn ze złożonością czasową O (2 ⁿ n) . Jednak zamiast iterować δ w szarej kolejności kodu, zamiast tego zamieniłem mnożenie przez wektor-wiersz pojedynczą operacją xor (mapowanie (1, -1) → (0, 1)). Sumę wektorową można również znaleźć w jednej operacji, biorąc n minus dwukrotność liczby popcount.

Rakieta 84 bajtów

Poniższa prosta funkcja działa dla mniejszych matryc, ale zawiesza się na mojej maszynie dla większych matryc:

(for/sum((p(permutations(range(length l)))))(for/product((k l)(c p))(list-ref k c)))

Nie golfowany:

(define (f ll) 
  (for/sum ((p (permutations (range (length ll))))) 
    (for/product ((l ll)(c p)) 
      (list-ref l c))))

Kod można łatwo zmodyfikować dla nierównej liczby wierszy i kolumn.

Testowanie:

(f '[[ 1 -1 -1  1]
     [-1 -1 -1  1]
     [-1  1 -1  1]
     [ 1 -1 -1  1]])

(f '[[ 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 -1  1 -1  1  1  1  1]])

Wydajność:

-4
192

Jak wspomniałem powyżej, testowanie polega na testowaniu następujących elementów:

(f '[[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 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 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 -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 -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 -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 -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 -1]])