Let $N(n)$ be a set of integers, which can be presented using first $n$ Egyptian fractions:

$$ N(n):=\{m\in\mathbb{Z}:\ \ m=\sum_{i=1}^n\frac{\epsilon_i}{i},\ \epsilon_i=0\ \text{or}\ 1\} $$

I want to write a code in Mathematica that gives $N(n)$, but I think A[e,n] does not define what I need above: for example



P.S. example

\begin{align} 0=&0\\ 1=&1\\ 2=&1+\frac12+\frac13+\frac04+\frac05+\frac16 \end{align} therefore $N(6)=\{0,1,2\}$

  • $\begingroup$ A[e,n] does not depend on e. So you can use A[n_]:=... Otherwise it looks good. $\endgroup$
    – Sumit
    Commented Nov 26, 2017 at 8:50
  • 1
    $\begingroup$ A solution that scales (very, very) poorly is Union[Select[Total /@ Subsets[1/Range[6]], IntegerQ]] (for n=6, for example). $\endgroup$
    – evanb
    Commented Nov 26, 2017 at 16:26
  • 2
    $\begingroup$ I changed the title to make the question to be more attractive. I hope, you don't mind. $\endgroup$
    – ybeltukov
    Commented Nov 26, 2017 at 21:44
  • 1
    $\begingroup$ Straight-forward method using Solve (and not suffering from exponential memory usage) would be something like this: n0 /. (With[{n = 50}, Solve[{Sum[a[i]/i, {i, n}] == n0, Array[a, n] \[Element] Cuboid[ConstantArray[0, n]]}, {n0}, Integers]] /. ConditionalExpression[v_, ___] :> v) (* {0, 1, 2, 3} *) ... this doesn't solve the exponential time complexity, though. $\endgroup$
    – kirma
    Commented Nov 27, 2017 at 6:39
  • 2
    $\begingroup$ This problem is essentially the same as computing the OEIS sequence A101877. With a little bit of cleverness one can search the values up to at least 300 (which would correspond with $N(300)=\{0,1,2,3,4,5\}$), just by using straight-forward Solve and constraints on the problem size based on solutions already found. Going up to the sum corresponding to 6 (at 469) might be infeasible using this method (probably takes at least several days of computing time, possibly a lot more). $\endgroup$
    – kirma
    Commented Nov 29, 2017 at 7:19

3 Answers 3


This solution acquires answers for couple hundred first values of $N$ - probably not much further in a reasonable amount of time. Much more efficient - and more explicit - code for finding larger values can be found in a link from OEIS A101877.

This pair of functions computes both $N$ and its companion function which returns a list of new integer solutions (value and contributing integers in sum of fractions) for any $N$. The construct is recursive, always searching only for solutions for sums which haven't been seen for lower values.

ClearAll[n, sn];
n[0] = {0};
n[x_] := n[x] = Union[n[x - 1], sn[x][[All, 1]]];
sn[0] = {{0, {}}};
sn[x_] := sn[x] = With[{r = Array[c, x]},
    {a /. #, Flatten@Position[r /. #, 1]} & /@
     Solve[Sum[c[i]/i, {i, x}] == a && c[x] == 1 && 
       Unequal[a, Sequence @@ n[x - 1]] && 
       And @@ (0 <= # <= 1 & /@ r), Append[r, a], Integers]];

Even to reach 5 it can take quite a while, over five hours of CPU time - but this is dramatically more efficient than exhaustive search, still:


{19020.6, {0, 1, 2, 3, 4, 5}}

At the same time this call saves solutions for minimum-$N$ integer solutions, for instance:

Grid[DeleteCases[Table[sn[i], {i, 0, 65}], {}][[All, 1]], 
 Alignment -> Left, Frame -> All]

enter image description here


Let indicator variables e_i encode a number in binary. That is let $b=\sum_{i} 2^i e_i$. Then we can loop over all values of $b$ and decode it into a vector of $e$s:

e[n_][b_] := ArrayPad[#, {n - Length[#], 0}] &@IntegerDigits[b, 2]

(We need to pass n to indicate how long the vector should be).

Now we can loop over all $2^n$ choices for $b$ and compute the sum by taking the dot product with a vector of $1/i$:

allSums[n_] := Map[(1/Range[n]).e[n][#] &, Range[2^n]]

and now filter / sort:

set[n_] := Union@Select[IntegerQ]@allSums[n]

I'm afraid that this will still take $\sim2^{1000}$ steps for $n=1000$, but at least it doesn't have an exponentially large memory requirement like the procedure mentioned in my comment :-) If you're really interested in $n=1000$ you'll need an algorithm that's better than the "obvious" one.

f[x_] := Sum[x[[i]]/i,{i,1,Length[x]}]
Select[IntegerQ][f /@ Tuples[{0,1},20]] 

gives still {0,1,2}, with the raspberry taking 5 min.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.