\\ A326489: Number of product-free subsets of {1..n} \\ ------------------------------------------------- \\ Basic Algorithm. \\ Recursively either exclude k or include k in the set. \\ k is processed in order 2..n \\ No checks are needed if excluding. \\ if including need to check that k is not a product of two existing inclusions. \\ This is ok up to about n = 30. A326489v1(n)={ my(lim=vector(n, k, sqrtint(k))); my(recurse(k, b)= k++; if(k > n, 1, my(t = self()(k, b)); for(i=2, lim[k], if(bittest(b,i) && !(k%i) && bittest(b,k/i), return(t))); t += self()(k, b + (1<= n/2 cannot affect future choices. \\ These values can then be processed independently of each other. \\ All information is known for k after completing k/2 since two is the smallest divisor we need to test for. \\ This is ok up to about n = 60. A326489v2(n)={ my(lim=vector(n, k, sqrtint(k))); my(accept(b,k) = for(i=2, lim[k], if(bittest(b,i) && !(k%i) && bittest(b,k/i), return(0))); 1); my(recurse(k, b)= my(m=1); for(j=max(2*k,n\2+1), min(2*k+1,n), if(accept(b,j), m*=2)); k++; if(k > n\2, m, my(t = m*self()(k, b)); if(accept(b,k), t += m*self()(k, b + (1< sqrt(n) cannot affect choices for values that don't have p as a prime divisor. \\ For example when n=100, the multiples of 11: 11, 22, 33, .. 99 can be processed after k reaches 9 and \\ this processing is independent of selections for other k > 9. \\ Furthermore we don't actually care what the prime factor is: only the number of multiples floor(n/p) is relevant. \\ For example when n=100, there are 10 primes p > 50. These immediately produce a multiplication factor \\ of 2^10 since they can be independently included or excluded in the set. \\ This was used to compute values up to n = 100. A326489v3(n)={ \\ initialise some data \\ lim[k] is the integer square root of k. \\ h[k] is the number of primes p > sqrt n with n\p = k. \\ hiset is the set of values with highest prime factor > sqrt n. my(lim=vector(n, k, sqrtint(k)), h = vector(n\nextprime(sqrtint(n)+1), k, primepi(n\k)-primepi(n\(k+1))), hiset = Set(select(i->my(f=factor(i)[,1]); f[#f]^2 > n, [2..n])) ); my(accept(b,k) = if(setsearch(hiset,k), 0, for(i=2, lim[k], if(bittest(b,i) && !(k%i) && bittest(b,k/i), return(0))); 1)); my(indrec(k, j, b)= j++; if(j > k, 1, my(t = self()(k, j, b)); for(i=2, j, if(bittest(b,i) && !(j%i) && bittest(b, j/i+k), return(t))); t += self()(k, j, b + (1<<(k+j))); ); ); my(recurse(k, b)= my(m=1); if(k <= #h && h[k], m*=indrec(k, 0, b)^h[k]); for(j=max(2*k,n\2+1), min(2*k+1,n), if(accept(b,j), m*=2)); k++; if(k > n\2, m, my(t = m*self()(k, b)); if(accept(b,k), t += m*self()(k, b + (1<