How do I ensure a sequence of quaternions from matrix.decompose() is continuous?
After doing a little more research I think I have found a technique that will work.  Examining the matrices that are computed from the quaternion led me to conclude that q and -q both represent the same orientation, so I compare both q and -q to the previous quaternion orientation and pick whichever one is closest in 4-space.  The following udpated code sample illustrates the technique: 

import bpy
from math import *
from mathutils import *

def mission(obj):
    res = 36
    qs = QuaternionStabilizer()
    for z in range(res):
        theta = 2*pi *z / res
        mat = Matrix( [ [ cos(theta), sin(theta), 0],
                        [ -sin(theta), cos(theta), 0],
                        [ 0,0,1]]).to_4x4()
        (loc,quat,scale) = mat.decompose()

        obj.rotation_quaternion = qs.stabilize(quat)
        for ai in range(len(quat)):
            obj.keyframe_insert(frame=z*5, data_path="rotation_quaternion", index=ai)


class QuaternionStabilizer:
    def __init__(self):
        self.old=None

    def stabilize(self, q):
        if self.old is None:
            rval = q
        else:
            d1 = (self.old-q).magnitude
            d2 = (self.old+q).magnitude
            if (d1&lt;d2):
                rval = q
            else:
                rval = -q
        self.old = rval
        return rval

mission(bpy.context.active_object)


And now my fcurves don't look discontinuous:


This will not give good results if your orientations are flailing about wildly, but that's a case of Garbage In Garbage Out.