A natural number $p$ is prime if and only if it has exactly two solution sets $\{x,y\}$ to the equation $$\frac1p = \frac1x + \frac1y, \text{ or equivalently, }p=\frac{xy}{x+y}.$$
Those solution sets are $\{p,p(p+1)\}$ and $\{2p\}$ (where $x=y=2p$).
Is this a known fact about prime numbers? If so, could somebody provide a reference?
A proof is given below.
Proof: $(\implies)$ Suppose $p$ is prime. Then if $$p=\frac{xy}{x+y},$$ we can express $x=p^ak$ and $y=p^bl$ for $k$ and $l$ coprime to $p$. Without loss of generality, suppose $0\leq a\leq b$ and $1\leq b$. Then we have $$p=p^b\frac{kl}{k+lp^{b-a}}.$$
Suppose that $b>a$. We then have that $p\not\mid k+lp^{b-a}$, and therefore $b=1$ and $a=0$. Therefore $kl=k+lp$, which implies $k(l-1)=lp$. Thus, $p\mid l-1$.
Consider any natural number $m\mid k$. Then $m\mid k(l-1)=lp$. Since $p\not\mid k$, we have $p\not\mid m$. Thus, $m\mid l$.
Similarly, for any natural $n\mid l$, we have $n\mid l(k-p)=k$.
So $l$ and $k$ have the exact same divisors, and thus $l=k$. Thus,
$$k^2=k+kp=k(p+1),$$ meaning $k=p+1$. So $x=p+1$, which implies $y=p(p+1)$. This is the first solution set.
Next, consider the case in which $b=a$. Then $$p=p^b\frac{kl}{k+l},$$ so $$p^{1-b}=\frac{kl}{k+l}.$$
Since $k$ and $l$ are coprime to $p$, it follows $b>0$. We also have $$\frac{1}{p^{b-1}}=\frac1k+\frac1l,$$ implying that $$\frac{1}{p^{a+b-1}}=\frac{1}{p^{2b-1}}=\frac1x+\frac1y,$$ so $2b-1=1$ and thus $a=b=1$. Finally, that implies $$p=p\frac{kl}{k+l},$$ so $kl=k+l$ and therefore $k=l=2$. This gives the solution $x=y=2p$. So if $p$ is prime then it has two solution sets.
$(\Longleftarrow)$ If $p$ is not prime, then $p=qr$ where $q$ is prime and $r\neq 1$. Then $$\frac1p=\frac1{p+r}+\frac1{q(p+r)}$$ gives us another solution set.