Why is my rotation method sometimes incorrect for parented objects

I wrote Python’s function to rotate a object ‘o‘ direction, step by step, to ‘directionW‘ direction.

When the o has no parent, the function works great. But, when the o has one or more parent, function sometimes works correctly, sometimes wrong. When it works wrong the direction of o don’t going to directionW by shortest way, but going in space by way (like 8-digit) around directionW.

Next is the function:

    def RotateToDirectionW( o, directionW, step, adaptivestep, maxunlock ):
    """o - object to step rotate to directionW"""
    """directionW - direction to achieve step by step (unit vector in world coordinates)"""
    """adaptivestep - if True, then rotation angle has minimal value when direction of o near to directionW,"""
    """maxunlock - maximal vector size to unlock, 1E-5 or 0"""

    matLtoW = o.matrix_world
    matWtoL = matLtoW.copy(); matWtoL.invert()
    if o.parent != None:
        matWtoP = o.parent.matrix_world.copy(); matWtoP.invert()
    else:
        matWtoP = Matrix.Identity(4)

    # direction of o ( axe z+ ) in world coordinates
    dirW = ( matLtoW @ Vector((0,0,1)) ) - ( matLtoW @ Vector((0,0,0)) )
    dirW.normalize()        # need to normalize angle

    # calculate axe, around what object will be rotated (in world c.)
    axeRotW = dirW.cross( directionW )
    
    power = axeRotW.length
    
    # unlock
    # if vector is too small to rotate, take any bigger
    if abs(power)<maxunlock:
        axeRotW = Vector( (maxunlock,maxunlock,maxunlock) )
        power = axeRotW.length

    # transform axe of rotation to parent base coordinates
    # (problem somewhere there)
    axeRotP = (matWtoP.to_quaternion() @ axeRotW

    if adaptivestep:
        angle = step * power
    else:
        angle = Vector.angle( dirW, directionW )
        if (angle > step): angle = step

    # rotation
    qutRot = Quaternion( axeRotP, angle )
    o.rotation_euler.rotate( qutRot )

I check every possible intermediate results, they are correct at all. The axe of rotation is correct in world space. I think, problem is when axeRotW converted to parent coordinate space. It look like there is transitional transformation between parent’s transformation and o own transformation. I think, i need convert axeRotW to this transitional coordinate system, but there is no information about it.

What is a way to solve the problem? I waste several day of my time to do it, but now i into deadlock.

P.S.: I try, also, to rotate o by bps.transform.rotate(), but this operation can’t work with custom rotation axe, XYZ-only.

P.P.S.: A founded 3 conditions for the correct result:

• 1. The object o has not parent
• 2. The object has parent object with identity (equal 1) ‘matrix_parent_inverse’.
• 3. The object has parent 3-verticies with same direction of coordinates as parent, and identity ‘matrix_parent_inverse’.

P.P.P.S.: To investigate matrix chain calculation, i run script from blender.stackexchange.com/a/160580/15543

import bpy
from mathutils import Matrix

context = bpy.context
ob = context.active_object
M = Matrix()
print("==== calc word matrix ====")
print(ob.matrix_world)
while ob.parent:
    M = (ob.matrix_parent_inverse @ ob.matrix_basis) @ M
    ob = ob.parent

M = ob.matrix_basis @ M

print(M)

If generic child object was selected, then script prints two equal matrices. But, if somewhere in parent-child chain exists 3-vertices parenting, then the printed matrices was different!!!. F.e.:

---- calc word matrix ----
<Matrix 4x4 (-0.0198,  0.1843,  0.4570, 15.1996)
            ( 0.0012,  0.4574, -0.1844, -4.9117)
            (-0.4928, -0.0063, -0.0188,  3.4083)
            ( 0.0000,  0.0000,  0.0000,  1.0000)>
<Matrix 4x4 (-0.0465, -0.2971,  1.1550, 36.2083)
            (-0.0120,  1.1559,  0.2969, 11.8968)
            (-1.1926,  0.0000, -0.0480,  1.8346)
            ( 0.0000,  0.0000,  0.0000,  1.0000)>

It is possible, there is one addition transformation in parent-child chain. And i need to take into account yellow colored block:

How i see, there is no python access to extra v3-parenting transformation matrix. I suppose, with help of some math we can calculate this extra matrix. But how?

matrix_world = parent.matrix_world * parent.matrix_x * matrix_parent_inverse * matrix_basis

# transform axe of rotation to parent base coordinates
matWtoP = o.matrix_basis @ o.matrix_world.inverted()
axeRotP = matWtoP.to_quaternion() @ axeRotW