How to represent integers using Egyptian fractions?



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\}$


Posted 2017-11-26T07:56:18.283

Reputation: 778

A[e,n] does not depend on e. So you can use A[n_]:=... Otherwise it looks good. – Sumit – 2017-11-26T08:50:34.533

Your A function takes e == 0 and runs through all the i-values, and then takes e == 1 and runs through all the i-values, giving 0/1 + 0/2 + ... + 0/n + 1/1 + 1/2 + ... + 1/n. Try Table[{e, i}, {e, 0, 1}, {i, 1, 6}] to see what the indices are doing. – aardvark2012 – 2017-11-26T10:49:57.820

Doesn't N(n) always contain 0? For example, for N(6), 0 = 0/1 + 0/2 + 0/3 + 0/4 + 0/5 + 0/6. – evanb – 2017-11-26T16:23:24.463

1A solution that scales (very, very) poorly is Union[Select[Total /@ Subsets[1/Range[6]], IntegerQ]] (for n=6, for example). – evanb – 2017-11-26T16:26:52.137

You can brute force it for $n \lesssim 25$. Will that work, or are you really interested in $n=1000$? – Carl Woll – 2017-11-26T16:31:42.050

2I changed the title to make the question to be more attractive. I hope, you don't mind. – ybeltukov – 2017-11-26T21:44:39.093

@evanb, thanks, I corrected. Also I like to have representation for each of them in output as you see in P.S. – asad – 2017-11-27T05:39:41.927

1Straight-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. – kirma – 2017-11-27T06:39:33.060


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).

– kirma – 2017-11-29T07:19:06.137



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


Posted 2017-11-26T07:56:18.283

Reputation: 13 550


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.


Posted 2017-11-26T07:56:18.283

Reputation: 4 517


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.


Posted 2017-11-26T07:56:18.283

Reputation: 31