# Summation over integers satisfying some conditions

$n$ is a fixed positive integer and $p$ is the largest prime $\le n$:

p = Prime[PrimePi[n]]


For each subset $L$ of positive composite integers less than or equal to $n$--i.e., for each element $L$ of:

Subsets[Select[Range[n], CompositeQ]]


I would like to sum $\sum \frac{1}{i}$ over all positive integers $i$ such that

(i) $i$ has no prime factor larger than $p$:

FactorInteger[i][[-1, 1]] <= p


and

(ii)

j/GCD(i,j)


is composite for each $j\in L$.

How can I incorporate (i) and (ii) into a sum given $L$?

Update: The fact that are are infinitely many $i$'s corresponding to each $L$ may cause difficulty. Perhaps I can find all $i$ satisfying (i) and (ii) where $i$ is bounded above by some $N$?

• What is the magnitude of n you want to test? The number of subsets grows quite rapidly... Moreover, are you sure that there will be any i satisfying your conditions? – Henrik Schumacher Dec 29 '17 at 14:11
• Btw.: := is not the way to set a variable (e.g., p) in Mathematica. Better use p = Prime[PrimePi[n]] (have a look at the documentation about SetDelayed and what it does). – Henrik Schumacher Dec 29 '17 at 14:12
• I would like to test many values of n, but perhaps I'd like to start with 28 (since I've resolved the problem for smaller n). Yes, there will always be i's; for example, we can always take i=1 since then j/GCD(1,j) = j is composite for each j in L. – The Substitute Dec 29 '17 at 23:10
• And can there be another, nontrivial i? Otherwise, this going to be a very expensive way to compute many ones. Maybe you should give a detailed example for some n in the question. This is a task in which one has to invest quite a lot of effort and people won't attack it without knowing that this will lead somewhere... – Henrik Schumacher Dec 30 '17 at 12:25
• Since composite numbers only have prime factors <= n/2, any i consisting only of primes greater than or equal to n/2 will always work. E.g., if n=12, then for any choice of L, any j in L has no prime factor larger than 5. So any i of the form i=7^a*11^b will work since the GCD will always be 1. This furnishes infinitely many values of i. For more examples of i, it will depend on the set L. – The Substitute Dec 30 '17 at 14:06

Updated to fix a bug pointed out by DanielLichtblau

The restriction $\frac{j}{\gcd(j, i)}$ is composite means that $j$ must contain at least two prime factors not in $i$. The following function encodes this restriction:

