Skip to content

Trovare il numero di matrici n per n (-1, 0, 1) con zero permanente nel più breve tempo possibile

Ti consigliamo di rivedere questa risposta in un ambiente controllato prima di inviarla alla produzione, saluti.

Soluzione:

Ruggine, $A(5) = 186481694371$ in 0,2 s, $A(6) = 19733690332538577$ in 580 s

(Tempi non ufficiali su un Core i7-10710U con 6 core/12 thread).

src/main.rs

use itertools::Itertools;
use num_integer::gcd;
use rayon::prelude::*;
use std::cmp::Ordering;
use std::collections::HashMap;
use std::convert::TryFrom;
use std::sync::Mutex;

type Count = u64;
type Coeff = i32;

fn sort_scols(n: u32, i: u32, sets: &[u32], state: &[Coeff]) -> Vec {
    let mut scol_ranks: Vec = (0..2 * n).map(|scol| (scol >> 1 > i) as u32).collect();
    let mut sorted_scols: Vec = (0..2 * n).collect();
    let mut scol_descs = vec![vec![]; 2 * n as usize];
    for (set_index, (&set, &coeff)) in sets.iter().zip(state).enumerate() {
        if coeff != 0 {
            for col in 0..n {
                if set & 1 << col != 0 {
                    scol_descs[(col as usize) << 1].push(set_index);
                    scol_descs[(col as usize) << 1 | 1].push(set_index);
                }
            }
        }
    }
    let mut sorted_set_indexes: Vec<_> = (0..sets.len()).collect();
    let mut set_ranks: Vec = vec![!0; sets.len()];
    let mut set_descs = vec![(0, vec![]); sets.len()];

    while {
        for (&coeff, set_desc) in state.iter().zip(&mut set_descs) {
            set_desc.0 = coeff;
            set_desc.1.clear();
        }
        let mut seen = 0;
        for &scol in &sorted_scols {
            if seen & 1 << (scol >> 1) == 0 {
                seen |= 1 << (scol >> 1);
                if scol & 1 == 0 {
                    for (&set, set_desc) in sets.iter().zip(&mut set_descs) {
                        if set_desc.0 != 0 && set & 1 << (scol >> 1) != 0 {
                            set_desc.1.push(scol_ranks[scol as usize]);
                        }
                    }
                } else {
                    for (&set, set_desc) in sets.iter().zip(&mut set_descs) {
                        if set_desc.0 != 0 && set & 1 << (scol >> 1) != 0 {
                            set_desc.0 = -set_desc.0;
                            set_desc.1.push(scol_ranks[scol as usize]);
                        }
                    }
                }
            }
        }

        let set_index_key =
            |&set_index: &usize| (set_descs[set_index].0.abs(), &set_descs[set_index].1);
        sorted_set_indexes.sort_by_key(set_index_key);
        if sets.len() != 0 {
            let mut next_set_rank = 0;
            set_ranks[sorted_set_indexes[0]] = next_set_rank;
            for (&set_index0, &set_index1) in sorted_set_indexes.iter().tuple_windows() {
                if set_index_key(&set_index0) < set_index_key(&set_index1) {
                    next_set_rank += 1;
                }
                set_ranks[set_index1] = next_set_rank;
            }
        }

        let scol_set_index_key = |scol: u32| {
            let sign = 1 - 2 * (scol & 1) as Coeff;
            let set_ranks = &set_ranks[..];
            let set_descs = &set_descs[..];
            move |&set_index: &usize| (set_ranks[set_index], set_descs[set_index].0 * sign)
        };
        for (scol, scol_desc) in scol_descs.iter_mut().enumerate() {
            scol_desc.sort_by_key(scol_set_index_key(scol as u32));
        }
        let scol_cmp = |&scol0: &u32, &scol1: &u32| {
            scol_ranks[scol0 as usize]
                .cmp(&scol_ranks[scol1 as usize])
                .then_with(|| {
                    scol_descs[scol0 as usize]
                        .iter()
                        .map(scol_set_index_key(scol0))
                        .cmp(
                            scol_descs[scol1 as usize]
                                .iter()
                                .map(scol_set_index_key(scol1)),
                        )
                })
        };
        sorted_scols.sort_by(scol_cmp);

        let mut new_scol_ranks = vec![!0; 2 * n as usize];
        let mut next_rank = 0;
        let mut changed = false;
        new_scol_ranks[sorted_scols[0] as usize] = next_rank;
        for (&scol0, &scol1) in sorted_scols.iter().tuple_windows() {
            if scol_ranks[scol0 as usize] < scol_ranks[scol1 as usize] {
                next_rank += 1;
            } else if scol_cmp(&scol0, &scol1) == Ordering::Less {
                changed = true;
                next_rank += 1;
            }
            new_scol_ranks[scol1 as usize] = next_rank;
        }
        scol_ranks = new_scol_ranks;

        changed
    } {}

    let mut seen = 0;
    sorted_scols.retain(|scol| {
        seen & 1 << (scol >> 1) == 0 && {
            seen |= 1 << (scol >> 1);
            true
        }
    });
    assert_eq!(seen, !(!0 << n));

    sorted_scols
}

