home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / ledar34.zip / leda-r-3_4_tar / LEDA-3.4 / src / plane / _rat_segment.c < prev    next >
C/C++ Source or Header  |  1996-09-03  |  18KB  |  590 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA 3.4
  4. +
  5. +  _rat_segment.c
  6. +
  7. +  This file is part of the LEDA research version (LEDA-R) that can be 
  8. +  used free of charge in academic research and teaching. Any commercial
  9. +  use of this software requires a license which is distributed by the
  10. +  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
  11. +  (fax +49 681 31104).
  12. +
  13. +  Copyright (c) 1991-1996  by  Max-Planck-Institut fuer Informatik
  14. +  Im Stadtwald, 66123 Saarbruecken, Germany     
  15. +  All rights reserved.
  16. + 
  17. *******************************************************************************/
  18.  
  19. //------------------------------------------------------------------------------// rat_segment : segments with rational coordinates
  20. //
  21. // by S. Naeher (1995,1996)
  22. //------------------------------------------------------------------------------
  23.  
  24. #include <LEDA/rat_segment.h>
  25. #include <math.h>
  26. #include <ctype.h>
  27.  
  28.  
  29.  
  30. #if defined(sparc)
  31. #define FABS(x) (*(unsigned long*)&x) &= 0x7FFFFFFF
  32. #else
  33. #define FABS(x) x=fabs(x)
  34. #endif
  35.  
  36.  
  37. mutex rat_segment_rep::mutex_id_counter;
  38.  
  39. unsigned long rat_segment_rep::id_counter = 0;
  40.  
  41.  
  42. rat_segment_rep::rat_segment_rep(const rat_point& p, const rat_point& q) : start(p),end(q)
  43. { dx = q.X()*p.W() - p.X()*q.W();
  44.   dy = q.Y()*p.W() - p.Y()*q.W();
  45.   dxd = dx.todouble();
  46.   dyd = dy.todouble();
  47.   mutex_id_counter.lock();
  48.   id  = id_counter++; 
  49.   mutex_id_counter.unlock();
  50.  }
  51.  
  52. rat_segment_rep::rat_segment_rep() : start(0,0,1), end(0,0,1)
  53. { dx  = 0;
  54.   dy  = 0;
  55.   dxd = 0;
  56.   dyd = 0;
  57.   mutex_id_counter.lock();
  58.   id  = id_counter++; 
  59.   mutex_id_counter.unlock();
  60.  }
  61.  
  62. rat_segment::rat_segment(const rat_point& x, const rat_vector& v) 
  63. { PTR = new rat_segment_rep(x,x+v); }
  64.   
  65.  
  66. // MOVED TO _pglobal.c
  67. /*
  68.  * int rat_segment::use_filter = true;
  69.  */
  70.  
  71.  
  72. segment rat_segment::to_segment() const
  73. { return segment(ptr()->start.to_point(),ptr()->end.to_point()); }
  74.  
  75.  
  76. rat_segment rat_segment::reversal() const
  77. { return rat_segment(ptr()->end,ptr()->start); }
  78.  
  79.  
  80. rat_segment rat_segment::translate(const integer& dx, const integer& dy, const integer& dw) const
  81. { rat_point p = ptr()->start.translate(dx,dy,dw);
  82.   rat_point q = ptr()->end.translate(dx,dy,dw);
  83.   return rat_segment(p,q);
  84.  }
  85.  
  86. rat_segment rat_segment::translate(const rational& dx, const rational& dy) const
  87. { rat_point p = ptr()->start.translate(dx,dy);
  88.   rat_point q = ptr()->end.translate(dx,dy);
  89.   return rat_segment(p,q);
  90. }
  91.  
  92. rat_segment rat_segment::translate(const rat_vector& v) const
  93. { rat_point p = ptr()->start.translate(v);
  94.   rat_point q = ptr()->end.translate(v);
  95.   return rat_segment(p,q);
  96.  }
  97.  
  98.  
  99. rat_segment rat_segment::rotate90() const
  100. { return rat_segment(start(),end().rotate90(start())); }
  101.  
  102. rat_segment rat_segment::rotate90(const rat_point& q) const
  103. {  return rat_segment(start().rotate90(q),end().rotate90(q)); }
  104.  
  105.  
  106. rat_segment rat_segment::reflect(const rat_point& p, const rat_point& q) const
  107. { return rat_segment(start().reflect(p,q),end().reflect(p,q)); }
  108.  
  109. rat_segment rat_segment::reflect(const rat_point& p) const
  110. { return rat_segment(start().reflect(p),end().reflect(p)); }
  111.  
  112.  
  113.  
  114. rational rat_segment::sqr_length() const
  115. { integer x = dx();
  116.   integer y = dy();
  117.   integer w = W1()*W2();
  118.   return rational(x*x+y*y,w*w);
  119.  }
  120.  
  121.  
  122. ostream& operator<<(ostream& out, const rat_segment& s) 
  123. { out << "[" << s.start() << "===" << s.end() << "]"; 
  124.   return out;
  125.  } 
  126.  
  127. istream& operator>>(istream& in, rat_segment& s) 
  128. { // syntax: {[} p {===} q {]}
  129.  
  130.   rat_point p,q; 
  131.   char c;
  132.  
  133.   do in.get(c); while (isspace(c));
  134.   if (c != '[') in.putback(c);
  135.  
  136.   in >> p;
  137.  
  138.   do in.get(c); while (isspace(c));
  139.   while (c== '=') in.get(c);
  140.   while (isspace(c)) in.get(c);
  141.   in.putback(c);
  142.  
  143.   in >> q; 
  144.  
  145.   do in.get(c); while (c == ' ');
  146.   if (c != ']') in.putback(c);
  147.  
  148.   s = rat_segment(p,q); 
  149.   return in; 
  150.  
  151.  } 
  152.  
  153.  
  154. bool rat_segment::contains(const rat_point& p) const
  155. { rat_point a = source();
  156.   rat_point b = target();
  157.   if (orientation(a,b,p) == 0)
  158.   { if ( a == b) return ( a == p? true: false);
  159.     rat_point c = b.rotate90(a);
  160.     rat_point d = a.rotate90(b);
  161.     return (orientation(a,c,p) <= 0)  && (orientation(b,d,p) <= 0);
  162.    }
  163.   return false;
  164.  }
  165.  
  166.  
  167. bool rat_segment::intersection(const rat_segment& t)
  168. { // decides whether this and |t| intersect
  169.   if (   is_trivial() ) return t.contains(  source());
  170.   if ( t.is_trivial() ) return   contains(t.source());
  171.    
  172.   int o1 = orientation(*this,t.start()); 
  173.   int o2 = orientation(*this,t.end());
  174.   int o3 = orientation(t,start());
  175.   int o4 = orientation(t,end());
  176.  
  177.   if ( o1 == 0 && o2 == 0 )
  178.      // same slope
  179.      return ( t.contains(  source()) || t.contains(  target()) ||
  180.                 contains(t.source()) ||   contains(t.target())    );
  181.   else 
  182.      return o1 != o2 && o3 != o4;
  183. }
  184.  
  185.  
  186.  
  187. bool rat_segment::intersection(const rat_segment& t, rat_point& I) const
  188. { if (   is_trivial() )
  189.   { if ( t.contains(  source()) )
  190.        { I = source(); return true; }
  191.     else
  192.        return false;
  193.   }
  194.   if ( t.is_trivial() ) 
  195.   { if ( contains(t.source()) )
  196.        { I = t.source(); return true; }
  197.     else 
  198.        return false;
  199.   }
  200.    
  201.   int o1 = orientation(*this,t.start()); 
  202.   int o2 = orientation(*this,t.end());
  203.   int o3 = orientation(t,start());
  204.   int o4 = orientation(t,end());
  205.  
  206.   if ( o1 == 0 && o2 == 0 )
  207.   { // same slope
  208.     if ( t.contains(  source()) )
  209.        {  I = source(); return true; }
  210.     if ( t.contains(  target()) )
  211.        { I = target(); return true;  }
  212.     if ( contains(t.source()) )
  213.        {  I = t.source(); return true; }
  214.     if ( contains(t.target()) )
  215.        { I = t.target(); return true; }
  216.     return false;
  217.   }
  218.   
  219.   if ( o1 != o2 && o3 != o4 )
  220.   { integer w  =   dy()*t.dx() - t.dy()*dx();
  221.     integer c1 =   X2()*  Y1() -   X1()*  Y2();
  222.     integer c2 = t.X2()*t.Y1() - t.X1()*t.Y2();
  223.  
  224.     I = rat_point(c2*dx()-c1*t.dx(), c2*dy()-c1*t.dy(), w);
  225.     return true;
  226.   }
  227.   return false;
  228. }
  229.  
  230. bool rat_segment::intersection(const rat_segment& t, rat_segment& I) const
  231. { if (   is_trivial() )
  232.   { if ( t.contains(  source()) )
  233.        { I = rat_segment(source(),source()); return true; }
  234.     else
  235.        return false;
  236.   }
  237.   if ( t.is_trivial() ) 
  238.   { if ( contains(t.source()) )
  239.        { I = rat_segment(t.source(),t.source()); return true; }
  240.     else 
  241.        return false;
  242.   }
  243.    
  244.   int o1 = orientation(*this,t.start()); 
  245.   int o2 = orientation(*this,t.end());
  246.   int o3 = orientation(t,start());
  247.   int o4 = orientation(t,end());
  248.  
  249.   if ( o1 == 0 && o2 == 0 )
  250.   { rat_point sa = source(); 
  251.     rat_point sb = target();
  252.     if ( compare (sa,sb) > 0 )
  253.        { rat_point h = sa; sa = sb; sb = sa; }
  254.  
  255.     rat_point ta = t.source(); rat_point tb = t.target();
  256.  
  257.     if ( compare (ta,tb) > 0 )
  258.        { rat_point h = ta; ta = tb; tb = ta; }
  259.  
  260.     rat_point a = ta;
  261.     rat_point b = tb; 
  262.  
  263.     if (compare(sa,a) < 0) a = sa;
  264.     if (compare(sb,b) < 0) b = sb;
  265.  
  266.     if ( compare(a,b) <= 0 )
  267.     { I = rat_segment(a,b);
  268.       return true;
  269.     }
  270.     return false;
  271.   }
  272.   
  273.   if ( o1 != o2 && o3 != o4 )
  274.   { integer w  =   dy()*t.dx() - t.dy()*dx();
  275.     integer c1 =   X2()*  Y1() -   X1()*  Y2();
  276.     integer c2 = t.X2()*t.Y1() - t.X1()*t.Y2();
  277.  
  278.     rat_point p(c2*dx()-c1*t.dx(), c2*dy()-c1*t.dy(), w);
  279.     I = rat_segment(p,p);
  280.     return true;
  281.   }
  282.   return false;
  283. }
  284.  
  285. bool rat_segment::intersection_of_lines(const rat_segment& t, rat_point& inter) const
  286. { 
  287.   /* decides whether the lines induced by |s| and |this| segment 
  288.      intersect and, if so, returns the intersection in |inter|. 
  289.      It is assumed that both segments have non-zero length
  290.    */
  291.  
  292.   integer w = dy()*t.dx() - dx()*t.dy();
  293.  
  294.   if (w == 0) //same slope
  295.   { inter = source();
  296.     return true;
  297.   }
  298.  
  299.   integer c1 =   X2()*  Y1() -   X1()*  Y2();
  300.   integer c2 = t.X2()*t.Y1() - t.X1()*t.Y2();
  301.  
  302.   inter = rat_point(c2*dx()-c1*t.dx(), c2*dy()-c1*t.dy(), w);
  303.  
  304.   return true;
  305. }
  306.  
  307.  
  308. const double eps0 = ldexp(1.0,-53);
  309. const double eps2 = ldexp(1.0,-47); // 64 * eps0
  310.  
  311. int cmp_slopes(const rat_segment& s1, const rat_segment& s2)
  312. { 
  313.   if (rat_segment::use_filter)
  314.   { if ( s2.dxD() == 0 || s1.dxD() == 0 )
  315.     { // one of the lines is vertical
  316.       if ( s2.dxD() == 0 && s1.dxD() == 0 ) return 0;
  317.       if ( s2.dxD() == 0) return -1;  // s1 has smaller slope
  318.       return +1;  // s2 has smaller slope
  319.     }
  320.  
  321.     double dy1dx2 = s1.dyD()*s2.dxD();
  322.     double dy2dx1 = s2.dyD()*s1.dxD();
  323.     double E = dy1dx2 - dy2dx1; 
  324.     if ( s1.dxD() * s2.dxD() < 0 ) E = -E;
  325.  
  326.     //-----------------------------------------------------------------------
  327.     //  ERROR BOUNDS
  328.     //-----------------------------------------------------------------------
  329.     // mes(E) = 2 * (mes(dy1*dx2) + mes(dy2*dx2))
  330.     //        = 2 * (mes(dy1)*mes(dx2) + mes(dy2)*mes(dx2))
  331.     //        = 2 * (4*fabs(dy1dx2) + 4*fabs(dy2dx2))
  332.     //        = 8 * (fabs(dy1dx2) + fabs(dy2dx2))
  333.     //
  334.     // ind(E) = (ind(dy1*dx2) + ind(dy2*dx2) + 1)/2
  335.     //        = (ind(dy1) + ind(dx2) + 0.5 + ind(dy2) + ind(dx2) + 0.5 + 1)/2
  336.     //        = (0.5 + 0.5 + 0.5 + 0.5 + 0.5 + 0.5 + 1)/2
  337.     //        = 2
  338.     //
  339.     // eps(E) = ind(E) * mes(E) * eps0
  340.     //        = 16 * (fabs(dy1dx2) + fabs(dy2dx2)) * eps0
  341.     //-----------------------------------------------------------------------
  342.   
  343.     FABS(dy1dx2);
  344.     FABS(dy2dx1);
  345.     double eps = 16 * (dy1dx2+dy2dx1) * eps0;
  346.   
  347.     if (E >  eps) return  1;
  348.     if (E < -eps) return -1;
  349.     if (eps < 1)  return 0;
  350.   }
  351.   
  352.   // use big integers
  353.   if ( s2.dx() == 0 || s1.dx() == 0 )
  354.   { // one of the lines is vertical
  355.     if ( s2.dx() == 0 && s1.dx() == 0 ) return 0;
  356.     if ( s2.dx() == 0) return -1;  // s1 has smaller slope
  357.     return +1;  // s2 has smaller slope
  358.   }
  359.   integer E = s1.dy() * s2.dx() - s2.dy() * s1.dx();
  360.   int s = sign(E);
  361.   return ( sign(s1.dx()) * sign(s2.dx()) < 0 ) ? -s : s;
  362. }
  363.  
  364. int rat_segment::cmp_slope(const rat_segment& s1) const
  365. { return cmp_slopes(*this,s1); }
  366.  
  367.  
  368. int orientation(const rat_segment& s, const rat_point& p)
  369. { 
  370.   rat_point a = s.start();
  371.  
  372.   if (rat_segment::use_filter)
  373.   { 
  374.     double dx = s.dxD();
  375.     double dy = s.dyD();
  376.     double axpw = a.XD() * p.WD();
  377.     double aypw = a.YD() * p.WD();
  378.     double pxaw = p.XD() * a.WD();
  379.     double pyaw = p.YD() * a.WD();
  380.   
  381.     double E =  dy * (axpw - pxaw) - dx * (aypw - pyaw); 
  382.  
  383.     //-----------------------------------------------------------------------
  384.     //  ERROR BOUNDS
  385.     //-----------------------------------------------------------------------
  386.     //
  387.     // mes(E) = 2 * (mes(dx) * 2*(mes(aypw)+mes(pyaw)) +
  388.     //               mes(dy) * 2*(mes(axpw)+mes(pxaw)))
  389.     //        = 4 * (mes(dx) * (mes(aypw)+mes(pyaw)) +
  390.     //               mes(dy) * (mes(axpw)+mes(pxaw)))
  391.     //        = 8 * (fabs(dx) * (2*fabs(aypw)+2*fabs(pyaw)) +
  392.     //               fabs(dy) * (2*fabs(axpw)+2*fabs(pxaw)))
  393.     //        =16 * (fabs(dx) * (fabs(aypw)+fabs(pyaw)) +
  394.     //               fabs(dy) * (fabs(axpw)+fabs(pxaw)))
  395.     //
  396.     // ind(E) =  (ind(dx) + ind(aypw - pyaw) + 0.5 +
  397.     //            ind(dy) + ind(axpw - pxaw) + 0.5 + 1)/2
  398.     //        =  (ind(dx) + (ind(aypw) + ind(pyaw) + 1)/2 + 0.5 +
  399.     //            ind(dy) + (ind(axpw) + ind(pxaw) + 1)/2 + 0.5 + 1)/2
  400.     //        =  (0.5 + (1.5 + 1.5 + 1)/2 + 0.5
  401.     //            0.5 + (1.5 + 1.5 + 1)/2 + 0.5 + 1)/2
  402.     //        =  (0.5 + 2 + 0.5 + 0.5 + 2 + 0.5 + 1)/2
  403.     //        =  3.5 
  404.     //
  405.     // mes(E) = ind(E) * mes(E) * eps0
  406.     //        = 56 * (fabs(dx) * (fabs(aypw)+fabs(pyaw)) +
  407.     //                fabs(dy) * (fabs(axpw)+fabs(pxaw))) * eps0
  408.     //-----------------------------------------------------------------------
  409.   
  410.     FABS(aypw);
  411.     FABS(pyaw);
  412.     FABS(axpw);
  413.     FABS(pxaw);
  414.     FABS(dx);
  415.     FABS(dy);
  416.   
  417.     double eps = 56 * ((aypw+pyaw)*dx + (axpw+pxaw)*dy) * eps0;
  418.   
  419.     if (E >  eps) return  1;
  420.     if (E < -eps) return -1;
  421.     if (eps < 1)  return 0;
  422.    }
  423.  
  424.   // big integers
  425.  
  426.   return sign( s.dy() * (a.X()*p.W() - p.X()*a.W()) -
  427.                s.dx() * (a.Y()*p.W() - p.Y()*a.W()) );
  428.  
  429.  }
  430.  
  431.  
  432.  
  433. int cmp_segments_at_xcoord(const rat_segment& s1, const rat_segment& s2, 
  434.                            const rat_point& r)
  435. { 
  436.   if (rat_segment::use_filter)
  437.     {
  438.       double s1x = s1.XD1();
  439.       double s1y = s1.YD1();
  440.       double s1w = s1.WD1();
  441.       
  442.       double s2x = s2.XD1();
  443.       double s2y = s2.YD1();
  444.       double s2w = s2.WD1();
  445.       
  446.       double rx = r.XD();
  447.       double ry = r.YD();
  448.       double rw = r.WD();
  449.       
  450.       double d1x = s1.dxD();
  451.       double d1y = s1.dyD();
  452.       
  453.       double d2x = s2.dxD();
  454.       double d2y = s2.dyD();
  455.       
  456.       double E = d2x * s2w * (s1y * d1x * rw + d1y * (rx * s1w - s1x * rw))
  457.           - d1x * s1w * (s2y * d2x * rw + d2y * (rx * s2w - s2x * rw));
  458.  
  459.  
  460.       /* ----------------------------------------------------------------------
  461.        *                             ERROR BOUNDS
  462.        * ----------------------------------------------------------------------
  463.        *
  464.        * mes(E) = mes(d2x*s2w *(s1y*d1x*rw+d1y*(rx*s1w-s1x*rw))
  465.        *              -d1x*s1w*(s2y*d2x*rw+d2y*(rx*s2w-s2x*rw)))
  466.        *
  467.        *        = 2*(mes(d2x)*mes(s2w)*mes(s1y*d1x*rw+d1y*(rx*s1w-s1x*rw))
  468.        *             + mes(d1x)*mes(s1w)*mes(s2y*d2x*rw+d2y*(rx*s2w-s2x*rw)))
  469.        *
  470.        *        = 4*(mes(d2x)*mes(s2w)*
  471.        *               (mes(s1y)*mes(d1x)*mes(rw)+
  472.        *                  mes(d1y)*2(mes(rx)*mes(s1w)+mes(s1x)*mes(rw)))
  473.        *           +mes(d1x)*mes(s1w)*
  474.        *              (mes(s2y)*mes(d2x)*mes(rw)+
  475.        *                  mes(d2y)*2(mes(rx)*mes(s2w)+mes(s2x)*mes(rw))))
  476.        *
  477.        *        = 4*(fabs(d2x)*fabs(s2w)*
  478.        *               (fabs(s1y)*2*fabs(d1x)*fabs(rw)+
  479.        *                  fabs(d1y)*2*(fabs(rx)*fabs(s1w)+fabs(s1x)*fabs(rw)))
  480.        *            +fabs(d1x)*fabs(s1w)*
  481.        *               (fabs(s2y)*2*fabs(d2x)*fabs(rw)+
  482.        *                  fabs(d2y)*2*(fabs(rx)*fabs(s1w)+fabs(s2x)*fabs(rw))))
  483.        *
  484.        * ----------------------------------------------------------------------
  485.        *
  486.        * ind(E) = ind(d2x*s2w *(s1y*d1x*rw+d1y*(rx*s1w-s1x*rw))
  487.        *              -d1x*s1w*(s2y*d2x*rw+d2y*(rx*s2w-s2x*rw)))
  488.        * 
  489.        *        = (ind(d2x*s2w *(s1y*d1x*rw+d1y*(rx*s1w-s1x*rw)))
  490.        *          +ind(d1x*s1w*(s2y*d2x*rw+d2y*(rx*s2w-s2x*rw))) + 1)/2
  491.        *
  492.        *        = (ind(d2x)+ind(s2w *(s1y*d1x*rw+d1y*(rx*s1w-s1x*rw)))+0.5
  493.        *          +ind(d1x)+ind(s1w*(s2y*d2x*rw+d2y*(rx*s2w-s2x*rw)))+0.5+1)/2
  494.        *
  495.        *        = (1 + ind(s1y*d1x*rw+d1y*(rx*s1w-s1x*rw)) + 1
  496.        *           + 1 ind(s2y*d2x*rw+d2y*(rx*s2w-s2x*rw)) + 1 + 1)/2
  497.        *
  498.        *        = (5 + (ind(s1y*d1x*rw) + ind(d1y*(rx*s1w-s1x*rw)) + 1)/2
  499.        *          + (ind(s2y*d2x*rw) + ind(d2y*(rx*s2w-s2x*rw)) + 1)/2 )/2
  500.        *
  501.        *        = (5 + (ind(s1y)+ind(d1x)+ind(rw) + 1 + 
  502.        *                ind(d1y*(rx*s1w-s1x*rw)) + 1)/2
  503.        *          + (ind(s2y)+ind(d2x)+ind(rw) + 1 +
  504.        *             ind(d2y*(rx*s2w-s2x*rw))+ 1)/2)/2
  505.        *
  506.        *        = (5 + (0.5 + 0.5 + 0.5 + 1 + ind(d1y*(rx*s1w-s1x*rw)) + 1)/2
  507.        *          + (0.5 + 0.5 + 0.5 + 1 + ind(d2y*(rx*s2w-s2x*rw)) + 1)/2 )/2
  508.        *
  509.        *        = (5 + (3.5 + ind(d1y) + ind(rx*s1w-s1x*rw) + 0.5)/2
  510.        *          + (3.5 + ind(d2y) + ind(rx*s2w-s2x*rw) + 0.5)/2 )/2
  511.        *
  512.        *        = (5 + (3.5 + 0.5 +(ind(rx*s1w) + ind(s1x*rw) + 1)/2 + 0.5)/2
  513.        *          + (3.5 + 0.5 +(ind(rx*s2w) + ind(s2x*rw) + 1)/2 + 0.5)/2 )/2 
  514.        *
  515.        *        = (5 + (3.5 + 0.5 + (1.5 + 1.5 + 1 )/2 + 0.5 )/2
  516.        *             + (3.5 + 0.5 + (1.5 + 1.5 + 1 )/2 + 0.5 )/2)/2
  517.        *
  518.        *        = (5 + (4 + 2 + 0.5)/2 + (4 + 2 + 0.5)/2)/2 
  519.        *
  520.        *        = (5 + 13/2)/2 = 23/4
  521.        *
  522.        * ----------------------------------------------------------------------
  523.        *
  524.        * eps(E) = ind(E) * mes(E) * eps0
  525.        *
  526.        * eps    = 23*(fabs(d2x)*fabs(s2w)*
  527.        *               (fabs(s1y)*fabs(d1x)*fabs(rw)+
  528.        *                  fabs(d1y)*2*(fabs(rx)*fabs(s1w)+fabs(s1x)*fabs(rw)))
  529.        *            +fabs(d1x)*fabs(s1w)*
  530.        *               (fabs(s2y)*fabs(d2x)*fabs(rw)+
  531.        *                  fabs(d2y)*2*(fabs(rx)*fabs(s1w)+fabs(s2x)*fabs(rw))))
  532.        *            *eps0
  533.        *
  534.        *
  535.        * ------------------------------------------------------------------- */
  536.  
  537.       FABS(s1x);
  538.       FABS(s1y);
  539.       FABS(s1w);
  540.       
  541.       FABS(s2x);
  542.       FABS(s2y);
  543.       FABS(s2w);
  544.       
  545.       FABS(rx);
  546.       FABS(ry);
  547.       FABS(rw);
  548.       
  549.       FABS(d1x);
  550.       FABS(d1y);
  551.       
  552.       FABS(d2x);
  553.       FABS(d2y);
  554.       
  555.       double eps = 23 * 
  556.                    (d2x * s2w * (s1y * d1x * rw + 2*d1y * 
  557.                                  (rx * s1w + s1x * rw))
  558.                    + d1x * s1w * (s2y * d2x * rw + 2*d2y * 
  559.                                   (rx * s2w + s2x * rw)))
  560.                    * eps0; 
  561.  
  562.       if(E > eps)  return  1;
  563.       if(E < -eps) return -1;
  564.       if (eps < 1) return 0;
  565.  
  566.     }
  567.   
  568.   return sign( s2.dx() * s2.W1() *
  569.                  ( s1.Y1() * s1.dx() * r.W() + 
  570.                    s1.dy() * (r.X() * s1.W1() - s1.X1() * r.W()))
  571.                -
  572.                s1.dx() * s1.W1() * 
  573.                  ( s2.Y1() * s2.dx() * r.W() + 
  574.                    s2.dy() * ( r.X() * s2.W1() - s2.X1() * r.W()))); 
  575. }
  576.  
  577.  
  578. rat_segment rat_segment::perpendicular(const rat_point& q) const
  579. { if ( is_trivial() )
  580.      error_handler(1,"rat_segment::perpendicular: this must be nontrivial");
  581.   rat_point p = source();
  582.   integer  dx = q.X()*p.W() - p.X()*q.W();
  583.   integer  dy = q.Y()*p.W() - p.Y()*q.W();
  584.   integer  dw = p.W()*q.W();
  585.   rat_point r;
  586.   intersection_of_lines(translate(dx,dy,dw).rotate90(), r);
  587.   return rat_segment(q,r);
  588.  }
  589.  
  590.