The current algorithm for the power function of p-box objects are not correct for cases that straddle 0.
Below are two simple examples to illustrate the problem. For the first Gaussian p-box, it seems an additional envelope to 0 will solve the issue. But this solution seems inadequate for the second example of a Uniform case.
The function verify_power_calculations_plot(a, 2) means that we are calculating $a^2$. reference suggests a discretisation-interval propagation-stacking interval solution; while direct suggests the current algorithm which naively takes takes the power of the edges.