fn count(n: u32) -> Count {
    static ERR: &str = "n is too large";
    Count::checked_pow(3, TryFrom::try_from(n * n).expect(ERR)).expect(ERR);
    (1..=Coeff::try_from(n).expect(ERR)).fold(1, |x: Coeff, y| x.checked_mul(y).expect(ERR));

    let mut cur_sets = vec![0];
    let mut cur_set_indexes = vec![0];

    let mut cur_states = HashMap::new();
    cur_states.insert(vec![1], 1);

    for k in 0..n {
        for i in 0..n {
            let keep_set = |&set: &u32| {
                if set.count_ones() == k {
                    !set & !1 << i & !(!0 << n) != 0
                } else {
                    set & 1 << i == 0
                }
            };
            let enlarge_set = |&set: &u32| set.count_ones() == k && set & 1 << i == 0;
            let new_sets: Vec = cur_sets
                .iter()
                .copied()
                .filter(keep_set)
                .chain(
                    cur_sets
                        .iter()
                        .copied()
                        .filter(enlarge_set)
                        .map(|set| set | 1 << i),
                )
                .collect();
            let mut new_set_indexes = vec![!0; 1 << n];
            for (set_index, &set) in new_sets.iter().enumerate() {
                new_set_indexes[set as usize] = set_index;
            }

            let new_states = Mutex::new(HashMap::new());
            cur_states.into_par_iter().for_each(|(cur_state, count)| {
                for elt in -1..=1 {
                    let mut new_state: Vec = cur_sets
                        .iter()
                        .copied()
                        .zip(&cur_state)
                        .filter(|(set, _)| keep_set(set))
                        .map(|(_, &coeff)| coeff)
                        .chain(
                            cur_sets
                                .iter()
                                .copied()
                                .zip(&cur_state)
                                .filter(|(set, _)| enlarge_set(set))
                                .map(|(set, &coeff)| {
                                    let mut new_coeff = coeff * elt;
                                    if set & !(!0 << i) != 0 {
                                        new_coeff +=
                                            cur_state[cur_set_indexes[set as usize | 1 << i]];
                                    }
                                    new_coeff
                                }),
                        )
                        .collect();

                    let factor = new_state.iter().copied().fold(0, gcd);
                    if factor > 1 {
                        for coeff in &mut new_state {
                            *coeff /= factor;
                        }
                    }

                    let sorted_scols = sort_scols(n, i, &new_sets, &new_state);

                    let new_state: Vec = new_sets
                        .iter()
                        .map(|&set| {
                            let mut permuted_set = 0;
                            let mut sign = 0;
                            for (col, &permuted_scol) in sorted_scols.iter().enumerate() {
                                sign ^= set >> col & permuted_scol & 1;
                                permuted_set |= (set >> col & 1) << (permuted_scol >> 1);
                            }
                            (1 - 2 * sign as Coeff)
                                * new_state[new_set_indexes[permuted_set as usize] as usize]
                        })
                        .collect();

                    *new_states.lock().unwrap().entry(new_state).or_insert(0) += count;
                }
            });

            cur_sets = new_sets;
            cur_set_indexes = new_set_indexes;
            cur_states = new_states.into_inner().unwrap();
        }
    }

    *cur_states.get(&vec![0]).unwrap_or(&0)
}

