addhelp(is, "Is n a term of A048196? I.e. Has Binomial(n, floor(n/2)) the as number of unitary and non-unitary divisors?"); is(n) = {my(f = binom(n, n \ 2)); sqf = 1<<#f; ndiv = prod(i = 1, #f, f[i][2] + 1); sqf << 1 == ndiv } addhelp(binom, "List of prime factors of binomial(n, k) with their exponents.") binom(n, k) = {my(res = List()); forprime(p = 2, n, e = val(n, p) - val(n - k, p) - val(k, p); if(e > 0, listput(res, [p, e]~) ) ); res } addhelp(val, "exponent of prime p in n!"); val(n, p) = my(r=0); while(n, r+=n\=p);r fasteris(n) = { testbinom(n, n \ 2) } testbinom(n, k) = {my(count3 = 0); forprime(p = 2, n, e = val(n, p) - val(n - k, p) - val(k, p); if(e > 1, if(e == 3, count3++; if(count3 > 1, return(0) ) , return(0); ) ) ); count3 == 1 }