See this and the following discussion, we should likely fix the double-double value for np.finfo(np.longdouble).smallest_normal. @wkschwartz, lets continue discussion here.

I might be missing some context here, so forgive me if I'm confusing matters. It appears that you are able to compute the smallest positive subnormal 𝜎, and a few lines below you define macheps, so I take it you are able to compute the machine epsilon 𝜀. From there, you can in fact compute the smallest positive normal 𝜔 as 𝜎 ⨸ 𝜀, where ⨸ means division in the relevant floating point system (which is what the / operator does in C and Python). In Python terms, you can do approximately the following, but you'll need to move up the definition of macheps from below.

    # Leave the same value for the smallest subnormal as double
    smallest_subnormal_dd = ld(nextafter(0., 1.))
    smallest_normal_dd = smallest_subnormal / macheps

Originally posted by @wkschwartz in #18536 (comment)