fn main() {
    for n in 1..=6 {
        println!("{} {}", n, count(n));
    }
}

Cargo.toml

[package]
name = "permanent"
version = "0.1.0"
authors = ["Anders Kaseorg <[email protected]>"]
edition = "2018"

[dependencies]
num-integer = "0.1.41"
rayon = "1.2.1"
itertools = "0.8.2"

Costruzione ed esecuzione

$ cargo build --release
$ time target/release/permanent
1 1
2 33
3 7555
4 13482049
5 186481694371
6 19733690332538577
5995.30user 282.67system 9:42.97elapsed 1076%CPU (0avgtext+0avgdata 6479028maxresident)k
0inputs+0outputs (0major+3259222minor)pagefaults 0swaps

Come funziona

Questo conta le matrici utilizzando la programmazione dinamica, quindi non è necessario enumerarle tutte direttamente. Fondamentalmente, lo stato di una matrice parzialmente riempita è rappresentato dai coefficienti di un'espressione per il permanente come un polinomio negli elementi rimanenti non riempiti. In particolare, corrisponde a un'esecuzione parziale del seguente algoritmo per il calcolo del permanente:

$operatorname{perm} A = c_{\{0, 1, dotsc, n - 1}}$, dove c_varnothing = 1, $c_S = sum_{i \\code(0144)} in S} c_{S mallsetminus {i}}a_{|S| - 1, i}$ per S ´ne ´varnothing´$.

Per esempio, ecco la formula $n = 4$ scritto esplicitamente:

c = 1

c0 = c * a00
c1 = c * a01
c2 = c * a02
c3 = c * a03

c01  = c1 * a10; c02  = c2 * a10; c03  = c3 * a10
c01 += c0 * a11; c12  = c2 * a11; c13  = c3 * a11
c02 += c0 * a12; c12 += c1 * a12; c23  = c3 * a12
c03 += c0 * a13; c13 += c1 * a13; c23 += c2 * a13

c012  = c12 * a20; c013  = c13 * a20; c023  = c23 * a20
c012 += c02 * a21; c013 += c03 * a21; c123  = c23 * a21
c012 += c01 * a22; c023 += c03 * a22; c123 += c13 * a22
c013 += c01 * a23; c023 += c02 * a23; c123 += c12 * a23

c0123  = c123 * a30
c0123 += c023 * a31
c0123 += c013 * a32
c0123 += c012 * a33

return c0123

Dopo l'esecuzione di ogni riga, possiamo dimenticare gli elementi compilati a… e rappresentare lo stato corrente con i soli elementi c… che sono ancora attive. Questa rappresentazione fa collassare molte matrici diverse parzialmente riempite nello stesso vettore di coefficienti. Ad esempio, scambiando due righe riempite si ottiene lo stesso vettore di coefficienti.

Il sort_scols comprime ulteriormente lo spazio degli stati utilizzando le simmetrie delle permutazioni delle colonne e il capovolgimento dei segni delle colonne per cercare di trovare una rappresentazione canonica per ogni stato. Non fa un lavoro perfetto, ma questa imperfezione ha un impatto solo sulle prestazioni, non sulla correttezza.

C# ($textrm{A330134}(5)=186481694371$ in meno di 10 secondi)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;

namespace Sandbox
{
    // See https://codegolf.stackexchange.com/q/196712
    static class Permanent
    {
        const int numThreads = 12;

        static void Main(string[] args)
        {
            var sw = new System.Diagnostics.Stopwatch();
            sw.Start();
            for (int n = 1; n < 8; n++)
            {
                Console.WriteLine($"{n}t{_PPCG196712_LaplaceMT(n)}t{sw.ElapsedMilliseconds}ms");
            }
        }

