Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

guar Namespace Reference


Compounds

struct  guar::Interval
class  guar::rational
struct  guar::rational::rational_rep
class  guar::ULong

Functions

void IntAdd (const unsigned int na, const unsigned long *a, const unsigned int nb, const unsigned long *b, unsigned int *nc, unsigned long *c)
int IntSub (const unsigned int na, const unsigned long *a, const unsigned int nb, const unsigned long *b, unsigned int *nc, unsigned long *c)
unsigned long IntMultShort (const unsigned int na, const unsigned long *a, const unsigned long b, unsigned long *c)
void IntMult (const unsigned int na, const unsigned long *a, const unsigned int nb, const unsigned long *b, unsigned int *nc, unsigned long *c)
unsigned long IntDivShort (const unsigned int na, const unsigned long *a, const unsigned long b, unsigned int *nq, unsigned long *q)
int IntDiv (const unsigned int na, const unsigned long *a, const unsigned int nb, const unsigned long *b, unsigned int *nq, unsigned long *q, unsigned int *nr, unsigned long *r)
int DoubleToInt (const double d, const unsigned int na, unsigned long *a, int *retLen)
Interval operator * (const Interval &a, const Interval &b)
GULAPI rational operator * (const rational &a, const rational &b)
GULAPI rational operator+ (const rational &A, const rational &B)
GULAPI rational operator/ (const rational &A, const rational &B)
GULAPI rational operator- (const rational &A, const rational &B)
template<class T> void IntTo (const unsigned int na, const unsigned long *a, T &t)
int IntSub (const int unsigned na, const unsigned long *a, const int unsigned nb, const unsigned long *b, unsigned int *nc, unsigned long *c)
int DoubleToInt (const double d, const unsigned int na, unsigned long *a, unsigned int *retLen)
double IntToDouble (const unsigned int na, const unsigned long *a)
Interval operator+ (const Interval &a, const Interval &b)
Interval operator- (const Interval &a, const Interval &b)
Interval operator/ (const Interval &a, const Interval &b)
Interval Sqrt (const Interval &a)
int Test (const Interval &a)
int Test (const Interval &a, const Interval &b)
double GetRatio (const Interval &a)
bool RegularIntersectLines (const point< rational > &A, const point< rational > &B, const point< rational > &a, const point< rational > &b, rational *lambda, rational *mu)
template<class T> GULAPI void BarycentricCoordinates (const point< T > &P1, const point< T > &P2, const point< T > &P3, const point< T > &Pin, T *w, T *u, T *v)
template GULAPI void BarycentricCoordinates (const point< rational > &P1, const point< rational > &P2, const point< rational > &P3, const point< rational > &Pin, rational *w, rational *u, rational *v)
template<class T> GULAPI void BarycentricCoordinates (const point< T > &P1, const point< T > &P2, const point< T > &P3, const point< T > &P12, const point< T > &P13, const point< T > &P23, const point< T > &N, const point< T > &Pin, T *w, T *u, T *v)
template GULAPI void BarycentricCoordinates (const point< rational > &P1, const point< rational > &P2, const point< rational > &P3, const point< rational > &P12, const point< rational > &P13, const point< rational > &P23, const point< rational > &N, const point< rational > &Pin, rational *w, rational *u, rational *v)
GULAPI int IntersectTriangles (const triangle< rational > &tri0, const triangle< rational > &tri1, Ptr< point< rational > > &retP)
GULAPI bool IntersectTrianglesUV (const triangle< rational > &tri1, const triangle2< rational > &tri1uv, const triangle< rational > &tri2, const triangle2< rational > &tri2uv, point2< rational > *S1uv, int *S1flag, point2< rational > *S2uv, int *S2flag, point< rational > *S)
template<class T> GULAPI void BarycentricCoordinates (const gul::point< T > &P1, const gul::point< T > &P2, const gul::point< T > &P3, const gul::point< T > &Pin, T *w, T *u, T *v)
template<class T> GULAPI void BarycentricCoordinates (const gul::point< T > &P1, const gul::point< T > &P2, const gul::point< T > &P3, const gul::point< T > &P12, const gul::point< T > &P13, const gul::point< T > &P23, const gul::point< T > &N, const gul::point< T > &Pin, T *w, T *u, T *v)
GULAPI int IntersectTriangles (const gul::triangle< rational > &tri0, const gul::triangle< rational > &tri1, gul::Ptr< gul::point< rational > > &retP)
GULAPI bool IntersectTrianglesUV (const gul::triangle< rational > &tri1, const gul::triangle2< rational > &tri1uv, const gul::triangle< rational > &tri2, const gul::triangle2< rational > &tri2uv, gul::point2< rational > *S1uv, int *S1flag, gul::point2< rational > *S2uv, int *S2flag, gul::point< rational > *S)

Function Documentation

template<class T>
GULAPI void BarycentricCoordinates (  const gul::point< T > &    P1,
const gul::point< T > &    P2,
const gul::point< T > &    P3,
const gul::point< T > &    P12,
const gul::point< T > &    P13,
const gul::point< T > &    P23,
const gul::point< T > &    N,
const gul::point< T > &    Pin,
T *    w,
T *    u,
T *    v
) 
 

template<class T>
GULAPI void BarycentricCoordinates (  const gul::point< T > &    P1,
const gul::point< T > &    P2,
const gul::point< T > &    P3,
const gul::point< T > &    Pin,
T *    w,
T *    u,
T *    v
) 
 

template GULAPI void BarycentricCoordinates (  const point< rational > &    P1,
const point< rational > &    P2,
const point< rational > &    P3,
const point< rational > &    P12,
const point< rational > &    P13,
const point< rational > &    P23,
const point< rational > &    N,
const point< rational > &    Pin,
rational *    w,
rational *    u,
rational *    v
) 
 

template<class T>
GULAPI void BarycentricCoordinates (  const point< T > &    P1,
const point< T > &    P2,
const point< T > &    P3,
const point< T > &    P12,
const point< T > &    P13,
const point< T > &    P23,
const point< T > &    N,
const point< T > &    Pin,
T *    w,
T *    u,
T *    v
) 
 

Definition at line 220 of file guar_intersect.cpp. 

00226 {
00227   point2<T> W,U,V,P;
00228   point2<T> WU,WV,UV,WP,UP;
00229   point<T> UVW,PVW,PWU,PUV;
00230   T aUVW,aPVW,aPWU,aPUV;
00231   rational one_neg(ULong(1),-1);
00232 
00233   // calculate barycentric coordinates in xy-plane
00234   if( test(N.z) )
00235   {
00236     W.x = P1.x;
00237     W.y = P1.y;
00238 
00239     U.x = P2.x;
00240     U.y = P2.y;
00241 
00242     V.x = P3.x;
00243     V.y = P3.y;
00244 
00245     WU.x = P12.x;
00246     WU.y = P12.y;
00247 
00248     WV.x = P13.x;
00249     WV.y = P13.y;
00250 
00251     UV.x = P23.x;
00252     UV.y = P23.y;
00253 
00254     P.x = Pin.x;
00255     P.y = Pin.y;
00256   }
00257   else if( test(N.y) ) // calculate barycentric coordinates in xz-plane
00258   {
00259     W.x = P1.x;
00260     W.y = P1.z;
00261 
00262     U.x = P2.x;
00263     U.y = P2.z;
00264 
00265     V.x = P3.x;
00266     V.y = P3.z;
00267 
00268     WU.x = P12.x;
00269     WU.y = P12.z;
00270 
00271     WV.x = P13.x;
00272     WV.y = P13.z;
00273 
00274     UV.x = P23.x;
00275     UV.y = P23.z;
00276 
00277     P.x = Pin.x;
00278     P.y = Pin.z;
00279   }
00280   else // calculate barycentric coordinates in yz-plane
00281   {
00282     gul::Assert<InternalError>( ndebug || test(N.x) );
00283 
00284     W.x = P1.y;
00285     W.y = P1.z;
00286 
00287     U.x = P2.y;
00288     U.y = P2.z;
00289 
00290     V.x = P3.y;
00291     V.y = P3.z;
00292 
00293     WU.x = P12.y;
00294     WU.y = P12.z;
00295 
00296     WV.x = P13.y;
00297     WV.y = P13.z;
00298 
00299     UV.x = P23.y;
00300     UV.y = P23.z;
00301 
00302     P.x = Pin.y;
00303     P.y = Pin.z;
00304   }
00305 
00306   /*
00307   cout << "N = (" << N << ")\n";
00308   cout << "P12 = (" << P12 << "), ";
00309   cout << "WU = (" << WU << ")\n";
00310   cout << "P13 = (" << P13 << "), ";
00311   cout << "WV = (" << WV << ")\n";
00312   */
00313 
00314   UVW = cross_product(WU,WV);
00315   aUVW = UVW.z;
00316   if( test(aUVW) < 0 ) aUVW = aUVW * one_neg;
00317 
00318   WP = P-W;
00319   PVW = cross_product(WP,WV);
00320   aPVW = PVW.z;
00321   if( test(aPVW) < 0 ) aPVW = aPVW * one_neg;
00322 
00323   *u = aPVW/aUVW;
00324 
00325   PWU = cross_product(WP,WU);
00326   aPWU = PWU.z;
00327   if( test(aPWU) < 0 ) aPWU = aPWU *one_neg;
00328    
00329   *v = aPWU/aUVW;
00330 
00331   UP = P-U;    
00332   PUV = cross_product(UV,UP);
00333   aPUV = PUV.z;
00334   if( test(aPUV) < 0 ) aPUV = aPUV * one_neg;
00335 
00336   *w = aPUV/aUVW;
00337 }

