@Gaia Clary (gaiaclary) do you feel up to adding some unit tests for this function? That way it's much easier to document certain edge cases, and to show that this change doesn't have a negative impact on existing situations. The test could include the examples from T23954, T31333, T27675, and T30438, to ensure that these situations actually remain fixed.

I have lots of questions about this particular bit of code, and I hope now that you've dived into it you can help answer them.

due to float precision errors, we can have nor = (0.0, 0.99999994, 0.0)...

Why is this bad? What is theta? There is no comment that explains it. It's first defined as 1+nor[1] (so in this example theta is almost 2.0), but later it's redefined as nor[0]*nor[0] + nor[2]*nor[2] (so in this example it is 0.0).

I have found other situations where nor == (x, 0.99999994, z) with x and/or z very close to 0, but not exactly 0. However the test (nor[0] || nor[2]) can be true for tiny numbers.

A zero check for floats is usually bad, so changing those to some epsilon-around-zero check would likely be better, indeed.

Because of this i believe it is better to not only check (x or z != 0) but also to make sure that (nor[0] * nor[0] + nor[2] * nor[2]) is above THETA_THRESHOLD_NEGY_CLOSE.

This would create a tight coupling between that check, and the redefinition of theta as nor[0]*nor[0] + nor[2]*nor[2]. I think having a separate variable, with a more describing name, would be better. That way if the new variable is too small, the approach that uses the old variable can still be used. This would also make it possible to declare both theta and the new var as const.

Remark: I was once able to create a three bone situation by hand which shows the described failure, but i could not reproduce the situation ever since then. I guess the problem is mathematical unprecision. However i see this issue very frequently when using a generated rig (with python) where the bones are combined with a rather complicated set of constraints. Depending on the generator parameters occasionally the rig implodes due to this issue.

This makes it even more important to write some automated tests. These will just use a vector and a float as input, instead of requiring an entire rig.

This seems to make sense yes, but indeed would be nice to have some regression tests on this now...

Below is a python script that generates a rig that suffers from this bug. Run the script, then switch from edit mode to pose mode and see how "bad bone" disappears. Because its pose bone matrix is totally broken. The problem originates in the super small bone roll values which finally lead to the bug:

import bpy

bone_data = [
    ["good bone", (0.0, 0.0, 1.2038891315460205), (0.0, 0.01536799967288971, 1.3546786308288574), -4.1202467175263267e-16],
    ["bad bone", (0.0, 0.01536799967288971, 1.3546786308288574), (0.0, 0.0, 1.2038891315460205), -1.570796257510665e-07],
]

def make_armature(ctx, name):
    data = bpy.data.armatures.new(name)
    obj  = bpy.data.objects.new(name, data)
    ctx.collection.objects.link(obj)
    ctx.view_layer.objects.active = obj
    bpy.ops.object.mode_set(mode='EDIT')
    parent = None
    for bname, head, tail, roll in bone_data:
        bone = obj.data.edit_bones.new(bname)
        if parent:
            bone.parent = parent
        parent=bone
        bone.head=head
        bone.tail=tail
        bone.roll = roll
    return obj

armobj=make_armature(bpy.context, "armature")

I added a simple workaround to our Addon such that it just limits the precision of the bone roll to only 6 digits. This is enough to keep Blender from doing bad stuff. The correct fix though is to apply this differential to master :)