        internal static BigInteger _PPCG196712_LaplaceMT(int n)
        {
            // If we have a zero row then the permanent is zero.
            // If the first zero row is number i then the previous i rows have (3^n)-1 possibilities each, and the remaining ones have 3^n.
            BigInteger result = 0;
            int threePowN = (int)BigInteger.Pow(3, n);
            for (int i = 0; i < n; i++)
            {
                result += BigInteger.Pow(threePowN - 1, i) * BigInteger.Pow(threePowN, n - i - 1);
            }

            // Let's call a row containing a permutation of 1^1 0^{n-1} a "singleton row".
            // Under Laplace expansion it simplifies to a single {n-1}x{n-1} permanent.
            // If we have two singleton rows with the same non-zero index, that {n-1}x{n-1} permanent has a zero row, so the nxn permanent must be zero.
            // If we have multiple distinct singleton rows, we still need to worry about double-counting.
            // But all this effort would give roughly a 1.5% improvement at n=6, so is probably not worth it.

            // Now the interesting rows. We cluster similar rows as a technique to optimise column permutations.
            var clusters = new List();
            for (int ones = 1; ones <= n; ones++)
                for (int negs = 0; negs <= ones && ones + negs <= n; negs++)
                    clusters.Add(new Cluster(n, ones, negs));

            // Consider partitions of n-1 by cluster.
            foreach (var partition in PartitionAssignments(clusters, n-1))
            {
                var subclusters = partition.OrderBy(tuple =>
                {
                    var (cluster, k) = tuple;
                    // What's the speedup factor of pivoting on this cluster?
                    return k < cluster.DistinctCounts.Count ? cluster.DistinctCounts[k] / (double)cluster.NaiveCounts[k] : 1;
                }).ToList();

                // We pivot on the first cluster, and take a Cartesian product.
                result += Enumerable.Range(0, numThreads).AsParallel().
                    Select(threadId => CountByPartition(n, subclusters, 0, new sbyte[n - 1][], 0, 1, threadId)).
                    Aggregate(BigInteger.Zero, (a, b) => a + b);
            }

            return result;
        }

        private static BigInteger CountByPartition(int n, IReadOnlyList<(Cluster cluster, int repeat)> clusters, int clusterIdx, sbyte[][] matrix, int rowIdx, BigInteger weight, int threadId)
        {
            if (clusterIdx == clusters.Count) return weight * Count_Laplace(matrix) << (n - 1);

            int k = clusters[clusterIdx].repeat;
            if (clusterIdx == 0 && clusters[0].repeat < clusters[0].cluster.DistinctCounts.Count)
            {
                // Pivot
                weight *= Binom(n - 1, k) * Factorial(n - 1 - k); // Choose the pivoted rows and then permute the others freely

                BigInteger sum = 0;
                var reps = clusters[0].cluster.Representatives[k];

                int delta = reps.Count >= numThreads || clusterIdx == clusters.Count - 1 ? numThreads : 1;
                int skip = delta == 1 ? 0 : threadId;

                for (int repIdx = skip; repIdx < reps.Count; repIdx += delta)
                {
                    var (prefixWeight, prefix) = reps[repIdx];
                    for (int i = 0; i < k; i++) matrix[rowIdx + i] = prefix[i];
                    sum += CountByPartition(n, clusters, clusterIdx + 1, matrix, rowIdx + k, weight * prefixWeight, delta == 1 ? threadId : -1);
                }
                return sum;
            }

            {
                if (clusterIdx == 0) weight = Factorial(n - 1); // Permute all rows freely

                var rows = clusters[clusterIdx].cluster.Rows;
                var indices = new int[k];
                BigInteger sum = 0;

                int skip = threadId > 0 ? threadId : 0;
                int delta = threadId >= 0 ? numThreads : 1;

                while (true)
                {
                    if (skip == 0)
                    {
                        for (int i = 0; i < k; i++)
                        {
                            var j = indices[i];
                            var inRow = rows[j];
                            matrix[rowIdx + i] = inRow;
                        }

                        // Adjust the weight to account for duplicates in this cluster.
                        var recurWeight = weight;
                        int currRun = 1;
                        for (int i = 1; i < k; i++)
                        {
                            if (indices[i] == indices[i - 1]) currRun++;
                            else { recurWeight /= Factorial(currRun); currRun = 1; }
                        }
                        recurWeight /= Factorial(currRun);

                        sum += CountByPartition(n, clusters, clusterIdx + 1, matrix, rowIdx + k, recurWeight, -1);
                    }

                    skip--;
                    if (skip < 0) skip += delta;

                    if (indices[0] == rows.Count - 1) break;
                    NextOrderedMultiset(indices, rows.Count);
                }

                return sum;
            }
        }

