Here is my second take on explaining my multiparent script.

What it does? Standard technique in rigging is linking. You link one object A to another object B so that its position and rotation of object A is fixed if we manipulate with object B. Quite often one would like to do converse, so manipulate with A and keep position and orientation of B fixed against A. Normally you have to destroy the original link(A to B) and link B to the A. Here I propose solution where you don't have to relink those objects. But everything comes with its price, so you loose direct control of those object but to that later.

Ok we start by developing such a link constraint between two objects A and B such that you can rotate center of A around center of B and vice versa. The orientation of A and B won't be affected by such rotation, this is necessary because you need some unconstrained controls.
To achieve that we have to introduce a helper object H to which will be centers of A and B fixed and rotation of A,B will somehow affect position and orientation of H.

Next we face difficulty that position of A,B at time $t$ does depend on rotation history of those objects. Have a look here for an example http://www.youtube.com/watch?v=lOfaFJq5Wqk as you can see rotation of A,B at time $t$ cannot fully describe the system at time $t$. You have to know the rotation history of A,B to fully determine position of those objects. So we will define exact positions of A,B  at some time $t_0$ and to find out position of A,B at time $t$ we will march frame by frame from time the $t_0$ to the time $t$.

In opposite to the previous post I won't use quaternions except when I will combine rotations but I will propose different approach with matrix exponential.

$$\begin{align} & x_A(t) \quad \text{position of object A} \\ & x_B(t) \quad \text{position of object H} \\ & x_H(t) \quad \text{position of object B} \\ & R_A^t \quad \text{rotation matrix of object A at time t} \\ & R_B^t \quad \text{rotation matrix of object B at time t} \\ & R_H^t \quad \text{rotation matrix of object H at time t} \end{align}$$

As already said position of centers of A and B are fixed against H, so
$$\begin{align} x_A(t) &= x_H(t) + R_H^tp_A \\ p_A &= (R_H^{t_0})^{-1}(x_A(t_0) - x_H(t_0))\\ x_B(t) &= x_H(t) + R_H^tp_B \\ p_B &= (R_H^{t_0})^{-1}(x_B(t_0) - x_H(t_0)) \end{align}$$

Now if we rotate A a little bit we want to move with B, but B is fixed against H so it is sufficient to move with H. So just for a moment forget about B and imagine that H is linked to the A. Than if we rotate with A position of H will change as
$$\begin{align} x_H( t + \Delta t )  = x_H( t ) + R^{t+\Delta t}_A (R^t_A)^{-1} (x_H(t)-x_A(t)) \end{align}$$
and the orientation of H will change as
$$\begin{align} R_H^{t + \Delta t } =  R^{t+\Delta t}_A (R^t_A)^{-1} R^t_H \end{align}$$
Now if we take account the B's contribution we get equation

$$\begin{align} x_H( t + \Delta t )  = x_H( t ) + R^{t+\Delta t}_A (R^t_A)^{-1} (x_H(t)-x_A(t))  + R^{t+\Delta t}_B (R^t_B)^{-1} (x_H(t)-x_B(t)) \end{align}$$
But there is trouble with orientation update  of H because matrix multiplication is not commutative so we don't know in which order to multiply those matrix. For now we just write something and we will  explain later what $\{ \cdot , \cdot \}$ means.
$$\begin{align} R_H^{t + \Delta t } =  \{R^{t+\Delta t}_A (R^t_A)^{-1}, R^{t+\Delta t}_A (R^t_A)^{-1} \}R^t_H \end{align}$$

What  $\{ \cdot , \cdot \}$ does is that it takes two rotations and produce another new rotation, which somehow captures those two rotations. From this we require from  $\{ \cdot , \cdot \}$ these identities
$$\begin{align} \{ R_1 , R_2 \}  &= \{ R_2 , R_1 \} \\ \{R_1, I \} &= R_1 \end{align}$$
where $I$ is identity(ie no rotation). First says that it does not depend on the order of the matrices and the second say that if we combine some rotation with identity(ie no rotation) we should get the original rotation.

Now to define  $\{ \cdot , \cdot \}$  without quaternions we need to know a little bit about rotations and matrix exponential.

If you have rotation around axis $n$ by angle $\omega$. Than its rotation matrix $R$ can be expressed as
$$R  = e^{\omega [n]_\times} = \sum_{k=0}^\infty \frac{ \omega^k}{k!} [n]_\times^k$$
where $[n]_\times$ is cross-product matrix .
Other way around if you have antisymmetric matrix $A$ that $e^A$ is rotation matrix.

We are ready to define  $\{ \cdot , \cdot \}$. Let's have two rotation matrices $R_1 = e^{A_1}, R_2 = e^{A_2}$. Than
$$\{ R_1, R_2 \} =  \{ e^{A_1}, e^{A_2} \}  = e^{A_1+A_2}$$
Observe that $A_1+A_2$ is again antisymmetric matrix so $e^{A_1+A_2}$ is rotation matrix. Next $A_1+A_2 = A_2+A_1$ therefore  $\{ e^{A_1}, e^{A_2} \} =  \{ e^{A_2}, e^{A_1} \}$. Lastly $I = e^0$ so  $\{ e^{A_1}, e^0 \} = e^{A_1 + 0} = e^{A_1}$. So  $\{ \cdot , \cdot \}$  satisfy all identities we wanted.

How to program this then.

What are the inputs?
We have to specify time $t_f$ at which we want to get positions and orientations of A,B,H. Next we have to know positions of A,B,H at time $t_0$ and orientation of H at time $t_0$. Than we have to know the whole history of rotation matrices of A,B from time $t_0$ to time $t_f$.
What is the output?
Position of A,B,H at time $t$ and orientation of H.

So the code would be something like this:

1. precalculate values  $p_A,p_B$

2. then use update equations, start at time $t_0$ and and at time $t_f$.
$$\begin{align} x_H( t + \Delta t )  &= x_H( t )  - R^{t+\Delta t}_A (R^t_A)^{-1} R_H^t p_A  - R^{t+\Delta t}_B (R^t_B)^{-1}R_H^t p_B\\ R_H^{t + \Delta t } &=  \{R^{t+\Delta t}_A (R^t_A)^{-1}, R^{t+\Delta t}_A (R^t_A)^{-1} \}R^t_H \end{align}$$

3. from valuse $x_H(t_f), R_H^{t_f}$ calculate $x_A(t_f),x_B(t_f)$