[Question]
  [
An "Egyptian fraction" is a list of *distinct* fractions with a numerator of \$1\$. For example:
\$
\frac 1 1+ \frac 1 2 + \frac 1 3 + \frac 1 6 
\$
The "size" of an Egyptian fraction is just the number of terms involved.
Your task is to take a positive integer \$n\$ and output the smallest Egyptian fraction that sums to \$n\$. In the case of ties (there are ties) you may output any of the tied fractions.
Since the numerator is always 1 you may output your answer as a list of denominators instead.
This is [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'") so the goal is to minimize the size of your source code within the bounds of your chosen language and algorithm with answers being scored in bytes.
Fast answers are encouraged but must take steps to achieve the goal such as removing unnecessary whitespace or giving the variables short names. Feel free to include **(fast)** in your answer title along with the language and byte count if you do give it a go for speed.
## Precision
While our default precision is that the algorithm needs to be *theoretically* correct, if you try to do this challenge with floating point numbers you will get results that are wrong almost immediately. It's also not at all difficult to solve this challenge using only integers. For this reason in addition to being theoretically correct answers must also give correct calculations for \$n < 4\$.
As an additional clarification *theoretically* correct means that you cannot for example assume some arbitrary bound on the size of the denominator to restrict your search space. Any bounds used to restrict the search space must be theoretically justified.
## Test cases
The sizes for the first 3 solutions next to a potential output are:
\$
1, \frac 1 1\\
4, \frac 1 1+ \frac 1 2 + \frac 1 3 + \frac 1 6 \\
13, \frac 1 1+\frac 1 2+\frac 1 3+\frac 1 4+\frac 1 5+\frac 1 6+\frac 1 7+\frac 1 8+\frac 1 9+\frac 1{15}+\frac 1{18}+\frac 1{40}+\frac 1{42} \\
\$
The first two I have calculated and verified by hand. The third one was verified by a computer program that ran for about 5 seconds. Since the search space is infinite, I am not 100% sure that my program doesn't miss some potential solution. But I have hand verified that there are no solutions up to 11. And I do have a solution of size 13.
My program doesn't find a solution for 4 in a reasonable amount of time but it is easy to check that the answer is at least 31 terms long, and my program confirms that it is at least 32 terms long.
  
[Answer]
# [Haskell](https://www.haskell.org/), 103 bytes
```
head.([1..]>>=).g 1 0
g b i a d=[[]|a<1]++[j:o|j<-[max(div(b-1)a)i+1..div(b*d)a],o<-g(b*j)j(a*j-b)$d-1]
```
[Try it online!](https://tio.run/##HYuxDoIwFEV3vuINDi3QxhonAyzORkaT0piHhdIKhYgxDny7Fd3uOTm3w/ne9H0Y0HrIYcDpdAVSeWAFTA/rn0B8Ci14SkEKzvcqateuCl2DmpOfUkWRU25AwDYyUIMFBJ1LqRbMhEoS6Q7j4jImB3wTbV@kZoIitcn6/WOsKap0zJhZt6OOYOxYTTeaCRXC59b2aObALseyDOy8@wI "Haskell – Try It Online")
This uses iterative deepening depth-first search, using the bound \$\frac a{bd} ≤ \frac 1j ≤ \frac ab\$ on the first term \$\frac1j\$ of a size-\$d\$ Egyptian fraction for \$\frac ab\$. Unclear if this qualifies as “fast”, but it does solve \$n = 4\$ in a few seconds:
\begin{split}
4 = {}&\frac11 + \frac12 + \frac13 + \frac14 + \frac15 + \frac16 + \frac17 + \frac18 + \frac19 + \frac1{10} + \frac1{11} + \frac1{12}
+ \frac1{13} \\ &+ \frac1{14} + \frac1{15} + \frac1{16} + \frac1{17} + \frac1{18} + \frac1{19} + \frac1{20} + \frac1{21} + \frac1{22} + \frac1{23} + \frac1{24} \\ &+ \frac1{25} + \frac1{26} + \frac1{27} + \frac1{28} + \frac1{30} + \frac1{34} + \frac1{100} + \frac1{11934} + \frac1{14536368}.
\end{split}
[Answer]
# [Python 3](https://docs.python.org/3/), 159 130 125 bytes
```
f=lambda x,s={1},i=1,p=1,l=[]:i<7**x**x**x and f(x,s^{i},1+i*(i in s),p*i,l+[s]*(sum(p//y for y in s)==x*p))or min(l,key=len)
```
[Try it online!](https://tio.run/##NY1BCoMwFET3PcWHbpL4QYILQZpLdCsW0hrrpzEGY6lBPLvVSmFmMzPM83Fse5eta6Os7u61hgmDmuWCpCT6zVaVVUGXXIjpEGhXQ8O23W2mBWVCghGQg8DRC0KblKESLLw75tM0QtMPEI9eqUl4zregI8csvkxU1ji@nuFqvNUPA/kfs5M@NLaQn/aHH2HQ7mmYxIwXJ/ADuZERNow4X78 "Python 3 – Try It Online")
Bruteforces trough all Egyptian fractions, in order of maximum denominator. A limit of \$7^{x^{x^x}}\$ is used to stop the program. Then the shortest fraction is picked.
Now the elephant in the room is the upper bound \$7^{x^{x^x}}\$. This corresponds to the maximum denominator in an optimal solution. Let's see how this bound is obtained.
Let's first make an upper bound for the length of the Egyptian fraction. To construct the egyptian fraction for `x`, we first greedily pick unit fractions, until the next one would get us over `x`. For example, if `x` is 2, then we would do \$\frac{1}{1}+\frac{1}{2}+\frac{1}{3}\$, stopping before \$\frac{1}{4}\$, since that would result in a sum greater then 2. Note that this initial tail has length less than \$e^x\$.
Now we have \$2-(\frac{1}{1}+\frac{1}{2}+\frac{1}{3})=\frac{1}{6}\$. Here we got lucky, and the difference is already a unit fraction. If it isn't, we can use the [greedy algorithm](https://en.wikipedia.org/wiki/Greedy_algorithm_for_Egyptian_fractions) to turn the difference into a unit fraction. The amount of fractions that the greedy algorithm introduces is at most the denominator of the difference. The denominator is less than the numerator which is less than \$e^x!\$.
Finally, if we have an integer egyptian fraction of length `l`, the numerator can be at most `l^l`. To see why, note that the numerator of the last fraction can't be greater than the product of the other numerators.
Now we have the upper bound \$(e^x!)^{e^x!}\$.
When \$x\ge 4\$, \$(e^x!)^{e^x!}<7^{x^{x^x}}\$. The seven takes care of values of \$x\$ less than 4.
Proof for \$(e^x!)^{e^x!}<x^{x^{x^x}}\$:
\$(e^x!)^{e^x!}<((e^x)^{e^x})^{(e^x){e^x}}=(e^{xe^x})^{e^{xe^x}}=e^{xe^xe^{xe^x}}=e^{xe^{x+xe^x}}\$
\$x+xe^x< x^x\$
\$xe^{x+xe^x} < x^{x^x}\$
\$e^{xe^{x+xe^x}}<7^{x^{x^x}}\$
[Answer]
# [Python 3](https://docs.python.org/3/), ~~180~~ 161 bytes
```
def f(n):
 u,*r=0,(0,1,n,1)
 for l,q,a,b,*x in r:
  if a==0:
   if(l>u>0)<1:u=l;o=x
  if u<1or(u-l)*b>=a*q>0:r+=(l+1,q+1,a*q-b,q*b,*x,q),(l,q+1,a,b,*x)
 return o
```
[Try it online!](https://tio.run/##JY7RCoQgEEXf@4p5VJtA6a1N/8XY2g1kzCGh/XpX62HgcrnMOcfv/EYaS3mvG2yC5NRBRsVWo9BokNDIDrbIEDChxwXVBTsB1x3sG3hrdYs1i@Cy03I2U7bhFe31LPJsIos8BKkWZ71KTk/cWxF6g6lebYYFk2qfMUkU4alvVGXzemYmiKVJ0M329FmFQRibLRy80ykIb31Z/g "Python 3 – Try It Online")
This does a kind of breadth-first search over sets of fractions. When it finds a first valid representation, it uses the length of that as an upper bound.  
With the upper bound we can cut search branches when the target number can't be reached with the remaining number of fractions.
This only uses integer arithmetic, given enough memory and time it should be able to calculate any value.
Local results:
```
$ time python3 ef.py
1 [1]
2 [1, 2, 3, 6]
3 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 230, 57960]
python3 ef.py  15.89s user 2.03s system 98% cpu 18.148 total
```
---
At the cost of ~60 bytes we can use sharper bounds for cutting the search, skip fractions that would lead to higher sums, and traverse the search tree in a depth-first order. This allows us to get \$n=4\$ on TIO:
```
def f(n):
 u,*r=0,(0,1,n,1)
 while r:
  l,q,a,b,*x=r.pop()
  if a*q==b:
   if(l+1>u>0)<1:u=l+1;o=x+[q]
  elif u<1or(u+~l)*b>=a*q>0:A=a*q-b;B=q*b;D,M=divmod(B,A);r+=(l,q+1,a,b,*x),(l+1,max(q+1,D+(M>0)),A,B,*x,q)
 return o
```
[Try it online!](https://tio.run/##LY9NboQwDIX3cwovHeKpiKpuhgkSaLZzgqoLEKGDlEmCRVrY9OrUVF355z3r80vb8ojh9Zy2tO374EYYMajLCTIVbEvCkgwFMuoE34/JO2DRwNNMHfVUrJZfUkwoMkwjdMVsbX84ZEKvTZ3rUl3NJVsZqmhX/T5/iOy8uPPVRMasf7wq@trKcV1emqOe@6q1c9FXN7rbYfp6xgFbalTF2qLAtfnHKzoo9OxWPJY3jXcBKmqoFZVm@YvdkjlA3MfIEGAKwF34dGgI3o6gwGD/QkubeAoLBgLvArIiYLX/Ag "Python 3 (PyPy) – Try It Online")
[Answer]
# Python3, 386 bytes:
```
def f(n):
 A,I=[1],[]
 while sum(1/i for i in A)+1/(A[-1]+1)<=n:A.append(A[-1]+1)
 q = [(A[:],n,0)]
 while 1:
  c,d,k = q.pop(0)
  if(s:=sum(1/i for i in c))==n:return c
  F=1
  if d:
   l= c[-1]
   while (s+1/(l:=l+1))>n or l in I:
    if l>10e7:I.append(A.pop());q=[(A[:],n,0)];F=0;break
   if F:q.append((c+[l],d-(k:=l>(c[-1]+1)),k))
  if k and F:
   q.append((c[:-1]+[c[-1]+1],d,k))
```
[Try it online!](https://tio.run/##ZZDNboMwEITvPMUevcJJcHNoZWokLkg8g8WB8qMgXGN@oqpPT9eQREi97Vrzzcza/S63wV4/3LSuddNCyyzKAFKeKy0KrosAfm6daWC@fzNx6aAdJuigs5BiKC4s1SdRhAI/lZXpuXSusfXrMYARFGjaZcEtj/DlJigDKl7zngTj2Q2ORSSHrmWzVP@yKkRFAVOz3CfaSJgpscmh9k5gFFQ@1M97Apt9PSOVoSKYWCAv473yDfCoSUTUvMv8VXvrgRiP6tg5zlQUf01N2Qc7l8nxibAq1Kbg9Yn1lJSw6nE58h73e6CH0tbEePjAaemV@gEU/isQVzd1dmEtI4fgOb8d5itp/gA)
Local results (Macbook Pro):
```
1 [1]
2 [1, 2, 3, 6]
3 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 230, 57960]
4 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 31, 32, 33, 34, 45, 7397, 82673165]
python3 egyptian_tests.py  29.14s user 0.15s system 99% cpu 29.480 total
```
]