template GULAPI void BarycentricCoordinates (  const point< rational > &    P1,
const point< rational > &    P2,
const point< rational > &    P3,
const point< rational > &    Pin,
rational *    w,
rational *    u,
rational *    v
) 
 

template<class T>
GULAPI void BarycentricCoordinates (  const point< T > &    P1,
const point< T > &    P2,
const point< T > &    P3,
const point< T > &    Pin,
T *    w,
T *    u,
T *    v
) 
 

Definition at line 104 of file guar_intersect.cpp. 

00108 {
00109   point2<T> W,U,V,P;
00110   point2<T> WU,WV,UV,WP,UP;
00111   point<T> UVW,PVW,PWU,PUV,P12,P13,N;
00112   T aUVW,aPVW,aPWU,aPUV;
00113   rational one_neg(ULong(1),-1);
00114 
00115   P12 = P2 - P1;
00116   P13 = P3 - P1;
00117   N = cross_product( P12, P13 );
00118 
00119   // calculate barycentric coordinates in xy-plane
00120   if( test(N.z) )
00121   {
00122     W.x = P1.x;
00123     W.y = P1.y;
00124 
00125     U.x = P2.x;
00126     U.y = P2.y;
00127 
00128     V.x = P3.x;
00129     V.y = P3.y;
00130 
00131     WU.x = P12.x;
00132     WU.y = P12.y;
00133 
00134     WV.x = P13.x;
00135     WV.y = P13.y;
00136 
00137     P.x = Pin.x;
00138     P.y = Pin.y;
00139   }
00140   else if( test(N.y) ) // calculate barycentric coordinates in xz-plane
00141   {
00142     W.x = P1.x;
00143     W.y = P1.z;
00144 
00145     U.x = P2.x;
00146     U.y = P2.z;
00147 
00148     V.x = P3.x;
00149     V.y = P3.z;
00150 
00151     WU.x = P12.x;
00152     WU.y = P12.z;
00153 
00154     WV.x = P13.x;
00155     WV.y = P13.z;
00156 
00157     P.x = Pin.x;
00158     P.y = Pin.z;
00159   }
00160   else // calculate barycentric coordinates in yz-plane
00161   {
00162     gul::Assert<InternalError>( ndebug || test(N.x) );
00163 
00164     W.x = P1.y;
00165     W.y = P1.z;
00166 
00167     U.x = P2.y;
00168     U.y = P2.z;
00169 
00170     V.x = P3.y;
00171     V.y = P3.z;
00172 
00173     WU.x = P12.y;
00174     WU.y = P12.z;
00175 
00176     WV.x = P13.y;
00177     WV.y = P13.z;
00178 
00179     P.x = Pin.y;
00180     P.y = Pin.z;
00181   }
00182 
00183   UVW = cross_product(WU,WV);
00184   aUVW = UVW.z;
00185   if( test(aUVW) < 0 ) aUVW = aUVW * one_neg;
00186 
00187   WP = P-W;
00188   PVW = cross_product(WP,WV);
00189   aPVW = PVW.z;
00190   if( test(aPVW) < 0 ) aPVW = aPVW * one_neg;
00191 
00192   *u = aPVW/aUVW;
00193 
00194   PWU = cross_product(WP,WU);
00195   aPWU = PWU.z;
00196   if( test(aPWU) < 0 ) aPWU = aPWU *one_neg;
00197    
00198   *v = aPWU/aUVW;
00199 
00200   UV = V-U;
00201   UP = P-U;    
00202 
00203   PUV = cross_product(UV,UP);
00204   aPUV = PUV.z;
00205   if( test(aPUV) < 0 ) aPUV = aPUV * one_neg;
00206 
00207   *w = aPUV/aUVW;
00208 }

int DoubleToInt (  const double    d,
const unsigned int    na,
unsigned long *    a,
unsigned int *    retLen
) 
 

int DoubleToInt (  const double    d,
const unsigned int    na,
unsigned long *    a,
int *    retLen
) 
 

Definition at line 349 of file guar_exact.cpp. 

00351 {
00352   double f,m;
00353   int i,sign,len;
00354   
00355   if( d < 0.0 )
00356   {
00357     sign = -1;
00358     m = fabs( d );
00359   }
00360   else
00361   {
00362     sign = 0;
00363     m = d;
00364   }
00365 
00366   len = 0;
00367   while( m >= 1.0 )
00368   {
00369     m = m / (double)(LL(0x100000000));
00370     len++;
00371   }
00372   for( i = len-1; i >= 0; i-- )
00373   {
00374     m = m*(double)(LL(0x100000000));
00375     f = floor(m); 
00376     a[i] = (unsigned long)f;
00377     m -= f;
00378   }
00379 
00380   *retLen = len;
00381   return(sign);
00382 }

double GetRatio (  const Interval &    a )  [inline]
 

Definition at line 471 of file guar_exact.h. 

00472 {
00473   double h,l;
00474 
00475   h = a.m_high;
00476   l = a.m_low;
00477 
00478   if( l < 0.0 )
00479   {
00480     if( h >= 0.0 ) return gul::rtr<double>::giant(); // relative error = inf
00481 
00482     h = -l;
00483     l = -h;
00484   }
00485    
00486   if( l == 0.0 )
00487   {
00488     if( h == 0.0 ) return 1.0;
00489     return gul::rtr<double>::giant();
00490   }
00491 
00492   return h/l;
00493 }

void guar::IntAdd (  const unsigned int    na,
const unsigned long *    a,
const unsigned int    nb,
const unsigned long *    b,
unsigned int *    nc,
unsigned long *    c
) 
 

Definition at line 37 of file guar_exact.cpp. 

00040 {
00041   uint64 l,z=0;
00042   unsigned int i,len;
00043   
00044   if( na < nb )
00045   {
00046     len = nb+1;
00047     
00048     for( i = 0; i < na; i++ )
00049     {
00050       l = (uint64)a[i] + (uint64)b[i] + z;
00051       c[i] = (unsigned long)(l & 0xffffffff);
00052       z = (l >> 32) & 0xffffffff;
00053     }
00054     for( i = na; i < nb; i++ )
00055     {
00056       l = (uint64)b[i] + z;
00057       c[i] = (unsigned long)(l & 0xffffffff);
00058       z = (l >> 32) & 0xffffffff;
00059     }
00060     if( (c[i] = (unsigned long)z) == 0 ) len -= 1;
00061   }
00062   else
00063   {
00064     len = na+1;
00065     
00066     for( i = 0; i < nb; i++ )
00067     {
00068       l = (uint64)a[i] + (uint64)b[i] + z;
00069       c[i] = (unsigned long)(l & 0xffffffff);
00070       z = (l >> 32) & 0xffffffff;
00071     }
00072     for( i = nb; i < na; i++ )
00073     {
00074       l = (uint64)a[i] + z;
00075       c[i] = (unsigned long)(l & 0xffffffff);
00076       z = (l >> 32) & 0xffffffff;
00077     }
00078     if( (c[i] = (unsigned long)z) == 0 ) len -= 1;
00079   }
00080 
00081   *nc = len;
00082 }

