np.linalg.norm overflows at 1e154 and underflows at 1e-200

NumPy 2.5.4 evidently computes vector 2-norms as sqrt(sum of squares) without scaling. np.linalg.norm([v, v, v]) was inf for v = 1e154 and 1e200, 0.0 for v = 1e-200, and 5.6e-6 relative off at 1e-160 (the squares are subnormal). On float32 it overflows already at v = 1e20. scipy.linalg.norm(x) dispatches to BLAS nrm2 for a 1-D float array with ord None/2 and no axis (scipy 1.18.1 source), and matched v*sqrt(3) to the printed 6 digits in all these cases. With axis= it returned inf like NumPy. math.hypot(*x) is safe for single vectors. For rows of a matrix, scale first: m = abs(X).max(axis=1, keepdims=True); m[m == 0] = 1; m[:, 0] * np.linalg.norm(X / m, axis=1). That gave <= 4.4e-16 relative error at scales 1e200, 1e-200 and 1 (100k x 50, vs math.hypot) at 2.7x the time of the plain norm; np.hypot.reduce(X, axis=1) was also accurate (<= 1e-15) but 8x slower.