        private static IEnumerable> PartitionAssignments(IReadOnlyList elts, int n)
        {
            IEnumerable> Inner(int m, int maxPart, int off, List<(T, int)> prefix, ISet indicesUsed)
            {
                if (m == 0)
                {
                    yield return new List<(T, int)>(prefix);
                    yield break;
                }

                for (int part = Math.Min(m, maxPart); part > 0; part--)
                {
                    // When we reduce the part, we re-enable skipped elts.
                    for (int i = part == maxPart ? off : 0; i < elts.Count; i++)
                    {
                        if (indicesUsed.Contains(i)) continue;

                        prefix.Add((elts[i], part));
                        indicesUsed.Add(i);

                        foreach (var soln in Inner(m - part, part, i + 1, prefix, indicesUsed)) yield return soln;

                        indicesUsed.Remove(i);
                        prefix.RemoveAt(prefix.Count - 1);
                    }
                }
            }

            return Inner(n, n, 0, new List<(T, int)>(), new HashSet());
        }

        private static void NextOrderedMultiset(int[] indices, int @base)
        {
            int j = indices.Length - 1;
            while (indices[j] == @base - 1) j--;
            indices[j]++;
            j++;
            while (j < indices.Length) { indices[j] = indices[j - 1]; j++; }
        }

        private static BigInteger Count_Laplace(sbyte[][] matrix)
        {
            int m = matrix.Length;
            int n = m + 1;
            BigInteger[] sums = new BigInteger[m];
            uint gray = 0;
            int sign = 1;
            BigInteger[] weights = new BigInteger[n];
            for (uint i = 0; i < 1u << n; i++)
            {
                BigInteger term = 1;
                foreach (var sum in sums) term *= sum;

                for (int j = 0; j < n; j++)
                {
                    if (((gray >> j) & 1) == 1) weights[j] += sign * term;
                }

                // Gray code update
                int flipPos = (i + 1).CountTrailingZeros();
                if (flipPos == n) break;

                uint bit = 1u << flipPos;
                if ((gray & bit) == bit)
                {
                    for (int j = 0; j < m; j++) sums[j] -= matrix[j][flipPos];
                }
                else
                {
                    for (int j = 0; j < m; j++) sums[j] += matrix[j][flipPos];
                }
                gray ^= bit;
                sign = -sign;
            }

            // Count weighted subsets which total zero, optimising with meet-in-the-middle
            var distribL = SubsetSums(weights, 0, n >> 1);
            var distribR = SubsetSums(weights, n >> 1, n - (n >> 1));
            BigInteger total = -1; // Discount the zero column solution
            foreach (var kvp in distribL)
            {
                if (distribR.TryGetValue(-kvp.Key, out var v2)) total += kvp.Value * v2;
            }

            return total;
        }

        private static Dictionary SubsetSums(BigInteger[] weights, int off, int len)
        {
            var distrib = new Dictionary { [0] = 1 };
            for (int i = 0; i < len; i++)
            {
                var weight = weights[off + i];
                var nextDistrib = new Dictionary(distrib);
                foreach (var kvp in distrib)
                {
                    nextDistrib.TryGetValue(kvp.Key - weight, out var tmp);
                    nextDistrib[kvp.Key - weight] = tmp + kvp.Value;

                    nextDistrib.TryGetValue(kvp.Key + weight, out tmp);
                    nextDistrib[kvp.Key + weight] = tmp + kvp.Value;
                }

                distrib = nextDistrib;
            }

            return distrib;
        }

