Hmm, it looks like an angle wrapping problem to me, although I can't figure out exactly where it's ending up backward. Maybe during the interpolation. But I would recommend using a fun technique, where angles are represented as a 16-bit value, with 65536 being a full revolution so it automatically wraps when it overflows. You can simply read the value as signed 16-bit if you need a -180 to +180 degree value rather than 0-360.
At the very least, it'll save you doing those wrapping if checks, and it will keep things wrapped during interpolation. You can subtract, then read as signed, and it naturally takes the shortest path positive or negative to get from one value to the other, as opposed to radians where you might end up trying to go from -180 to +179 degrees by going all the way around, rather than just one degree backward.
There's really not much to it, but here's my code. Scalar is just a typedef of float or double depending on what system I'm coding for, and s16/u16 are signed short/unsigned short. The ScalarConstant macro either does nothing, or adds an f, like 65536.0f. Oh, and sometimes I make Scalar a fixed-point value, although the multiplies in this particular file wouldn't work right unless using C++ with operator overloading...
I don't think it's a wrapping problem because it doesn't matter what combination of roll values the nodes have - it will still twist in that situation no matter what. I'm starting to think there's no simple fix for this and the calculations needed to do this are much more complicated than the basic ones I'm doing.