What is the idea behind expm1 to avoid cancellation error? It is well known that when x is close to 0, computing exp(x) - 1 introduces significant cancelation errors. As such, we have expm1 implemented in c99 and python. 
My question is how expm1 avoids cancellation error? Can anyone give me a general idea without too many nitty-gritty implementation details?
 A: The reason for making the expm1 function is to avoid precision loss when computing $e^x-1$ when $x$ is close to $0$. For large $x$ or very negative $x$ the expm1 function will give the same result as computing the usual exp and then subtracting $1$ so we will only focus on $x\approx 0$ where differences are found. But I can mention that a common way of computing exp (and a similar thing works for for expm1) for large $|x|$ is to write $e^{x} = (e^{x/2})^2 = ((e^{x/4})^2)^2$ and so on (this trick is explained more in the other answer). This allows us to compute exp for a much smaller argument, for which a Taylor series-like approach works well, and then do a few multiplications to get back the value we want.
Now to explain what happens close to $x = 0$ lets use the simplest way of computing $e^x$ for very small $x$, namely ${\rm exp}(x) = 1 + x$ as an example. Likewise our implementation of expm1 will be ${\rm expm1}(x) = x$. This is not far away from how an actual implementation would do it for the smallest $x$.
To understand why expm1$(x)$ is superior to exp$(x)-1$ close to $x = 0$ let's quickly review how floating point numbers works. They take the form $a\cdot 2^b$ where we have a certain number of bits to represent $a$ and $b$ (here $a$ is some number between $1$ and $2$ and $b$ some number between $-N$ and $N$). For the case of double precision numbers $a$ has a precision of $10^{-16}$ and $N\sim 1023$. This is very simplified, but it's all we need here.
The smallest number we can represent above $1$ will be $\approx 1 + 10^{-16}$ (the exponent $b=0$ so the precision is the precision of $a$). On the other hand the smallest number next to zero we can represent is $2^{-1023} \approx 10^{-308}$ (here $a=1$ and the precision is determined by $b$). This is the key point to understand: numbers we can represent aren't spaced evenly (like $0,\epsilon,2\epsilon,3\epsilon,\ldots, 1-\epsilon,1,1+\epsilon, \ldots$) but the size of the spacing between two numbers we can represent depends on which number you are looking at and it's much smaller close to $0$ than close to $1$.
Because of this if $x$ is less than $10^{-16}$ then $1+x$ simply evaluates to $1$ and $(1+x) - 1$ evaluates to $0$. If we want to compute $\frac{e^x-1}{x}$ using the $e^x$ function this will therefore give the very wrong result $\frac{0}{x} = 0$. If we on the other hand use the expm1 function then we have $\frac{(e^x-1)}{x} = \frac{x}{x} = 1$ which is correct (of course if $|x|<10^{-308}$ then this will also give $0$ as we cannot represent smaller numbers).
This was just a simple example, but it explains the main problem: we get huge cancellation errors from adding $1$ to something small and then later subtracting $1$ when computing $e^x-1$. To read more about how it's actually implemented see this note and the Google Go implementation (these links were given in comments by David K and Bungo).
A: Around $x=0$ the exponential is computed as some version of $1+xp(x)$ or $\frac{1+xg(x^2)}{1-xg(x^2)}$ or similar with some polynomials that approximate the exact term behind them in an uniform fashion on some interval. 
For other values the logarithm laws $e^x=(e^{x/2^m})^{2^m}$ or $e^x=2^k3^me^{x-k\ln2-m\ln3}$ are used to reduce the argument to some small value.
Then it is easy to modify those approximating expressions to return the expm1(x) value, in the mentioned cases $xp(x)$ or $\frac{2xg(x^2)}{1-xg(x^2)}$.