restriction[n_] := With[{fi=FactorInteger[n]},
Total[Ramp[#2-f[#1]]& @@@ fi] >= 2
]


For example:

restriction[12]


Ramp[2 - f[2]] + Ramp[1 - f[3]] >= 2

The above says that the prime factorization of $i$ can have no factor of 2 (and any power of 3), or it can have 1 factor of 2 and no factor of 3. Applying this restriction for each element of $L$, and also specifying that the exponents must be nonnegative yields:

allowed[set_List] := With[{restrict = restriction /@ set},
With[{ff = Cases[restrict, _f, Infinity] //Union},
Reduce[And @@ restrict && And @@ Thread[ff >= 0], ff, Integers]
]
]


For example:

allowed[{4, 12, 21}]


f[2] == 0 && f[3] == 0 && f[7] == 0

This says that the exponents of 2, 3 and 7 must be 0, but the remaining exponents are not restricted. So, the set of elements $i$ consists of:

5^a 11^b 13^c 17^d 19^e 23^f


for all nonnegative integers $a, b, c, d, e,$ and $f$. Finally, using:

$$\sum _{a_1=0}^{\infty } \sum _{a_2=0}^{\infty } \cdots \sum _{a_k=0}^{\infty }\frac{1}{p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}}=\prod _{i=1}^k \frac{1}{1-\frac{1}{p_i}}$$

we obtain:

unrestrictedSum[p_List] := Times @@ (1/(1 - 1/p))

unrestrictedSum[{5, 11, 13, 17, 19, 23}]


96577/55296

Now, the allowed set can sometimes contain a parameter. For example:

allowed[{12}]


(f[2] == 1 && f[3] == 0) || (f[2] == 0 && f[3] == 1) || (f[2] == 0 && f[3] == 0) || (C[1] ∈ Integers && C[1] >= 2 && f[2] == 0 && f[3] == C[1])

Here is a function that sums up the reciprocal contributions of each Or element:

reciprocalSums[a_Or]:=With[{logical = LogicalExpand[a]},
Total[reciprocalSum /@ List @@ logical]
]
reciprocalSums[a_]:=reciprocalSum[a]

reciprocalSum[e_]:=Module[{fs = Cases[e, _f, Infinity] //Union, res, cs, n},
res = m /. First @ Solve[Reduce[m == Times @@ Replace[fs, f[n_]:>n^-f[n], {1}] && e, m], m, fs];
cs = Sequence @@ Cases[e, C[i_] >= l_ :> {n[i], l, Infinity}, Infinity];
res /. ConditionalExpression[m_, conds_] :> Sum @@ {m /. C[i_]:>n[i], cs}
]


For the above example we have:

reciprocalSums[allowed[{12}]]


2

Or, breaking it down to pieces:

reciprocalSum[f[2] == 1 && f[3] == 0]
reciprocalSum[f[2] == 0 && f[3] == 1]
reciprocalSum[f[2] == 0 && f[3] == 0]
reciprocalSum[C[1] ∈ Integers && C[1] >= 2 && f[2] == 0 && f[3] == C[1]]


1/2

1/3

1

1/6

Finally, packaging up the above as a function:

restrictedSum[set_, p_] := Module[{a = allowed[set], r, u},
r = Cases[a, f[q_]->q, Infinity]//Union;
u = Complement[Prime[Range[PrimePi[p]]], r];

reciprocalSums[a] unrestrictedSum[u]
]


Some examples:

restrictedSum[{8}, 23]
restrictedSum[{12}, 23]
restrictedSum[{25}, 23]
restrictedSum[{12 25}, 23]
restrictedSum[{12, 25}, 23]
restrictedSum[{12, 15}, 23]


676039/147456

676039/165888

676039/138240

49350847/8294400

676039/207360

676039/276480

This function should work for any set $L$. @DanielLichtblau's answer only works when the set $L$ is coprime, and in that case, my answer and his answer agree.

• Take the subset {30}. The restrictions above are {{f[2]->0,f[3]->0,f[5]->0},{f[2]->0,f[3]->0,f[5]->1},{f[2]->0,f[3]->1,f[5]->0},{f[2]->1,f[3]->0,f[5]->0}. This does not allow a factor such as 1/2^n for n>1. This is why it is giving results that are too low when given values in L that have three or more prime factors (counting multiplicity). – Daniel Lichtblau Jan 3 '18 at 17:18
• @DanielLichtblau Thanks, I see where I went wrong now. I will patch it up shortly. – Carl Woll Jan 3 '18 at 17:39

Straightforward implementation. First, some definitions:

n = 4;
p = Prime[PrimePi[n]];
imax = 10;
L = Subsets[Select[Range[n], CompositeQ]];

set = Select[Range[imax], (FactorInteger[#][[-1, 1]] <= p) && (And @@
Table[CompositeQ[j/GCD[#, j]], {j, Flatten@L}]) &]

(* {1, 3, 9} *)


Now, the sum is trivial to implement:

Sum[1/i, {i, set}]
(* 13/9 *)


Is this the expected output?

• Impressive, but very clear approach. Wow. – Riccardo Cazzin Jan 2 '18 at 23:35
• From what I understand, the L defined by the OP is actually the elements of what you define L to be. The OP wants a sum for each of these elements. From what I can tell, this method doesn't consider these elements separately due to the Flatten function. – BioPhysicist Jan 3 '18 at 14:58
• @AaronStevens You do have a point. The OP is slightly ambiguous IMHO. I left a comment asking for clarification. – AccidentalFourierTransform Jan 3 '18 at 21:14

I have something, but I am unsure if it is what you are looking for or if it is even fast enough for larger values of n. Using n = 28:

p = Prime[PrimePi[28]]
23

sets = Subsets[Select[Range[28], CompositeQ]];


set = sets[[1000]]
{4, 6, 8, 25}


The first condition for i is already given by you. For the second condition:

gcdTest[i_, subset_] := And @@ (CompositeQ[#/GCD[i, #]] & /@ subset)


I am basically just testing each element in the subset L with a given i, then using the logical And to pick out only the i values where this is true for all elements in L.

The 2 tests can then be used to pick out the desired i values for some given cutoff value (to avoid infinite possible i values):

cutoff = 1000;

iList = Select[
Range[cutoff], (FactorInteger[#][[-1, 1]] <= p &&
gcdTest[#, set]) &]
{1, 7, 11, 13, 17, 19, 23, 49, 77, 91, 119, 121, 133, 143, 161, 169,
187, 209, 221, 247, 253, 289, 299, 323, 343, 361, 391, 437, 529, 539,
637, 833, 847, 931}


From here we can determine the desired total:

N@Total[1/iList]
1.60518


This process could then be applied to all subsets.

I am far from an expert in mathematica, so there is probably a way better way to do this, but I thought I would take a shot.

UPDATE: Here is a more compact version that will give a list of the specified sums for all of the subsets L rather than just one L as I have demonstrated above. Speed wise, this will take about 50 minutes to find the sums for the 262143 possible sets L when n = 28 (it takes about 3 seconds to get through about 260 of the sets).

gcdTest[i_, subset_] := And @@ (CompositeQ[#/GCD[i, #]] & /@ subset);

sums[n_, iMax_] :=
Module[

{p = Prime[PrimePi[n]],
sets = Subsets[Select[Range[28], CompositeQ]][[2 ;;]]
},

Table[Total[1/#] &
[Select[
Range[iMax], (FactorInteger[#][[-1, 1]] <= p &&
gcdTest[#, sets[[a]]]) &]]
, {a, Length@sets}]
]

• We did almost the same thing. I'm quite confident this is the way to go, even if I'm not an expert at all. I'm totally upvoting you. – Riccardo Cazzin Jan 2 '18 at 18:21

There is a fairly simple-to-code way to handle the (very) special case of subsets of one element, and a straightforward extension for when the set is all pairwise coprime.

First the singleton case. Suppose plist is the set of primes up to the cutoff. The unrestricted sum of all produc-of-reciprocal-prime-powers is Product[1/(1 - 1/pj), {pj, plist}]. Given another integer n, from this unrestricted sum certain summands need to be removed in a way that involves the prime factorization of n (this is to support requirement #2). It amounts to carving out infinite rectangular regions in the integer lattice of integer powers e.g. if a term in the unrestricted sum is 1/p1^e1*1/p2^e2*...*1/pk^ek then the lattice point is {e1,e2,...,ek}. I won't try to explain in detail why this works, but it involves removing first a rectangle that is infinite in all dimensions involving prime factors ofn, then removing pieces that are each infinite in one less dimension (each corresponding to one particular prime factor).

sum[n_, plist_] := Module[{fax = FactorInteger[n], flist, prod},
prod = Product[1/(1 - 1/pj), {pj, plist}];
flist = fax[[All, 1]];
(1 - 1/n)*prod - Sum[p/n*prod*(1 - 1/p), {p, flist}]]


Example: We have the subset {12} and we are working with maximum prime of 23.

s12 = sum[12, Prime[Range[9]]]

(* Out[2126]= 676039/165888 *)

N[s12]

(* Out[2127]= 4.07527367863 *)


Here are a few more examples. They might be compared to some of the more general methods already presented.

Map[sum[#, Prime[Range[9]]] &, {8, 12, 14, 25, 12*25, 12*35}] // N

(* Out[2133]= {4.58468288845, 4.07527367863, 2.6198187934, \
4.89032841435, 5.94989957079, 5.90914683401} *)


When the set is pairwise coprime one can basically iterate the method of carving away some semiinfinite rectangles.

sumCoprimes[nlist_, plist_] /; Apply[CoprimeQ, nlist] :=
Module[{res, nj, fax, flist},
res = Product[1/(1 - 1/pj), {pj, plist}];
Do[nj = nlist[[j]];
fax = FactorInteger[nj];
flist = fax[[All, 1]];
res = ((1 - 1/nj) - Sum[p/nj*(1 - 1/p), {p, flist}])*res;
, {j, Length[nlist]}];
res
]


Examples:

sumCoprimes[{12, 25}, Prime[Range[9]]] // N

(* Out[2150]= 3.2602189429 *)

sumCoprimes[{12, 25, 7^2*11}, Prime[Range[9]]] // N

(* Out[2151]= 3.15739200036 *)
`

If ever I figure out how to lift the coprimality assumption and extend to the general case, e.g. by inclusion-exclusion, I'll post an update. I thought even without that it might be useful to have an exact method for this special case.