I want to do this with as few copies as possible

This should not be a concern while your function is incorrect (your "sign error"). But also: should this -

vec_a = (1, 0)
vec_b = (0, 1)

really have the opposite angle sign from this?

vec_a = (1, 0)
vec_b = (0, 2)

A vector's argument is invariant under uniform scaling, so I think the answer must be "no".

Other bits -

Don't use two calls to asarray. If you're making a library-like function that does array coercion, then you're better off issuing a single call to broadcast_vectors.

I really don't understand why your code is prepending a new axis. It shouldn't need to do that. I also don't understand why it's comparing vector norms. That doesn't make geometric sense to me.

Given the non-confidence you had in your code's accuracy, and me not really being able to follow half of it, I gave up in the attempt to reverse-engineer it and produced something that is probably slower, but that I understand and think is more likely to be correct.

import numpy as np


def angle_between(vec_a: np.typing.ArrayLike, vec_b: np.typing.ArrayLike, *, axis: int = -1) -> np.ndarray:
    """
    Computes the angle in signed radians from vec_a to vec_b in a right-handed frame.

    In all cases the magnitude of the angle follows simple cosine similarity.

    In the two-dimensional case, the angle sign is based on a simple vector-argument comparison.

    In the three-dimensional case, consider the plane shared by vec_a and vec_b. The sign of
    the plane's normal is chosen so that it minimises the angular deflection from an
    all-ones vector to disambiguate chirality. The angle sign is then the same as the 2D case
    for the two vectors projected to the plane oriented by that normal.

    Signed rotation demands a well-defined plane of rotation to disambiguate chirality.
    Since two vectors cannot fix a hyperplane in four-dimensional or higher space, those
    higher cases are all returned with a positive sign.
    """

    vec_a, vec_b = np.broadcast_arrays(vec_a, vec_b)

    # Simple cosine similarity
    abs_angle = np.acos(
        np.vecdot(vec_a, vec_b, axis=axis)
        /np.linalg.norm(vec_a, axis=axis)
        /np.linalg.norm(vec_b, axis=axis)
    )

    dim = vec_a.shape[axis]
    if 2 <= dim <= 3:
        if dim == 3:
            # normal vector orienting the plane intersecting a and b
            normal = np.linalg.cross(vec_a, vec_b, axis=axis)

            # coerce the normal to conventional orientation
            ones_basis = np.full(shape=normal.shape, fill_value=np.sqrt(1 / dim))  # use the normal closest to this
            normal_to_ones = np.acos(
                np.vecdot(normal, ones_basis, axis=axis)
                / np.linalg.norm(normal, axis=axis)
                # the ones basis is already a unit vector, so no second norm in denominator
            )
            normal *= 1 - 2*(normal_to_ones > 0.5*np.pi)  # -1 if > 180  # now the normal has conventional orientation

            # project to the plane
            u_basis = vec_a
            v_basis = np.linalg.cross(normal, vec_a, axis=axis)
            ax = np.vecdot(vec_a, u_basis, axis=axis)  # u-dimension in first vector
            bx = np.vecdot(vec_b, u_basis, axis=axis)  # u-dimension in second vector
            ay = np.vecdot(vec_a, v_basis, axis=axis)  # v-dimension in first vector
            by = np.vecdot(vec_b, v_basis, axis=axis)  # v-dimension in second vector
        else:
            indexer = [slice(None)] * len(vec_a.shape)  # to start, select all axes

            xidx = list(indexer)  # copy
            xidx[axis] = 0        # select first (x) position in axis of interest
            xidx = tuple(xidx)    # multi-axis indexer must be tuple
            ax = vec_a[xidx]      # x-dimension in first vector
            bx = vec_b[xidx]      # x-dimension in second vector

            yidx = list(indexer)  # copy
            yidx[axis] = 1        # select second (y) position in axis of interest
            yidx = tuple(yidx)    # multi-axis indexer must be tuple
            ay = vec_a[yidx]      # y-dimension in first vector
            by = vec_b[yidx]      # y-dimension in second vector

        ta = np.atan2(ay, ax)  # argument, first vector
        tb = np.atan2(by, bx)  # argument, second vector
        sign = np.ones_like(ta)      # sign vector, preserving shape and dtype
        sign[tb < ta] = -1
    else:
        return abs_angle  # implied sign of 1

    return sign*abs_angle



def unit_test() -> None:
    vec_a = (1, 0, 0)
    vec_b = (0, 1, 0)
    result = angle_between(vec_a, vec_b)
    assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
    result = angle_between(vec_b, vec_a)
    assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)

    vec_a = (-0.5, 0, 0)
    vec_b = ( 0,   3, 0)
    result = angle_between(vec_a, vec_b)
    assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
    result = angle_between(vec_b, vec_a)
    assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)

    vec_a = (1, 0)
    vec_b = (0, 1)
    result = angle_between(vec_a, vec_b)
    assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
    result = angle_between(vec_b, vec_a)
    assert np.isclose(result, -0.5*np.pi)

    vec_a = (1, 0)
    vec_b = (0, 2)
    result = angle_between(vec_a, vec_b)
    assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
    result = angle_between(vec_b, vec_a)
    assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)

    vec_a = (
        (1, 0),
        (2, 0),
        (3, 0),
    )
    vec_b = (
        (0, 2),
        (0, 1),
        (0, 0.5),
    )
    result = angle_between(vec_a, vec_b)
    assert np.allclose(result, (0.5*np.pi, 0.5*np.pi, 0.5*np.pi), atol=0, rtol=1e-14)
    result = angle_between(vec_b, vec_a)
    assert np.allclose(result, (-0.5*np.pi, -0.5*np.pi, -0.5*np.pi), atol=0, rtol=1e-14)


if __name__ == '__main__':
    unit_test()