code 41558;
real procedure BIVANORM(H, K, RHO); value H, K, RHO;
    real H, K, RHO;
    begin real B;
    
      real procedure V(H, K, EPS); value H, K, EPS;
          real H, K, EPS;
          if H = 0 ∨ K = 0 then V:= 0 else
          if ABS(H) < ABS(K) then
          V:= (PHI(H) - .5) × (PHI(K) - .5) - V(K, H, EPS)
          else
          if ABS(K)> 8 then
              V:=.15915 49430 9189 × ARCTAN (K/H)
          else
      begin real M, L, L2, S, R, T, SS, TSN; integer N;
          L:= K / H; M:= H × H / 2; L2:= L × L; R:= EXP(-M);
          S:= 1 - R; T:= L; SS:= T × S;
          for N:= 1, N + 1 while ABS(TSN) ≥ EPS do
          begin R:= R × M / N; S:= S - R; T:= -T × L2;
              TSN:=S × T / (2 × N + 1);
              SS:= SS + TSN
          end;
          V:= SS × .15915 49430 9189
      end V;
    
    if H < -8 ∨ K < -8 then B:=0 else
    if H > 8 ∧ K > 8 then B:=1 else
    B:= if ABS(RHO) > 1 then
            STATAL3 ERROR(“BIVANORM”, 3, RHO)
        else if ABS(RHO) = 1 then
           (if RHO = 1 then (if K ≤ H then PHI(K)
                                              else PHI(H))
                          else (if H ≤ -K then 0
                                 else PHI(K) - PHI(H)))
        else V(H,(K - RHO × H)/ SQRT(1 - RHO ⭡ 2), 10-14)
             + V(K,(H - RHO × K)/ SQRT(1 - RHO ⭡ 2), 10-14)
             + .5 × (PHI(H) + PHI(K))
             - .15915 49430 9189 × ARCCOS(RHO);
    if B < 0 then BIVANORM:=0 else BIVANORM:=B;
end BIVANORM;
eop