[Question]
  [
## Background
A classical [logic gate](https://en.wikipedia.org/wiki/Logic_gate) is an idealized electronic device implementing a Boolean function, i.e. one that takes a certain number of Boolean inputs and outputs a Boolean. We only consider two-input logic gates in this challenge.
There are \$2^{2 \times 2} = 16\$ distinct two-input logic gates (2×2 input combinations and 2 possible outputs for each input combination). Out of the 16, only two, namely NAND and NOR gates, are [**universal**](https://en.wikipedia.org/wiki/Logic_gate#Universal_logic_gates), in the sense that you can build arbitrary circuits (any Boolean function) by composing the single type of gates. CGCCers might recognize NAND gates because we had a few [challenges](https://codegolf.stackexchange.com/questions/tagged/atomic-code-golf+logic-gates) to build arbitrary circuits out of them.
```
NAND | 0   1      NOR | 0   1
-----+-------    -----+-------
  0  | 1   1       0  | 1   0
  1  | 1   0       1  | 0   0
```
We can generalize the concept of logic gates to n-ary inputs and outputs instead of Booleans, for example ternary for \$n=3\$. (Physical three-state circuits are also being studied.)
For ternary logic gates, the allowed states are often written as `- 0 +`. There are \$3^{3 \times 3} = 19683\$ distinct ternary logic gates in total, and Craig Gidney managed to [enumerate 3774 universal ternary logic gates](https://web.archive.org/web/20221102160829/http://twistedoakstudios.com/blog/Post7878_exploring-universal-ternary-gates) out of the 19683. ([Full C# code](https://github.com/Strilanc/UniversalTernaryGates))
Typos in the blog post and an off-by-one error corrected thanks to @AnttiP
One such logic gate looks like this (named `tand` in the blog post):
```
tand | -   0   +
-----+-----------
  -  | 0   +   +
  0  | +   +   +
  +  | +   +   -
```
Another looks like this: (X can take any value out of `-0+`. All three gates represented by this pattern are universal.)
```
  ?  | -   0   +
-----+-----------
  -  | +   -   0
  0  | X   -   -
  +  | +   +   -
```
On the other hand, [Malbolge's Crazy operator](https://en.wikipedia.org/wiki/Malbolge#Crazy_operation) is NOT universal.
```
 crz | -   0   +
-----+-----------
  -  | 0   -   -
  0  | 0   -   +
  +  | +   +   0
```
For \$n=1\$, there is only one possible logic gate (constant), and it is trivially universal because the one logic gate can build every possible (the only one) logic gate.
## Challenge
Given the value of \$n \ge 1\$, compute the number of distinct universal n-ary logic gates.
Standard [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'") rules apply. The shortest code in bytes wins. Note that your code must theoretically calculate the answer for any valid input, when run in an idealized machine (integers and/or floats have unlimited precision) with enough but finite time and memory.
The known values are:
```
1 -> 1
2 -> 2
3 -> 3774
4 -> 942897552 (first calculated by @AnttiP, needs cross-check)
```
  
[Answer]
# [Python 3.8 (pre-release)](https://docs.python.org/3.8/), 287 282 242 bytes
-5 bytes thanks to @Dialfrost
```
f=lambda n,g=0,c=0:g<n**len(e:=range(n*n))!=(s:={g,sum(i//n*n**i for i in e),sum(i%n*n**i for i in e)})and f(n,g+1,c+([s:=s|{sum(g//n**(l//n**i%n*n+r//n**i%n)%n*n**i for i in e)}for _ in"a"*n**n**2for l in s for r in s]!=len(s)==n**n**2))or c
```
[Try it online!](https://tio.run/##bY/BCsIwEER/JRWE3TTQVi9S3C8RkVjbGGi3JakHqX57TaqeFAI7zJtMssN9vPa83Q1unhtqdXe@aMHKUK4qykuzZynbmqEuyWk2NbBkxITAlzQZ5W8d2CwLppRWNL0TVlgWNb7J@hc8UfNFNBDeSAtVpXAITf4xxbyJTRLaZSyXU/fV@K8r6lPQK72KMJxNtNqI/ZJ0izwmFJfwSPSJIQZYzYOzPEIDBarwp@DOLw "Python 3.8 (pre-release) – Try It Online")
Uses a simple integer encoding of the gates, and composes them to depth `n**(n*n)` which should be sufficient. It's able to calculate up to `n=2`.
Full disclosure, this python golf is just to qualify this answer as code golf. The rest of this answer will explain how to calculate `n=4`
# Explanation of the algorithm used to calculate n=4
The following algorithm works for all `n>=2`. In this explanation I refer to gates (also) as functions. `f(a,b)` denotes the value of gate `f` with inputs `a` and `b`. Gates can be composed just like functions. To keep every function binary, I often use the special functions `l(a,b)=a` and `r(a,b)=b` (which just pass the left/right argument).
## Core algorithm
### Unary completeness
A gate `f` is said to be unary complete if it's possible to create every unary function using `f`. Obviously, `f` has to be unary complete to be functionally complete.
We first check if `f` is unary complete. If not, then it's not functionally complete. Then we move to the next section.
### Functional completeness
Assuming that `f` is unary complete, we can check if `f` is functionally complete with the two following tests.
### Lonelyness-test
We want to find some symbols `a, b, c, d` where `a≠b` and `c≠d`, so that `f(a, c)` is distinct from `f(a, d)`, `f(b, c)` and `f(b, d)`. This means that `f` is "lonely".
Using these values, and unary functions, we can simulate all of boolean logic.
If there are no such `a, b, c, d` I claim that `f` can't be functionally complete.
To see why, let's look at how non-lonely functions compose. Note that the functions `l(x,y)=x` and `r(x,y)=y` are not lonely. Assume `L` and `R` are non-lonely functions. `g(a,b)=f(L(a, b), R(a,b))`. In general, `L`, `R` and `f` have the one of the three following "truth" tables, for the relevant inputs:
```
. a b
c x x
d y y
. a b
c x y
d x y
. a b
c x y
d y x
```
It should be clear that composing these does not make the resulting function lonely. This is a finitary problem and can be easily proved with brute-force. It is the same reason `xor` and `not` are not functionally complete.
In other words, composing non-lonely functions makes non-lonely functions, and since some gates are lonely, a non-lonely gate cannot be universal.
### Constructible-test
Again, we assume `f` is unary complete, and also lonely.
A function `f` is n+1-constructible if there are some subsets of the alphabet `L` and `R` where `|L|=|R|=n`, `|f[L,R]|≥n+1`. A function `f` is constructible iff it is n-constructible for all `n≥3` up to the length of the alphabet.
If `f` is constructible, then it is definitely functionally complete. You can just convert your input to base-2, do whatever calculations you want, and then convert back to the full alphabet.
Now is it possible to not be constructible and still be functionally complete? No.
Let's assume that there is some `n` so that `f` is not n-constructible. Note that the functions `l(x,y)=x` and `r(x,y)=y` are never n-constructible. Assume `L` and `R` are non-n-constructible functions. Then `g(a,b)=f(L(a,b),R(a,b))` is also non-n-constructible. Because the image of `L` and `R` has cardinality `≤n`, so does `g`. Therefore `g` is also non-n-constructible
Therefore `f` has to be constructible to be functionally complete, because again, non-n-constructible gates compose to create non-n-constructible gates and some gates are n-constructible.
## Optimizations
### Data Format
A duadic function `f` is represented as a `n*n` list, so that `f(a, b) = list[a*n+b]`.
This is a bit slow, so for `n=4` I use a 32-bit integer.
Every function has a unique index. The indices start from 0 and have no gaps.
Every monadic function `f(a,b)=g(b)` has a unique index. It also starts from 0 and has no gaps.
### Transposition
A transposition `fᵀ` is defined as `fᵀ(a,b)=f(b,a)`
### Permutation
`p` is a permutation of the alphabet. `p(x)` applies the permutation. `p⁻¹(x)` applies the inverse. `p(p⁻¹(x))=p⁻¹(p(x))=x`.
`fₚ` is defined as `fₚ(a,b)=p⁻¹(f(p(a),p(b)))`
### Theorem 1
Claim: `f` is functionally complete, iff `fᵀ` also is.
Proof: You can just swap the arguments lol
### Theorem 2
Claim: `f` is functionally complete, iff `fₚ` also is (for all permutations `p`)
Proof: When you do function composition with `fₚ`, the outer and inner permutations cancel out. So the "meat" of the function doesn't change, just that the alphabet is relabled.
### Minimal representation
`f` is said to be minimal if `f≤fₚ` and `f≤fₚᵀ` for all permutations `p` (here `≤` refers to some total order). By theorem 1 and 2, we only need to consider minimal elements.
# Rust code for n=4
## main.rs
```
#![feature(let_else)]
#![feature(map_first_last)]
#![feature(adt_const_params)]
#![feature(generic_const_exprs)]
#![feature(label_break_value)]
use indicatif::{ProgressBar, ProgressStyle};
use itertools::Itertools;
use rand::Rng;
use rayon::prelude::*;
use std::sync::atomic::AtomicBool;
use std::sync::atomic::AtomicUsize;
use std::sync::atomic::Ordering::*;
use std::thread;
#[cfg(not(debug_assertions))]
macro_rules! get_unsafe {
    [$a:expr, $i:expr] => {
        *unsafe {$a.get_unchecked($i)}
    };
}
#[cfg(debug_assertions)]
macro_rules! get_unsafe {
    [$a:expr, $i:expr] => {
        *$a.get($i).unwrap()
    };
}
#[cfg(not(debug_assertions))]
macro_rules! get_mut_unsafe {
    [$a:expr, $i:expr] => {
        *unsafe {$a.get_unchecked_mut($i)}
    };
}
#[cfg(debug_assertions)]
macro_rules! get_mut_unsafe {
    [$a:expr, $i:expr] => {
        *$a.get_mut($i).unwrap()
    };
}
trait Function: std::hash::Hash + Clone + std::cmp::Eq + std::fmt::Debug + std::cmp::Ord {
    const N: usize;
    fn impl_eval(&self, a: usize, b: usize) -> usize;
    #[cfg(debug_assertions)]
    fn eval(&self, a: usize, b: usize) -> usize {
        assert!(
            a < Self::N && b < Self::N,
            "Called eval with a={}, b={} but N={}",
            a,
            b,
            Self::N
        );
        self.impl_eval(a, b)
    }
    #[cfg(not(debug_assertions))]
    fn eval(&self, a: usize, b: usize) -> usize {
        use std::hint::unreachable_unchecked;
        if a >= Self::N || b >= Self::N {
            // I am speed
            unsafe {
                unreachable_unchecked();
            }
        }
        self.impl_eval(a, b)
    }
    fn pass_left() -> Self;
    fn pass_right() -> Self;
    fn compose(&self, left: &Self, right: &Self) -> Self;
    fn unary_compose(&self, left: usize, right: usize) -> usize {
        self.compose(
            &Self::from_unary_index(left),
            &Self::from_unary_index(right),
        )
        .unary_index()
    }
    // Used by fuzzers
    fn random<R: Rng>(rng: &mut R) -> Self;
    // Returns true if should be discarded
    fn low_effort_discard(&self) -> bool {
        for n in 0..Self::N {
            if self.eval(n, n) == n {
                return true;
            }
        }
        false
    }
    // Unique, no holes, for r-unary ops
    fn unary_index(&self) -> usize {
        let mut ret = 0;
        for a in (0..Self::N).rev() {
            ret *= Self::N;
            ret += self.eval(a, a);
        }
        ret
    }
    fn from_unary_index(i: usize) -> Self;
    fn is_unary_complete(&self) -> bool
    where
        [(); Self::N.pow(Self::N as u32)]:,
    {
        // Bruteforce-checks if self is a unary-complete function
        // Keeps track of all the visited nodes and has a queue of functions to add
        // A function is popped from the queue and composed with every previous function, plus itself
        // The new compositions are added to the queue
        // visited is a "Linked list" containing the already visited unary functions
        // 0 - not visited
        // n - visited, next visited in n steps (n is possibly negative)
        // Hack - this list can store BOTH the queue AND the already visited nodes
        let mut visited = [0isize; Self::N.pow(Self::N as u32)];
        let vl = visited.len() as isize;
        let mut start = vl;
        let mut queue;
        let pri = Self::pass_right().unary_index();
        get_mut_unsafe![visited, pri] = vl - pri as isize;
        queue = pri as isize;
        while queue != vl {
            let ui = queue as usize;
            // Pop from queue
            queue += get_unsafe![visited, ui];
            // Redundant, visited[ui] will get overwritten anyway
            // visited[ui] = 0;
            // Start iterating
            let mut ptr = start;
            while ptr != vl {
                let add = [
                    self.unary_compose(ptr as usize, ui as usize),
                    self.unary_compose(ui as usize, ptr as usize),
                ];
                // Push to queue
                for x in add {
                    if get_unsafe![visited, x] == 0 {
                        get_mut_unsafe![visited, x] = queue - x as isize;
                        queue = x as isize;
                    }
                }
                // Next visited function
                ptr += get_unsafe![visited, ptr as usize];
            }
            // Self-composition
            let x = self.unary_compose(ui as usize, ui as usize);
            // Push
            if get_unsafe![visited, x] == 0 {
                get_mut_unsafe![visited, x] = queue - x as isize;
                queue = x as isize;
            }
            // Add self
            get_mut_unsafe![visited, ui] = start - ui as isize;
            start = ui as isize;
        }
        visited.iter().all(|&x| x != 0)
    }
    fn is_functionally_complete(&self) -> bool
    where
        [(); Self::N.pow(Self::N as u32)]:,
    {
        let n = Self::N;
        let unary = self.is_unary_complete();
        if !unary {
            return false;
        }
        // Loneliness
        let mut pass = false;
        'a: for a in 0..n {
            for b in 0..a {
                for c in 0..n {
                    for d in 0..c {
                        let ac = self.eval(a, c);
                        let bc = self.eval(b, c);
                        let ad = self.eval(a, d);
                        let bd = self.eval(b, d);
                        let mut arr = [ac, bc, ad, bd];
                        arr.sort_unstable();
                        if arr[0] != arr[1] || arr[2] != arr[3] {
                            pass = true;
                            break 'a;
                        }
                    }
                }
            }
        }
        if !pass {
            return false;
        }
        // k+1 - completeness
        self.is_k_complete()
    }
    fn is_k_complete(&self) -> bool {
        let n = Self::N;
        let mut counter = vec![0; n];
        let mut k = 2;
        'b: while k < n {
            let cmb: Vec<_> = (0..n).combinations(k).collect();
            for l in &cmb {
                for r in &cmb {
                    for &a in l {
                        for &b in r {
                            counter[self.eval(a, b)] = 1;
                        }
                    }
                    let s = counter.iter().sum();
                    counter.fill(0);
                    if s > k {
                        k = s;
                        continue 'b;
                    }
                }
            }
            return false;
        }
        true
    }
    fn from_index(i: usize) -> Self;
    fn to_index(&self) -> usize;
    // If this function is minimal in it's equivalence class, returns the number of functions in the equivalence class
    // Else returns zero
    fn cifminelse0(&self) -> usize;
    // Pretty prints this gate
    fn pretty_print(&self) {
        println!("Gate i={}", self.to_index());
        let pad = " ";
        let space = " ";
        print!("{}?", pad);
        for b in 0..Self::N {
            print!("{}{}",space,b);
        }
        println!();
        for a in 0..Self::N {
            print!("{}{}", pad, a);
            for b in 0..Self::N {
                print!("{}{}",space,self.eval(a,b));
            }
            println!();
        }
    }
}
// Four bit functions
impl Function for u32 {
    const N: usize = 4;
    fn impl_eval(&self, a: usize, b: usize) -> usize {
        *self as usize >> (a * 4 + b) * 2 & 3
    }
    fn pass_left() -> u32 {
        0xFFAA5500
    }
    fn pass_right() -> u32 {
        0xE4E4E4E4
    }
    fn compose(&self, l: &Self, r: &Self) -> Self {
        let mut ret = 0;
        for a in (0..4).rev() {
            for b in (0..4).rev() {
                ret <<= 2;
                ret |= self.eval(l.eval(a, b), r.eval(a, b));
            }
        }
        ret as u32
    }
    fn random<R: Rng>(rng: &mut R) -> Self {
        rng.gen()
    }
    fn to_index(&self) -> usize {
        *self as usize
    }
    fn from_index(i: usize) -> Self {
        i as u32
    }
    fn unary_index(&self) -> usize {
        *self as usize & 255
    }
    fn from_unary_index(i: usize) -> Self {
        let r = i | (i << 8);
        let r = r | (r << 16);
        r as u32
    }
    fn cifminelse0(&self) -> usize {
        const PERM4: [[usize; 4]; 24] = [
            [0, 1, 2, 3],
            [0, 1, 3, 2],
            [0, 2, 1, 3],
            [0, 2, 3, 1],
            [0, 3, 1, 2],
            [0, 3, 2, 1],
            [1, 0, 2, 3],
            [1, 0, 3, 2],
            [1, 2, 0, 3],
            [1, 2, 3, 0],
            [1, 3, 0, 2],
            [1, 3, 2, 0],
            [2, 0, 1, 3],
            [2, 0, 3, 1],
            [2, 1, 0, 3],
            [2, 1, 3, 0],
            [2, 3, 0, 1],
            [2, 3, 1, 0],
            [3, 0, 1, 2],
            [3, 0, 2, 1],
            [3, 1, 0, 2],
            [3, 1, 2, 0],
            [3, 2, 0, 1],
            [3, 2, 1, 0],
        ];
        fn flipped(f: u32) -> u32 {
            f & 0xc0_30_0c_03
                | f << 6 & 0x30_0c_03_00
                | f << 12 & 0x0c_03_00_00
                | f << 18 & 0x03_00_00_00
                | f >> 6 & 0x00_c0_30_0c
                | f >> 12 & 0x00_00_c0_30
                | f >> 18 & 0x00_00_00_c0
        }
        fn permuted(f: u32, p: &[usize; 4]) -> u32 {
            let mut invp = [0, 0, 0, 0];
            for i in 0..4 {
                invp[p[i]] = i as u32;
            }
            let mut r = 0;
            for a in (0..4).rev() {
                for b in (0..4).rev() {
                    r <<= 2;
                    r |= invp[f.eval(p[a], p[b])];
                }
            }
            r
        }
        let mut m = vec![*self; 48];
        m[1] = flipped(*self);
        if m[1] < *self {
            return 0;
        }
        for (i, p) in PERM4.iter().enumerate().skip(1) {
            m[i * 2] = permuted(*self, p);
            m[i * 2 + 1] = flipped(m[i * 2]);
            if m[i * 2] < *self || m[i * 2 + 1] < *self {
                return 0;
            }
        }
        m.sort_unstable();
        m.dedup();
        m.len()
    }
}
impl<const N: usize> Function for [[usize; N]; N] {
    const N: usize = N;
    fn impl_eval(&self, a: usize, b: usize) -> usize {
        self[a][b]
    }
    fn pass_left() -> Self {
        let mut r = [[0; N]; N];
        for a in 0..N {
            r[a] = [a; N];
        }
        r
    }
    fn pass_right() -> Self {
        let mut d = [0; N];
        for b in 0..N {
            d[b] = b;
        }
        [d; N]
    }
    fn compose(&self, l: &Self, r: &Self) -> Self {
        let mut ret = [[0; N]; N];
        for a in 0..N {
            for b in 0..N {
                ret[a][b] = self.eval(l.eval(a, b), r.eval(a, b));
            }
        }
        ret
    }
    fn random<R: Rng>(rng: &mut R) -> Self {
        let mut r = [[0; N]; N];
        for a in 0..N {
            for b in 0..N {
                r[a][b] = rng.gen_range(0..N);
            }
        }
        r
    }
    fn from_index(mut i: usize) -> Self {
        let mut r = [[0; N]; N];
        for a in 0..N {
            for b in 0..N {
                r[a][b] = i % N;
                i /= N;
            }
        }
        r
    }
    fn to_index(&self) -> usize {
        let mut r = 0;
        for a in (0..N).rev() {
            for b in (0..N).rev() {
                r *= N;
                r += self[a][b];
            }
        }
        r
    }
    fn from_unary_index(mut i: usize) -> Self {
        let mut x = [0; N];
        for n in 0..N {
            x[n] = i % N;
            i /= N;
        }
        [x; N]
    }
    fn cifminelse0(&self) -> usize {
        let flip = |i: &Self| i.compose(&Self::pass_right(), &Self::pass_left());
        let permute_self = |p: &[usize]| {
            let mut inv = [0; N];
            for i in 0..N {
                inv[p[i]] = i;
            }
            let mut ret = [[0; N]; N];
            for a in 0..N {
                for b in 0..N {
                    ret[a][b] = inv[self.eval(p[a], p[b])];
                }
            }
            ret
        };
        let mut m = vec![*self; (1..=N).product::<usize>() * 2];
        m[1] = flip(self);
        if m[1] < *self {
            return 0;
        }
        for (i, p) in (0..N).permutations(N).enumerate().skip(1) {
            m[i * 2] = permute_self(&p);
            m[i * 2 + 1] = flip(&m[i * 2]);
            if m[i * 2] < *self || m[i * 2 + 1] < *self {
                return 0;
            }
        }
        m.sort_unstable();
        m.dedup();
        m.len()
    }
}
#[test]
fn test_pass() {
    fn test_pass_g<F: Function + std::fmt::Debug>() {
        let l = F::pass_left();
        let r = F::pass_right();
        for a in 0..F::N {
            for b in 0..F::N {
                assert_eq!(l.eval(a, b), a, "a is {}, b is {}, l is {:?}", a, b, l);
                assert_eq!(r.eval(a, b), b, "a is {}, b is {}, r is {:?}", a, b, r);
            }
        }
    }
    test_pass_g::<u32>();
    test_pass_g::<[[usize; 3]; 3]>();
    test_pass_g::<[[usize; 4]; 4]>();
}
#[test]
fn test_compose() {
    use rand::rngs::SmallRng;
    use rand::SeedableRng;
    fn test_compose_g<F: Function + std::fmt::Debug>() {
        let mut rng = SmallRng::seed_from_u64(42);
        let repeats = 1_000_000;
        for _ in 0..repeats {
            let f = F::random(&mut rng);
            let l = F::random(&mut rng);
            let r = F::random(&mut rng);
            let c = f.compose(&l, &r);
            let a = rng.gen_range(0..F::N);
            let b = rng.gen_range(0..F::N);
            assert_eq!(
                c.eval(a, b),
                f.eval(l.eval(a, b), r.eval(a, b)),
                "a:{} b:{} f:{:?} l:{:?} r:{:?} c:{:?}",
                a,
                b,
                f,
                l,
                r,
                c
            );
        }
    }
    test_compose_g::<u32>();
    test_compose_g::<[[usize; 3]; 3]>();
    test_compose_g::<[[usize; 4]; 4]>();
}
#[test]
fn test_low_effort_discard() {
    use rand::rngs::SmallRng;
    use rand::SeedableRng;
    fn test_low_effort_discard_g<F: Function + std::fmt::Debug>()
    where
        [(); F::N.pow(F::N as u32)]:,
    {
        let mut rng = SmallRng::seed_from_u64(42);
        let repeats = 1_000;
        for _ in 0..repeats {
            let f = F::random(&mut rng);
            assert!(
                !f.low_effort_discard() || !f.is_functionally_complete(),
                "Discarded f:{:?}",
                f
            );
        }
    }
    test_low_effort_discard_g::<u32>();
    test_low_effort_discard_g::<[[usize; 3]; 3]>();
    test_low_effort_discard_g::<[[usize; 4]; 4]>();
}
#[test]
fn test_specifics() {
    // Using the examples from the blog post + = 0, 0 = 1, 1 = 2
    let triplets = [
        ([[0; 3]; 3], false),
        ([[1; 3]; 3], false),
        ([[2; 3]; 3], false),
        ([[0, 1, 2]; 3], false),
        ([[0; 3], [1; 3], [2; 3]], false),
        ([[0, 1, 2], [1, 1, 2], [2, 2, 2]], false), // Min
        ([[0, 1, 2], [0, 1, 1], [0, 0, 0]], false), // Imp
        ([[0; 3], [1; 3], [0; 3]], false),          // Imp composition
        ([[2, 0, 0], [0, 0, 0], [0, 0, 1]], true),  // Tand
        ([[2, 1, 1], [1, 0, 1], [1, 1, 1]], true),  // Modified Tand
        ([[2, 2, 2], [2, 0, 2], [2, 2, 1]], true),  // Modified Tand
        ([[1, 0, 0], [0, 2, 0], [0, 0, 0]], true),  // Modified Tand
        ([[1, 1, 1], [1, 2, 1], [1, 1, 0]], true),  // Modified Tand
        ([[1, 2, 2], [2, 2, 2], [2, 2, 0]], true),  // Modified Tand
        ([[2, 0, 1], [0, 0, 0], [1, 0, 1]], true),  // Pointy Tand
        // Experimentally found
        ([[0, 2, 0], [0, 0, 0], [0, 0, 0]], false),
        ([[1, 2, 0], [0, 0, 0], [0, 0, 0]], true),
        // From the post
        ([[2, 0, 1], [0, 0, 0], [2, 2, 0]], true),
        ([[2, 0, 1], [1, 0, 0], [2, 2, 0]], true),
        ([[2, 0, 1], [2, 0, 0], [2, 2, 0]], true),
        ([[1, 0, 0], [1, 0, 2], [2, 2, 1]], false),
    ];
    for (f, r) in triplets {
        println!("{:?} -> {}", f, r);
        assert_eq!(f.is_functionally_complete(), r);
    }
}
fn smart_method<F: Function>()
where
    [(); F::N.pow(F::N as u32)]:,
{
    let amount = F::N.pow(F::N as u32 * F::N as u32);
    (0..amount).into_par_iter().for_each(|i| {
        let n = F::from_index(i);
        let m = n.cifminelse0();
        if m == 0 {
            return;
        }
        let r = n.low_effort_discard() || !n.is_functionally_complete();
        if !r {
            TC.fetch_add(m, Relaxed);
        } else {
            FC.fetch_add(m, Relaxed);
        }
    });
    CONT.store(false, Release);
}
static TC: AtomicUsize = AtomicUsize::new(0);
static FC: AtomicUsize = AtomicUsize::new(0);
static CONT: AtomicBool = AtomicBool::new(true);
// Change this to [[usize;3];3] for n=3
type T = u32;
fn main() {
    // Fluff
    let start = std::time::Instant::now();
    let amount = T::N.pow(T::N as u32 * T::N as u32);
    let pb = ProgressBar::new(amount as u64);
    pb.set_style(
        ProgressStyle::default_bar().template(
            "{wide_bar:.green/red}\n{pos}/{len} - {percent}% - {per_sec} - {eta} - {msg}",
        ),
    );
    let t = thread::spawn(move || {
        while CONT.load(Acquire) {
            thread::sleep(std::time::Duration::from_millis(100));
            let t = TC.load(Relaxed);
            let f = FC.load(Relaxed);
            pb.set_position(t as u64 + f as u64);
            pb.set_message(format!(
                "Ratio at {}/{} = {:.5}",
                t,
                t + f,
                t as f64 / (t as f64 + f as f64)
            ));
            pb.tick();
        }
    });
    // Call
    smart_method::<T>();
    // Fluff
    match t.join() {
        Err(_e) => {
            println!("Failed to join ui-thread :/");
        }
        _ => {}
    };
    let t = TC.load(Relaxed);
    let f = FC.load(Relaxed);
    println!(
        "Ended with {} universal gates and {} non-universal gates. Ratio is {}/{} = {:.5}",
        t,
        f,
        t,
        t + f,
        t as f64 / (t as f64 + f as f64)
    );
    println!("Took {:?}", start.elapsed());
    // Remember to reset TC, FC and CONT if you want to call again
}
```
## cargo.toml
```
[dependencies]
itertools = "0.10"
rand = {version="0.8", features=["small_rng"]}
rayon = "1.5"
indicatif = {version = "0.16", features = ["rayon"]}
```
Takes around 20 minutes on my machine™. The actual code is not that interesting, maybe with the exception of the unary completeness check, which uses a buffer with two non-overlapping linked lists, to form a very efficient deduplicated queue + set data structure.
# n=5?
For `n=5` you basically need a new idea (or a lot of computing power). Checking each gate individually becomes impractical (if you could do it in one cycle, it would still take two years). If someone manages to calculate `n=5`, I'm happy to forward Bubblers +500 bounty to them.
[Answer]
# [Python](https://www.python.org), 1648 bytes
```
lambda n:(K:=G(n))and n**(n*n)-sum(U([O({((x,y),):C(K)if x!=y else[(B[x],)]for x,y in W(K,K)},j)for j in[{((x,y),):C(S)for x,y in W(S,S)}for j in G(1,n)for S in combinations(K,j)]+[{Z(W(X,V)):sum(([*W(*[S]*len(X)*len(V))]for S in P if not S==X==V),[])for X,V in W(P,P)}for P in R(K)if 1<len(P)<n]+sum(([{Z(zip(X,V[i:]+V[:i])):[S[i:]+S[:i]for S in P for i in G(p)]for X,V,i in W(P,P,G(p))}for P in R(K)if set(map(len,P))=={p}]for p in G(2,n+1)if n%p<1),[])])for B in W(*[K]*n))or 1
from itertools import*
import math
from collections import*
W=product
Z=tuple
H=Counter
G=range
def T(l):
	d={(i,):C(c)for i in G(len(l[0]))if{*H(x[:i]+x[i+1:]for x in l).values()}=={len((c:=H(x[i]for x in l)))}!=2>len({*c.values()})};r=Z({*G(len(l[0]))}-{*sum(d,())})
	if r:d[r]=[*{Z(x[i]for i in r)for x in l}]
	return d
E=lambda s:math.prod(map(len,s.values()))
C=lambda l:[(i,)for i in l]
def O(s,o):
	q={i:[dict(zip(I,l))for l in s[I]]for I in s for i in I};D={}
	for I in o:
		d=[dict(zip(I,l))for l in o[I]];S=set()
		for i in I:S|={*q[i][0]};d=[h for a,b in W(d,q[i])for h in[a|{}]if all(b[j]==h.setdefault(j,b[j])for j in b)]
		for i in S:q[i]=d
	for d in q.values():
		if[]==d:return{Z(q):d}
		t=Z(sorted(d[0]));D[t]=[Z(l[j]for j in t)for l in d]
	for i in[*D]:
		if d:=T(D[i]):D.pop(i);D|={Z(i[k]for k in j):d[j]for j in d}
	return D
def M(l):
	if l:
		for a,b in M(l[1:]):yield(l[0],)+a,b;yield a,(l[0],)+b
	else:yield(),()
R=lambda l:[(a+(l[0],),*t)for(a,b)in M(l[1:])for t in R(b)]if l else[()]
def U(L):
	L=[s for s in L if E(s)]
	if[]==L:return 0
	x=L[0];L2=[];L1=[]
	for i in L:
		if i!=x!=E((j:=O(i,x)))<E(i):L1+=[i];L2+=[j]
	return E(x)+U(L1)-U(L2)
```
[Attempt This Online! (only up to n=4)](https://ato.pxeger.com/run?1=dVXBbuM2EL0VsL6CORQhZSZYu1lklw4LdOMkG8RFglU2CcLyIFtyQq8sKpJcONXqS3oJULTf0F9pv6YzoiK7KHqRzOHjezNvRvSvf2TP5aNNX347lT_9virne-_-_iZLwuU0Ckkq6IWQZzRlLEwjkvo-Tf2U7RWrJf1M1SWtKF3zZ8aZOKYXzMzJekc-kzgpYkU_qLXmTM9tTgBDTEpu6QW_YDVfMAwuIKS2CQL2L2zAA1a_AskZHfC0AQS4nNnl1KRhaWxaAOuC6b6q7uktveM3jAlMkCr_lvoq0H4Sp_SONS_Y1B3JFYGMU1uSQMo7KW8YV7qRABKXwxW_cjlc4fqTq3FwhExX7CjVfScEyr-YDLWVEbp_o4TRkIUKmmWAyy1R_GlcTZnLBg5y0ylyjP9XtohLugwzCuKQFZOyyurmdOa4hjztDxCYfpsdDZpaXDUfHLOvLjQ0j0Fk4M1zuySmjPPS2qQgZpnZvPQ99ybLsHx0kJlNknjW2NyBbmWW22g1K717Wa6yJPY-ymO7SoHNO5N5mD7EXhTPyTVNmPB6kayoaRo8Y1ulo4eJegM-mXnlf6RrdKm_VqY_EG5oEJew_Z_DZBUXlNVQMR6iMyERbrZR4NeOHH6P-5U_25xh9SiX9xDbFqz3Kh8bF3EKC-b1wLRcRCrXUvnQy1fuJtOcbWRq7fXyuFzlKYm8E9l-JYVAv_bRlK5BRZcCY97xKzIRCq3ouBPdGHVJC27RqidZGaEiMyubeTrnUBiCEwQX6lw3aZ03q80Yndejsaxqr9dtWuAC3_-PySLTKJA4UVB9b8Mkgq-y8p-gfjCqHgHFY6MT8qkboojjZkP1iN9v-LWqNdgXJgmdqoWW8nEfaKGqcJWUdMEx2H3uZMr0tl4gkE1GLvUIQ0-dcViDmSugjIQzHVrzxEQElfZKaGoB0xhHNGp6OhqrEtp3Dy1e6E6v3BQdaa8TVv5YO3oSCXlNx1iTGO9nNqMGqMCEe2rUl4boC55egO42MSbRTsK46eGPbtiBMRFtia1psKNgpJl4NnESNSPIWR82R00AYK-xqdfDq7MFMphO79P26IT9Fsn9pjAKJGxLAUVLd2OA0ZhKexUzN2ef6QRznEjlpqdA7AQvwRNaYGec3ZPWbvLG663lBBRHk6FU8BzAc-MimbQemh0J9_4JpQshL2HA1zD0RyfgpJgM-hK8hfPwXmy-nhO6Zn1IZ8D24Dlk7o_nrz8frI2IJNf5KvZQJuUPKIQjrAZ8yL_jB5q7X4eHB_z9wfDd-8O3b4ca67Jw8hT-q5ouWLIjyQMmmOUmLel895RWac0AU9kaN6uHehenv9U8DcGqtgHdKcs84EIExF6JfkgSUsZFWZAsLIo42tlt8395ce9_AA)
Golfed version of my code for `n=5`.
# Explanation
This algorithm relies almost entirely on the fundamental theorem from the paper mentioned in the comments: [On n-Valued Sheffer Functions by
Roy O. Davies](https://doi.org/10.1002/malq.19790251903).
The fundamental theorem lists three conditions that are both sufficient and necessary for a given n-ary logic gate to be complete.
Let \$U=\{1,2,...,n\}\$
A function \$f:U\times U\rightarrow U\$ is complete if and only if none of the following hold:
* There exists a nonempty set \$S\subsetneq U\$ such that for all \$x,y\in S\$, \$f(x,y)\in S\$ (e.g. if `f(0,0)==0` then `f` must be incomplete)
* There exists a nontrivial equivalence relation on \$U\$ where if \$x\_1\sim x\_2\$ and \$y\_1\sim y\_2\$, then \$f(x\_1,y\_1)\sim f(x\_2,y\_2)\$ (that is, you can view \$f\$ as acting on the equivalence classes). (Here an equivalence relation is trivial if \$a\$ and \$b\$ are either always related, or related only when \$a=b\$).
* There exists a permutation \$\pi\$ on \$U\$ consisting of \$\frac{n}{p}\$ cycles where \$p\$ is a prime factor of \$n\$ such that for all \$x,y\in U\$: \$\pi(f(x,y))=f(\pi(x),\pi(y))\$. (For example, if \$(a,b,c)\$ is a cycle and \$(x,y,z)\$ is a cycle, then \$(f(a,x), f(b,y), f(c,z))\$ is a cycle)
So we can count the number of complete gates by counting the number of incomplete gates and subtracting that from the total number of gates. The number of incomplete gates can be calculated by considering each reason a gate could be incomplete individually and then combining the results using inclusion-exclusion.
For example, consider the set of gates that fail the first test where \$n=4\$ and \$S=\{0,1\}\$. It's pretty easy to see that \$f(0,0), f(0,1), f(1,0), f(1,1)\$ each have 2 choices (0 and 1), and the other 12 outputs have 4 choices each.
The other two cases are a little more complicated because values are no longer independent, so the representation I settled on is to partition \$U\times U\$ into pairwise disjoint sets \$S\_i\$, and for each such set list the outputs that the whole set can take.
The set of gates where \$f(0,0)=f(1,1)\$ would be represented by
* \$(f(0,0),f(1,1))\in\{(0,0),(1,1),(2,2),(3,3)\}\$
* \$f(0,1)\in U\$
* \$f(0,2)\in U\$
* ...
With this representation, calculating the size of the set is easy, finding the intersection of two sets is relatively easy (the number of options can easily blow up, so the code tries to compress the representation when possible), and the rest of the cases can be represented.
There is one other major optimization as well:
If we fix \$f(i,i)\$ for \$i\in U\$ then the number of cases to consider is dramatically reduced (which is helpful because PIE is exponential in the number of cases) at the cost of having \$n^n\$ disjoint cases. Moreover, if you consider the directed graph induced by \$f(i,i)\$, isomorphic graphs correspond to the same number of complete gates, reducing the amount of work for \$n=5\$ by nearly 250 times! In fact, for \$n=5\$ it actually also fixes \$f(0,1)\$ which fixes one very slow case. (this optimization (among others) is removed in the golfed code).
Full code:
```
import itertools
from collections import Counter
def length(x):
    # counts the number of functions in the set
    z = 1
    for i in x.values():
        z *= len(i)
    return z
def intersect(A, B):
    # finds the intersection of two sets
    q = {}
    for I, L in A.items():
        d = [dict(zip(I,l)) for l in L]
        for i in I:
            q[i] = d
    for I, L in B.items():
        d = [dict(zip(I,l)) for l in L]
        seen = set()
        for i in I:
            d1 = q[i]
            if len(d1) == 0: return {tuple(q):[]}
            if min(d1[0]) in seen: continue
            seen.update(d1[0])
            d2 = []
            for a in d:
                for b in d1:
                    h = a.copy()
                    for j, k in b.items():
                        if h.setdefault(j, k) != k:
                            break
                    else:
                        d2.append(h)
            d = d2
            if len(d) == 0: return {tuple(q):[]}
        for i in seen:
            q[i] = d
    D = {}
    for i, d in q.items():
        if len(d) == 0:
            return {tuple(q):[]}
        elif i == min(d[0]):
            t = tuple(sorted(d[0]))
            D[t] = [tuple(l[j] for j in t) for l in d]
    return simplify(D)
def simplify(self):
    # simplifies the set representation
    if length(self) == 0:
        self = {(j,):[] for i in self for j in i}
    else:
        l0 = length(self)
        for i in list(self):
            l = self[i]
            d = split_up(l)
            if len(d) > 1:
                for j, o in d.items():
                    self[tuple(i[k] for k in j)] = o
                self.pop(i)
        assert l0 == length(self)
    return self
def split_up(l):
    # helps to simplify a list of assignments
    # e.g. [(0,0),(0,1),(0,2),(1,0),(1,1),(1,2)] can be expressed as (0,1)x(0,1,2)
    d = {}
    for i in range(len(l[0])):
        c = Counter(x[i] for x in l)
        if len(set(c.values())) == 1:
            q = Counter(x[:i]+x[i+1:] for x in l)
            if set(q.values()) == {len(c)}:
                d[i,] = [(j,) for j in c]
    remaining = tuple(set(range(len(l[0])))-set(sum(d,())))
    if remaining:
        d[remaining] = list({tuple(x[i] for i in remaining) for x in l})
    return d
def basic(l):
    # returns the representation for the set of functions f where f(i,i) = l[i]
    r = range(len(l))
    a = [(i,) for i in r]
    return {((x,y),):a if x != y else [(l[x],)] for x in r for y in r}
def nonempty_proper_subsets(l):
    for i in range(1, len(l)):
        yield from itertools.combinations(l, i)
def splits(l):
    if len(l) == 0:
        yield (), ()
        return
    x, *l1 = l
    for a,b in splits(l1):
        yield (x,)+a, b
        yield a, (x,)+b
def partitions(l):
    if len(l) == 0:
        yield ()
        return
    x, *l1 = l
    for a, b in splits(l1):
        a += (x,)
        for t in partitions(b):
            yield (a,) + t
def special_permutations(n, p):
    for P in partitions(list(range(n))):
        if set(map(len, P)) == {p}:
            d = {}
            options = [S[i:]+S[:i] for S in P for i in range(p)]
            for S1 in P:
                for S2_ in P:
                    for i in range(p):
                        S2 = S2_[i:] + S2_[:i]
                        d[tuple(zip(S1, S2))] = options
            yield d
def factors(n):
    p = 2
    s = set()
    while n > 1:
        if n%p == 0:
            s.add(p)
            n //= p
            while n%p == 0:
                n //= p
        p += 1 + p%2
    return s
def canonize_(l, d, r1, r2):
    options = []
    len0 = len(d)
    r = r1 or r2
    if len(r) == 0:
        return ()
    for i in r:
        while i not in d:
            d.append(i)
            i = l[i]
        p = (len(d)-len0, d.index(i))
        options.append(p + canonize_(l, d,
            [j for j in r1 if j not in d],
            [j for j in r2 if j not in d]
        ))
        while len(d) > len0: d.pop()
    return max(options)
def canonize(l):
    r2 = set(l)
    r1 = set(range(len(l))) - r2
    return canonize_(l, [], r1, r2)
def union_length(L):
    L = [s for s in L if length(s) > 0]
    if len(L) == 0: return 0
    x = max(L, key=length)
    L2 = []
    L1 = []
    for i in L:
        if i is not x:
            ix = intersect(i,x)
            if length(ix) < length(i):
                L1.append(i)
                L2.append(ix)
    return length(x) + union_length(L1) - union_length(L2)
def magic(n):
    if n == 1: return 1
    l = list(range(n))
    Z = []
    for S in nonempty_proper_subsets(l):
        allowed = [(i,) for i in S]
        Z.append({
            ((x, y),): allowed
                for x in S
                for y in S
        })
    for P in partitions(l):
        if len(P) == n: continue
        if len(P) == 1: continue
        options = {}
        d = {}
        for S1 in P:
            for S2 in P:
                if S1 is S2:
                    option = []
                    for S in P:
                        if S is S1: continue
                        option.extend(itertools.product(S, repeat=len(S1)*len(S1)))
                    d[tuple(itertools.product(S1, S1))] = option
                else:
                    ab = len(S1) * len(S2)
                    if ab not in options:
                        option = []
                        for S in P:
                            option.extend(itertools.product(S, repeat=ab))
                        options[ab] = option
                    d[tuple(itertools.product(S1, S2))] = options[ab]
        Z.append(d)
    for p in factors(n):
        Z.extend(special_permutations(n, p))
    if True:
        cache = {}
        out = 0
        for ll in itertools.product(range(n), repeat=n):
            if any(i==j for i,j in enumerate(ll)): continue
            c = canonize(ll)
            if c in cache:
                out += cache[c]
            else:  
                standard = basic(ll)
                q = 0
                rs = max(n-4, 0)
                for t in itertools.product(range(n), repeat=rs):
                    x = standard
                    if rs:
                        x = intersect(x, {
                            ((0,i+1),):[(t[i],)]
                            for i in range(rs)
                        })
                    q += length(x) - union_length([intersect(x,i) for i in Z])
                cache[c] = q
                out += q
        return out
for i in [1,2,3,4,5]:
    print(magic(i))
```
[Attempt This Online! (takes about a minute for n=5)](https://ato.pxeger.com/run?1=nVlLb9tGEL77H_S2RRGAG9OKqKZAIVQFkuYSQIcA6imEYKzElbU2taT5qKkYvvZP9JJL-6P6azqzL3L5cIzyIEu7M7PfvGfpv_7Jz9Uxk1-__l1Xh6uf__3uT3HKs6IiouJFlWVpeXEoshPZZ2nK95XIZEkMxW9ZLYHo4iLhB5JyeVMdg4YuLwg8PwAD7JakOnIi69OOFyQ7kEMtrQyptkpeKfovZEUi9e2QFUTgdjP7g6U1LwMjUpO9XuFRgaBqreBVXUjyRWMQCKcElMG7kLx3SA5CJhqIIwAICKd6yBBBqQjvAcLjk8PwMSRrhPFuBpY4eSgSoIwTAed8EXnwMUwpVSwp0q-3js6p8rHlVSfFYgsiksFZ7___WSXnEkhBm4B-E0ASASmi8FaF8mKQRJSsVmS-tNZ9rOo85cE9Xcbbpz7HSSBHPN9SPAdRLMH1shKy5h4tbs3qPGEVNww-pAUq6gNC-AzFJj58u7dTe9FwE58jyGOzfZafAzpKgBJuQ3KHQnZDw_cf0PU4A_NCoLE6rQJkpeT7Fbmb5sFnV3B2N0rB05JP8yaLGctzLpPg2LMUhs5i1HMvcpwLCuWs6cD84OeDCOFk4Lofmqp3vCfyWSg8BVaBXCqMMCh87gowaM4SKg5PNI1vjw9xhaBjTZfGt1vtWlViOsmSbLslo4QiBqefgw9UFw-3UPL04GqHWRW8tOUKBOQFL7msGJaRi9YAWP8Uc88MuIa2hJBB9bsOgA2HVWi7-FGRzsnKkz10YyrKqgvasapykB76aY7xU4JS1XWdBymdCKRfyUhemZTJlDmfTxl1snaJiO-00irVbil6KxswIcMsz3Jb2vFhZQlNSBlhxArWk7BgXNhqZf135GkOrsuce6GgoMGw-oN0cSNP4MnSUPPZzYzEwTyc0xA-I_W5gM9IrURqJYKVLdkzqBqc8AaDoeQJSCOKpcFPILmwtu7mEBqgYPIG4hTMnKpgbu23B2rTVIMGUxF5GuVj2k82LPR71yGpirmex-49cUuxvQShl9FyXK6RjXLvW7ko9hHP29OnoaOTWIQq9zC020je20w7MSGFvGmTGIT31adXuFrWpyAJURFqM8qxd5ph7BbxWBX5prI4e2kbWzLa0fXJC5tEx8yOlWLfCRi9q7Pdz3QlyRYBb5Y5kIcjLzg5BCIUFIHZnCvgR0dhoxxTJhPGZBqwV5weg6AJzxTKBUNLNNhmzqoyAF8aN9uQdpxYqK9n9fVJayUzyU95db7OiyznxXVZ73DOcXr2gjEKicHX2voseJoQNfu5URDa6WknpDIHCAuJoJ3Ma8WbGE37lVDLDGhIOj1Z66x-NiF5neJskjqULFRd3h4QDRCCpeglC8mutw5LamunEeasqITB_UKUL0ZIJiEycrlSMLyyXSF5B9CuV0MNAAYBckkqa2G-Fyy9Bm-e6sp4QIYk73j0U0-uyg_tYUmp37Ax6U4sx7gMySeT6Hkvxzvlyz5ZrmMeIngTi-X2coOVRR2_weM_9WMrp8OhbhMp0vEGs1lcT-xaCk_69AS1wZESpCFMMCR-W_Z6oaesaVc4Z28gIzYLqluV1njEQ6aEHNi-ygrwhsGSA5Oez0pvIn84ihQuQ35rBVfIV_nI4FTOWJKAft6iJG_erEjurRmxo0LGeHKMyQgMkr9aeH1UawONLZPiC7_G_E5CUoApioVRreN9bUgIHzOhwNDQlryIgJuKRTfNin6amWONcVq_thRaMwHVrBq5BSR2Pha9JtYtwNYfgUZ4hXhDnF1kwhvgbFmNalYoWKlvCu-U-Lbtd6AuKHnrcG6fIV30SB1lB4rW2w1iiHkJmHE48lrYiTWBgU1957kaVyxMCJpOX0Tmt9eVKLmy3jKiPdXjrQsDfUwt4cxrM5CtzVFrjIpSqaqu9-vuWIx6zLfdeFj3LitzXV9BCKq1htsVP680u4a-7twQ11H73UXO2ssqWCqVmRs_bAQe0b4pEGEzNgIjaNFQ8ov7MVJo1tFECGq4brPxvObelECE-YaM0A_-krX4id3AmCLbziX1wGeF6ncnqR2KXNFXy599Y6k6_a0JAR-WptkDT4bTyqYN289Wy0fPAjjAEDXBWCmjtV7NL5vRrbO_9dSWiUGbG15FP6ngGnsV4RFEIwRtiet0vl4jnOxhun9NtC84GplKoBhvWvro4XsQT_ozvdEeos4Y0238wBlvKhWnbsyDmEhqSI5NiDMwZxUmIjRF-tr8peOvVGwLHZGEDTXqNtSBgOnXIWxnWgxIIK_1t8U4AlAfqE11Na6cNta3DI7PS4zeynqJLdluwnytmDJmu2dMhc83bO0PLyhumLJJm1M5atgfZDSx0Wh6AHU3tt-LuuPBPdvDZcnLm6zGFzpzL49S9XJmqIUtYc5ssleB0dXyHIjV6ta8oFIdlsv6xAt8y5jifWY8CfCm3fbK4S14r66xqMDQ6ajE5Urvxns_cFQMEzLgKSsmE1ZgGTEXznQYAveebexTlKYlyqu3IZkP2dyN4gU2LMqJeRmbogU5lVjFM5nkN1Wo_Y_PJksQzENxGWF_iIMKhrWQTmegVbIz9oMek-RP41v36LS29_YabdwFLzqt7vN2KM66Hl-iT4VHu2E6NKxfXDixcRQuwh_Dt-FPW3NpKABBoNs8jqX63zHmvzL2vzP_AQ)
As an aside (also from the paper), another sufficient and necessary condition for completeness is if for each distinct \$a,b,c\in U\$, \$f\$ generates a function \$g:U\rightarrow U\$ such that \$f(a)=a\$ and \$f(b)=c\$, which means that unary completeness implies completeness.
[Answer]
# [JavaScript (Node.js)](https://nodejs.org), 194 bytes
```
f=(n,c=[],k=0)=>k-n?1/c[n*n]?(F=>g=s=>s[n**n*8]?eval(`try{for(A=n;A--;)for(B=n;B--;)eval(s)+~(A<B?A:B)%n&&G;1}catch{}`):[...'AB()F'].some(_=>g(s+_)))(t=>u=>c[t*n+u])``:f(n,c,k+1)+f(n,[...c,k]):0
```
[Try it online!](https://tio.run/##Hc5BboMwEAXQq3TTZAaDG3YVdIzsRXoIZAXkQtpCxxV2IlVRenWKu5sn/a8/n/21D275@I4F@7dhXUcCzh21Np/ogKSmgpvyybWcsW3gSOpMgVTYnHH2bJvh2s/QxeXnNvoFNHGti6LGBLPBJPxnAopf0C@m0ZXBR97tXuvy7vro3m/3DqtWSrnXBvC4tzL4rwFO2xYEcUJEiKQupFwbMxYXi11XjenNfBIlinSm@kaL1WGtH5zn4OdBzv4MI5SI6x8 "JavaScript (Node.js) – Try It Online")
If a function can express `f(x,y)=max(x,y)+1`, then it's universal:
* `f(x,x)=x+1`, therefore we can get `x+c` for any `c`
* `max(x,x+1,...,x+n-2)` is `n-1`, unless `x` is zero where the result is `n-2`
* `max(x,max(a,a+1,...,a+n-2,b,b+1,...,b+n-2)+1)` is `x` if either `a` or `b` is not 0, and `n+1` if `a` and `b` are both 0
* Set the value using this
The `eval` never meet a function that can't fit the exact length:
* The expression built with `AB()F` is always of odd length (`F(A)(B)`)
* One can write `A` as `((((A))))` to use even spaces
]