Q:

Is there a relatively quick way to solve my integer-valued equation?

This is a very simple problem. It can be solved by nested Forloops that exhaustively test every possible case, throwing out all non-solutions.

Find all the distinct positive integer solutions $0<a_i \leqslant 100$, such that:
$$\sum\limits_{i=1}^{10}\dfrac{1}{a_i}=1$$

There are 69014 solutions in total. But it will take a long time to find them all by nested Forloops. Is there an efficient method for this problem, i.e, one that will find all the solutions in less than 20 minutes?

A:

Abstract
This is not exactly the requested solution, but it is parametrized, and I think it will give some insight and can possibly be improved. 
With this code (modified 03.09.14 11:30) it is possible to determine all integer solutions to $1=\sum _{i=1}^n \frac{1}{a(i)}$ without the $a(i)$ being bounded by $amax$, by letting $amax=\infty$.
In this sense the code is complementary to the elegant solution of kirma which requires $amax$ to be given explicitly beforehand.
For $amax = 100$ the code runs in "reasonable" time (some minutes) up to $n = 8$. For $n = 9$ time and memory requirements explode.
For $amax=\infty$ the code runs in "reasonable" time (some seconds) up to $n = 6$. For $n = 7$ time and memory requirements explode.
1. Statement of the problem
calculate all solutions of the equation 
$\sum _{i=1}^n \frac{1}{a(i)}$
under the conditions
$a(i)\in \mathbb{Z}$ and $0<a(1)<a(2)<...<a(n)\leq \text{amax}$
for $n=10$ and $amax=100$ in "reasonable" time (~ 20 Minutes) with Mathematica.
Although I didn't meet the goal with these parameters I have developed a code which works "in principle" and in reasonable time for "appropriate" parameters.
2. Reduction of possible values of $a(i)$
The main idea is to determine the range of the possible numbers $a(i)$.
Consider the explanatory example $n=4$ and $\text{amax}=\infty$ (i.e. no restriction)
Let us determine the interval of possible values for the first number a.
Because $b>a,c>a,d>a$ we have the inequality
$1=\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d}<\frac{4}{a}$
or
$a < 4$
But
$1-\frac{1}{a}=\frac{1}{b}+\frac{1}{c}+\frac{1}{d}$
which is evidently $> 0$. Hence
$a > 1$
so that the range for a becomes
$a\in [2,3]$
In this first step we have thus obtained the following set of two partial solutions
v[1] = {{2}, {3}}

Now we take the first partial solution {2} and calculate the possible range of the next element $b$.
We have with $a = 2$:
$1-\frac{1}{2}=\frac{1}{2}=\frac{1}{b}+\frac{1}{c}+\frac{1}{d}<\frac{3}{b}\to b<6$
and
$\frac{1}{2}-\frac{1}{b}>0\to b>2$
Hence 
$b \in [3,4,5]$
and the set of partial solutions starting with 2 becomes
v22 = {{2,3}, {2,4}, {2,5}}

For $a = 3$ we find
$1-\frac{1}{3}=\frac{2}{3}<\frac{3}{b}\to b<\frac{9}{2}=4.5$
and
$\frac{2}{3}-\frac{1}{b}>0 \to b>\frac{3}{2} = 1.5$
or
$b \in [2,3,4]$
But $b$ must also be $> a (=3)$ hence only $b = 4$ remains,
and the set of partial solutions starting with 3 is
v23 = {{3,4}}

The total set of partial solutions of length 2 is therefore
v[2] = {{2,3}, {2,4}, {2,5}, {3,4}}

Now the same procedure for each of the four partial solutions
$\{2,3\}\to 1-\left(\frac{1}{2}+\frac{1}{3}\right)=\frac{1}{2}-\frac{1}{3}=\frac{1}{6}=\frac{1}{c}+\frac{1}{d}<\frac{2}{c}\to c<12$
and
$\frac{1}{6}-\frac{1}{c}>0\to c>6$
$\to c \in [7,8,9,10,11]$
v31 = {{2,3,7},{2,3,8},{2,3,9},{2,3,10},{2,3,11}}

$\{2,4\}\to 1-\left(\frac{1}{2}+\frac{1}{4}\right)=\frac{1}{4}=\frac{1}{c}+\frac{1}{d}<\frac{2}{c}\to c<8$
and
$c>4 \to c \in [5,6,7]$
v32 = {{2,4,5},{2,4,6},{2,4,7}}

$\{2,5\}\to 1-\left(\frac{1}{2}+\frac{1}{5}\right)=\frac{3}{10}=\frac{1}{c}+\frac{1}{d}<\frac{2}{c}\to c<\frac{20}{3}=6.666$
Because of $c > b (=5)$ we find $c = 6$, and
v33 = {{2,5,6}}

$\{3,4\}\to 1-\left(\frac{1}{3}+\frac{1}{4}\right)=\frac{5}{12}=\frac{1}{c}+\frac{1}{d}<\frac{2}{c}\to c<\frac{24}{5}=4.8\to c\leq 4$
but this contradicts the condition $c > b (=4)$. Hence there is no extension of the partial solution in this case.
Collecting things gives for the partial solutions of length 3
v[3] = {{2,3,7},{2,3,8},{2,3,9},{2,3,10},{2,3,11},{2,4,5},{2,4,6},{2,4,7},{2,5,6}}

Now the last step. Here we simply calculate d for each partial solution from
$\frac{1}{d}=1-\left(\frac{1}{a}+\frac{1}{b}+\frac{1}{c}\right)$
Using the function
is[x_] := Plus @@ (1/x)

