Intro
One of the basic building blocks of impicit procedural modeling (such when building a distance field for raymarching based on basic primitives) is the union operator.
float opU( float d1, float d2 ) { return min( d1, d2 ); }This operator works great, but has the problem that the resulting shape has discontinuities in its derivatives. Or in other words, the resulting surface of unifying two smooth objects is not a smooth surface anymore. This is often inconvenient from a looks perspective, such as when trying to model organic shapes.
Regular min() based primitive union 
smoothmin() based primitive union 
The images above are a zoom into a procedurally modelled creature (which you can see in realtime online here: https://www.shadertoy.com/view/Mss3zM), where the legs are made of (deformed) cylinders and the body of a (deformed) sphere. When using the regular min() based primitive union, shown in the image to the let, the intersection between the legs and the body becomes apparent. In the right side image, the smoothmin based union explained in this article is used instead, resulting in a much more organic and visually pleasing connection between the parts.
Several implementations
The way to smoothly blend the shapes is to get rid of the discontinuity of the min() function, of course. But we want our smoothmin function to behave quite like min() when one of the two primitives is way further that the other. It's only in the area where the two values get similar that we want to apply the smoothness.
// exponential smooth min (k = 32);
float smin( float a, float b, float k )
{
float res = exp2( k*a ) + exp2( k*b );
return log2( res )/k;
}

// polynomial smooth min (k = 0.1);
float smin( float a, float b, float k )
{
float h = clamp( 0.5+0.5*(ba)/k, 0.0, 1.0 );
return mix( b, a, h )  k*h*(1.0h);
}

// power smooth min (k = 8);
float smin( float a, float b, float k )
{
a = pow( a, k ); b = pow( b, k );
return pow( (a*b)/(a+b), 1.0/k );
}

As a curiosity, it might be worth noting that the exponential and power based smoothmin functions, both generalize to more than two distances, so they are probably better suited for computing minimun distances to big sets of points beyond 2, for example when you want to compute smooth voronoi patterns or interpolate pointclouds. In the case of the power based smoothmin function, the expression a*b/(a+b) generalizes with the same formula as when computing the global resistance of N parallel resistors: 1/ ( 1/a + 1/b + 1/c + ... ), or in other words. For example, for three distances, you get a*b*c / (b*c + c+a + a*b).
Besides accepting an aribtrary number of points, another advantage of the exponential smin() over the polynomial smin() is that the exponential smin() is commutative when taken in pairs, while the polynomial is not. So, smin(a,smin(b,c)) is equal to smin(b,smin(a,c)) for the exponential, but not for the polynomial. That means that the exponential smin() allows one to process long lists of distances in any arbitrary order and slowly compute the smin, while the polynomial is ordering dependent.
EDIT 5 years later  I have learnt that Media Molecule used the polynomial smin() for their game "Dreams" (credited to Dave Smith), although their rewrote it to look like this:
// polynomial smooth min (k = 0.1);
float smin( float a, float b, float k )
{
float h = max( kabs(ab), 0.0 )/k;
return min( a, b )  h*h*k*(1.0/4.0);
}
As noted by Shadertoy user TinyTexel, this can be generalized to higher levels of continuity than the quadratic polynomail offers (C1), which might be important for preventing lighting artifacts. Moving on to a cubic curve gives us C2 continuity, and doesn't get a lot more expensive than the quadratic one anyways:
// polynomial smooth min (k = 0.1);
float sminCubic( float a, float b, float k )
{
float h = max( kabs(ab), 0.0 )/k;
return min( a, b )  h*h*h*k*(1.0/6.0);
}
Two functions 
The polynomial smooth min of the two funcitons 
Derivation
Deriving the polynomial smin() fuction is not difficult. The easiest is probably to start with the simplified version, and opearte in 1 dimension to makes things easier (the generalization to 2 or more dimensions or trivial):
smin( f(x), g(x), k ) = min( f(x), g(x) )  w(x,k)
If x=a is the point at which f(x)g(x) = k, and x=b is the point where f(x)g(x) = k, and a middle point x=c is where f(x)=g(x), then we know that
w(a)=0, w(b)=0, w(c)=s
with s being the maximum value. Since we want w(x) to be a smooth function that connects nicely to the curves f(x) and g(x) at the points x=a and x=b, we can choose
w(x) = s·hⁿ(x)
with h(x) = 1 ± [f(x)g(x)]/k
meaning h(x) is a power curve of degree n, with range between 0 and 1, and with the sings being positive if x < c and negative if x > c. This gives two versions of w(x), which we can call wl(x) and wr(x) for "left" and "right".
Since we want continuity also at x=c, we need to make sure that the derivatives of smin match when coming from either left or right of x=c. Therefore,
f'(c) + wl'(c) = g'(c) + wr'(c)
which means that
f'(c) + [f'(c)g'(c)]nS/k = g'(c)  [f'(c)g'(c)]·n·s/k
This can only be solved if 1+2·n·s/k = 0, which gives
s = k/2n
which is what we used in the quadratic and cubic implementations above.
The quadratic smin() though doesn't have second derivatives for w(x), meaning tha the only way the condition
f''(c) + wl''(c) = g''(c) + wr''(c)
can be met is making w(x) a cubic. Fortunately, since h(c)=1, the same condition 1+2·n·s/k = 0 needs to be met for ensuring continuity of the second (or any higher order) derivative.
Sumary
Properties of all the most useful smoothminimum functions:
Continuity  Commutative  Generalizes to N values  Expression  
Quadratic  C1  No  No  smin(a,b,k) = min(a,b)  h²/(4k), h = max( kabs(ab), 0 ) 
Cubic  C2  No  No  smin(a,b,k) = min(a,b)  h³/(6k²), h = max( kabs(ab), 0) 
Exponential  Cinf  Yes  Yes  smin(a,b,k) = ln( e^(k⋅a) + e^(k⋅b) )/k; 
Results
In general, the polynomial smoothmin function works very well, predictably and fast. It can be used to make surfaces connect, such as snow and archnitecture in the images below.
Regular min() based primitive union 
Polynomial smoothmin() based primitive union 
It obviously becomes very handy for connecting the different pieces of one same character, such as the fingers, hands, arms and body (wich in this case are made of cylinder/cones for the fingers and arms and a sphere for the body). The image bellow is interactive, and you can click and move the mouse over it to see the effect of the polynomial smoothmin function compared to the regular min function.
And finally, this is the result under animation of the body parts that are blended (note that texture coordinates are not fixed so the textures swim, but that's a different problem unrelated to the smoothmin, which works beautifully also with moving objects):