int guar::IntDiv (  const unsigned int    na,
const unsigned long *    a,
const unsigned int    nb,
const unsigned long *    b,
unsigned int *    nq,
unsigned long *    q,
unsigned int *    nr,
unsigned long *    r
) 
 

Definition at line 235 of file guar_exact.cpp.

Referenced by guar::rational::div_mod(). 

00239 {
00240   unsigned long *ad, *bd;
00241   uint64 d, z, l, b1, b2, m, ml, qh;
00242   unsigned int i,j;
00243   int sign;
00244 
00245   if( na == 0 )
00246   {
00247     *nq = 0;
00248     *nr = 0;
00249     return(-1);
00250   }
00251 
00252   if( nb < 2 )
00253   {
00254     if( nb == 0 ) return(0);
00255     
00256     r[0] = IntDivShort( na, a, b[0], nq, q );
00257     if( r[0] != 0 ) *nr = 1; else *nr = 0;
00258     return(-1);
00259   }
00260 
00261   if( na < nb )
00262   {
00263     memcpy( r, a, na*sizeof(unsigned long) );
00264     *nr = na;
00265     *nq = 0;
00266     return(-1);
00267   }
00268   
00269   ad = (unsigned long *)alloca( sizeof(unsigned long)*(na+1) );
00270   if( ad == NULL ) return(0);
00271   bd = (unsigned long *)alloca( sizeof(unsigned long)*nb );
00272   if( bd == NULL ) return(0);
00273 
00274   d = LL(0x100000000)/(((uint64)b[nb-1]) + LL(1));
00275   
00276   ad[na] = IntMultShort( na, a, (unsigned long)d, ad );
00277   IntMultShort( nb, b, (unsigned long)d, bd );
00278   
00279   b1 = bd[nb-1];
00280   b2 = bd[nb-2];
00281 
00282   for( i = na; i >= nb; i-- )
00283   {
00284     l = (((uint64)ad[i])<<32) + (uint64)ad[i-1];
00285 
00286     if( ad[i] == bd[nb-1] )
00287       qh = 0xffffffff;
00288     else
00289       qh =  l/b1;
00290       
00291     while( b2*qh > ((l - b1*qh)<<32) + (uint64)ad[i-2] ) qh--;
00292         
00293     z = 1;    
00294     sign = 0;
00295     m = 0;
00296     for( j = 0; j < nb; j++ )
00297     {
00298       m = qh*(uint64)bd[j] + m;
00299       ml = m & 0xffffffff;
00300       m = m >> 32;
00301       
00302       l = 0xffffffff + (uint64)ad[i-nb+j] - ml + z;
00303       ad[i-nb+j] = (unsigned long)(l & 0xffffffff);
00304       z = (l >> 32) & 0xffffffff;
00305     }        
00306     l = 0xffffffff + (uint64)ad[i] - m + z;
00307     ad[i] = (unsigned long)(l & 0xffffffff);
00308     z = (l >> 32) & 0xffffffff;
00309     sign = (z - LL(1) != 0 ? -1 : 0);
00310     
00311     if( sign )        /* add borrow back */
00312     {
00313       qh--;
00314 
00315       z = 0;
00316       for( j = 0; j < nb; j++ )
00317       {
00318         l = (uint64)ad[i-nb+j] + (uint64)bd[j] + z;
00319         ad[i-nb+j] = (unsigned long)(l & 0xffffffff);
00320         z = (l >> 32) & 0xffffffff;
00321       }        
00322       l = (uint64)ad[i] + z;
00323       ad[i] = (unsigned long)(l & 0xffffffff);
00324     }
00325 
00326     q[i-nb] = (unsigned long)qh;
00327   }
00328   
00329   for( i = nb-1; (i > 0) && (ad[i]==0); i-- );
00330 
00331   if( (i > 0) || (ad[0] != 0) )
00332     // IntDivShort( nb, ad, (unsigned long)d, nr, r );
00333     IntDivShort( i+1, ad, (unsigned long)d, nr, r );
00334   else
00335     *nr = 0;
00336 
00337   if( q[na-nb] == 0 )
00338     *nq = na-nb;
00339   else
00340     *nq = na-nb+1; 
00341 
00342   return(-1);
00343 }

unsigned long guar::IntDivShort (  const unsigned int    na,
const unsigned long *    a,
const unsigned long    b,
unsigned int *    nq,
unsigned long *    q
) 
 

Definition at line 209 of file guar_exact.cpp. 

00212 {
00213   uint64 l, z, d = b;
00214   int i;
00215   
00216   z = 0;
00217   for( i = na-1; i >= 0; i-- )
00218   {
00219     l = (z<<32) + ((uint64)a[i]);
00220     q[i] = (unsigned long)(l / d);
00221     z = l % d;
00222   }
00223 
00224   if( q[na-1] == 0 )
00225     *nq = na-1;
00226   else
00227     *nq = na;
00228 
00229   return((unsigned long)z);
00230 }

GULAPI int IntersectTriangles (  const gul::triangle< rational > &    tri0,
const gul::triangle< rational > &    tri1,
gul::Ptr< gul::point< rational > > &    retP
) 
 

GULAPI int IntersectTriangles (  const triangle< rational > &    tri0,
const triangle< rational > &    tri1,
Ptr< point< rational > > &    retP
) 
 

Definition at line 355 of file guar_intersect.cpp. 