        class Cluster
        {
            internal readonly IReadOnlyList Rows;

            internal readonly IReadOnlyList NaiveCounts;
            internal readonly IReadOnlyList DistinctCounts;
            internal readonly IReadOnlyList prefix)>> Representatives;

            internal Cluster(int n, int ones, int negs)
            {
                if (ones == 0) throw new ArgumentOutOfRangeException(nameof(ones));
                if (negs > ones) throw new ArgumentOutOfRangeException(nameof(negs));

                Rows = BuildRows(n, ones, negs);

                // Build pairs and triples, and calculate their advantage factors.
                var naiveCounts = new List { 1, Rows.Count };
                var distinctCounts = new List { 1, 1 };
                var representatives = new List)>>
                {
                    new List<(int, IReadOnlyList)>{ (1, new List()) },
                    new List<(int, IReadOnlyList)>{ (Rows.Count, new List { Rows.First() }) }
                };

                for (int k = 2; k < 4; k++)
                {
                    var counts = new Dictionary();
                    var reps = new Dictionary>();
                    int uncollapsed = 0;

                    var indices = new int[k];
                    while (true)
                    {
                        var subset = new sbyte[k][];
                        for (int i = 0; i < k; i++) subset[i] = Rows[indices[i]];

                        BigInteger count = Factorial(k);
                        int currRun = 1;
                        for (int i = 1; i < k; i++)
                        {
                            if (indices[i] == indices[i - 1]) currRun++;
                            else { count /= Factorial(currRun); currRun = 1; }
                        }
                        count /= Factorial(currRun);

                        uncollapsed += (int)count;
                        var encoding = CanonicalEncoding(subset);
                        counts.TryGetValue(encoding, out var oldCount);
                        counts[encoding] = oldCount + (int)count;
                        if (oldCount == 0) reps[encoding] = (sbyte[][])subset.Clone();

                        if (indices[0] == Rows.Count - 1) break;
                        NextOrderedMultiset(indices, Rows.Count);
                    }

                    naiveCounts.Add(uncollapsed);
                    distinctCounts.Add(reps.Count);
                    representatives.Add(counts.Select(kvp => (kvp.Value, reps[kvp.Key])).ToList());
                }

                NaiveCounts = naiveCounts;
                DistinctCounts = distinctCounts;
                Representatives = representatives;
            }

            private static List BuildRows(int n, int ones, int negs)
            {
                var rows = new List();
                void Append(sbyte[] row)
                {
                    // Ensure that we don't have both a row and its negative by testing canonicity.
                    if (ones > negs || Array.IndexOf(row, (sbyte)1) < Array.IndexOf(row, (sbyte)-1)) rows.Add((sbyte[])row.Clone());
                }

                var perm = new sbyte[n];
                for (int i = 0; i < negs; i++) perm[i] = -1;
                for (int i = n - ones; i < n; i++) perm[i] = 1;

                while (true)
                {
                    Append(perm);

                    // Next permutation
                    int k;
                    for (k = n - 2; k >= 0; k--)
                    {
                        if (perm[k] < perm[k + 1]) break;
                    }
                    if (k == -1) return rows;

                    int l;
                    for (l = n - 1; l > k; l--)
                    {
                        if (perm[k] < perm[l]) break;
                    }

                    var tmp = perm[k]; perm[k] = perm[l]; perm[l] = tmp;

                    for (int i = k + 1, j = n - 1; i < j; i++, j--)
                    {
                        tmp = perm[i]; perm[i] = perm[j]; perm[j] = tmp;
                    }
                }
            }

            private static BigInteger CanonicalEncoding(sbyte[][] rows)
            {
                return rows.Permutations().Select(_CanonicalEncodingInner).Min();
            }

