ENH: `broadcast_arrays()` with option to not repeat arrays

Proposed new feature or change:

Consider a function $f(r, \theta) = [r \cos \theta, r \sin \theta]$ and when r.shape == (3, 1) and theta.shape == (4,) thus the shape of the result is supposed to be (2, 3, 4) (or (3, 4, 2)) by broadcasting.

import numpy as np


def ferror(r, theta): # (3, 1) and (4,)
	cossintheta = np.stack([np.cos(theta), np.sin(theta)]) # (2, 4) (or (4, 2) if axis=-1)
	return r * cossintheta # (3, 1) and (2, 4) (or (4, 2)) -> error

def fslow(r, theta): # (3, 1) and (4,)
	r, theta = np.broadcast_arrays(r, theta) # (3, 4) and (3, 4)
	cossintheta = np.stack([np.cos(theta), np.sin(theta)]) # (2, 3, 4)
	return r * cossintheta

def ffast(r, theta): # (3, 1) and (4,)
	r, theta = np.broadcast_arrays(r, theta, no_repeat=True) # (3, 1) and (1, 4)
	cossintheta = np.stack([np.cos(theta), np.sin(theta)]) # (2, 1, 4)
	return r * cossintheta

Where broadcast_arrays(*arrays, no_repeat=True) does

• check if the arrays are broadcastable
• for each array, add np.newaxis to the head until the dimension gets max([a.ndim for a in arrays])

Apparently the latter function ffast is less computationally demanding than fslow. Hence, it would be useful to have such an option.

import numpy as np

def broadcast_without_repeating(*arrays):
    np.broadcast_shapes([a.shape for a in arrays])    
    max_dim = max(a.ndim for a in arrays)
    return tuple(array[(None,) * (max_dim - array.ndim), ...] for array in arrays)

Start by reviewing the existing broadcast_arrays() behavior and the proposed np.broadcast_shapes() check. Determine how an option such as no_repeat should preserve broadcastability while only prepending dimensions, then compare the result shapes with the r and theta examples. Done means the API decision, implementation, and coverage for compatible and incompatible shapes are settled.