Calculate a term in the expansion of the fraction of configurations that don't cause a black hole to form; see https://math.stackexchange.com/questions/3620483.

Question3620483.java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.function.Consumer;
import java.util.stream.Collectors;

import BallsInBins; // https://gist.github.com/joriki/d324726fbc3fe4ccdf33628aa02b790c
import BigRational; // https://gist.github.com/joriki/3097452fbb9e1983daf1035f57e72722
import Binomials;   // https://gist.github.com/joriki/5a5c03d2c286effa69584e77a064a6e3

public class Question3620483 {
    final static int k = 2;
    
    static int l,n;

    static List<Disjunction> conditions = new ArrayList<>();

    static BigRational [] deltaCounts = new BigRational [4 * k + 1];
    
    static class Disjunction {
        List<Long> clauses = new ArrayList<>();
        
        public Disjunction () {}
        
        public Disjunction (Disjunction d1,Disjunction d2) {
            for (Long c1 : d1.clauses)
                for (Long c2 : d2.clauses)
                    addConstraint(c1,c2);
        }

        void addConstraint(long c1,long c2) {
            if ((c1 | c2) < 0)
                throw new Error(Long.toHexString(c1) + " / " + Long.toHexString(c2));
            if (((c1 ^ c2) & (6 * (c1 & c2 & 0x9249249249249249L))) == 0)
                clauses.add(c1 | c2);
        }
        
        public String toString() {
            return clauses.stream().map(constraint -> {
                    long [] [] constraints = new long [2] [2 * k];
                    for (int i = 0;i < 2;i++)
                        for (int j = 0;j < 2 * k;j++) {
                            long c = constraint & 7;
                            constraint >>= 3;
                            constraints [i] [j] = (c & 1) == 0 ? -1 : c >> 1;
                        }
                    return Arrays.deepToString(constraints);
            }).collect(Collectors.joining(" ∨ "));
        }
    }
    
    public static void main(String [] args) {
        for (l = 1;l <= k;l++)
            for (int i = 0;i <= l;i++) {
                final int x = i;
                for (int j = 0;j <= l;j++) {
                    final int y = j;
                    final Disjunction d = new Disjunction();
                    for (int nballs = 0;nballs < l;nballs++)
                        BallsInBins.traverseConfigurations(nballs,2 * l,new Consumer<int[]>() {
                            public void accept(int [] t) {
                                long constraint = 0;
                                for (int a = 0;a < l;a++) {
                                    int shift = a;
                                    if (shift >= x)
                                        shift += 2 * k - l;
                                    constraint |= ((t [a] << 1) + 1L) << (3 * shift);
                                }
                                for (int b = 0;b < l;b++) {
                                    int shift = b;
                                    if (shift >= y)
                                        shift += 2 * k - l;
                                    constraint |= ((t [b + l] << 1) + 1L) << (3 * (shift + 2 * k));
                                }
                                if (constraint < 0)
                                    throw new Error();
                                d.clauses.add(constraint);
                            }
                        });
                    conditions.add(d);
                }
            }
        
        for (Disjunction d : conditions)
            System.out.println(d);
        
        n = conditions.size();
        System.out.println(n + " conditions");
        
        Arrays.fill(deltaCounts, BigRational.ZERO);

        Disjunction emptyDisjunction = new Disjunction();
        emptyDisjunction.clauses.add(0L);

        recurse (emptyDisjunction,0,1);
        
        double c = 0;
        for (int i = 0;i < deltaCounts.length;i++) {
            System.out.println(i + " : " + deltaCounts [i]);
            c += deltaCounts [i].doubleValue() * Math.exp(-i);
        }
        System.out.println();
        System.out.println(c);
    }

    static int count;

    static void recurse (Disjunction disjunction,int index,int sign) {
        if (index == n) {
            if ((++count & 0xfff) == 0)
                System.out.println(Integer.toHexString(count));
            
            for (long constraint : disjunction.clauses) {
                int delta = 0;
                long factor = 1;
                while (constraint != 0) {
                    int c = (int) (constraint & 7);
                    constraint >>= 3;
                    if ((c & 1) != 0) {
                        delta++;
                        factor *= Binomials.factorial(c >> 1).intValue();
                    }
                }
                deltaCounts [delta] = BigRational.sum(deltaCounts [delta],new BigRational(sign,factor));
            }
        }
        else {
            recurse (disjunction,index + 1,sign);
            recurse (new Disjunction(disjunction,conditions.get(index)),index + 1,-sign);
        }
    }
}