Compute the probability that a random permutation of the 54 stickers on a Rubik's cube includes no adjacent squares of the same colour; see https://math.stackexchange.com/questions/4618350.

Question4618350.java

import java.math.BigInteger;
import java.util.Arrays;
import java.util.HashMap;
import java.util.Iterator;
import java.util.Map;
import java.util.function.BiConsumer;
import java.util.stream.Collectors;
import java.util.stream.IntStream;

public class Question4618350 {
    final static int d = 3;
    final static int ncolors = 6;
    final static BigInteger size = BigInteger.valueOf(d * d);
    
    static class Monomial {
        // These are negative so we can easily sort them in descending order so the zeros are at the end
        BigInteger [] exponents;
        
        public Monomial () {}
        
        public Monomial (int [] exponents) {
            this.exponents = IntStream.of(exponents).mapToObj(exponent -> BigInteger.valueOf(-exponent)).toArray(BigInteger[]::new);
            Arrays.sort(this.exponents);
        }
        
        public Monomial (BigInteger [] exponents) {
            this.exponents = exponents.clone();
            Arrays.sort(this.exponents);
        }
        
        public boolean equals (Object o) {
            return Arrays.equals(((Monomial) o).exponents, exponents);
        }
        
        public int hashCode () {
            return Arrays.hashCode(exponents);
        }
        
        public String toString() {
            return Arrays.asList(exponents).stream().map(exponent -> exponent.negate().toString()).collect(Collectors.joining(","));
        }
        
        BigInteger weight() {
            int weight = 1;
            for (int i = 0;i < ncolors && exponents [i].signum() != 0;i++)
                weight *= ncolors - i;
            
            return BigInteger.valueOf (weight);
        }
    }
    
    public static class Polynomial {
        Map<Monomial,BigInteger> terms = new HashMap<>();
        
        public void add(Monomial monomial,BigInteger coefficient) {
            BigInteger old = terms.get(monomial);
            terms.put(monomial,old == null ? coefficient : coefficient.add(old));            
        }
        
        static Polynomial product (Polynomial p1,Polynomial p2) {
            Polynomial product = new Polynomial();
            for (Monomial term : p1.terms.keySet()) {
                BigInteger c = p1.terms.get(term);
                symmetrize (term,0,p2,(monomial,coefficient) -> {
                    product.add(monomial,coefficient.multiply(c)); });
            }
            for (Monomial term : product.terms.keySet()) {
                BigInteger coefficient = product.terms.get(term);
                BigInteger [] div = coefficient.divideAndRemainder(term.weight());
                if (div [1].signum() != 0)
                    throw new Error();
                product.terms.put(term,div [0]);
            }
            return product;
        }
        
        static void symmetrize (Monomial term,int i,Polynomial p2,BiConsumer<Monomial,BigInteger> consumer) {
            if (i == ncolors || term.exponents [i].signum() == 0)
                for (Monomial term2 : p2.terms.keySet()) {
                    Monomial newTerm = new Monomial ();
                    newTerm.exponents = term2.exponents.clone();
                    for (int j = 0;j < ncolors;j++)
                        if (map [j] != -1)
                            newTerm.exponents [j] = newTerm.exponents [j].add(term.exponents [map [j]]);
                    Arrays.sort(newTerm.exponents);
                    consumer.accept(newTerm,p2.terms.get(term2).multiply(term2.weight()));
            }
            else
                for (int j = 0;j < ncolors;j++)
                    if (map [j] == -1) {
                        map [j] = i;
                        symmetrize (term,i + 1,p2,consumer);
                        map [j] = -1;
                    }
        }
        
        public String toString() {
            return terms.size () + " / " + sum () + " : " + terms.keySet().stream().map(term -> terms.get(term) + " * " + term).collect(Collectors.joining(" + "));
        }
        
        public BigInteger sum () {
            BigInteger sum = BigInteger.ZERO;
            for (Monomial term : terms.keySet())
                sum = sum.add(terms.get(term).multiply(term.weight()));
            return sum;
        }
    }
    
    final static int [] map = new int [ncolors];
    static int [] [] colors = new int [d] [d];
    static int [] counts = new int [ncolors];
    static Polynomial side = new Polynomial();
    
    static void countAdmissibleConfigurations (int x,int y,int max) {
        if (++x == d) {
            if (++y == d) {
                side.add(new Monomial (counts),BigInteger.ONE);
                return;
            }
            x = 0;
        }
        for (int color = 0;color <= max;color++) {
            if (x != 0 && colors [x - 1] [y] == color)
                continue;
            if (y != 0 && colors [x] [y - 1] == color)
                continue;
            colors [x] [y] = color;
            counts [color]++;
            countAdmissibleConfigurations (x,y,color != max || max == ncolors - 1 ? max : max + 1);
            counts [color]--;
        }
    }

    public static void main(String [] args) {
        Arrays.fill(map,-1);

        countAdmissibleConfigurations (d - 1,-1,0);
        
        Polynomial power = side;
        
        for (int i = 1;i < ncolors;i++) {
            power = Polynomial.product (power,side);
            Iterator<Monomial> terms = power.terms.keySet().iterator();
            while (terms.hasNext())
                if (terms.next().exponents [0].negate().compareTo(size) > 0)
                    terms.remove();
            System.out.println(power);
        }        
    }
}