DGtal 1.4.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
42template <typename TInteger>
43inline
44DGtal::IntegerComputer<TInteger>::~IntegerComputer()
45{}
46//-----------------------------------------------------------------------------
47template <typename TInteger>
48inline
49DGtal::IntegerComputer<TInteger>::IntegerComputer()
50{}
51//-----------------------------------------------------------------------------
52template <typename TInteger>
53inline
54DGtal::IntegerComputer<TInteger>::IntegerComputer( const Self & /*other*/ )
55{}
56//-----------------------------------------------------------------------------
57template <typename TInteger>
58inline
59DGtal::IntegerComputer<TInteger> &
60DGtal::IntegerComputer<TInteger>::operator=( const Self & /*other*/ )
61{
62 return *this;
63}
64//-----------------------------------------------------------------------------
65// ----------------------- Integer services ------------------------------
66//-----------------------------------------------------------------------------
67template <typename TInteger>
68inline
69bool
70DGtal::IntegerComputer<TInteger>::
71isZero( IntegerParamType a )
72{
73 return a == NumberTraits<Integer>::ZERO;
74}
75//-----------------------------------------------------------------------------
76template <typename TInteger>
77inline
78bool
79DGtal::IntegerComputer<TInteger>::
80isNotZero( IntegerParamType a )
81{
82 return a != NumberTraits<Integer>::ZERO;
83}
84//-----------------------------------------------------------------------------
85template <typename TInteger>
86inline
87bool
88DGtal::IntegerComputer<TInteger>::
89isPositive( IntegerParamType a )
90{
91 return a > NumberTraits<Integer>::ZERO;
92}
93//-----------------------------------------------------------------------------
94template <typename TInteger>
95inline
96bool
97DGtal::IntegerComputer<TInteger>::
98isNegative( IntegerParamType a )
99{
100 return a < NumberTraits<Integer>::ZERO;
101}
102//-----------------------------------------------------------------------------
103template <typename TInteger>
104inline
105bool
106DGtal::IntegerComputer<TInteger>::
107isPositiveOrZero( IntegerParamType a )
108{
109 return a >= NumberTraits<Integer>::ZERO;
110}
111//-----------------------------------------------------------------------------
112template <typename TInteger>
113inline
114bool
115DGtal::IntegerComputer<TInteger>::
116isNegativeOrZero( IntegerParamType a )
117{
118 return a <= NumberTraits<Integer>::ZERO;
119}
120//-----------------------------------------------------------------------------
121template <typename TInteger>
122inline
123typename DGtal::IntegerComputer<TInteger>::Integer
124DGtal::IntegerComputer<TInteger>::
125abs( IntegerParamType a )
126{
127 if ( isPositiveOrZero( a ) )
128 return a;
129 else
130 return -a;
131}
132//-----------------------------------------------------------------------------
133template <typename TInteger>
134inline
135typename DGtal::IntegerComputer<TInteger>::Integer
136DGtal::IntegerComputer<TInteger>::
137max( IntegerParamType a, IntegerParamType b )
138{
139 return ( a >= b ) ? a : b;
140}
141//-----------------------------------------------------------------------------
142template <typename TInteger>
143inline
144typename DGtal::IntegerComputer<TInteger>::Integer
145DGtal::IntegerComputer<TInteger>::
146max( IntegerParamType a, IntegerParamType b, IntegerParamType c )
147{
148 return max( max(a,b), c );
149}
150//-----------------------------------------------------------------------------
151template <typename TInteger>
152inline
153typename DGtal::IntegerComputer<TInteger>::Integer
154DGtal::IntegerComputer<TInteger>::
155min( IntegerParamType a, IntegerParamType b )
156{
157 return ( a <= b ) ? a : b;
158}
159//-----------------------------------------------------------------------------
160template <typename TInteger>
161inline
162typename DGtal::IntegerComputer<TInteger>::Integer
163DGtal::IntegerComputer<TInteger>::
164min( IntegerParamType a, IntegerParamType b, IntegerParamType c )
165{
166 return min( min(a,b), c );
167}
168//-----------------------------------------------------------------------------
169template <typename TInteger>
170inline
171void
172DGtal::IntegerComputer<TInteger>::
173getEuclideanDiv( Integer & q, Integer & r,
174 IntegerParamType a, IntegerParamType b ) const
175{
176 q = a / b;
177 r = a % b;
178}
179//-----------------------------------------------------------------------------
180template <typename TInteger>
181inline
182typename DGtal::IntegerComputer<TInteger>::Integer
183DGtal::IntegerComputer<TInteger>::
184floorDiv( 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//-----------------------------------------------------------------------------
209template <typename TInteger>
210inline
211typename DGtal::IntegerComputer<TInteger>::Integer
212DGtal::IntegerComputer<TInteger>::
213ceilDiv( 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//-----------------------------------------------------------------------------
238template <typename TInteger>
239inline
240void
241DGtal::IntegerComputer<TInteger>::
242getFloorCeilDiv( 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//-----------------------------------------------------------------------------
260template <typename TInteger>
261inline
262typename DGtal::IntegerComputer<TInteger>::Integer
263DGtal::IntegerComputer<TInteger>::
264staticGcd( 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//-----------------------------------------------------------------------------
280template <typename TInteger>
281inline
282typename DGtal::IntegerComputer<TInteger>::Integer
283DGtal::IntegerComputer<TInteger>::
284gcd( 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//-----------------------------------------------------------------------------
299template <typename TInteger>
300inline
301void
302DGtal::IntegerComputer<TInteger>::
303getGcd( 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//-----------------------------------------------------------------------------
319template <typename TInteger>
320inline
321typename DGtal::IntegerComputer<TInteger>::Integer
322DGtal::IntegerComputer<TInteger>::
323getCFrac( 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//-----------------------------------------------------------------------------
339template <typename TInteger>
340template <typename OutputIterator>
341inline
342typename DGtal::IntegerComputer<TInteger>::Integer
343DGtal::IntegerComputer<TInteger>::
344getCFrac( 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//-----------------------------------------------------------------------------
361template <typename TInteger>
362inline
363typename DGtal::IntegerComputer<TInteger>::Point2I
364DGtal::IntegerComputer<TInteger>::
365convergent( 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//-----------------------------------------------------------------------------
390template <typename TInteger>
391inline
392void
393DGtal::IntegerComputer<TInteger>::
394reduce( 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//-----------------------------------------------------------------------------
401template <typename TInteger>
402inline
403typename DGtal::IntegerComputer<TInteger>::Integer
404DGtal::IntegerComputer<TInteger>::
405crossProduct( 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//-----------------------------------------------------------------------------
412template <typename TInteger>
413inline
414void
415DGtal::IntegerComputer<TInteger>::
416getCrossProduct( 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//-----------------------------------------------------------------------------
423template <typename TInteger>
424inline
425typename DGtal::IntegerComputer<TInteger>::Integer
426DGtal::IntegerComputer<TInteger>::
427dotProduct( 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//-----------------------------------------------------------------------------
434template <typename TInteger>
435inline
436void
437DGtal::IntegerComputer<TInteger>::
438getDotProduct( 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//-----------------------------------------------------------------------------
445template <typename TInteger>
446typename DGtal::IntegerComputer<TInteger>::Vector2I
447DGtal::IntegerComputer<TInteger>::
448extendedEuclid( 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//-----------------------------------------------------------------------------
503template <typename TInteger>
504inline
505void
506DGtal::IntegerComputer<TInteger>::
507getCoefficientIntersection( 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//-----------------------------------------------------------------------------
521template <typename TInteger>
522inline
523void
524DGtal::IntegerComputer<TInteger>::
525getValidBezout ( 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//-----------------------------------------------------------------------------
552template <typename TInteger>
553inline
554void
555DGtal::IntegerComputer<TInteger>::
556reduce( 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//-----------------------------------------------------------------------------
563template <typename TInteger>
564inline
565typename DGtal::IntegerComputer<TInteger>::Integer
566DGtal::IntegerComputer<TInteger>::
567dotProduct( 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//-----------------------------------------------------------------------------
575template <typename TInteger>
576inline
577void
578DGtal::IntegerComputer<TInteger>::
579getDotProduct( 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 */
595template <typename TInteger>
596inline
597void
598DGtal::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 */
607template <typename TInteger>
608inline
609bool
610DGtal::IntegerComputer<TInteger>::isValid() const
611{
612 return true;
613}
614
615
616
617///////////////////////////////////////////////////////////////////////////////
618// Implementation of inline functions //
619
620template <typename TInteger>
621inline
622std::ostream&
623DGtal::operator<< ( std::ostream & out,
624 const IntegerComputer<TInteger> & object )
625{
626 object.selfDisplay( out );
627 return out;
628}
629
630// //
631///////////////////////////////////////////////////////////////////////////////
632
633