home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch2_7 / len4.c next >
Text File  |  1995-03-06  |  2KB  |  39 lines
  1. /*
  2.  * len4.c - find length of the 4-vector v having components v = [a b c d]
  3.  * programmed by: Alan Paeth, 17-Oct-1994
  4.  *
  5.  * This Euclidean distance estimator
  6.  *
  7.  *   1) can only overestimate, making it useful for containment heuristics,
  8.  *   2) is exact for vectors [1 0 0 0] (len4 = 1) and [1 1 1 1] (len4 = 2) 
  9.  *      under all permutations, sign changes and (integer) scaling (*),
  10.  *   3) is nearly exact for vectors [1 1 0 0] and [1 1 1 0] (len4=irrational)
  11.  *      subject to the same conditions as in (2) above,
  12.  *   4) offers best 3D fit using the nested form len3(a,b,c) = len4(a,b,c,0)
  13.  *   5) may be easily reprogrammed for floating point operation by rewriting
  14.  *      all occurences of the string "int" with "double" or "float" and by
  15.  *      removing the "a++" rounding fudge immediately before the return.
  16.  *  (*) - exact estimates hold only in the presence of the mod: one line up.
  17.  *
  18.  * Note: worst-case error occurs for vectors of the form [60c 25c 19c 16c]
  19.  *       and is worsened (because of increased a++ rounding imprecision)
  20.  *       with vector [2 1 1 1]. Otherwise (c large), the maximum error in
  21.  *       overestimation is always bounded by 16% [1.1597+ = sqrt(538)/20].
  22.  */
  23.   
  24. #define absv(x)    if (x < 0) x = -x
  25. #define inorder(x,y)    {int t; if ((t = a - b) < 0) {a -= t; b += t; } }
  26.  
  27. int len4(a, b, c, d)
  28.     {
  29.     absv(a); absv(b);            /* get the absolute values */
  30.     absv(c); absv(d);            /* (component magnitudes) */
  31.     inorder(a, b); inorder(c, d);    /* everyone has a chance to play */
  32.     inorder(a, c); inorder(b, d);    /* (a,d) are big (winner, loser) */
  33.     inorder(b, c);            /* playoff for 2nd and 3rd slots */
  34.     a += (25*b + 19*c + 16*d)/60;    /* compute 4D approximate length */
  35. /*  a +=  (5*b +  4*c +  3*d)/12;    .. only .1% worse; easy to eval  */
  36.     a++;                /* Roundoff -> underestimation   */
  37.     return(a);                /* omit the above one bit jitter */
  38.     }
  39.