#! /usr/bin/gawk -f
# Last edited on 2015-02-24 18:05:58 by stolfilocal

# ----------------------------------------------------------------------
    

            # Now assume {pav_new} is correct and adjust {vbt,vcr} so that {pav_new*vbt - vcr = 0}.
            # Convert {vbt,vcr} to scale where the errors have same deviation:
            U_avg = vbt/ulp_vbt;
            V_avg = vcr/ulp_vcr;
            
            # The equation {pav_new*vbt - vcr = 0} becomes {cU*U + cV*V = 0}, where:
            cU = pav_new*ulp_vbt;
            cV = -ulp_vcr;
            
            # Project {U_avg,V_avg} onto the line {cU*U + cV*V = 0}:
            dd = (cU*U_avg + cV*V_avg)/(cU*cU + cV*cV);
            U_sol = U_avg - dd*cU;
            V_sol = V_avg - dd*cV;
            
            # Convert back to original scale:
            vbt_new = U_sol*ulp_vbt;
            vcr_new = V_sol*ulp_vcr;

# ----------------------------------------------------------------------

function check_prices(pop,phi,plo,pcl,vbt,vcr,pav,pmd,  ploc,phic,tol_pav)
  {
    # Checks consistency of prices: {pop,phi,plo,pcl}, {pav} (nominal mean), and {pmd} (computed mean).
    # Uses global parameters {ulp_vbt,ulp_vcr,ulp_pav,ulp_phl}.
    
    if ((vbt == 0) != (vcr == 0))
      { data_error(("inconsistent vbt = \"" vbt "\" vcr = \"" vcr "\"")); }
    
    if ((vbt < min_vbt) || (vcr < min_vcr))
      { 
        ...
        
        # Checks whether the computed average price is within {plo,phi}, apart from rounding errors:
        # Compute range {[ploc _ phic]} of true average price, assuming worst rounding: 
        ploc = (vcr - 0.50001*ulp_vcr)/(vbt + 0.50001*ulp_vbt);
        phic = (vcr + 0.50001*ulp_vcr)/(vbt - 0.50001*ulp_vbt);

        if ((phic < plo - 0.500001*ulp_phl) || (ploc > phi + 0.500001*ulp_phl)) 
          { 
            printf "%s:%s: !! price range", FILENAME, FNR > "/dev/stderr";
            printf " is [ %.5f _ %.5f ]", plo, phi > "/dev/stderr";
            printf " should intersect [ %.5f _ %.5f ]\n", ploc, phic > "/dev/stderr";
            data_error(("excessive rounding error"));
          }

        # Check whether computed average price differs too much from given one:
        tol_pav = (phic - ploc);
        if (tol_pav < ulp_pav) { tol_pav = ulp_pav; }
        if ((pmd - pav) > 0.500001*tol_pav)
          { printf "%s:%s: !! average price inconsistent", FILENAME, FNR > "/dev/stderr";
            printf " is %.5f should be %.5f\n", pav, pmd > "/dev/stderr";
            # Check whether the difference can be ascribed to bad rounding:
            if ((pav > phic + 0.500001*tol_pav) || (pav < ploc - 0.500001*tol_pav))
              { data_error(("excessive rounding error")); }
          }
          
        ...
      }
   }