            private static BigInteger _CanonicalEncodingInner(sbyte[][] rows)
            {
                BigInteger @base = BigInteger.Pow(3, rows.Length);
                BigInteger[] digits = new BigInteger[rows[0].Length];
                for (int i = 0; i < digits.Length; i++)
                {
                    for (int j = 0; j < rows.Length; j++)
                    {
                        digits[i] = digits[i] * 3 + rows[j][i] + 1;
                    }
                }
                return digits.OrderBy(x => x).Aggregate(BigInteger.Zero, (accum, digit) => accum * @base + digit);
            }
        }

        private static readonly int[] _DeBruijn32 = new int[]
        {
            0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
            31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
        };

        // Returns 32 for input 0.
        public static int CountTrailingZeros(this uint i)
        {
            if (i == 0) return 32;

            int j = (int)i;
            i = (uint)(j & -j);
            var idx = (i * 0x077CB531U) >> 27;
            return _DeBruijn32[idx];
        }

        public static IEnumerable Permutations(this IEnumerable elts)
        {
            T[] arr = elts.ToArray();
            if (arr.Length == 0)
            {
                yield return arr;
                yield break;
            }

            int[] indices = Enumerable.Range(0, arr.Length).ToArray();
            while (true)
            {
                yield return arr.ToArray();
                // Next permutation
                int i = indices.Length - 1;
                while (i > 0 && indices[i - 1] > indices[i]) i--;

                if (i == 0) yield break;

                int j = indices.Length - 1;
                while (indices[j] < indices[i - 1]) j--;

                _DoubleSwap(arr, indices, i - 1, j);
                for (j = indices.Length - 1; i < j; i++, j--) _DoubleSwap(arr, indices, i, j);
            }
        }

        private static void _DoubleSwap(T[] elts, int[] indices, int i, int j)
        {
            var t1 = elts[i];
            elts[i] = elts[j];
            elts[j] = t1;

            var t2 = indices[i];
            indices[i] = indices[j];
            indices[j] = t2;
        }

        private static BigInteger Factorial(int n) => n < 2 ? BigInteger.One : n * Factorial(n - 1);

        private static BigInteger Binom(int n, int k)
        {
            BigInteger result = 1;
            for (int i = 0; i < k; i++) result = result * (n - i) / (i + 1);
            return result;
        }
    }
}

Demo online

Questo utilizza i seguenti trucchi:

  • Se c'è una riga zero, il permanente è zero. Questi possono essere contati molto velocemente e quindi si può assumere che non ci sia una riga zero.
  • Moltiplicare una riga per $-1$ inverte il segno del permanente. Quindi possiamo assumere un segno canonico per ogni riga e poi moltiplicare per $2^n$.
  • La permutazione delle righe lascia invariato il permanente, quindi possiamo assumere una permutazione canonica delle righe e poi moltiplicare per un fattore appropriato.
  • Il permanente può essere calcolato in un tempo migliore della forza bruta utilizzando la formula di Ryser con enumerazione in codice Gray di sottoinsiemi di colonne.
  • L'espansione di Laplace dell'ultima riga trasforma un fattore di $frac{3^n - 1}{2}$ in un'istanza del problema della somma di sottoinsiemi, che può essere affrontato utilizzando la programmazione dinamica con il meet-in-the-middle.
  • Possiamo anche permutare le colonne, il che è quanto segue Cluster serve a questo. Credo di aver spinto il Cluster al massimo, ma vorrei trovare un modo efficiente per sostituirlo, perché se prendiamo la prima riga [1, 1, 0, -1, -1, -1] (con peso 60) sarebbe bello sfruttare la simmetria delle prime due colonne e delle ultime tre colonne nel cluster successivo. Il problema è che fare questo in modo ingenuo aggiunge un fattore di $n!$ che è più costoso del risparmio.
  • Le permutazioni sono indipendenti, quindi dopo un po' di pre-elaborazione abbiamo un problema imbarazzantemente parallelo.

Ho compilato con VS2017 Community puntando su .Net Framework 4.7.1 e compilando per "Release" con le impostazioni predefinite. Si potrebbe preferire l'uso di .Net Core.

$n=6$ è ancora tecnicamente fuori portata - l'ho calcolato in poco meno di 3 ore su un AMD a 4 core con numThreads = 8.

Commenti e valutazioni



Utilizzate il nostro motore di ricerca

Ricerca
Generic filters

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.