Need Help Point to Point Implementation

[sukrit1234]

I'm just undergraduate student and I'm learning about rigid dynamics
and now I try to implement constraint of multibody begin with point to point constraint(P2P)

my P2P implement
first I set the P2P point and transform it to both rigidbody coordinate
p0 and p1

then I assume p2p has 3 constraint impulse axis X , Y, and Z axis
so i compute 3 jacobian entries for each axis like that

b0,b1 is rigid body
r0,r1 is vector from center of mass to constraint point
a is unit vector that impulse will apply

I compute jacobain for each axis a by

terms0 = b0->invMass + b1->invMass
terms1 = (b0->invIWorld.mul(r0.cross(axis))).cross(r0).dot(axis)
terms2 = (b1->invIWorld.mul(r1.cross(axis))).cross(r1).dot(axis)

so jac = term0 + terms1 + terms2

then I transform both P2P point (p0 and p1) to world coordinate gp0 and gp1
and I compute velocity at both point with vp = v + (wxr)

so i get relative velocity by v_rel = vp0 - vp1 (b0 relative to b1)
and i get position error by p_err = gp0 - gp1
compute bias(position correction data) by

bias = (p_err * bias_factor)/dt

bias_factor i use 0.2f (from Box2D Lite)

now I compute impulse
apply to both rigid bodies
(compute all three axis X Y and Z so lamda is vector of impulse

vector3 impulse(0,0,0);
for(int i = 3 ; i < 0 ; ++i)
{
impulse = (bias - v_rel)/jacc;
/*now I apply impulse to both rigid body like this:*/
b0->linear_v += b0->invMass * impulse
b1->linear_v -= b0->invMass * impulse
b0->ang_v += b0->invIWorld * r0.cross(impulse);
b1->ang_v -= b1->invIWorld * r1.cross(impulse);
impulse = 0;
}
and my simuation loop is

1. apply external forces to rigid body
2. update verocity with these forces
3. compute jacobian and bias
4. compute impulse and apply to rigid body velocity
with iterative method (i get this idea from GDC 2008 slide)
5 update position

ok now i test my implementation
first i initial rigid body position as bellow picture

[attachment=2]initial%20position.JPG[/attachment]

and then i take the force (0,0,1000) to lowest rigid body
ok it look great that i expect

[attachment=1]pic1.JPG[/attachment]

and then
i take tourque (0,100,0) and now problem occur like this

[attachment=0]ploblem.JPG[/attachment]

need to know what's happen but
but i think (I try to debug two or htree days ) that position correction impulse
(from bias/jac ) and velocity error impulse (-v_rel/jac) is opposite each other
when angular velocity has rotaion more than one axis but I don't know how to fix it
help me please ><