Repository: scipy/scipy
File: scipy/stats/_stats_pythran.py (Python)

Code before:
    elif dh == -math_inf:
        if dk == -math_inf:
            p = 1.
        else:
            p = phid(-dk)
    # elseif dk == -inf, p = phid(-dh);
    elif dk == -math_inf:
        p = phid(-dh)
    # elseif r == 0, p = phid(-dh)*phid(-dk);
    elif r == 0:
        p = phid(-dh) * phid(-dk)
    # else, tp = 2*pi; h = dh; k = dk; hk = h*k; bvn = 0;
    else:
        tp = 2*math_pi
        h = dh
        k = dk
        hk = h*k
        bvn = 0.
        # if abs(r) < 0.3      % Gauss Legendre points and weights, n =  6
        #     w(1:3) = [0.1713244923791705 0.3607615730481384 0.4679139345726904];
        #     x(1:3) = [0.9324695142031522 0.6612093864662647 0.2386191860831970];
        if abs(r) < 0.3:
            w = [0.1713244923791705, 0.3607615730481384, 0.4679139345726904]
            x = [0.9324695142031522, 0.6612093864662647, 0.2386191860831970]
        # elseif abs(r) < 0.75 % Gauss Legendre points and weights, n = 12
        #     w(1:3) = [.04717533638651177 0.1069393259953183 0.1600783285433464];
        #     w(4:6) = [0.2031674267230659 0.2334925365383547 0.2491470458134029];
        #     x(1:3) = [0.9815606342467191 0.9041172563704750 0.7699026741943050];
        #     x(4:6) = [0.5873179542866171 0.3678314989981802 0.1252334085114692];
        elif abs(r) < 0.75:
            w = [.04717533638651177, 0.1069393259953183, 0.1600783285433464,
                 0.2031674267230659, 0.2334925365383547, 0.2491470458134029]
            x = [0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
                 0.5873179542866171, 0.3678314989981802, 0.1252334085114692]
        # else,                % Gauss Legendre points and weights, n = 20
        #     w(1:3) = [.01761400713915212 .04060142980038694 .06267204833410906];
        #     w(4:6) = [.08327674157670475 0.1019301198172404 0.1181945319615184];
        #     w(7:9) = [0.1316886384491766 0.1420961093183821 0.1491729864726037];
        #     w(10) =   0.1527533871307259;
        #     x(1:3) = [0.9931285991850949 0.9639719272779138 0.9122344282513259];
        #     x(4:6) = [0.8391169718222188 0.7463319064601508 0.6360536807265150];
        #     x(7:9) = [0.5108670019508271 0.3737060887154196 0.2277858511416451];
        #     x(10) =   0.07652652113349733;
        else:
            w = [.01761400713915212, .04060142980038694, .06267204833410906,
                 .08327674157670475, 0.1019301198172404, 0.1181945319615184,
                 0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
                 0.1527533871307259]
            x = [0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
                 0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
                 0.5108670019508271, 0.3737060887154196, 0.2277858511416451,

Review comment (suggestion):
I still think it'd make sense to use `roots_legenendre`,  https://github.com/scipy/scipy/issues/23469#issuecomment-3353487433

I believe the whole extent of the change is 
https://github.com/scipy/scipy/compare/main...ev-br:scipy:mvn_2D?expand=1#diff-b5afc2abf94ff5b2f535a0f4d6a2f09b5838effa4f4686913c61637b624abde0R617-R624

This will also remove `np.concatenate` and its related build warning.

Code after:
        p = 0.
    # elseif dh == -inf, if dk == -inf, p = 1; else p = phid(-dk); end
    elif dh == -math_inf:
        if dk == -math_inf:
            p = 1.
        else:
            p = phid(-dk)
    # elseif dk == -inf, p = phid(-dh);
    elif dk == -math_inf:
        p = phid(-dh)
    # elseif r == 0, p = phid(-dh)*phid(-dk);
    elif r == 0:
        p = phid(-dh) * phid(-dk)
    # else, tp = 2*pi; h = dh; k = dk; hk = h*k; bvn = 0;
    else:
        tp = 2*math_pi
        h = dh
        k = dk
        hk = h*k
        bvn = 0.
        # if abs(r) < 0.3      % Gauss Legendre points and weights, n =  6
        #     w(1:3) = [0.1713244923791705 0.3607615730481384 0.4679139345726904];
        #     x(1:3) = [0.9324695142031522 0.6612093864662647 0.2386191860831970];
        if abs(r) < 0.3:
            w = [0.1713244923791705, 0.3607615730481384, 0.4679139345726904]
            x = [0.9324695142031522, 0.6612093864662647, 0.2386191860831970]
        # elseif abs(r) < 0.75 % Gauss Legendre points and weights, n = 12
        #     w(1:3) = [.04717533638651177 0.1069393259953183 0.1600783285433464];
        #     w(4:6) = [0.2031674267230659 0.2334925365383547 0.2491470458134029];
        #     x(1:3) = [0.9815606342467191 0.9041172563704750 0.7699026741943050];
        #     x(4:6) = [0.5873179542866171 0.3678314989981802 0.1252334085114692];
        elif abs(r) < 0.75:
            w = [.04717533638651177, 0.1069393259953183, 0.1600783285433464,
                 0.2031674267230659, 0.2334925365383547, 0.2491470458134029]
            x = [0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
                 0.5873179542866171, 0.3678314989981802, 0.1252334085114692]
        # else,                % Gauss Legendre points and weights, n = 20
        #     w(1:3) = [.01761400713915212 .04060142980038694 .06267204833410906];
        #     w(4:6) = [.08327674157670475 0.1019301198172404 0.1181945319615184];
        #     w(7:9) = [0.1316886384491766 0.1420961093183821 0.1491729864726037];
        #     w(10) =   0.1527533871307259;
        #     x(1:3) = [0.9931285991850949 0.9639719272779138 0.9122344282513259];
        #     x(4:6) = [0.8391169718222188 0.7463319064601508 0.6360536807265150];
        #     x(7:9) = [0.5108670019508271 0.3737060887154196 0.2277858511416451];
        #     x(10) =   0.07652652113349733;
        else:
            w = [.01761400713915212, .04060142980038694, .06267204833410906,
                 .08327674157670475, 0.1019301198172404, 0.1181945319615184,
                 0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
                 0.1527533871307259]
            x = [0.9931285991850949, 0.9639719272779138, 0.9122344282513259,