we find
(1 - is[#]) & /@ {{2, 3, 7}, {2, 3, 8}, {2, 3, 9}, {2, 3, 10}, {2, 3, 11}, {2, 4, 5}, {2, 4, 6}, {2, 4, 7}, {2, 5, 6}}

$\left\{\frac{1}{42},\frac{1}{24},\frac{1}{18},\frac{1}{15},\frac{5}{66},\frac{1}{20},\frac{1}{12},\frac{3}{28},\frac{2}{15}\right\}$
We observe that there are cases in which the numerator is different from 1. These cases must be ruled out, of couse.
This leaves us the 6 solutions for n = 4.
v[4] = {{2, 3, 7, 42}, {2, 3, 8, 24}, {2, 3, 9, 18}, {2, 3, 10, 15}, {2, 4, 5,
    20}, {2, 4, 6, 12}}

We could check the result using Reduce:
Reduce[1 == 1/a + 1/b + 1/c + 1/d && 0 < a < b < c < d, Integers]

(* (a == 2 && b == 3 && c == 7 && d == 42) || (a == 2 && b == 3 && c == 8 && 
   d == 24) || (a == 2 && b == 3 && c == 9 && d == 18) || (a == 2 && b == 3 &&
    c == 10 && d == 15) || (a == 2 && b == 4 && c == 5 && 
   d == 20) || (a == 2 && b == 4 && c == 6 && d == 12) *)
We can improve the estimation and increase the lower bound by taking into account that condition $a <= amax$, Taking the example $n=4$ again, we have
$1-\frac{1}{a}=\frac{1}{b}+\frac{1}{c}+\frac{1}{d}$
But because $d <= amax$ the rhs. is $> \frac{3}{amax}$, resulting in
$1-\frac{3}{\text{amax}}>\frac{1}{a}\to a>\frac{1}{1-\frac{3}{\text{amax}}}$
This again can be sharpened a bit using 
$\frac{1}{\text{amax}}+\frac{1}{\text{amax}-1}+\frac{1}{\text{amax}-2}$
or, generally, the function
sm[n_, k_] := Sum[1/(amax - j), {j, 0, n - k - 1}]

This long example containes all ingredients of the code presented in the next section.
3. The code
The code is not particilarly elegant, and can be improved doublessly.
The solution vectors v[r] are global.
I have also left the loop variable k global in order to display its values during the calculation using Dynamic.
First we dfine the two auxiliary functions introduced earlier:
is[x_] := Plus @@ (1/x)

sm[n_, k_] := Sum[1/(amax - j), {j, 0, n - k - 1}]

calcV6 is the main routine. It calculates the partial solution v[r] assuming the previous one, v[r-1], given.
The normal usage of calcV6 is therefore to execute the routine for $r = 2$ first, then for $r = 3$, and so on up to $r = n$.
The routine has some additional "security" measures not mentioned in the example but easily recognizable by the attentive reader.
calcV6[n_, r_] := 
 Module[{rr, smr, ss, ss1, z, z1, ymin, ymax, u, us, us1},
  If[amax < \[Infinity],
   v[1] = Table[{i}, {i, 2, Min[n - 1, amax]}],
   v[1] = Table[{i}, {i, 2, n - 1}]];
  If[r == n, Goto[mrk]];
  v0[r] = {};
  If[amax < \[Infinity], smr = sm[n, r], smr = 0]; (* calcV6 *)
  For[k = 1, k <= Length[v[r - 1]], k++,
   x = v[r - 1][[k]];
   ss1 = 1 - is[x];
   If[ss1 > 0, z1 = 1/ss1, Goto[cont]];
   ss = 1 - is[x] - smr;
   z = If[ss > 0, 1/ss, z1];
   ymin = Max[1 + x[[-1]], Floor[z]];
   If[amax < \[Infinity],
    ymax = Min[Ceiling[(n - r + 1) z1] - 1, amax],
    ymax = Ceiling[(n - r + 1) z1] - 1];
   rr = Range[ymin, ymax];
   If[ymin <= ymax, AppendTo[v0[r], Join[x, {#}] & /@ rr]];
   Label[cont];
   ];
  v[r] = Flatten[v0[r], 1];
  Return[];
  Label[mrk];
  (* Last step: calculate v[n] *)
  u = v[r - 1];
  us = Select[u, 1 == Numerator[1 - is @@ {#}] &];
  us1 = Join[#, {Denominator[1 - is[#]]}] & /@ us ;
  If[amax < \[Infinity], v[n] = Select[us1, #[[-1]] <= amax &], 
   v[n] = us1]; (* calcV6 *)
  ]

4. Applications
Now we are going to apply the code in some example cases
Example 1: $n = 4, amax = \infty$ (i.e. no restriction)
amax = \[Infinity];
Clear[v];
calcV6[4, 2]
v[2]
calcV6[4, 3]
v[3]
calcV6[4, 4]
v[4]

(* {{2, 3}, {2, 4}, {2, 5}, {3, 4}} *)

(* {{2, 3, 6}, {2, 3, 7}, {2, 3, 8}, {2, 3, 9}, {2, 3, 10}, {2, 3, 11}, {2, 4, 
  5}, {2, 4, 6}, {2, 4, 7}, {2, 5, 6}} *)

(* {{2, 3, 7, 42}, {2, 3, 8, 24}, {2, 3, 9, 18}, {2, 3, 10, 15}, {2, 4, 5, 
  20}, {2, 4, 6, 12}} *)

Example 2 : $n = 4, amax = 13$ (i.e. restriction operating)
amax = 13;
Clear[v];
calcV6[4, 2]
v[2]
calcV6[4, 3]
v[3]
calcV6[4, 4]
v[4]

(* {{2, 3}, {2, 4}, {2, 5}, {3, 4}} *)

(* {{2, 3, 11}, {2, 4, 5}, {2, 4, 6}, {2, 4, 7}, {2, 5, 6}} *)

(* {{2, 4, 6, 12}} *)

Example 3 : $n = 10, amax = 30$ (the required n but much smaller amax)
The output is:
column 1    n
column 2    amax
column 3    table: for each r = 2 .. n : 
{r, Timing, number of partial soultions}, 
number of partial solutions,
{shortest (partial) solution vector, longest (partial) solution vector; according to Euclidian norm}
amax = 30;
nmax = 10;
Clear[v];
t = Table[{n, amax,
    Table[{r, Timing[calcV6[n, r]][[1]], Length[v[r]]}, {r, 2, n}], 
    v[n] // Length, {mi = Min[#.# & /@ v[n]]; 
     Select[v[n], #.# == mi &][[1]], mx = Max[#.# & /@ v[n]]; 
     Select[v[n], #.# == mx &][[1]]}}, {n, 3, nmax}];
t // TableForm

$\begin{array}{lllll}
 3 & 30 & 
\begin{array}{lll}
 2 & 0. & 1 \\
 3 & 0. & 1 \\
\end{array}
 & 1 & 
\begin{array}{lll}
 2 & 3 & 6 \\
 2 & 3 & 6 \\
\end{array}
 \\
 4 & 30 & 
\begin{array}{lll}
 2 & 0. & 4 \\
 3 & 0. & 9 \\
 4 & 0. & 5 \\
\end{array}
 & 5 & 
\begin{array}{llll}
 2 & 4 & 6 & 12 \\
 2 & 3 & 8 & 24 \\
\end{array}
 \\
 5 & 30 & 
\begin{array}{lll}
 2 & 0. & 8 \\
 3 & 0. & 26 \\
 4 & 0. & 75 \\
 5 & 0. & 11 \\
\end{array}
 & 11 & 
\begin{array}{lllll}
 3 & 4 & 5 & 6 & 20 \\
 2 & 3 & 12 & 20 & 30 \\
\end{array}
 \\
 6 & 30 & 
\begin{array}{lll}
 2 & 0. & 14 \\
 3 & 0. & 50 \\
 4 & 0. & 169 \\
 5 & 0.016 & 458 \\
 6 & 0. & 26 \\
\end{array}
 & 26 & 
\begin{array}{llllll}
 3 & 4 & 6 & 10 & 12 & 15 \\
 2 & 3 & 20 & 21 & 28 & 30 \\
\end{array}
 \\
 7 & 30 & 
\begin{array}{lll}
 2 & 0. & 20 \\
 3 & 0. & 73 \\
 4 & 0. & 292 \\
 5 & 0.015 & 959 \\
 6 & 0.047 & 2241 \\
 7 & 0.016 & 36 \\
\end{array}
 & 36 & 
\begin{array}{lllllll}
 3 & 4 & 9 & 10 & 12 & 15 & 18 \\
 2 & 4 & 12 & 20 & 21 & 28 & 30 \\
\end{array}
 \\
 8 & 30 & 
\begin{array}{lll}
 2 & 0. & 29 \\
 3 & 0. & 138 \\
 4 & 0.015 & 693 \\
 5 & 0.016 & 2322 \\
 6 & 0.093 & 4688 \\
 7 & 0.234 & 8341 \\
 8 & 0.094 & 43 \\
\end{array}
 & 43 & 
\begin{array}{llllllll}
 4 & 5 & 6 & 9 & 10 & 15 & 18 & 20 \\
 2 & 6 & 8 & 20 & 21 & 24 & 28 & 30 \\
\end{array}
 \\
 9 & 30 & 
\begin{array}{lll}
 2 & 0. & 36 \\
 3 & 0. & 158 \\
 4 & 0.016 & 642 \\
 5 & 0.015 & 2256 \\
 6 & 0.109 & 6437 \\
 7 & 0.359 & 14356 \\
 8 & 1.077 & 22835 \\
 9 & 0.265 & 39 \\
\end{array}
 & 39 & 
\begin{array}{lllllllll}
 4 & 6 & 8 & 9 & 10 & 12 & 15 & 18 & 24 \\
 2 & 8 & 9 & 18 & 20 & 21 & 24 & 28 & 30 \\
\end{array}
 \\
 10 & 30 & 
\begin{array}{lll}
 2 & 0. & 46 \\
 3 & 0. & 212 \\
 4 & 0.015 & 850 \\
 5 & 0.032 & 2929 \\
 6 & 0.14 & 8261 \\
 7 & 0.515 & 18998 \\
 8 & 1.607 & 34984 \\
 9 & 4.196 & 47403 \\
 10 & 0.577 & 24 \\
\end{array}
 & 24 & 
\begin{array}{llllllllll}
 5 & 6 & 8 & 9 & 10 & 12 & 15 & 18 & 20 & 24 \\
 3 & 6 & 8 & 9 & 18 & 20 & 21 & 24 & 28 & 30 \\
\end{array}
 \\
\end{array}$
Example 4 : $amax = 100$ (the required amax, but n only up to n = 8)
amax = 100;
nmax = 8;
Clear[v];
t = Table[{n, amax,
    Table[{r, Timing[calcV6[n, r]][[1]], Length[v[r]]}, {r, 2, n}], 
    v[n] // Length, {mi = Min[#.# & /@ v[n]]; Select[v[n], #.# == mi &][[1]], 
     mx = Max[#.# & /@ v[n]]; Select[v[n], #.# == mx &][[1]]}}, {n, 3, nmax}];
t // TableForm

$\begin{array}{lllll}
 3 & 100 & 
\begin{array}{lll}
 2 & 0. & 1 \\
 3 & 0. & 1 \\
\end{array}
 & 1 & 
\begin{array}{lll}
 2 & 3 & 6 \\
 2 & 3 & 6 \\
\end{array}
 \\
 4 & 100 & 
\begin{array}{lll}
 2 & 0. & 4 \\
 3 & 0. & 10 \\
 4 & 0. & 6 \\
\end{array}
 & 6 & 
\begin{array}{llll}
 2 & 4 & 6 & 12 \\
 2 & 3 & 7 & 42 \\
\end{array}
 \\
 5 & 100 & 
\begin{array}{lll}
 2 & 0. & 8 \\
 3 & 0. & 30 \\
 4 & 0. & 179 \\
 5 & 0. & 42 \\
\end{array}
 & 42 & 
\begin{array}{lllll}
 3 & 4 & 5 & 6 & 20 \\
 2 & 3 & 7 & 78 & 91 \\
\end{array}
 \\
 6 & 100 & 
\begin{array}{lll}
 2 & 0.015 & 14 \\
 3 & 0. & 59 \\
 4 & 0. & 446 \\
 5 & 0.016 & 4045 \\
 6 & 0.047 & 229 \\
\end{array}
 & 229 & 
\begin{array}{llllll}
 3 & 4 & 6 & 10 & 12 & 15 \\
 2 & 3 & 10 & 22 & 90 & 99 \\
\end{array}
 \\
 7 & 100 & 
\begin{array}{lll}
 2 & 0. & 20 \\
 3 & 0. & 100 \\
 4 & 0. & 893 \\
 5 & 0.047 & 9437 \\
 6 & 0.811 & 83669 \\
 7 & 0.811 & 1131 \\
\end{array}
 & 1131 & 
\begin{array}{lllllll}
 3 & 4 & 9 & 10 & 12 & 15 & 18 \\
 2 & 4 & 5 & 70 & 78 & 84 & 91 \\
\end{array}
 \\
 8 & 100 & 
\begin{array}{lll}
 2 & 0. & 29 \\
 3 & 0. & 152 \\
 4 & 0. & 1330 \\
 5 & 0.078 & 15974 \\
 6 & 2.043 & 185821 \\
 7 & 495.475 & 1523993 \\
 8 & 15.959 & 4934 \\
\end{array}
 & 4934 & 
\begin{array}{llllllll}
 4 & 5 & 6 & 9 & 10 & 15 & 18 & 20 \\
 2 & 3 & 10 & 55 & 70 & 77 & 90 & 99 \\
\end{array}
 \\
\end{array}$
Example 5: all Solutions for $n = 6$
amax = \[Infinity];
nmax = 6;
Clear[v];
t = Table[{n, amax,
    Table[{r, Timing[calcV6[n, r]][[1]], Length[v[r]]}, {r, 2, n}], 
    v[n] // Length, {mi = Min[#.# & /@ v[n]]; 
     Select[v[n], #.# == mi &][[1]], mx = Max[#.# & /@ v[n]]; 
     Select[v[n], #.# == mx &][[1]]}}, {n, 3, nmax}];
t // ColumnForm;
t // TableForm

$\begin{array}{lllll}
 3 & \infty  & 
\begin{array}{lll}
 2 & 0. & 1 \\
 3 & 0. & 1 \\
\end{array}
 & 1 & 
\begin{array}{lll}
 2 & 3 & 6 \\
 2 & 3 & 6 \\
\end{array}
 \\
 4 & \infty  & 
\begin{array}{lll}
 2 & 0. & 4 \\
 3 & 0. & 10 \\
 4 & 0. & 6 \\
\end{array}
 & 6 & 
\begin{array}{llll}
 2 & 4 & 6 & 12 \\
 2 & 3 & 7 & 42 \\
\end{array}
 \\
 5 & \infty  & 
\begin{array}{lll}
 2 & 0. & 8 \\
 3 & 0. & 30 \\
 4 & 0. & 231 \\
 5 & 0. & 72 \\
\end{array}
 & 72 & 
\begin{array}{lllll}
 3 & 4 & 5 & 6 & 20 \\
 2 & 3 & 7 & 43 & 1806 \\
\end{array}
 \\
 6 & \infty  & 
\begin{array}{lll}
 2 & 0. & 14 \\
 3 & 0. & 60 \\
 4 & 0. & 593 \\
 5 & 0.046 & 24690 \\
 6 & 0.234 & 2320 \\
\end{array}
 & 2320 & 
\begin{array}{llllll}
 3 & 4 & 6 & 10 & 12 & 15 \\
 2 & 3 & 7 & 43 & 1807 & 3263442 \\
\end{array}
 \\
\end{array}$
Conclusions
We have provided a code which solves the problem for apropriate parameters $n$ and $amax$.
All solutions can be obtained by letting $amax = \infty$.
The reader can play with the code in order to detect the boundaries of reasonable calculation time and memory usage, respectively. A tradeoff between $n$ and $amax$ must be made.
Regards,
Wolfgang

A:

I shoot for reasonable clarity instead of performance - in my opinion observed <1.5 hour runtime is acceptable, considering clarity - at least to my eyes - of the solution:
With[{b = LCM @@ Range@100}, 
  Select[b/IntegerPartitions[b, {10}, b/Range@100 // Reverse], Unequal@@#&]]

Here all numbers are multiplied by b, least common multiplier for all needed fractions. b, now representing total sum of 1, is partitioned to ten integer partitions, which have to be the size of listed fractions (multiplied by b). Values of original contributing integers are recovered by dividing b by the results, and finally, values where every integer is distinct are selected from this result.
Only non-obvious part in this is Reverse, which would appear to speed up the partitions search. Otherwise it has no effect, except on the order of results.

A:

One might use integer linear programming (one would perhaps be foolish to do so, but one might all the same). I am posting this not as a competitive choice for this problem, but to give some idea of the possibilities. And I do not by any means rule out that there might be significant improvements to be had. As it stands, this is around 20x slower than what was requested.
So here goes. We begin with 99 variables that are intended as taking values of 0 or 1, with exactly ten of them equal to 1. We impose the constraint that their dot product with the reciprocals of {1/2, 1/3, ..., 1/100} equal 1. Multiply by the lcm of these denominators to get an integer equation.
n = 100;
lcm = LCM @@ Range[n];
vals = Rest[lcm/Range[n]];
vars = Array[a, n - 1, 2];
c0 = vars.vals == lcm;
c1 = Map[0 <= # <= 1 &, vars];
c2 = Total[vars] == 10;

Ignoring the constraints that exactly 10 of them contribute exactly one time, form an obvious solution.
soln = 2*UnitVector[n - 1, 1];

Now we reformulate as an integer problem wherein we want not just one solution but also generators for the null space.
lat = Transpose[Prepend[IdentityMatrix[n - 1], vals]];
hnf = HermiteDecomposition[lat];
nulls = hnf[[2, 2 ;; -1, 2 ;; -1]];
newvars = Array[b, Length[nulls]];
SetSystemOptions[
  "LatticeReduceOptions" -> {"LatticeReduceRatioParameter" -> .999}];
rednulls = LatticeReduce[nulls];
c3 = Thread[newvars.rednulls + soln == vars];
constraints = Flatten[{c1, c2, c3}];
allvars = Join[newvars, vars];

Here is a crucial point. The original variables must be an integer linear combination of the solution vector plus integer multiples of the null vectors.
newpolys = newvars.rednulls + soln - vars;

Now we will use some linear programming with machine number coefficients. This next setting gives a modest speed boost but is not something I would admit to having shown.
SetOptions[LinearProgramming, Method -> "CLP"];

We will get extremal values for the new variables, that is, the multipliers of the null vectors. I will repeat this process one time, accounting for prior extrema on the second round. A brief experiment will show that this stabilizes so more repetitions give no improvement.
ctmp = {};
eps = 1/10^5;
Do[mins = 
   Table[Quiet[
     val = NMinimize[{allvars[[j]], Join[constraints, ctmp]}, allvars];
     If[Head[val] === NMinimize || ! FreeQ[val, Indeterminate], 
      bad = True; Break[]];
     val = First[val]; val, NMinimize::"nsol"], {j, 
     Length[allvars]}];
  maxs = Table[
    Quiet[val = 
      NMaximize[{allvars[[j]], Join[constraints, ctmp]}, allvars];
     If[Head[val] === NMaximize || ! FreeQ[val, Indeterminate], 
      bad = True; Break[]];
     val = First[val]; val, NMaximize::"nsol"], {j, 
     Length[allvars]}];
  ctmp = Join[Thread[LessEqual[allvars, Floor[Chop[maxs] + eps]]], 
    Thread[GreaterEqual[allvars, Ceiling[Chop[mins] - eps]]]];
  , {2}];

Now we see if we can remove some of the variables, either the 0-1 ones or the null space multipliers. It turns out that several of the 0-1 variables are forced to be zero.
aconstraints = Cases[ctmp, (GreaterEqual | LessEqual)[_a, _]];
bconstraints = DeleteCases[ctmp, Apply[Alternatives, aconstraints]];
aone = Cases[aconstraints, _a >= 1]
azero = Cases[aconstraints, _a <= 0]

(* Out[288]= {}

Out[289]= {a[23] <= 0, a[29] <= 0, a[31] <= 0, a[37] <= 0, a[41] <= 0,
  a[43] <= 0, a[46] <= 0, a[47] <= 0, a[49] <= 0, a[51] <= 0, 
 a[53] <= 0, a[58] <= 0, a[59] <= 0, a[61] <= 0, a[62] <= 0, 
 a[64] <= 0, a[67] <= 0, a[68] <= 0, a[69] <= 0, a[71] <= 0, 
 a[73] <= 0, a[74] <= 0, a[79] <= 0, a[81] <= 0, a[82] <= 0, 
 a[83] <= 0, a[86] <= 0, a[87] <= 0, a[89] <= 0, a[92] <= 0, 
 a[93] <= 0, a[94] <= 0, a[97] <= 0, a[98] <= 0} *)

Important remark: This means we now have not Binomial[99,10] but "only" Binomial[65,10] possibilities to account for. THis almost brings brute force into the realm of tractability. But wait...
We now remove these variables and begin anew. Some of this, I'm sure, is unnecessary code. But we'll cut to the chase soon enough.
new = Complement[Range[2, n], azero[[All, 1, 1]]]
newn = Length[new];
lcm = LCM @@ new;
vals = lcm/new;
vars = Array[a, newn];
c0 = vars.vals == lcm;
c1 = Map[0 <= # <= 1 &, vars];
c2 = Total[vars] == 10;
soln = 2*UnitVector[newn, 1];

(* Out[429]= {2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, \
18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 34, 35, 36, 38, \
39, 40, 42, 44, 45, 48, 50, 52, 54, 55, 56, 57, 60, 63, 65, 66, 70, \
72, 75, 76, 77, 78, 80, 84, 85, 88, 90, 91, 95, 96, 99, 100} *)

We now rinse and repeat. Again, pardon my confusion, this may be unnecessary.
lat = Transpose[Prepend[IdentityMatrix[newn], vals]];
hnf = HermiteDecomposition[lat];
nulls = hnf[[2, 2 ;; -1, 2 ;; -1]];
newvars = Array[b, Length[nulls]];
SetSystemOptions[
  "LatticeReduceOptions" -> {"LatticeReduceRatioParameter" -> .999}];
rednulls = LatticeReduce[nulls];
c3 = Thread[newvars.rednulls + soln == vars];
constraints = Flatten[{c1, c2, c3}];
allvars = Join[newvars, vars];

ctmp = {};
eps = 1/10^5;
Do[mins = 
   Table[Quiet[
     val = NMinimize[{allvars[[j]], Join[constraints, ctmp]}, allvars];
     If[Head[val] === NMinimize || ! FreeQ[val, Indeterminate], 
      bad = True; Break[]];
     val = First[val]; val, NMinimize::"nsol"], {j, 
     Length[allvars]}];
  maxs = Table[
    Quiet[val = 
      NMaximize[{allvars[[j]], Join[constraints, ctmp]}, allvars];
     If[Head[val] === NMaximize || ! FreeQ[val, Indeterminate], 
      bad = True; Break[]];
     val = First[val]; val, NMaximize::"nsol"], {j, 
     Length[allvars]}];
  ctmp = Join[Thread[LessEqual[allvars, Floor[Chop[maxs] + eps]]], 
    Thread[GreaterEqual[allvars, Ceiling[Chop[mins] - eps]]]];
  , {2}];

We now recheck size ranges for those null vector multipliers.
bconstraints = DeleteCases[ctmp, Apply[Alternatives, aconstraints]];
vranges = 
 Transpose[Reverse[Partition[bconstraints[[All, 2]], Length[newvars]]]]
diffs = -Subtract @@@ vranges

(* Out[454]= {{-3, 5}, {-3, 6}, {-4, 3}, {0, 1}, {-4, 5}, {-5, 5}, {-6, 
  3}, {-3, -1}, {-4, 4}, {0, 3}, {-4, 1}, {-2, 4}, {-4, 1}, {0, 
  4}, {0, 1}, {-1, 2}, {-2, 1}, {0, 1}, {-1, 1}, {-2, 0}, {-1, 
  2}, {-1, 0}, {-3, 0}, {-1, 2}, {-1, 2}, {0, 1}, {0, 1}, {-1, 
  1}, {-1, 1}, {-2, -1}, {0, 2}, {0, 2}, {-1, 1}, {0, 2}, {0, 1}, {-1,
   0}, {0, 1}, {0, 1}, {-1, 1}, {0, 1}, {0, 1}, {-1, 0}, {0, 0}, {0, 
  2}, {0, 1}, {0, 1}, {0, 1}, {-1, 0}, {0, 1}, {0, 0}, {0, 1}, {0, 
  1}, {-1, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {-1, 0}, {-1, 0}, {-1, 
  0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}}

Out[455]= {8, 9, 7, 1, 9, 10, 9, 2, 8, 3, 5, 6, 5, 4, 1, 3, 3, 1, 2, \
2, 3, 1, 3, 3, 3, 1, 1, 2, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, \
0, 2, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0} *)

Some have ranges of length zero, that is, they are exactly determined. So we can remove them from further consideration. Moreover we can reorder so that those with smaller ranges come earlier. This will be important for the branch-and-prune loop that comes later. Also we can remove more variables that occur in simple equations of the form v1==v2 or v1==-v2. Actually we could solve the linear system and remove more. I leave it as an exercise as to whether this will give a speed improvement (I did not try it).
ordr = Ordering[diffs];
newvarsordered = Pick[newvars[[ordr]], Sign[diffs[[ordr]]], 1];
restvars = DeleteCases[newvars, Alternatives @@ newvarsordered];
rule3 = Map[# -> 0 &, restvars];
constraints = (Union[constraints, ctmp] /. rule3) /. 
   True :> Sequence[];
pairs = Cases[constraints, b[_] == a[_]]
rule4 = Apply[Rule, pairs, {1}];
constraints3 = constraints /. rule4 /. True :> Sequence[];
removebs = pairs[[All, 1]];
newvarsordered2 = 
  DeleteCases[newvarsordered, Alternatives @@ removebs];
pairs2 = Cases[constraints, -b[_] == a[_]]

(* Out[461]= {b[15] == a[60], b[18] == a[8], b[26] == a[36], 
 b[27] == a[6], b[45] == a[28], b[45] == a[63], b[4] == a[65], 
 b[46] == a[52], b[4] == a[23], b[46] == a[40], b[47] == a[25], 
 b[47] == a[42]}

Out[466]= {-b[22] == a[13], -b[58] == a[18], -b[59] == 
  a[33], -b[59] == a[62], -b[58] == a[45], -b[60] == a[16], -b[60] == 
  a[30], -b[60] == a[58]} *)

We get rid of those null space multiplier vectors that can be removed by these obvious equalities.
rule5 = Map[-#[[1]] -> -#[[2]] &, pairs2];
removebs2 = -pairs2[[All, 1]];
constraints4 = constraints3 /. rule5 /. True :> Sequence[];
newvarsordered3 = 
 DeleteCases[newvarsordered2, Alternatives @@ removebs2]

(* Out[331]= {b[30], b[35], b[36], b[37], b[38], b[40], b[41], b[42], 
 b[48], b[49], b[51], b[52], b[53], b[8], b[19], b[20], b[28], b[29], 
 b[31], b[32], b[33], b[34], b[39], b[44], b[10], b[16], b[17], b[21],
  b[23], b[24], b[25], b[14], b[11], b[13], b[12], b[3], b[1], b[9], 
 b[2], b[5], b[7], b[6]} *)

Finally we seek out and remove those 0-1 variables that are equated to one another.
allvars = Join[newvarsordered3, vars];
constraints = constraints4;
eq = Cases[constraints4, Equal[a_, b_]];
aeq = Cases[constraints4, _a == _]

(* Out[335]= {a[28] == a[63], a[65] == a[23], a[52] == a[40], 
 a[25] == a[42], a[33] == a[62], a[18] == a[45], a[16] == a[30], 
 a[16] == a[58]} *)

rule6 = Map[#[[1]] -> #[[2]] &, aeq];
constraints5 = constraints4 /. rule6 /. True :> Sequence[];
vars2 = Complement[vars, aeq[[All, 1]]];
aeq = Cases[constraints5, _a == _];
rule7 = Map[#[[1]] -> #[[2]] &, aeq];
constraints6 = constraints5 /. rule7 /. True :> Sequence[];
vars3 = Complement[vars2, aeq[[All, 1]]];
allvars = Join[newvarsordered3, vars3];
constraints = Complement[Union[constraints6], c1];

On the off chance that I did not omit any critical lines, here is what we have.
In[513]:= allvars

(* Out[513]= {b[30], b[35], b[36], b[37], b[38], b[40], b[41], b[42], 
 b[48], b[49], b[51], b[52], b[53], b[8], b[19], b[20], b[28], b[29], 
 b[31], b[32], b[33], b[34], b[39], b[44], b[10], b[16], b[17], b[21],
  b[23], b[24], b[25], b[14], b[11], b[13], b[12], b[3], b[1], b[9], 
 b[2], b[5], b[7], b[6], a[1], a[2], a[3], a[4], a[5], a[6], a[7], 
 a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15], a[17], a[19], 
 a[20], a[21], a[22], a[24], a[26], a[27], a[29], a[30], a[31], a[32],
  a[34], a[35], a[36], a[37], a[38], a[39], a[41], a[42], a[43], 
 a[44], a[45], a[46], a[47], a[48], a[49], a[50], a[51], a[52], a[53],
  a[54], a[55], a[56], a[57], a[59], a[60], a[61], a[62], a[63], 
 a[64], a[65]} *)

In[514]:= constraints

(* Out[514]= {a[45] + a[62] == a[53], 
 a[1] + a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[8] + a[9] + 
   a[10] + a[11] + a[12] + a[13] + a[14] + a[15] + a[17] + a[19] + 
   a[20] + a[21] + a[22] + a[24] + a[26] + a[27] + a[29] + 3 a[30] + 
   a[31] + a[32] + a[34] + a[35] + a[36] + a[37] + a[38] + a[39] + 
   a[41] + 2 a[42] + a[43] + a[44] + 2 a[45] + a[46] + a[47] + a[48] +
    a[49] + a[50] + a[51] + 2 a[52] + a[53] + a[54] + a[55] + a[56] + 
   a[57] + a[59] + a[60] + a[61] + 2 a[62] + 2 a[63] + a[64] + 
   2 a[65] == 10, b[2] + b[5] + b[8] == a[4], -b[8] - b[10] == a[5], 
 b[7] + b[9] + b[10] == a[14], -b[12] - b[13] == 
  a[7], -b[11] + b[13] == a[35], -b[17] - b[19] == a[17], 
 b[16] + b[19] == a[32], -a[6] - b[20] == a[31], 
 b[25] - b[28] == a[44], b[21] + b[28] == a[47], 
 b[24] - b[29] == a[26], 
 b[1] - b[2] - b[30] == a[3], -b[1] - b[5] - b[30] == a[2], 
 2 + b[30] == a[1], b[33] + b[34] == a[54], 
 a[8] + a[60] - b[14] + b[16] - b[19] + b[21] + b[28] + b[35] == 
  a[51], b[31] - b[33] - b[35] - b[36] == 
  a[49], -a[42] + b[14] + b[17] - b[32] - b[33] - b[37] == a[38], 
 b[32] + b[33] - b[35] + b[37] == a[64], b[34] - b[38] == a[37], 
 b[37] + b[38] == a[10], b[35] - b[39] == a[59], 
 b[32] + b[36] - b[37] - b[39] == a[21], -b[36] + b[39] == a[29], 
 b[31] - b[38] + b[39] == a[43], 
 b[40] + b[41] == a[39], -b[40] - b[42] == a[56], -b[41] - b[42] == 
  a[15], -a[8] - a[60] - a[63] - b[11] + b[12] - b[14] - b[16] + 
   b[19] - b[21] + b[25] + b[39] - b[40] + b[41] + b[42] == a[22], 
 a[36] + b[23] + b[24] - b[25] + b[29] + b[33] + b[34] + b[44] - 
   b[48] == 
  a[50], -a[42] - a[52] - a[60] + b[3] + b[5] - b[6] - b[8] + b[14] + 
   b[17] + b[23] - b[25] + b[29] - b[31] - b[32] - b[37] - b[39] + 
   b[40] + b[42] + b[48] == a[27], b[44] + b[49] - b[51] == a[55], 
 b[48] + b[51] == a[41], 
 b[51] - b[52] == a[12], -a[6] - b[6] - b[7] + b[9] + b[11] + b[52] ==
   a[46], -a[6] - a[36] + a[60] - a[62] - a[65] + b[2] + b[3] + b[6] -
    b[7] + b[11] - b[24] - b[34] - b[38] - b[44] + b[52] == a[19], 
 b[49] + b[52] == 
  a[34], -a[30] + b[5] + b[10] + b[13] + b[20] - b[53] == 
  a[9], -b[49] - b[53] == a[24], -b[48] - b[52] - b[53] == 
  a[48], -a[6] - b[23] + b[53] == a[20], 
 b[44] - b[48] + b[53] == a[61], 
 a[13] + a[36] + b[20] + b[21] - b[25] - b[28] - b[29] + b[44] - 
   b[48] + b[53] == 
  a[57], -a[8] - a[13] - a[45] + b[1] - b[3] - b[9] + b[12] - b[20] + 
   b[30] - b[51] + b[53] == a[11], a[1] >= 0, a[2] >= 0, a[3] >= 0, 
 a[4] >= 0, a[5] >= 0, a[6] >= 0, a[7] >= 0, a[8] >= 0, a[9] >= 0, 
 a[10] >= 0, a[11] >= 0, a[12] >= 0, -a[13] >= -1, a[13] >= 0, 
 a[14] >= 0, a[15] >= 0, a[17] >= 0, a[19] >= 0, a[20] >= 0, 
 a[21] >= 0, a[22] >= 0, a[24] >= 0, a[26] >= 0, a[27] >= 0, 
 a[29] >= 0, -a[30] >= -1, a[30] >= 0, a[31] >= 0, a[32] >= 0, 
 a[34] >= 0, a[35] >= 0, a[36] >= 0, a[37] >= 0, a[38] >= 0, 
 a[39] >= 0, a[41] >= 0, a[42] >= 0, a[43] >= 0, 
 a[44] >= 0, -a[45] >= -1, a[45] >= 0, a[46] >= 0, a[47] >= 0, 
 a[48] >= 0, a[49] >= 0, a[50] >= 0, a[51] >= 0, a[52] >= 0, 
 a[53] >= 0, a[54] >= 0, a[55] >= 0, a[56] >= 0, a[57] >= 0, 
 a[59] >= 0, a[60] >= 0, a[61] >= 0, -a[62] >= -1, a[62] >= 0, 
 a[63] >= 0, a[64] >= 0, a[65] >= 0, b[1] >= -3, b[2] >= -3, 
 b[3] >= -4, b[5] >= -4, b[6] >= -5, b[7] >= -6, b[8] >= -3, 
 b[9] >= -4, b[10] >= 0, b[11] >= -4, b[12] >= -2, b[13] >= -4, 
 b[14] >= 0, b[16] >= -1, b[17] >= -2, b[19] >= -1, b[20] >= -2, 
 b[21] >= -1, b[23] >= -3, b[24] >= -1, b[25] >= -1, b[28] >= -1, 
 b[29] >= -1, b[30] >= -2, b[31] >= 0, b[32] >= 0, b[33] >= -1, 
 b[34] >= 0, b[35] >= 0, b[36] >= -1, b[37] >= 0, b[38] >= 0, 
 b[39] >= -1, b[40] >= 0, b[41] >= 0, b[42] >= -1, b[44] >= 0, 
 b[48] >= -1, b[49] >= 0, b[51] >= 0, b[52] >= 0, b[53] >= -1, 
 a[1] <= 1, a[2] <= 1, a[3] <= 1, a[4] <= 1, a[5] <= 1, a[6] <= 1, 
 a[7] <= 1, a[8] <= 1, a[9] <= 1, a[10] <= 1, a[11] <= 1, 
 a[12] <= 1, -a[13] <= 0, a[13] <= 1, a[14] <= 1, a[15] <= 1, 
 a[17] <= 1, a[19] <= 1, a[20] <= 1, a[21] <= 1, a[22] <= 1, 
 a[24] <= 1, a[26] <= 1, a[27] <= 1, a[29] <= 1, -a[30] <= 0, 
 a[30] <= 1, a[31] <= 1, a[32] <= 1, a[34] <= 1, a[35] <= 1, 
 a[36] <= 1, a[37] <= 1, a[38] <= 1, a[39] <= 1, a[41] <= 1, 
 a[42] <= 1, a[43] <= 1, a[44] <= 1, -a[45] <= 0, a[45] <= 1, 
 a[46] <= 1, a[47] <= 1, a[48] <= 1, a[49] <= 1, a[50] <= 1, 
 a[51] <= 1, a[52] <= 1, a[53] <= 1, a[54] <= 1, a[55] <= 1, 
 a[56] <= 1, a[57] <= 1, a[59] <= 1, a[60] <= 1, 
 a[61] <= 1, -a[62] <= 0, a[62] <= 1, a[63] <= 1, a[64] <= 1, 
 a[65] <= 1, b[1] <= 5, b[2] <= 6, b[3] <= 3, b[5] <= 5, b[6] <= 5, 
 b[7] <= 3, b[8] <= -1, b[9] <= 4, b[10] <= 3, b[11] <= 1, b[12] <= 4,
  b[13] <= 1, b[14] <= 4, b[16] <= 2, b[17] <= 1, b[19] <= 1, 
 b[20] <= 0, b[21] <= 2, b[23] <= 0, b[24] <= 2, b[25] <= 2, 
 b[28] <= 1, b[29] <= 1, b[30] <= -1, b[31] <= 2, b[32] <= 2, 
 b[33] <= 1, b[34] <= 2, b[35] <= 1, b[36] <= 0, b[37] <= 1, 
 b[38] <= 1, b[39] <= 1, b[40] <= 1, b[41] <= 1, b[42] <= 0, 
 b[44] <= 2, b[48] <= 0, b[49] <= 1, b[51] <= 1, b[52] <= 1, 
 b[53] <= 0} *)

Here is what we do next. With linear programming we can find all values for the a abd b variables (the 0-1 original variables and the null vector multipliers respectively) that give viable solutions. My belief is that when we solve relaxed LPs, that is, when inequalities but not integrality is enforced, it is best to branch on those null space multipler variables that have the smallest ranges. Empirically there is some evidence to suggest that this is the best way to branch (okay, one article on this is mine, but there are others as well). So we look for the first noninteger value from amongst the variables that we had reordered exactly for this purpose.
When a solution is found, we can now make new constraints that enforce that further solutions not be the same, and reenter the branch-and-prune loop.
Timing[sol = 
   Reap[Module[{vlen = Length[allvars], program, stack, badvar, vvals,
        val, counter = 0, var, len, eps = 1/10^5, soln = {}, cnstrnt, 
       newc, newv, pairdiff, vv, kk = 0, oldj, oneposns, newcons},
      program = constraints;
      stack = {program, {}};
      While[stack =!= {},
       counter++;
       program = stack[[1]];
       stack = stack[[2]];
       program = {1, program};
       Quiet[vvals = NMinimize[program, allvars]];
       If[Head[vals] == NMinimize, Continue[]];
       vvals = vvals[[2]];
       If[! FreeQ[vvals, Indeterminate], Continue[]];
       cnstrnt = program[[2]];
       badvar = 
        Position[allvars /. vvals, 
         a_ /; Abs[a - Round[a]] >= eps, {1}, 1, Heads -> False];
       If[badvar == {},
        kk++;
        If[Mod[kk, 100] == 0, Print["have soln ", {counter, kk}]];
        oneposns = 
         Cases[vvals, HoldPattern[a[_] -> aa_] /; Round[aa] == 1];
        oneposns = oneposns[[All, 1]];
        Do[
         newcons = 
          Join[Table[
            oneposns[[j]] == 1, {j, k - 1}], {oneposns[[k]] == 0}];
         stack = {Join[cnstrnt, newcons], stack};
         , {k, Length[oneposns]}];
        Sow[vals[[oneposns[[All, 1]]]]/lcm];
        Continue[];];
       badvar = badvar[[1, 1]];
       var = allvars[[badvar]];
       val = var /. vvals;
       stack = {Join[cnstrnt, {var <= Floor[val]}], stack};
       stack = {Join[cnstrnt, {var >= Ceiling[val]}], stack};];]][[2, 
    1]];]

After many many print statements, wherein we get around 160 solutions per minute, this finally ends.
(* {26562.945818, Null} *)

How long was that?
In[181]:= 26563./3600                                                           

(* Out[181]= 7.37861 *)

That's 7+ hours. Did we get the expected number of solutions?
Length[sol]                                                           

(* Out[182]= 69014 *)

Well, that much is good. Further work, using the equalities amongst variables that we found along the way, would be needed to recover all the actual sets of values.
So what are the salient features of this approach? Let me give my best effort at a summary.
(1) It is not an especially compelling method for the problem at hand. Moreover, as @Dr.Wolfgang Hintze observes in a different response, this would not scale at all well were we to allow reciprocals of integers arbitrarily large. I doubt it would scale were we even to allow them to be modestly larger. That said, we managed to recover a correct result in less than glacial time.
(2) We could have avoided the last code and just called Solve, specifying Integers for the domain over which to solve. Or could we? The problem is that Solve and the like will use exact linear programming. This is slow, but guaranteed to be correct. One cannot say the same for the approximate LP method shown above (though with the smallish integers involved, realistically it will not fail). One might ask why Solve does not do similarly. The answer is that we do not have a feel for what precision in the LP solver would be sufficient to guarantee correctness of the result.
(3) It can be a really good method when there are variables for which ranges are not very small (when we are not doing, say, 0-1 programming).
(4) One improvement I did not try would be to reduce to LinearProgramming in lieu of NMinimize. My past experience is that can reduce a run time by a factor of two or more (my past experience also indicates that this can increase the pain of programming by a factor of two or more).
(5) There could be ways to improve upon this, deep in the LP code, by "warm restarts" or use of cutting planes. This would not be trivial to implement though. Also it goes a ways outside of my own limited understanding of the topic.
(6) There are perhaps a few other improvements to be had that I simply never considered.

