function _rotationMatrix3x3 (theta, v, format) {
    const normV = norm(v)
    if (normV === 0) {
      throw new RangeError('Rotation around zero vector')
    }

    const Big = isBigNumber(theta) ? BigNumber : null

    const one = Big ? new Big(1) : 1
    const minusOne = Big ? new Big(-1) : -1
    const vx = Big ? new Big(v.get([0]) / normV) : v.get([0]) / normV
    const vy = Big ? new Big(v.get([1]) / normV) : v.get([1]) / normV
    const vz = Big ? new Big(v.get([2]) / normV) : v.get([2]) / normV
    const c = cos(theta)
    const oneMinusC = addScalar(one, unaryMinus(c))
    const s = sin(theta)

    const r11 = addScalar(c, _mul([vx, vx, oneMinusC]))
    const r12 = addScalar(_mul([vx, vy, oneMinusC]), _mul([minusOne, vz, s]))
    const r13 = addScalar(_mul([vx, vz, oneMinusC]), _mul([vy, s]))

    const r21 = addScalar(_mul([vx, vy, oneMinusC]), _mul([vz, s]))
    const r22 = addScalar(c, _mul([vy, vy, oneMinusC]))
    const r23 = addScalar(_mul([vy, vz, oneMinusC]), _mul([minusOne, vx, s]))

    const r31 = addScalar(_mul([vx, vz, oneMinusC]), _mul([minusOne, vy, s]))
    const r32 = addScalar(_mul([vy, vz, oneMinusC]), _mul([vx, s]))
    const r33 = addScalar(c, _mul([vz, vz, oneMinusC]))

    const data = [[r11, r12, r13], [r21, r22, r23], [r31, r32, r33]]

    return _convertToFormat(data, format)
  }
Returns a 3x3 matrix of rotation of angle theta around vector v

@param {number | BigNumber | Complex | Unit} theta The rotation angle
@param {Matrix} v                                  The rotation axis vector
@param {string} format                             The storage format of the resulting matrix
@returns {Matrix}
@private