<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def axis_rotation(axis, angle, vectors, axis_shift=(0, 0, 0)):
 """Rotate a vector or an array of vectors around an axis in 3d. The rotation is computed by `Rodrigues' rotation formula`_. Parameters axis : `array-like`, shape ``(3,)`` Rotation axis, assumed to be a unit vector. angle : float Angle of the counter-clockwise rotation. vectors : `array-like`, shape ``(3,)`` or ``(N, 3)`` The vector(s) to be rotated. axis_shift : `array_like`, shape ``(3,)``, optional Shift the rotation center by this vector. Note that only shifts perpendicular to ``axis`` matter. Returns ------- rot_vec : `numpy.ndarray` The rotated vector(s). References .. _Rodrigues' rotation formula: https://en.wikipedia.org/wiki/Rodrigues'_rotation_formula Examples -------- Rotating around the third coordinate axis by and angle of 90 degrees: True True The rotation can be performed with shifted rotation center. A shift along the axis does not matter: True The distance between the rotation center and the vector to be rotated determines the radius of the rotation circle: True True Rotation of multiple vectors can be done in bulk: True True """

    rot_matrix = axis_rotation_matrix(axis, angle)
    vectors = np.asarray(vectors, dtype=float)
    if vectors.shape == (3,):
        vectors = vectors[None, :]
    elif vectors.ndim == 2 and vectors.shape[1] == 3:
        pass
    else:
        raise ValueError('`vectors` must have shape (3,) or (N, 3), got array '
                         'with shape {}'.format(vectors.shape))

    # Get `axis_shift` part that is perpendicular to `axis`
    axis_shift = np.asarray(axis_shift, dtype=float)
    axis = np.asarray(axis, dtype=float)
    axis_shift = axis_shift - axis.dot(axis_shift) * axis

    # Shift vectors with the negative of the axis shift to move the rotation
    # center to the origin. Then rotate and shift back.
    centered_vecs = vectors - axis_shift[None, :]
    # Need to transpose the vectors to make the axis of length 3 come first
    rot_vecs = rot_matrix.dot(centered_vecs.T).T
    return axis_shift[None, :] + rot_vecs