00359 {
00360   point<rational>  v11,v12,v13,n,v21,v22,v23,g,t;
00361   point<rational>  P1,P2,P3;
00362   rational A,B,C,R,a,b,c,r,N;
00363   rational one(ULong(1)),zero;
00364   rational lambda,mu,l1,r1,l2,r2,l;
00365   bool result,init;
00366   int dim;
00367 
00368   // calculate plane equation for 1.triangle
00369   // (Ax + By + Cz = R)
00370 
00371   v11 = tri0.P2 - tri0.P1;       // 1. direction vector
00372   v12 = tri0.P3 - tri0.P1;       // 2. direction vector
00373   v13 = tri0.P3 - tri0.P2;       // (needed later)
00374   n = cross_product( v11, v12 ); // normal vector of the plane
00375   A = n.x;  B = n.y;  C = n.z;
00376   R = n * tri0.P1;               // right side of plane equation
00377 
00378   // calculate plane equation for 2.triangle
00379   // (ax + by + cz = r)
00380   v21 = tri1.P2 - tri1.P1;       // 1. direction vector
00381   v22 = tri1.P3 - tri1.P1;       // 2. direction vector
00382   v23 = tri1.P3 - tri1.P2;       // (needed later)
00383   n = cross_product( v21, v22 ); // normal vector of the plane
00384   a = n.x;  b = n.y;  c = n.z;
00385   r = n * tri1.P1;               // right side of plane equation
00386 
00387   // calculate intersection line
00388   // (G : x = g + lambda * t)
00389 
00390   // if intersection line is not parallel to the X-axis,
00391   // set g.x = 0 and t.x = 1
00392   // ==> g.y, g.z and t.x, t.z
00393    
00394   N = B*c - b*C;    // determinant
00395   dim = 1;          // assume intersection has dimension 1
00396 
00397   if( test(N) )     // if not parallel to the X-axis 
00398   {
00399     g.x = zero;
00400     g.y = (R*c - r*C)/N;
00401     g.z = (B*r - b*R)/N;
00402     t.x = one;
00403     t.y = (C*a - c*A)/N;
00404     t.z = (A*b - a*B)/N;    
00405   }
00406   else
00407   {
00408     N = A*c - a*C;        
00409     
00410     if( test(N) )   // if not parallel to the Y-axis
00411     {
00412       g.x = (R*c - r*C) / N;
00413       g.y = zero;
00414       g.z = (A*r - a*R) / N;
00415       t.x = (C*b - c*B) / N;
00416       t.y = one;
00417       t.z = (B*a - b*A) / N;
00418     }
00419     else
00420     {
00421       N = A*b - a*B;
00422           
00423       if( test(N) )   // if not parallel to the Z-axis
00424       {
00425         g.x = (R*b - r*B) / N;
00426         g.y = (A*r - a*R) / N;
00427         g.z = zero;
00428         t.x = (B*c - b*C) / N;
00429         t.y = (C*a - c*A) / N;
00430         t.z = one;      
00431       }
00432       else        // normal vectors of the two planes are linear dependent
00433       {
00434         // check if one point of the second triangle lies within the plane
00435         // of the first triangle
00436         if( test(A*tri1.P1.x + B*tri1.P1.y + C*tri1.P1.z - R) )
00437           return 0; // no intersection
00438         else
00439           dim = 2;
00440       }  
00441     }
00442   }
00443 
00444   if( dim == 1 )
00445   {
00446     // intersect intersection line of the planes with the border of the first
00447     // triangle. the result is an interval [l1,r1] in the parametric domain
00448     // of the intersection line.
00449 
00450     // convenient notation
00451     P1 = tri0.P1;
00452     P2 = tri0.P2;
00453     P3 = tri0.P3;
00454 
00455     // intersection line: X = g + mu * t
00456     // line through triangle edge P1P2:  X = P1 + lambda * v11
00457 
00458     init = false;
00459     result = RegularIntersectLines( P1, v11, g, t, &lambda, &mu );
00460     if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00461     {
00462       init = true;
00463       l1 = mu;
00464       r1 = mu;
00465     }      
00466 
00467     // line through triangle edge P1P3:  X = P1 + lambda * v12
00468 
00469     result = RegularIntersectLines( P1, v12, g, t, &lambda, &mu );
00470     if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00471     {
00472       if( init )
00473       {
00474         if( compare(mu,l1) < 0 ) l1 = mu;
00475         else if( compare(mu,r1) > 0 ) r1 = mu;
00476       }
00477       else
00478       {
00479         init = true;
00480         l1 = mu;
00481         r1 = mu;
00482       }      
00483     } 
00484 
00485     // line through triangle edge P2P3:  X = P2 + lambda * v13
00486 
00487     result = RegularIntersectLines( P2, v13, g, t, &lambda, &mu );
00488     if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00489     {
00490       if( init )
00491       {
00492         if( compare(mu,l1) < 0 ) l1 = mu;
00493         else if( compare(mu,r1) > 0 ) r1 = mu;
00494       }
00495       else
00496       {
00497         init = true;
00498         l1 = mu;
00499         r1 = mu;
00500       }      
00501     } 
00502 
00503     if( !init )
00504       return 0;
00505     cout << "l1 = " << l1.dump() << ", r1 = " << r1.dump() << "\n";
00506 
00507     // intersect intersection line of the planes with the border of the second
00508     // triangle. the result is an interval [l2,r2] in the parametric domain
00509     // of the intersection line.
00510 
00511     // convenient notation
00512     P1 = tri1.P1;
00513     P2 = tri1.P2;
00514     P3 = tri1.P3;
00515 
00516     // intersection line: X = g + mu * t
00517     // line through triangle edge P1P2:  X = P1 + lambda * v21
00518 
00519     init = false;
00520     result = RegularIntersectLines( P1, v21, g, t, &lambda, &mu );
00521     if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00522     {
00523       init = true;
00524       l2 = mu;
00525       r2 = mu;
00526     }      
00527   
00528     // line through triangle edge P1P3:  X = P1 + lambda * v22
00529 
00530     result = RegularIntersectLines( P1, v22, g, t, &lambda, &mu );
00531     if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00532     {
00533       if( init )
00534       {
00535         if( compare(mu,l2) < 0 ) l2 = mu;
00536         else if( compare(mu,r2) > 0 ) r2 = mu;
00537       }
00538       else
00539       {
00540         init = true;
00541         l2 = mu;
00542         r2 = mu;
00543       }      
00544     } 
00545 
00546     // line through triangle edge P2P3:  X = P2 + lambda * v23
00547 
00548     result = RegularIntersectLines( P2, v23, g, t, &lambda, &mu );
00549     if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
00550     {
00551       if( init )
00552       {
00553         if( compare(mu,l2) < 0 ) l2 = mu;
00554         else if( compare(mu,r2) > 0 ) r2 = mu;
00555       }
00556       else
00557       {
00558         init = true;
00559         l2 = mu;
00560         r2 = mu;
00561       }      
00562     } 
00563 
00564     if( !init )
00565       return 0;
00566     cout << "l2 = " << l2.dump() << ", r2 = " << r2.dump() << "\n";
00567 
00568     // intersection of the intervals [l1,r1] and [l2,r2]  delivers the 
00569     // resulting segment of the intersection line 
00570 
00571     if(compare(l1,l2)>0) l = l1; else l = l2;
00572     if(compare(r1,r2)<0) r = r1; else r = r2;
00573     if( compare(r,l) < 0 )
00574       return 0;
00575 
00576     // calculate the endpoints of the intersection segment
00577     retP[0] = g + l * t;
00578     retP[1] = g + r * t;
00579     return 2;
00580   }
00581 
00582   return 0;
00583 }

GULAPI bool IntersectTrianglesUV (  const gul::triangle< rational > &    tri1,
const gul::triangle2< rational > &    tri1uv,
const gul::triangle< rational > &    tri2,
const gul::triangle2< rational > &    tri2uv,
gul::point2< rational > *    S1uv,
int *    S1flag,
gul::point2< rational > *    S2uv,
int *    S2flag,
gul::point< rational > *    S
) 
 

GULAPI bool IntersectTrianglesUV (  const triangle< rational > &    tri1,
const triangle2< rational > &    tri1uv,
const triangle< rational > &    tri2,
const triangle2< rational > &    tri2uv,
point2< rational > *    S1uv,
int *    S1flag,
point2< rational > *    S2uv,
int *    S2flag,
point< rational > *    S
) 
 

Definition at line 593 of file guar_intersect.cpp. 

