DGtal  1.5.beta
IntegerComputer.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file IntegerComputer.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2012/03/05
23  *
24  * Implementation of inline methods defined in IntegerComputer.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 ///////////////////////////////////////////////////////////////////////////////
35 // IMPLEMENTATION of inline methods.
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
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  if ( k % 2 != 0 )
484  { // odd case
485  _m_v[ 0 ] = abs( b ) + _m_bezout[ 2 ][ k ];
486  _m_v[ 1 ] = -abs( a ) + _m_bezout[ 3 ][ k ];
487  }
488  // choose sgn(a) = sgn(x) when x != 0, iff c > 0
489  // |x| <= |bc|, |y| < |ac|
490  _m_v *= c / _m_bezout[ 0 ][ k ]; // c / gcd(a,b)
491  if ( isNegative( a ) ) _m_v[ 0 ] = - _m_v[ 0 ];
492  if ( isNegative( b ) ) _m_v[ 1 ] = - _m_v[ 1 ];
493  // std::cout << "a * x + b * y == c, "
494  // << a << " * " << _m_v[ 0 ] << " + " << b << " * " << _m_v[ 1 ] << " == " << c << std::endl;
495  ASSERT( (a*_m_v[ 0 ]+b*_m_v[ 1 ]) == c );
496  ASSERT( isNegative( c ) || ( ( isPositive( a ) == isPositive( _m_v[ 0 ] ) ) || isZero( _m_v[ 0 ] ) ) );
497  ASSERT( isPositive( c ) || ( ( isPositive( a ) == isNegative( _m_v[ 0 ] ) ) || isZero( _m_v[ 0 ] ) ) );
498  ASSERT( abs( _m_v[ 0 ] ) <= abs( b*c ) );
499  ASSERT( abs( _m_v[ 1 ] ) < abs( a*c ) );
500  return _m_v;
501 }
502 //-----------------------------------------------------------------------------
503 template <typename TInteger>
504 inline
505 void
506 DGtal::IntegerComputer<TInteger>::
507 getCoefficientIntersection( Integer & fl, Integer & ce,
508  const Vector2I & p,
509  const Vector2I & u,
510  const Vector2I & N,
511  IntegerParamType c ) const
512 {
513  getDotProduct( _m_c0, p, N );
514  getDotProduct( _m_c1, u, N );
515  _m_c2 = c - _m_c0;
516  fl = floorDiv( _m_c2, _m_c1 );
517  ce = ceilDiv( _m_c2, _m_c1 );
518  // getFloorCeilDiv( fl, ce, _m_c2, _m_c1 );
519 }
520 //-----------------------------------------------------------------------------
521 template <typename TInteger>
522 inline
523 void
524 DGtal::IntegerComputer<TInteger>::
525 getValidBezout ( Vector2I & v,
526  const Point2I & A, const Vector2I & u,
527  const Vector2I & N, IntegerParamType c,
528  const Vector2I & N2, IntegerParamType c2,
529  bool compute_v ) const
530 {
531  if ( compute_v )
532  {
533  v = extendedEuclid( -u[ 1 ], u[ 0 ], NumberTraits<Integer>::ONE );
534  _m_v0 = A + v;
535  getDotProduct( _m_c0, _m_v0, N );
536  if ( _m_c0 > c )
537  {
538  v[ 0 ] = -v[ 0 ];
539  v[ 1 ] = -v[ 1 ];
540  _m_v0 = A + v;
541  }
542  }
543  else _m_v0 = A + v;
544  getCoefficientIntersection( _m_a0, _m_a1,
545  _m_v0, u, N2, c2 );
546  _m_v1 = u * _m_a0; // floor value
547  v += _m_v1;
548  ASSERT( N2.dot( A + v ) <= c2 );
549  ASSERT( N2.dot( A + v + u ) > c2 );
550 }
551 //-----------------------------------------------------------------------------
552 template <typename TInteger>
553 inline
554 void
555 DGtal::IntegerComputer<TInteger>::
556 reduce( Vector3I & p ) const
557 {
558  _m_a = gcd( p[ 0 ], p[ 1 ] );
559  _m_r = gcd( _m_a, p[ 2 ] );
560  p /= _m_r;
561 }
562 //-----------------------------------------------------------------------------
563 template <typename TInteger>
564 inline
565 typename DGtal::IntegerComputer<TInteger>::Integer
566 DGtal::IntegerComputer<TInteger>::
567 dotProduct( const Vector3I & u, const Vector3I & v) const
568 {
569  _m_a0 = u[ 0 ] * v[ 0 ];
570  _m_a0 += u[ 1 ] * v[ 1 ];
571  _m_a0 += u[ 2 ] * v[ 2 ];
572  return _m_a0;
573 }
574 //-----------------------------------------------------------------------------
575 template <typename TInteger>
576 inline
577 void
578 DGtal::IntegerComputer<TInteger>::
579 getDotProduct( Integer & dp,
580  const Vector3I & u, const Vector3I & v) const
581 {
582  dp = u[ 0 ] * v[ 0 ];
583  dp += u[ 1 ] * v[ 1 ];
584  dp += u[ 2 ] * v[ 2 ];
585 }
586 
587 
588 ///////////////////////////////////////////////////////////////////////////////
589 // Interface - public :
590 
591 /**
592  * Writes/Displays the object on an output stream.
593  * @param out the output stream where the object is written.
594  */
595 template <typename TInteger>
596 inline
597 void
598 DGtal::IntegerComputer<TInteger>::selfDisplay ( std::ostream & out ) const
599 {
600  out << "[IntegerComputer]";
601 }
602 
603 /**
604  * Checks the validity/consistency of the object.
605  * @return 'true' if the object is valid, 'false' otherwise.
606  */
607 template <typename TInteger>
608 inline
609 bool
610 DGtal::IntegerComputer<TInteger>::isValid() const
611 {
612  return true;
613 }
614 
615 
616 
617 ///////////////////////////////////////////////////////////////////////////////
618 // Implementation of inline functions //
619 
620 template <typename TInteger>
621 inline
622 std::ostream&
623 DGtal::operator<< ( std::ostream & out,
624  const IntegerComputer<TInteger> & object )
625 {
626  object.selfDisplay( out );
627  return out;
628 }
629 
630 // //
631 ///////////////////////////////////////////////////////////////////////////////
632 
633