Summation over integers satisfying some conditions

Question

$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$?

Answer 1

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.

Answer 2

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]];

For testing I picked just one subset to start with:

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}]
  ]

Answer 3

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?

Answer 4

[Not a complete answer, sorry.]

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.