00599 {
00600   point<rational>  v11,v12,v13,n,N,v21,v22,v23,g,t;
00601   point<rational>  P1,P2,P3;
00602   rational A,B,C,R,a,b,c,r,T;
00603   rational one(ULong(1)),zero;
00604   rational lambda,mu,l1,r1,l2,r2,l;
00605   bool result,init;
00606   int dim;
00607   bool inter1_12, inter1_13, inter1_23;
00608   bool inter2_12, inter2_13, inter2_23;
00609   int  mid1,mid2,t0,t1,sign,side;
00610   point<rational> Pc1,Pc2;
00611   rational u,v,w;
00612   point2<rational> Pc1uv, Pc2uv;
00613 
00614   // calculate plane equation for 1.triangle
00615   // (Ax + By + Cz = R)
00616 
00617   v11 = tri1.P2 - tri1.P1;       // 1. direction vector
00618   v12 = tri1.P3 - tri1.P1;       // 2. direction vector
00619   v13 = tri1.P3 - tri1.P2;       // (needed later)
00620   N = cross_product( v11, v12 ); // normal vector of the plane
00621   A = N.x;  B = N.y;  C = N.z;
00622   R = N * tri1.P1;               // right side of plane equation
00623   // cout << "N = (" << N << ")\n";
00624 
00625   // calculate plane equation for 2.triangle
00626   // (ax + by + cz = r)
00627   v21 = tri2.P2 - tri2.P1;       // 1. direction vector
00628   v22 = tri2.P3 - tri2.P1;       // 2. direction vector
00629   v23 = tri2.P3 - tri2.P2;       // (needed later)
00630   n = cross_product( v21, v22 ); // normal vector of the plane
00631   a = n.x;  b = n.y;  c = n.z;
00632   r = n * tri2.P1;               // right side of plane equation
00633 
00634   // calculate intersection line
00635   // (G : x = g + lambda * t)
00636 
00637   // if intersection line is not parallel to the X-axis,
00638   // set g.x = 0 and t.x = 1
00639   // ==> g.y, g.z and t.x, t.z
00640    
00641   T = B*c - b*C;    // determinant
00642   dim = 1;          // assume intersection has dimension 1
00643 
00644   if( test(T) )     // if not parallel to the X-axis 
00645   {
00646     g.x = zero;
00647     g.y = (R*c - r*C)/T;
00648     g.z = (B*r - b*R)/T;
00649     t.x = one;
00650     t.y = (C*a - c*A)/T;
00651     t.z = (A*b - a*B)/T;    
00652   }
00653   else
00654   {
00655     T = A*c - a*C;        
00656     
00657     if( test(T) )   // if not parallel to the Y-axis
00658     {
00659       g.x = (R*c - r*C) / T;
00660       g.y = zero;
00661       g.z = (A*r - a*R) / T;
00662       t.x = (C*b - c*B) / T;
00663       t.y = one;
00664       t.z = (B*a - b*A) / T;
00665     }
00666     else
00667     {
00668       T = A*b - a*B;
00669           
00670       if( test(T) )   // if not parallel to the Z-axis
00671       {
00672         g.x = (R*b - r*B) / T;
00673         g.y = (A*r - a*R) / T;
00674         g.z = zero;
00675         t.x = (B*c - b*C) / T;
00676         t.y = (C*a - c*A) / T;
00677         t.z = one;      
00678       }
00679       else        // normal vectors of the two planes are linear dependent
00680       {
00681         // check if one point of the second triangle lies within the plane
00682         // of the first triangle
00683         if( test(A*tri2.P1.x + B*tri2.P1.y + C*tri2.P1.z - R) )
00684           return false; // no intersection
00685         else
00686           dim = 2;
00687       }  
00688     }
00689   }
00690 
00691   if( dim == 1 )
00692   {
00693     // intersect intersection line of the planes with the border of the first
00694     // triangle. the result is an interval [l1,r1] in the parametric domain
00695     // of the intersection line.
00696 
00697     mid1 = false; 
00698     inter1_12 = false;
00699     inter1_13 = false;
00700     inter1_23 = false;
00701 
00702     // convenient notation
00703     P1 = tri1.P1;
00704     P2 = tri1.P2;
00705     P3 = tri1.P3;
00706 
00707     // intersection line: X = g + mu * t
00708     // line through triangle edge P1P2:  X = P1 + lambda * v11
00709 
00710     init = false;
00711     result = RegularIntersectLines( P1, v11, g, t, &lambda, &mu );
00712     if( result )
00713     {
00714       t0 = test(lambda);
00715       if( t0 >= 0 )
00716       {
00717         t1 = compare(lambda,one);
00718         if( t1 <= 0 )
00719         {
00720           inter1_12 = true;
00721 
00722           init = true;
00723           l1 = mu;
00724           r1 = mu;
00725 
00726           if( (t0 > 0) && (t1 < 0) )
00727           {
00728             mid1 = true;
00729             Pc1 = tri1.P1;  // will be used to classify the two halves of
00730                             // tri 0 left and right of the intersection line
00731                             // as lying above or below tri2
00732           }
00733         }
00734       }
00735     }
00736 
00737     // line through triangle edge P1P3:  X = P1 + lambda * v12
00738 
00739     result = RegularIntersectLines( P1, v12, g, t, &lambda, &mu );
00740     if( result )
00741     {
00742       t0 = test(lambda);
00743       if( t0 >= 0 )
00744       {
00745         t1 = compare(lambda,one);
00746         if( t1 <= 0 )
00747         {
00748           inter1_13 = true;
00749 
00750           if( init )
00751           {
00752             if( compare(mu,l1) < 0 ) l1 = mu;
00753             else if( compare(mu,r1) > 0 ) r1 = mu;
00754           }
00755           else
00756           {
00757             init = true;
00758             l1 = mu;
00759             r1 = mu;
00760           }
00761 
00762           if( (t0 > 0) && (t1 < 0) )
00763           { 
00764             mid1 = true;
00765             Pc1 = tri1.P1;  // will be used to classify the two halves of
00766                             // tri 0 left and right of the intersection line
00767                             // as lying above or below tri2
00768           }
00769         }
00770       }
00771     }
00772 
00773     // line through triangle edge P2P3:  X = P2 + lambda * v13
00774 
00775     result = RegularIntersectLines( P2, v13, g, t, &lambda, &mu );
00776     if( result )
00777     {
00778       t0 = test(lambda);
00779       if( t0 >= 0 )
00780       {
00781         t1 = compare(lambda,one);
00782         if( t1 <= 0 )
00783         {
00784           inter1_23 = true;
00785 
00786           if( init )
00787           {
00788             if( compare(mu,l1) < 0 ) l1 = mu;
00789             else if( compare(mu,r1) > 0 ) r1 = mu;
00790           }
00791           else
00792           {
00793             init = true;
00794             l1 = mu;
00795             r1 = mu;
00796           }
00797 
00798           if( (t0 > 0) && (t1 < 0) )
00799           { 
00800             mid1 = true;
00801             Pc1 = tri1.P2;  // will be used to classify the two halves of
00802                             // tri 0 left and right of the intersection line
00803                             // as lying above or below tri2
00804           }
00805         }
00806       }
00807     }
00808 
00809     if( !init )
00810       return 0;
00811     // cout << "l1 = " << l1.dump() << ", r1 = " << r1.dump() << "\n";
00812 
00813     // intersect intersection line of the planes with the border of the second
00814     // triangle. the result is an interval [l2,r2] in the parametric domain
00815     // of the intersection line.
00816 
00817     mid2 = false;
00818     inter2_12 = false;
00819     inter2_13 = false;
00820     inter2_23 = false;
00821 
00822     // convenient notation
00823     P1 = tri2.P1;
00824     P2 = tri2.P2;
00825     P3 = tri2.P3;
00826 
00827     // intersection line: X = g + mu * t
00828     // line through triangle edge P1P2:  X = P1 + lambda * v21
00829 
00830     init = false;
00831     result = RegularIntersectLines( P1, v21, g, t, &lambda, &mu );
00832     if( result )
00833     {
00834       t0 = test(lambda);
00835       if( t0 >= 0 )
00836       {
00837         t1 = compare(lambda,one);
00838         if( t1 <= 0 )
00839         {
00840           inter2_12 = true;
00841 
00842           init = true;
00843           l2 = mu;
00844           r2 = mu;
00845 
00846           if( (t0 > 0) && (t1 < 0) )
00847           {
00848             mid2 = true;
00849             Pc2 = tri2.P1;   // will be used to classify the two halves of
00850                             // tri 0 left and right of the intersection line
00851                             // as lying above or below tri2
00852           }
00853         }
00854       }
00855     }
00856   
00857     // line through triangle edge P1P3:  X = P1 + lambda * v22
00858 
00859     result = RegularIntersectLines( P1, v22, g, t, &lambda, &mu );
00860     if( result )
00861     {
00862       t0 = test(lambda);
00863       if( t0 >= 0 )
00864       {
00865         t1 = compare(lambda,one);
00866         if( t1 <= 0 )
00867         {
00868           inter2_13 = true;
00869 
00870           if( init )
00871           {
00872             if( compare(mu,l2) < 0 ) l2 = mu;
00873             else if( compare(mu,r2) > 0 ) r2 = mu;
00874           }
00875           else
00876           {
00877             init = true;
00878             l2 = mu;
00879             r2 = mu;
00880           }
00881 
00882           if( (t0 > 0) && (t1 < 0) )
00883           {
00884             mid2 = true;
00885             Pc2 = tri2.P1;  // will be used to classify the two halves of
00886                             // tri 0 left and right of the intersection line
00887                             // as lying above or below tri2
00888           }
00889         }
00890       }
00891     }
00892 
00893     // line through triangle edge P2P3:  X = P2 + lambda * v23
00894 
00895     result = RegularIntersectLines( P2, v23, g, t, &lambda, &mu );
00896     if( result )
00897     {
00898       t0 = test(lambda);
00899       if( t0 >= 0 )
00900       {
00901         t1 = compare(lambda,one);
00902         if( t1 <= 0 )
00903         {
00904           inter2_23 = true;
00905 
00906           if( init )
00907           {
00908             if( compare(mu,l2) < 0 ) l2 = mu;
00909             else if( compare(mu,r2) > 0 ) r2 = mu;
00910           }
00911           else
00912           {
00913             init = true;
00914             l2 = mu;
00915             r2 = mu;
00916           }
00917 
00918           if( (t0 > 0) && (t1 < 0) )
00919           {
00920             mid2 = true;
00921             Pc2 = tri2.P2;  // will be used to classify the two halves of
00922                             // tri 0 left and right of the intersection line
00923                             // as lying above or below tri2
00924           }
00925         }
00926       }
00927     }
00928 
00929     if( !init )
00930       return 0;
00931     // cout << "l2 = " << l2.dump() << ", r2 = " << r2.dump() << "\n";
00932 
00933     // intersection of the intervals [l1,r1] and [l2,r2]  delivers the 
00934     // resulting segment of the intersection line 
00935 
00936     if(compare(l1,l2)>0) l = l1; else l = l2;
00937     if(compare(r1,r2)<0) r = r1; else r = r2;
00938     t1 = compare(r,l);
00939     if( t1 < 0 ) return false;
00940     else if( t1 == 0 ) // single point intersection
00941     {
00942       // retP[0] = g + l * t;
00943       return false;
00944     }
00945 
00946     // calculate the endpoints of the intersection segment
00947     S[0] = g + l * t;
00948     S[1] = g + r * t;
00949 
00950     // in triangle coordinates of triangle 1
00951     BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
00952                             S[0], &u, &v, &w );
00953     S1uv[0] = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00954 
00955     BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
00956                             S[1], &u, &v, &w );
00957     S1uv[1] = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00958 
00959     // calculate the edge/side pairs of the intersection segment in tri1
00960     if( mid1 )
00961     {
00962       sign = test(n * Pc1 - r);
00963 
00964       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
00965       // parametric domain are known
00966       BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
00967                               Pc1, &u, &v, &w );
00968       Pc1uv = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00969 
00970       side = test( cross_product(S1uv[1]-S1uv[0],Pc1uv-S1uv[0]).z );
00971 
00972       gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
00973       if( side < 0 ) side = 0; 
00974       else if( side > 0 ) side = 1;
00975       S1flag[side] = sign;
00976       S1flag[!side] = -sign;
00977     }
00978     else
00979     {
00980       if( inter1_12 && inter1_13 )
00981         Pc1 = tri1.P1;
00982       else if( inter1_12 && inter1_23 )
00983         Pc1 = tri1.P2;
00984       else
00985       {
00986         gul::Assert<InternalError>( ndebug || (inter1_13 && inter1_23));
00987         Pc1 = tri1.P3;
00988       }
00989 
00990       sign = test(n * Pc1 - r);
00991 
00992       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
00993       // parametric domain are known
00994       BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
00995                               Pc1, &u, &v, &w );
00996       Pc1uv = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;
00997 
00998       side = test( cross_product(S1uv[1]-S1uv[0],Pc1uv-S1uv[0]).z );
00999 
01000       gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
01001       S1flag[side] = sign;
01002       S1flag[!side] = 0;  // undetermined
01003     }
01004 
01005     // in triangle coordinates of triangle 2
01006     BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
01007                             S[0], &u, &v, &w );
01008     S2uv[0] = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01009 
01010     BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
01011                             S[1], &u, &v, &w );
01012     S2uv[1] = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01013 
01014     // calculate the edge/side pairs of the intersection segment in tri2
01015     if( mid2 )
01016     {
01017       sign = test(N * Pc2 - R);
01018 
01019       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
01020       // parametric domain are known
01021       BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
01022                               Pc2, &u, &v, &w );
01023       Pc2uv = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01024 
01025       side = test( cross_product(S2uv[1]-S2uv[0],Pc2uv-S2uv[0]).z );
01026 
01027       gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
01028       if( side < 0 ) side = 0; 
01029       else if( side > 0 ) side = 1;
01030       S2flag[side] = sign;
01031       S2flag[!side] = -sign;
01032     }
01033     else
01034     {
01035       if( inter2_12 && inter2_13 )
01036         Pc2 = tri2.P1;
01037       else if( inter2_12 && inter2_23 )
01038         Pc2 = tri2.P2;
01039       else
01040       {
01041         gul::Assert<InternalError>( ndebug || (inter2_13 && inter2_23));
01042         Pc2 = tri1.P3;
01043       }
01044 
01045       sign = test(N * Pc2 - R);
01046 
01047       // this is a bit of a overkill since the coordinates of 'Pc1' in the 
01048       // parametric domain are known
01049       BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
01050                               Pc2, &u, &v, &w );
01051       Pc2uv = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;
01052 
01053       side = test( cross_product(S2uv[1]-S2uv[0],Pc2uv-S2uv[0]).z );
01054 
01055       gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
01056       S2flag[side] = sign;
01057       S2flag[!side] = 0;  // undetermined
01058     }
01059 
01060     // two intersection points
01061     return true;
01062   }
01063 
01064   return false;
01065 }

