  We derive an exact expression for the probability density of work done on a
particle that diffuses in a parabolic potential with a stiffness varying by an
arbitrary piecewise constant protocol. Based on this result, the work
distribution for time-continuous protocols of the stiffness can be determined
up to any degree of accuracy. This is achieved by replacing the continuous
driving by a piecewise constant one with a number $n$ of positive or negative
steps of increasing or decreasing stiffness. With increasing $n$, the work
distributions for the piecewise protocols approach that for the continuous
protocol. The moment generating function of the work is given by the inverse
square root of a polynomial of degree $n$, whose coefficients are efficiently
calculated from a recurrence relation. The roots of the polynomials are real
and positive (negative) steps of the protocol are associated with negative
(positive) roots. Using these properties the inverse Laplace transform of the
moment generating function is carried out explicitly. Fluctuation theorems are
used to derive further properties of the polynomials and their roots.
