DGtal  0.6.devel
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
IntegerComputer.ih
1 
30 
31 #include <cstdlib>
33 
35 // IMPLEMENTATION of inline methods.
37 
39 // ----------------------- Standard services ------------------------------
40 
41 //-----------------------------------------------------------------------------
42 template <typename TInteger>
43 inline
44 DGtal::IntegerComputer<TInteger>::~IntegerComputer()
45 {}
46 //-----------------------------------------------------------------------------
47 template <typename TInteger>
48 inline
49 DGtal::IntegerComputer<TInteger>::IntegerComputer()
50 {}
51 //-----------------------------------------------------------------------------
52 template <typename TInteger>
53 inline
54 DGtal::IntegerComputer<TInteger>::IntegerComputer( const Self & /*other*/ )
55 {}
56 //-----------------------------------------------------------------------------
57 template <typename TInteger>
58 inline
59 DGtal::IntegerComputer<TInteger> &
60 DGtal::IntegerComputer<TInteger>::operator=( const Self & /*other*/ )
61 {
62  return *this;
63 }
64 //-----------------------------------------------------------------------------
65 // ----------------------- Integer services ------------------------------
66 //-----------------------------------------------------------------------------
67 template <typename TInteger>
68 inline
69 bool
70 DGtal::IntegerComputer<TInteger>::
71 isZero( IntegerParamType a )
72 {
73  return a == NumberTraits<Integer>::ZERO;
74 }
75 //-----------------------------------------------------------------------------
76 template <typename TInteger>
77 inline
78 bool
79 DGtal::IntegerComputer<TInteger>::
80 isNotZero( IntegerParamType a )
81 {
82  return a != NumberTraits<Integer>::ZERO;
83 }
84 //-----------------------------------------------------------------------------
85 template <typename TInteger>
86 inline
87 bool
88 DGtal::IntegerComputer<TInteger>::
89 isPositive( IntegerParamType a )
90 {
91  return a > NumberTraits<Integer>::ZERO;
92 }
93 //-----------------------------------------------------------------------------
94 template <typename TInteger>
95 inline
96 bool
97 DGtal::IntegerComputer<TInteger>::
98 isNegative( IntegerParamType a )
99 {
100  return a < NumberTraits<Integer>::ZERO;
101 }
102 //-----------------------------------------------------------------------------
103 template <typename TInteger>
104 inline
105 bool
106 DGtal::IntegerComputer<TInteger>::
107 isPositiveOrZero( IntegerParamType a )
108 {
109  return a >= NumberTraits<Integer>::ZERO;
110 }
111 //-----------------------------------------------------------------------------
112 template <typename TInteger>
113 inline
114 bool
115 DGtal::IntegerComputer<TInteger>::
116 isNegativeOrZero( IntegerParamType a )
117 {
118  return a <= NumberTraits<Integer>::ZERO;
119 }
120 //-----------------------------------------------------------------------------
121 template <typename TInteger>
122 inline
123 typename DGtal::IntegerComputer<TInteger>::Integer
124 DGtal::IntegerComputer<TInteger>::
125 abs( IntegerParamType a )
126 {
127  if ( isPositiveOrZero( a ) )
128  return a;
129  else
130  return -a;
131 }
132 //-----------------------------------------------------------------------------
133 template <typename TInteger>
134 inline
135 typename DGtal::IntegerComputer<TInteger>::Integer
136 DGtal::IntegerComputer<TInteger>::
137 max( IntegerParamType a, IntegerParamType b )
138 {
139  return ( a >= b ) ? a : b;
140 }
141 //-----------------------------------------------------------------------------
142 template <typename TInteger>
143 inline
144 typename DGtal::IntegerComputer<TInteger>::Integer
145 DGtal::IntegerComputer<TInteger>::
146 max( IntegerParamType a, IntegerParamType b, IntegerParamType c )
147 {
148  return max( max(a,b), c );
149 }
150 //-----------------------------------------------------------------------------
151 template <typename TInteger>
152 inline
153 typename DGtal::IntegerComputer<TInteger>::Integer
154 DGtal::IntegerComputer<TInteger>::
155 min( IntegerParamType a, IntegerParamType b )
156 {
157  return ( a <= b ) ? a : b;
158 }
159 //-----------------------------------------------------------------------------
160 template <typename TInteger>
161 inline
162 typename DGtal::IntegerComputer<TInteger>::Integer
163 DGtal::IntegerComputer<TInteger>::
164 min( IntegerParamType a, IntegerParamType b, IntegerParamType c )
165 {
166  return min( min(a,b), c );
167 }
168 //-----------------------------------------------------------------------------
169 template <typename TInteger>
170 inline
171 void
172 DGtal::IntegerComputer<TInteger>::
173 getEuclideanDiv( Integer & q, Integer & r,
174  IntegerParamType a, IntegerParamType b ) const
175 {
176  q = a / b;
177  r = a % b;
178 }
179 //-----------------------------------------------------------------------------
180 template <typename TInteger>
181 inline
182 typename DGtal::IntegerComputer<TInteger>::Integer
183 DGtal::IntegerComputer<TInteger>::
184 floorDiv( IntegerParamType na, IntegerParamType nb ) const
185 {
186  _m_a = na;
187  _m_b = nb;
188  if ( isNegative( _m_b ) )
189  {
190  _m_a=-_m_a;
191  _m_b=-_m_b;
192  }
193  // if ( isNegative( _m_a ) && isNegative( _m_b ) )
194  // {
195  // _m_a=-_m_a;
196  // _m_b=-_m_b;
197  // }
198  // else if ( isNegative( _m_b ) )
199  // {
200  // _m_a=-_m_a;
201  // _m_b=-_m_b;
202  // }
203  if ( isPositive( _m_a ) || isZero( _m_a % _m_b ) )
204  return _m_a/_m_b;
205  else
206  return _m_a/_m_b - NumberTraits<Integer>::ONE;
207 }
208 //-----------------------------------------------------------------------------
209 template <typename TInteger>
210 inline
211 typename DGtal::IntegerComputer<TInteger>::Integer
212 DGtal::IntegerComputer<TInteger>::
213 ceilDiv( IntegerParamType na, IntegerParamType nb ) const
214 {
215  _m_a = na;
216  _m_b = nb;
217  if ( isNegative( _m_b ) )
218  {
219  _m_a=-_m_a;
220  _m_b=-_m_b;
221  }
222  // if ( isNegative( _m_a ) && isNegative( _m_b ) )
223  // {
224  // _m_a=-_m_a;
225  // _m_b=-_m_b;
226  // }
227  // else if ( isNegative( _m_b ) )
228  // {
229  // _m_a=-_m_a;
230  // _m_b=-_m_b;
231  // }
232  if ( isNegative( _m_a ) || isZero( _m_a % _m_b ) )
233  return _m_a/_m_b;
234  else
235  return _m_a/_m_b + NumberTraits<Integer>::ONE;
236 }
237 //-----------------------------------------------------------------------------
238 template <typename TInteger>
239 inline
240 void
241 DGtal::IntegerComputer<TInteger>::
242 getFloorCeilDiv( Integer & fl, Integer & ce,
243  IntegerParamType na, IntegerParamType nb ) const
244 {
245  _m_a = na;
246  _m_b = nb;
247  if ( isNegative( _m_b ) )
248  {
249  _m_a=-_m_a;
250  _m_b=-_m_b;
251  }
252  fl = ce = _m_a/_m_b;
253  if ( isNotZero( _m_a % _m_b ) )
254  {
255  if ( isNegativeOrZero( _m_a ) ) --fl;
256  if ( isPositiveOrZero( _m_a ) ) ++ce;
257  }
258 }
259 //-----------------------------------------------------------------------------
260 template <typename TInteger>
261 inline
262 typename DGtal::IntegerComputer<TInteger>::Integer
263 DGtal::IntegerComputer<TInteger>::
264 staticGcd( IntegerParamType a, IntegerParamType b )
265 {
266  Integer _m_a = abs( a );
267  Integer _m_b = abs( b );
268  Integer _m_a0 = max( _m_a, _m_b );
269  Integer _m_a1 = min( _m_a, _m_b );
270  Integer _m_r;
271  while ( isNotZero( _m_a1 ) )
272  {
273  _m_r = _m_a0 % _m_a1;
274  _m_a0 = _m_a1;
275  _m_a1 = _m_r;
276  }
277  return _m_a0;
278 }
279 //-----------------------------------------------------------------------------
280 template <typename TInteger>
281 inline
282 typename DGtal::IntegerComputer<TInteger>::Integer
283 DGtal::IntegerComputer<TInteger>::
284 gcd( IntegerParamType a, IntegerParamType b ) const
285 {
286  _m_a = abs( a );
287  _m_b = abs( b );
288  _m_a0 = max( _m_a, _m_b );
289  _m_a1 = min( _m_a, _m_b );
290  while ( isNotZero( _m_a1 ) )
291  {
292  _m_r = _m_a0 % _m_a1;
293  _m_a0 = _m_a1;
294  _m_a1 = _m_r;
295  }
296  return _m_a0;
297 }
298 //-----------------------------------------------------------------------------
299 template <typename TInteger>
300 inline
301 void
302 DGtal::IntegerComputer<TInteger>::
303 getGcd( Integer & g, IntegerParamType a, IntegerParamType b ) const
304 {
305  // std::cerr << "gcd(" << a << ", " << b << ")=";
306  _m_a = abs( a );
307  _m_b = abs( b );
308  _m_a0 = max( _m_a, _m_b );
309  _m_a1 = min( _m_a, _m_b );
310  while ( isNotZero( _m_a1 ) )
311  {
312  _m_r = _m_a0 % _m_a1;
313  _m_a0 = _m_a1;
314  _m_a1 = _m_r;
315  }
316  g = _m_a0;
317 }
318 //-----------------------------------------------------------------------------
319 template <typename TInteger>
320 inline
321 typename DGtal::IntegerComputer<TInteger>::Integer
322 DGtal::IntegerComputer<TInteger>::
323 getCFrac( std::vector<Integer> & quotients,
324  IntegerParamType a, IntegerParamType b ) const
325 {
326  ASSERT( isPositiveOrZero( a ) && isPositiveOrZero( b ) );
327  _m_a0 = a;
328  _m_a1 = b;
329  while ( isNotZero( _m_a1 ) )
330  {
331  getEuclideanDiv( _m_q, _m_r, _m_a0, _m_a1 );
332  quotients.push_back( _m_q );
333  _m_a0 = _m_a1;
334  _m_a1 = _m_r;
335  }
336  return _m_a0;
337 }
338 //-----------------------------------------------------------------------------
339 template <typename TInteger>
340 template <typename OutputIterator>
341 inline
342 typename DGtal::IntegerComputer<TInteger>::Integer
343 DGtal::IntegerComputer<TInteger>::
344 getCFrac( OutputIterator outIt,
345  IntegerParamType a, IntegerParamType b ) const
346 {
347  BOOST_CONCEPT_ASSERT(( boost::OutputIterator< OutputIterator, Integer > ));
348  ASSERT( isPositiveOrZero( a ) && isPositiveOrZero( b ) );
349  _m_a0 = a;
350  _m_a1 = b;
351  while ( isNotZero( _m_a1 ) )
352  {
353  getEuclideanDiv( _m_q, _m_r, _m_a0, _m_a1 );
354  *outIt++ = _m_q;
355  _m_a0 = _m_a1;
356  _m_a1 = _m_r;
357  }
358  return _m_a0;
359 }
360 //-----------------------------------------------------------------------------
361 template <typename TInteger>
362 inline
363 typename DGtal::IntegerComputer<TInteger>::Point2I
364 DGtal::IntegerComputer<TInteger>::
365 convergent( const std::vector<Integer> & quotients,
366  unsigned int k ) const
367 {
368  _m_bezout[ 0 ].clear(); // p
369  _m_bezout[ 0 ].push_back( NumberTraits<Integer>::ZERO );
370  _m_bezout[ 0 ].push_back( NumberTraits<Integer>::ONE );
371  _m_bezout[ 1 ].clear(); // q
372  _m_bezout[ 1 ].push_back( NumberTraits<Integer>::ONE );
373  _m_bezout[ 1 ].push_back( NumberTraits<Integer>::ZERO );
374  if ( k >= quotients.size() )
375  k = (quotients.size() - 1);
376  for ( unsigned int i = 0; i <= k; ++i )
377  {
378  _m_bezout[ 0 ].push_back( quotients[ i ] * _m_bezout[ 0 ][ i + 1 ]
379  + _m_bezout[ 0 ][ i ] );
380  _m_bezout[ 1 ].push_back( quotients[ i ] * _m_bezout[ 1 ][ i + 1 ]
381  + _m_bezout[ 1 ][ i ] );
382  }
383  _m_v[ 0 ] = _m_bezout[ 0 ].back();
384  _m_v[ 1 ] = _m_bezout[ 1 ].back();
385  return _m_v;
386 }
387 //-----------------------------------------------------------------------------
388 // ----------------------- Point2I services ------------------------------
389 //-----------------------------------------------------------------------------
390 template <typename TInteger>
391 inline
392 void
393 DGtal::IntegerComputer<TInteger>::
394 reduce( Vector2I & p ) const
395 {
396  _m_a = gcd( p[ 0 ], p[ 1 ] );
397  if ( ( _m_a != NumberTraits<Integer>::ONE ) && ( isNotZero( _m_a ) ) )
398  p /= _m_a;
399 }
400 //-----------------------------------------------------------------------------
401 template <typename TInteger>
402 inline
403 typename DGtal::IntegerComputer<TInteger>::Integer
404 DGtal::IntegerComputer<TInteger>::
405 crossProduct( const Vector2I & u, const Vector2I & v) const
406 {
407  _m_a0 = u[ 0 ] * v[ 1 ];
408  _m_a1 = u[ 1 ] * v[ 0 ];
409  return _m_a0 - _m_a1;
410 }
411 //-----------------------------------------------------------------------------
412 template <typename TInteger>
413 inline
414 void
415 DGtal::IntegerComputer<TInteger>::
416 getCrossProduct( Integer & cp,
417  const Vector2I & u, const Vector2I & v) const
418 {
419  cp = u[ 0 ] * v[ 1 ];
420  cp -= u[ 1 ] * v[ 0 ];
421 }
422 //-----------------------------------------------------------------------------
423 template <typename TInteger>
424 inline
425 typename DGtal::IntegerComputer<TInteger>::Integer
426 DGtal::IntegerComputer<TInteger>::
427 dotProduct( const Vector2I & u, const Vector2I & v ) const
428 {
429  _m_a0 = u[ 0 ] * v[ 0 ];
430  _m_a1 = u[ 1 ] * v[ 1 ];
431  return _m_a0 + _m_a1;
432 }
433 //-----------------------------------------------------------------------------
434 template <typename TInteger>
435 inline
436 void
437 DGtal::IntegerComputer<TInteger>::
438 getDotProduct( Integer & dp,
439  const Vector2I & u, const Vector2I & v) const
440 {
441  dp = u[ 0 ] * v[ 0 ];
442  dp += u[ 1 ] * v[ 1 ];
443 }
444 //-----------------------------------------------------------------------------
445 template <typename TInteger>
446 typename DGtal::IntegerComputer<TInteger>::Vector2I
447 DGtal::IntegerComputer<TInteger>::
448 extendedEuclid( IntegerParamType a, IntegerParamType b,
449  IntegerParamType c ) const
450 {
451  if( isZero( a ) ) return Point2I( NumberTraits<Integer>::ZERO, b * c );
452  if( isZero( b ) ) return Point2I( a * c, NumberTraits<Integer>::ZERO );
453 
454  for ( unsigned int i = 0; i < 4; ++i )
455  _m_bezout[ i ].clear();
456 
457  _m_bezout[ 0 ].push_back( abs( a ) );
458  _m_bezout[ 0 ].push_back( abs( b ) );
459  _m_bezout[ 1 ].push_back( NumberTraits<Integer>::ZERO );
460  _m_bezout[ 2 ].push_back( NumberTraits<Integer>::ONE );
461  _m_bezout[ 2 ].push_back( NumberTraits<Integer>::ZERO );
462  _m_bezout[ 3 ].push_back( NumberTraits<Integer>::ZERO );
463  _m_bezout[ 3 ].push_back( NumberTraits<Integer>::ONE );
464 
465  unsigned int k = 0; // index of the iteration during the computation.
466  while( isNotZero( _m_bezout[ 0 ][ k+1 ] ) )
467  {
468  _m_bezout[ 1 ].push_back( _m_bezout[ 0 ][ k ]
469  / _m_bezout[ 0 ][ k+1 ] );
470  _m_bezout[ 0 ].push_back( _m_bezout[ 0 ][ k ]
471  % _m_bezout[ 0 ][ k+1 ] );
472  _m_bezout[ 2 ].push_back( _m_bezout[ 2 ][ k ]
473  - _m_bezout[ 1 ][ k+1 ]
474  * _m_bezout[ 2 ][ k+1 ] );
475  _m_bezout[ 3 ].push_back( _m_bezout[ 3 ][ k ]
476  - _m_bezout[ 1 ][ k+1 ]
477  *_m_bezout[ 3 ][ k+1 ] );
478  ++k;
479  }
480 
481  _m_v[ 0 ] = _m_bezout[ 2 ][ k ];
482  _m_v[ 1 ] = _m_bezout[ 3 ][ k ];
483 
484  if ( isNegative( a ) ) _m_v[ 0 ] = -_m_v[ 0 ];
485  if ( isNegative( b ) ) _m_v[ 1 ] = -_m_v[ 1 ];
486  //_m_v *= c;
487  return _m_v;
488 }
489 //-----------------------------------------------------------------------------
490 template <typename TInteger>
491 inline
492 void
493 DGtal::IntegerComputer<TInteger>::
494 getCoefficientIntersection( Integer & fl, Integer & ce,
495  const Vector2I & p,
496  const Vector2I & u,
497  const Vector2I & N,
498  IntegerParamType c ) const
499 {
500  getDotProduct( _m_c0, p, N );
501  getDotProduct( _m_c1, u, N );
502  _m_c2 = c - _m_c0;
503  fl = floorDiv( _m_c2, _m_c1 );
504  ce = ceilDiv( _m_c2, _m_c1 );
505  // getFloorCeilDiv( fl, ce, _m_c2, _m_c1 );
506 }
507 //-----------------------------------------------------------------------------
508 template <typename TInteger>
509 inline
510 void
511 DGtal::IntegerComputer<TInteger>::
512 getValidBezout ( Vector2I & v,
513  const Point2I & A, const Vector2I & u,
514  const Vector2I & N, IntegerParamType c,
515  const Vector2I & N2, IntegerParamType c2,
516  bool compute_v ) const
517 {
518  if ( compute_v )
519  {
520  v = extendedEuclid( -u[ 1 ], u[ 0 ], NumberTraits<Integer>::ONE );
521  _m_v0 = A + v;
522  getDotProduct( _m_c0, _m_v0, N );
523  if ( _m_c0 > c )
524  {
525  v[ 0 ] = -v[ 0 ];
526  v[ 1 ] = -v[ 1 ];
527  _m_v0 = A + v;
528  }
529  }
530  else _m_v0 = A + v;
531  getCoefficientIntersection( _m_a0, _m_a1,
532  _m_v0, u, N2, c2 );
533  _m_v1 = u * _m_a0; // floor value
534  v += _m_v1;
535  ASSERT( N2.dot( A + v ) <= c2 );
536  ASSERT( N2.dot( A + v + u ) > c2 );
537 }
538 //-----------------------------------------------------------------------------
539 template <typename TInteger>
540 inline
541 void
542 DGtal::IntegerComputer<TInteger>::
543 reduce( Vector3I & p ) const
544 {
545  _m_a = gcd( p[ 0 ], p[ 1 ] );
546  _m_r = gcd( _m_a, p[ 2 ] );
547  p /= _m_r;
548 }
549 //-----------------------------------------------------------------------------
550 template <typename TInteger>
551 inline
552 typename DGtal::IntegerComputer<TInteger>::Integer
553 DGtal::IntegerComputer<TInteger>::
554 dotProduct( const Vector3I & u, const Vector3I & v) const
555 {
556  _m_a0 = u[ 0 ] * v[ 0 ];
557  _m_a0 += u[ 1 ] * v[ 1 ];
558  _m_a0 += u[ 2 ] * v[ 2 ];
559  return _m_a0;
560 }
561 //-----------------------------------------------------------------------------
562 template <typename TInteger>
563 inline
564 void
565 DGtal::IntegerComputer<TInteger>::
566 getDotProduct( Integer & dp,
567  const Vector3I & u, const Vector3I & v) const
568 {
569  dp = u[ 0 ] * v[ 0 ];
570  dp += u[ 1 ] * v[ 1 ];
571  dp += u[ 2 ] * v[ 2 ];
572 }
573 
574 
576 // Interface - public :
577 
582 template <typename TInteger>
583 inline
584 void
585 DGtal::IntegerComputer<TInteger>::selfDisplay ( std::ostream & out ) const
586 {
587  out << "[IntegerComputer]";
588 }
589 
594 template <typename TInteger>
595 inline
596 bool
597 DGtal::IntegerComputer<TInteger>::isValid() const
598 {
599  return true;
600 }
601 
602 
603 
605 // Implementation of inline functions //
606 
607 template <typename TInteger>
608 inline
609 std::ostream&
610 DGtal::operator<< ( std::ostream & out,
611  const IntegerComputer<TInteger> & object )
612 {
613  object.selfDisplay( out );
614  return out;
615 }
616 
617 // //
619 
620 
Generated on Wed Dec 19 2012 19:10:27 for DGtal by   doxygen 1.8.1.1