void guar::IntMult (  const unsigned int    na,
const unsigned long *    a,
const unsigned int    nb,
const unsigned long *    b,
unsigned int *    nc,
unsigned long *    c
) 
 

Definition at line 176 of file guar_exact.cpp. 

00179 {
00180   uint64 l,z;
00181   unsigned int i,j;
00182 
00183   if( (na == 0) || (nb == 0) )
00184   { *nc = 0; return; }
00185 
00186   for( i = 0; i < na; i++ )
00187     c[i] = 0;
00188 
00189   for( j = 0; j < nb; j++ )
00190   {
00191     z = 0;
00192     for( i = 0; i < na; i++ )
00193     {
00194       l = (uint64)a[i]*(uint64)b[j] + 
00195           (uint64)c[i+j] + z;
00196       c[i+j] = (unsigned long)(l & 0xffffffff);
00197       z = (l >> 32) & 0xffffffff;
00198     }
00199     c[i+j] = (unsigned long)z;
00200   }
00201 
00202   if( c[na+nb-1] == 0 ) *nc = na+nb-1;
00203   else *nc = na+nb;
00204 }

unsigned long guar::IntMultShort (  const unsigned int    na,
const unsigned long *    a,
const unsigned long    b,
unsigned long *    c
) 
 

Definition at line 153 of file guar_exact.cpp. 

00156 {
00157   uint64 z, l, m;
00158   unsigned int i;
00159   
00160   z = 0;
00161   m = (uint64)b;
00162 
00163   for( i = 0; i < na; i++ )
00164   {
00165     l = (uint64)a[i]*m + z;
00166     c[i] = (unsigned long)(l & 0xffffffff);
00167     z = (l >> 32) & 0xffffffff;
00168   }
00169   return((unsigned long)z);
00170 }

int IntSub (  const int unsigned    na,
const unsigned long *    a,
const int unsigned    nb,
const unsigned long *    b,
unsigned int *    nc,
unsigned long *    c
) 
 

int IntSub (  const unsigned int    na,
const unsigned long *    a,
const unsigned int    nb,
const unsigned long *    b,
unsigned int *    nc,
unsigned long *    c
) 
 

Definition at line 87 of file guar_exact.cpp. 

00090 {
00091   uint64 l,z=1;
00092   unsigned int i,len;
00093   int sign;
00094 
00095   /* test if b > a */
00096   
00097   if( nb > na )
00098     sign = -1;
00099   else if( na > nb )
00100     sign = 0;
00101   else        /* na == nb */
00102   {
00103     i = na-1;
00104     while( (a[i] == b[i]) && (i > 0) )
00105       i--;
00106     sign = b[i] > a[i] ? -1 : 0;
00107   }
00108  
00109   if( sign )       /* calc b - a, na <= nb */
00110   {
00111     for( i = 0; i < na; i++ )
00112     {
00113       l = 0xffffffff + (uint64)b[i] - (uint64)a[i] + z;
00114       c[i] = (unsigned long)(l & 0xffffffff);
00115       z = (l >> 32) & 0xffffffff;
00116     }
00117     for( i = na; i < nb; i++ )
00118     {
00119       l = 0xffffffff + (uint64)b[i] + z;
00120       c[i] = (unsigned long)(l & 0xffffffff);
00121       z = (l >> 32) & 0xffffffff;
00122     }   
00123     len = nb;
00124   }
00125   else             /* calc a - b, na >= nb */
00126   {
00127     for( i = 0; i < nb; i++ )
00128     {
00129       l = 0xffffffff + (uint64)a[i] - (uint64)b[i] + z;
00130       c[i] = (unsigned long)(l & 0xffffffff);
00131       z = (l >> 32) & 0xffffffff;
00132     }
00133     for( i = nb; i < na; i++ )
00134     {
00135       l = 0xffffffff + (uint64)a[i] + z;
00136       c[i] = (unsigned long)(l & 0xffffffff);
00137       z = (l >> 32) & 0xffffffff;
00138     }    
00139     len = na;
00140   }
00141 
00142   /* strip leading zeros */
00143 
00144   for( ; (len > 0) && (c[len-1] == 0); len-- );
00145 
00146   *nc = len;
00147   return( sign );
00148 }

template<class T>
void IntTo (  const unsigned int    na,
const unsigned long *    a,
T &    t
) 
 

Definition at line 38 of file guar_exact.h. 

00039 {
00040   int i;
00041 
00042   if( na == 0 )
00043   {
00044     t = (T)0.0;
00045     return;
00046   }
00047   t = (T)a[na-1];
00048   
00049   for( i = na-2; i >= 0; i-- )
00050   {
00051     t = t*(T)(LL(0x100000000)) + (T)a[i];
00052   }  
00053 }

double IntToDouble (  const unsigned int    na,
const unsigned long *    a
) 
 

GULAPI rational operator * (  const rational &    a,
const rational &    b
) 
 

Definition at line 580 of file guar_exact.cpp. 

00581 {
00582   if( (a.m->m_na == 0) || (b.m->m_na == 0) )  /* a = 0 or b = 0 */
00583   {
00584     rational c;
00585     return(c);
00586   }
00587   else
00588   {
00589     rational c( a.m->m_sign ^ b.m->m_sign, 
00590                       a.m->m_na + b.m->m_na, a.m->m_nb + b.m->m_nb );
00591 
00592     IntMult( a.m->m_na, a.m->m_a, b.m->m_na, b.m->m_a, 
00593              &c.m->m_na, c.m->m_a );
00594 
00595     if( (a.m->m_nb != 0) && (b.m->m_nb != 0) )
00596     {
00597       IntMult( a.m->m_nb, a.m->m_b, b.m->m_nb, b.m->m_b, 
00598                     &c.m->m_nb, c.m->m_b );
00599     }
00600     else
00601     {
00602       if( b.m->m_nb != 0 )
00603       {
00604         memcpy( c.m->m_b, b.m->m_b, b.m->m_nb*sizeof(unsigned long) );
00605         c.m->m_nb = b.m->m_nb; 
00606       }       
00607       else if( a.m->m_nb != 0 )
00608       {
00609         memcpy( c.m->m_b, a.m->m_b, a.m->m_nb*sizeof(unsigned long) );
00610         c.m->m_nb = a.m->m_nb;
00611       }
00612     }
00613       
00614     return(c);
00615   }
00616 }

Interval operator * (  const Interval &    a,
const Interval &    b
) 
 

Definition at line 387 of file guar_exact.cpp. 

00388 {
00389   Interval c;
00390   double f1,f2;
00391   
00392   if( a.m_low >= 0.0 )       /* 0 <= alo <= ahi                   */
00393   {
00394     if( b.m_low >= 0.0 )     /* 0 <= alo <= ahi,  0 <= blo <= bhi */ 
00395     {
00396       c.m_low = a.m_low * b.m_low;
00397       c.m_high = a.m_high * b.m_high;
00398     }
00399     else                     /* blo < 0 */
00400     {
00401       if( b.m_high >= 0 )    /*  0 <= alo <= ahi,  blo < 0 <= bhi */
00402       {
00403         c.m_low = a.m_high * b.m_low;
00404         c.m_high = a.m_high * b.m_high;
00405       }
00406       else                   /*  0 <= alo <= ahi,  blo <= bhi < 0 */
00407       {
00408         c.m_low = a.m_high * b.m_low;
00409         c.m_high = a.m_low * b.m_high;
00410       }
00411     }
00412   }
00413   else                       /* alo < 0 */
00414   {
00415     if( a.m_high < 0.0 )     /* alo <= ahi < 0 */
00416     {
00417       if( b.m_low >= 0 )     /* alo <= ahi < 0,  0 <= blo <= bhi */
00418       {
00419         c.m_low = a.m_low * b.m_high;
00420         c.m_high = a.m_high * b.m_low;
00421       }
00422       else                   /* alo <= ahi < 0,  blo < 0 */
00423       {
00424         if( b.m_high < 0 )   /* alo <= ahi < 0,  blo <= bhi < 0 */
00425         {
00426           c.m_low = a.m_high * b.m_high;
00427           c.m_high = a.m_low * b.m_low;
00428         }
00429         else     /* alo <= ahi < 0,  blo < 0 <= bhi */
00430         {
00431           c.m_low = a.m_low * b.m_high;
00432           c.m_high = a.m_high * b.m_low;
00433         }
00434       }
00435     }
00436     else    /* alo < 0 <= ahi */
00437     {
00438       if( b.m_low >= 0 )   /* alo < 0 <= ahi,  0 <= blo <= bhi */
00439       {
00440         c.m_low = a.m_low * b.m_high;
00441         c.m_high = a.m_high * b.m_high;
00442       }
00443       else                 /* alo < 0 <= ahi,  blo < 0 */
00444       {
00445         if( b.m_high < 0 ) /* alo < 0 <= ahi,  blo <= bhi < 0 */
00446         { 
00447           c.m_low = a.m_high * b.m_low;
00448           c.m_high = a.m_low * b.m_high;
00449         }
00450         else              /* alo < 0 <= ahi,  blo < 0 <= bhi */
00451         {
00452           f1 = a.m_low * b.m_high;
00453           f2 = a.m_high * b.m_low;
00454           c.m_low = f1 < f2 ? f1 : f2;
00455 
00456           f1 = a.m_low * b.m_low;
00457           f2 = a.m_high * b.m_high;
00458           c.m_high = f1 > f2 ? f1 : f2;
00459         }
00460       }
00461     }
00462   }
00463       
00464   c.m_low = c.m_low - DBL_EPSILON*fabs(c.m_low);
00465   c.m_high = c.m_high + DBL_EPSILON*fabs(c.m_high);
00466     
00467   return( c );
00468 }

Interval operator+ (  const Interval &    a,
const Interval &    b
)  [inline]
 

Definition at line 413 of file guar_exact.h. 

00414 {
00415   Interval c;
00416   
00417   c.m_low = a.m_low + b.m_low;
00418   c.m_low = c.m_low - DBL_EPSILON*fabs(c.m_low);
00419 
00420   c.m_high = a.m_high + b.m_high;
00421   c.m_high = c.m_high + DBL_EPSILON*fabs(c.m_high);
00422     
00423   return( c );
00424 }

GULAPI rational operator+ (  const rational &    A,
const rational &    B
) 
 

Definition at line 618 of file guar_exact.cpp. 

00619 {
00620   unsigned long *a1b2,*a2b1,*a1b2_a2b1,*b1b2;
00621   unsigned int na1b2,na2b1,n_a1b2_a2b1,nb1b2;
00622   size_t size_a1b2,size_a2b1,size_a1b2_a2b1,size_b1b2;
00623   int sign,sign2;
00624   
00625   if( A.m->m_na == 0 ) return(B);
00626   if( B.m->m_na == 0 ) return(A);
00627   
00628   rational c;
00629    
00630   if( B.m->m_nb == 0 )
00631   {
00632     a1b2 = A.m->m_a;
00633     na1b2 = A.m->m_na;
00634     size_a1b2 = 0;
00635   }
00636   else
00637   {
00638     a1b2 = (unsigned long *)gust::PoolAlloc( 
00639                      (A.m->m_na+B.m->m_nb)*sizeof(unsigned long), &size_a1b2 );
00640     if( a1b2 == NULL ) throw gul::PoolAllocError();
00641     IntMult( A.m->m_na, A.m->m_a, B.m->m_nb, B.m->m_b, 
00642              &na1b2, a1b2 );
00643   }
00644 
00645   if( A.m->m_nb == 0 )
00646   {
00647     a2b1 = B.m->m_a;
00648     na2b1 = B.m->m_na;
00649     size_a2b1 = 0;
00650   }
00651   else
00652   {
00653     a2b1 = (unsigned long *)gust::PoolAlloc( 
00654                      (B.m->m_na+A.m->m_nb)*sizeof(unsigned long), &size_a2b1 );
00655     if( a2b1 == NULL ) throw gul::PoolAllocError();
00656     IntMult( B.m->m_na, B.m->m_a, A.m->m_nb, A.m->m_b, 
00657              &na2b1, a2b1 );
00658   }
00659 
00660   n_a1b2_a2b1 = na1b2 > na2b1 ? na1b2 : na2b1;
00661   if( n_a1b2_a2b1 == 0 )
00662   {
00663     if( size_a1b2 ) gust::PoolFree( a1b2, size_a1b2 );
00664     if( size_a2b1 ) gust::PoolFree( a2b1, size_a2b1 );
00665     return(c);
00666   }
00667 
00668   sign = A.m->m_sign ^ B.m->m_sign;
00669   if( !sign ) n_a1b2_a2b1++; 
00670 
00671   a1b2_a2b1 = (unsigned long *)gust::PoolAlloc( 
00672                  n_a1b2_a2b1*sizeof(unsigned long), &size_a1b2_a2b1 );
00673   if( a1b2_a2b1 == NULL ) throw gul::PoolAllocError();
00674 
00675   if( sign )
00676   {
00677     sign2 = IntSub( na1b2, a1b2, na2b1, a2b1, &n_a1b2_a2b1, a1b2_a2b1 );
00678     if( n_a1b2_a2b1 == 0 ) 
00679     {
00680       gust::PoolFree( a1b2_a2b1, size_a1b2_a2b1 );
00681       a1b2_a2b1 = NULL;
00682     }
00683   }
00684   else
00685   {
00686     sign2 = 0;
00687     IntAdd( na1b2, a1b2, na2b1, a2b1, &n_a1b2_a2b1, a1b2_a2b1 );
00688   }
00689 
00690   if( size_a1b2 ) gust::PoolFree( a1b2, size_a1b2 );
00691   if( size_a2b1 ) gust::PoolFree( a2b1, size_a2b1 );
00692   
00693   if( n_a1b2_a2b1 == 0 )  /* numerator = 0 */
00694   {
00695     sign = 0;
00696     b1b2 = NULL;
00697     nb1b2 = 0;
00698     size_b1b2 = 0;
00699   }
00700   else                    /* numerator != 0 */
00701   {  
00702     if( A.m->m_sign )
00703       sign = sign2 ^ -1;
00704     else
00705       sign = sign2;
00706      
00707     if( A.m->m_nb == 0 )
00708     {
00709       if( B.m->m_nb == 0 )
00710       {
00711         b1b2 = NULL;
00712         nb1b2 = 0;
00713         size_b1b2 = 0;
00714       }
00715       else
00716       {
00717         nb1b2 = B.m->m_nb;
00718         b1b2 = (unsigned long *)gust::PoolAlloc( 
00719                                     nb1b2*sizeof(unsigned long), &size_b1b2 );
00720         if( b1b2 == NULL ) throw gul::PoolAllocError();
00721         memcpy( b1b2, B.m->m_b, nb1b2*sizeof(unsigned long) );
00722       }
00723     }
00724     else if( B.m->m_nb == 0 )
00725     {
00726       nb1b2 = A.m->m_nb;
00727       b1b2 = (unsigned long *)gust::PoolAlloc( 
00728                                     nb1b2*sizeof(unsigned long), &size_b1b2 );
00729       if( b1b2 == NULL ) throw gul::PoolAllocError();
00730       memcpy( b1b2, A.m->m_b, nb1b2*sizeof(unsigned long) );
00731     }
00732     else
00733     {
00734       b1b2 = (unsigned long *)gust::PoolAlloc( 
00735                 (A.m->m_nb+B.m->m_nb)*sizeof(unsigned long), &size_b1b2 );
00736       if( b1b2 == NULL ) throw gul::PoolAllocError();
00737   
00738       IntMult( A.m->m_nb, A.m->m_b, B.m->m_nb, B.m->m_b,
00739                &nb1b2, b1b2 );
00740     }
00741   }
00742   
00743   c.m->m_sign = sign;
00744   c.m->m_na = n_a1b2_a2b1;
00745   c.m->m_size_a = size_a1b2_a2b1;
00746   c.m->m_a = a1b2_a2b1;
00747   c.m->m_nb = nb1b2;
00748   c.m->m_size_b = size_b1b2;
00749   c.m->m_b = b1b2;
00750   return(c);
00751 }

Interval operator- (  const Interval &    a,
const Interval &    b
)  [inline]
 

Definition at line 426 of file guar_exact.h. 

00427 {
00428   Interval c(-b.m_high,-b.m_low);
00429   
00430   return(a+c);
00431 }

GULAPI rational operator- (  const rational &    A,
const rational &    B
) 
 

Definition at line 767 of file guar_exact.cpp. 

00768 {
00769   rational c;  // (= 0)
00770 
00771   if( A.m == B.m )
00772     return c;
00773 
00774   if( A.m->m_na == 0 )   /* if A=0, A+B returns B !!! */
00775   {
00776     if( B.m->m_na != 0 )
00777     {
00778       c = B.get_copy();
00779       c.m->m_sign = B.m->m_sign ^ -1;
00780     }
00781     return c;
00782   }
00783     
00784   B.m->negative();
00785   c = A + B;
00786   B.m->negative();
00787 
00788   return(c);
00789 }

Interval operator/ (  const Interval &    a,
const Interval &    b
)  [inline]
 

Definition at line 433 of file guar_exact.h. 

00434 {
00435   if( (b.m_low < 0.0) && (b.m_high > 0.0) )
00436     throw gul::IntervalDivZeroError();
00437 
00438   Interval c(1.0/b.m_high,1.0/b.m_low);
00439   
00440   return(a*c);
00441 }

GULAPI rational operator/ (  const rational &    A,
const rational &    B
) 
 

Definition at line 753 of file guar_exact.cpp. 

00754 {
00755   rational c;
00756   
00757   if( A.m == B.m )
00758     return rational(ULong(1));
00759   
00760   B.m->reciprocal();
00761   c = A * B;
00762   B.m->reciprocal();
00763 
00764   return(c);
00765 }

bool RegularIntersectLines (  const point< rational > &    A,
const point< rational > &    B,
const point< rational > &    a,
const point< rational > &    b,
rational *    lambda,
rational *    mu
) 
 

Definition at line 46 of file guar_intersect.cpp. 

00054 {  
00055   rational A0,A1,A2,B0,B1,B2;
00056   rational a0,a1,a2,b0,b1,b2;
00057   rational N1,N2,N3;
00058   
00059   // more convenient notation
00060   A0 = A.x;  A1 = A.y;  A2 = A.z;
00061   B0 = B.x;  B1 = B.y;  B2 = B.z;
00062   a0 = a.x;  a1 = a.y;  a2 = a.z;
00063   b0 = b.x;  b1 = b.y;  b2 = b.z;
00064   
00065   N1 = B1*b0 - b1*B0;
00066 
00067   if( test(N1) ) // calculation without z (for exact arithmetic the
00068   {              // goodness of N1 doesn't matter)
00069     *lambda = (b0*(a1-A1) - b1*(a0-A0)) / N1; 
00070     *mu = (B0*(a1-A1) - B1*(a0-A0)) / N1; 
00071   }
00072   else
00073   {
00074     N2 = B2*b0 - b2*B0; 
00075 
00076     if( test(N2) )  // calculation without y (for exact arithmetic the
00077     {               // goodness of N2 doesn't matter)
00078       *lambda = (b0*(a2-A2) - b2*(a0-A0)) / N2;
00079       *mu = (B0*(a2-A2) - B2*(a0-A0)) / N2;
00080     }
00081     else
00082     {
00083       N3 = B2*b1 - b2*B1;
00084 
00085       if( test(N3) ) // calculation without y (for exact arithmetic the 
00086       {              // goodness of N3 doesn't matter)
00087         *lambda = (b1*(a2-A2) - b2*(a1-A1)) / N3;
00088         *mu = (B1*(a2-A2) - B2*(a1-A1)) / N3;
00089       }
00090       else
00091       {
00092         return false;  // no intersection point or lines colinear
00093       }      
00094     }
00095   }  
00096   return true;
00097 }

Interval Sqrt (  const Interval &    a )  [inline]
 

Definition at line 444 of file guar_exact.h. 

00445 {
00446   Interval c;
00447   
00448   c.m_low = sqrt(a.m_low);
00449   c.m_low = c.m_low - DBL_EPSILON*fabs(c.m_low);
00450 
00451   c.m_high = sqrt(a.m_high);
00452   c.m_high = c.m_high + DBL_EPSILON*fabs(c.m_high);
00453     
00454   return( c );
00455 }

int Test (  const Interval &    a,
const Interval &    b
)  [inline]
 

Definition at line 465 of file guar_exact.h. 

00466 {
00467   Interval c = a-b;
00468   return Test(c);
00469 }

int Test (  const Interval &    a )  [inline]
 

Definition at line 457 of file guar_exact.h. 

00458 {
00459   if( a.m_high < 0.0 ) return(-1);
00460   if( a.m_low > 0.0 ) return(+1);
00461 
00462   return(0);   /* dunno :) */
00463 }

Generated on Mon Jan 21 04:17:48 2002 for GUL 0.6 - Geometry Utility Library by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001