DGtal  0.9.2
SternBrocot.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 SternBrocot.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  * @author Xavier Proven├žal (\c xavier.provencal@univ-savoie.fr )
22  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
23  *
24  * @date 2012/03/07
25  *
26  * Implementation of inline methods defined in SternBrocot.h
27  *
28  * This file is part of the DGtal library.
29  */
30 
31 
32 //////////////////////////////////////////////////////////////////////////////
33 #include <cstdlib>
34 #include "DGtal/arithmetic/IntegerComputer.h"
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
43 
44 ///////////////////////////////////////////////////////////////////////////////
45 // DGtal::SternBrocot<TInteger, TQuotient>::Node
46 //-----------------------------------------------------------------------------
47 template <typename TInteger, typename TQuotient>
48 inline
49 DGtal::SternBrocot<TInteger, TQuotient>::Node::
50 Node( Integer p1, Integer q1, Quotient u1, Quotient k1,
51  Node* ascendant_left1, Node* ascendant_right1,
52  Node* descendant_left1, Node* descendant_right1,
53  Node* inverse1 )
54  : p( p1 ), q( q1 ), u( u1 ), k( k1 ),
55  ascendantLeft( ascendant_left1 ),
56  ascendantRight( ascendant_right1 ),
57  descendantLeft( descendant_left1 ),
58  descendantRight( descendant_right1 ),
59  inverse( inverse1 )
60 {
61  // std::cerr << "(" << p1 << "/" << q1 << "," << u1 << "," << k1 << ")";
62 }
63 
64 //////////////////////////////////////////////////////////////////////////////
65 // DGtal::SternBrocot<TInteger, TQuotient>::Fraction
66 //-----------------------------------------------------------------------------
67 template <typename TInteger, typename TQuotient>
68 inline
69 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
70 Fraction( Integer aP, Integer aQ, Fraction ancestor )
71 {
72  this->operator=( SternBrocotTree::fraction( aP, aQ, ancestor ) );
73 }
74 //-----------------------------------------------------------------------------
75 template <typename TInteger, typename TQuotient>
76 inline
77 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
78 Fraction( Node* sb_node )
79  : myNode( sb_node )
80 {
81 }
82 //-----------------------------------------------------------------------------
83 template <typename TInteger, typename TQuotient>
84 inline
85 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
86 Fraction( const Self & other )
87  : myNode( other.myNode )
88 {
89 }
90 //-----------------------------------------------------------------------------
91 template <typename TInteger, typename TQuotient>
92 inline
93 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction &
94 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
95 operator=( const Self & other )
96 {
97  if ( this != &other )
98  {
99  myNode = other.myNode;
100  }
101  return *this;
102 }
103 //-----------------------------------------------------------------------------
104 template <typename TInteger, typename TQuotient>
105 inline
106 bool
107 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
108 null() const
109 {
110  return myNode == 0;
111 }
112 //-----------------------------------------------------------------------------
113 template <typename TInteger, typename TQuotient>
114 inline
115 typename DGtal::SternBrocot<TInteger, TQuotient>::Integer
116 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
117 p() const
118 {
119  return myNode ? myNode->p : 0;
120 }
121 //-----------------------------------------------------------------------------
122 template <typename TInteger, typename TQuotient>
123 inline
124 typename DGtal::SternBrocot<TInteger, TQuotient>::Integer
125 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
126 q() const
127 {
128  return myNode ? myNode->q : 0;
129 }
130 //-----------------------------------------------------------------------------
131 template <typename TInteger, typename TQuotient>
132 inline
133 typename DGtal::SternBrocot<TInteger, TQuotient>::Quotient
134 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
135 u() const
136 {
137  ASSERT( myNode != 0 );
138  return myNode->u;
139 }
140 //-----------------------------------------------------------------------------
141 template <typename TInteger, typename TQuotient>
142 inline
143 typename DGtal::SternBrocot<TInteger, TQuotient>::Quotient
144 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
145 k() const
146 {
147  ASSERT( myNode != 0 );
148  return myNode->k;
149 }
150 //-----------------------------------------------------------------------------
151 template <typename TInteger, typename TQuotient>
152 inline
153 bool
154 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
155 equals( Integer p1, Integer q1 ) const
156 {
157  return ( this->p() == p1 ) && ( this->q() == q1 );
158 }
159 //-----------------------------------------------------------------------------
160 template <typename TInteger, typename TQuotient>
161 inline
162 bool
163 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
164 lessThan( Integer p1, Integer q1 ) const
165 {
166  Integer d = p() * q1 - q() * p1;
167  return d < NumberTraits<Integer>::ZERO;
168 }
169 //-----------------------------------------------------------------------------
170 template <typename TInteger, typename TQuotient>
171 inline
172 bool
173 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
174 moreThan( Integer p1, Integer q1 ) const
175 {
176  Integer d = p() * q1 - q() * p1;
177  return d > NumberTraits<Integer>::ZERO;
178 }
179 //-----------------------------------------------------------------------------
180 template <typename TInteger, typename TQuotient>
181 inline
182 bool
183 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
184 operator==( const Fraction & other ) const
185 {
186  return myNode == other.myNode;
187 }
188 //-----------------------------------------------------------------------------
189 template <typename TInteger, typename TQuotient>
190 inline
191 bool
192 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
193 operator!=( const Fraction & other ) const
194 {
195  return myNode != other.myNode;
196 }
197 //-----------------------------------------------------------------------------
198 template <typename TInteger, typename TQuotient>
199 inline
200 bool
201 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
202 operator<( const Fraction & other ) const
203 {
204  return this->lessThan( other.p(), other.q() );
205 }
206 //-----------------------------------------------------------------------------
207 template <typename TInteger, typename TQuotient>
208 inline
209 bool
210 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
211 operator>( const Fraction & other ) const
212 {
213  return this->moreThan( other.p(), other.q() );
214 }
215 //-----------------------------------------------------------------------------
216 template <typename TInteger, typename TQuotient>
217 inline
218 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
219 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
220 left() const
221 {
222  if ( myNode->descendantLeft == 0 )
223  {
224  Node* pleft = myNode->ascendantLeft;
225  Node* n = new Node( p() + pleft->p,
226  q() + pleft->q,
227  odd() ? u() + 1 : (Quotient) 2,
228  odd() ? k() : k() + 1,
229  pleft, myNode,
230  0, 0, 0 );
231  Fraction inv = Fraction( myNode->inverse );
232  Node* invpright = inv.myNode->ascendantRight;
233  Node* invn = new Node( inv.p() + invpright->p,
234  inv.q() + invpright->q,
235  inv.even() ? inv.u() + 1 : (Quotient) 2,
236  inv.even() ? inv.k() : inv.k() + 1,
237  myNode->inverse, invpright,
238  0, 0, n );
239  n->inverse = invn;
240  myNode->inverse->descendantRight = invn;
241  myNode->descendantLeft = n;
242  instance().nbFractions += 2;
243  }
244  return Fraction( myNode->descendantLeft );
245 }
246 //-----------------------------------------------------------------------------
247 template <typename TInteger, typename TQuotient>
248 inline
249 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
250 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
251 right() const
252 {
253  if ( myNode->descendantRight == 0 )
254  {
255  Fraction inv( myNode->inverse );
256  inv.left();
257  ASSERT( myNode->descendantRight != 0 );
258  }
259  return Fraction( myNode->descendantRight );
260 }
261 //-----------------------------------------------------------------------------
262 template <typename TInteger, typename TQuotient>
263 inline
264 bool
265 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
266 even() const
267 {
268  return ( k() & 1 ) == 0;
269 }
270 //-----------------------------------------------------------------------------
271 template <typename TInteger, typename TQuotient>
272 inline
273 bool
274 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
275 odd() const
276 {
277  return ( k() & 1 ) != 0;
278 }
279 //-----------------------------------------------------------------------------
280 template <typename TInteger, typename TQuotient>
281 inline
282 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
283 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
284 father() const
285 {
286  return Fraction( odd() ? myNode->ascendantRight : myNode->ascendantLeft );
287 }
288 //-----------------------------------------------------------------------------
289 template <typename TInteger, typename TQuotient>
290 inline
291 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
292 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
293 father( Quotient m ) const
294 {
295  if ( m > NumberTraits<Quotient>::ONE ) // > 1
296  {
297  Node* n = myNode;
298  while ( n->u > m )
299  n = odd() ? n->ascendantRight : n->ascendantLeft;
300  return Fraction( n );
301  }
302  else if ( m != NumberTraits<Quotient>::ZERO ) // == 1
303  {
304  return odd() ? previousPartial().right() : previousPartial().left();
305  }
306  else // == 0
307  return reduced( 2 ); //previousPartial().previousPartial();
308 }
309 //-----------------------------------------------------------------------------
310 template <typename TInteger, typename TQuotient>
311 inline
312 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
313 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
314 previousPartial() const
315 {
316  return Fraction( odd() ? myNode->ascendantLeft : myNode->ascendantRight );
317  // return Fraction( odd()
318  // ? myNode->ascendantLeft->descendantRight
319  // : myNode->ascendantRight->descendantLeft );
320 }
321 //-----------------------------------------------------------------------------
322 template <typename TInteger, typename TQuotient>
323 inline
324 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
325 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
326 inverse() const
327 {
328  return Fraction( myNode->inverse );
329 }
330 //-----------------------------------------------------------------------------
331 template <typename TInteger, typename TQuotient>
332 inline
333 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
334 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
335 partial( Quotient kp ) const
336 {
337  ASSERT( ( ((Quotient)-2) <= kp ) && ( kp <= k() ) );
338  return reduced( k() - kp );
339 }
340 //-----------------------------------------------------------------------------
341 template <typename TInteger, typename TQuotient>
342 inline
343 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
344 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
345 reduced( Quotient i ) const
346 {
347  ASSERT( ( ((Quotient)0) <= i ) && ( i <= ( k()+((Quotient)2) ) ) );
348  Node* n = this->myNode;
349 
350  bool bleft = ( n->k & NumberTraits<Quotient>::ONE )
351  != NumberTraits<Quotient>::ZERO;
352  while ( i-- > NumberTraits<Quotient>::ZERO )
353  {
354  n = bleft ? n->ascendantLeft : n->ascendantRight;
355  bleft = ! bleft;
356  }
357  return Fraction( n );
358 }
359 //-----------------------------------------------------------------------------
360 template <typename TInteger, typename TQuotient>
361 inline
362 void
363 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
364 push_back( const std::pair<Quotient, Quotient> & quotient )
365 {
366  pushBack( quotient );
367 }
368 //-----------------------------------------------------------------------------
369 template <typename TInteger, typename TQuotient>
370 inline
371 void
372 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
373 pushBack( const std::pair<Quotient, Quotient> & quotient )
374 {
375  // std::vector<Quotient> quots;
376  // if ( ! null() )
377  // {
378  // this->getCFrac( quots );
379  // std::cerr << "F[";
380  // for ( unsigned int i = 0; i < quots.size(); ++i )
381  // std::cerr << " " << quots[ i ];
382  // }
383  // else std::cerr << "[";
384  // std::cerr << "] + " << "(" << quotient.first
385  // << "," << quotient.second << ")";
386  if ( null() )
387  {
388  ASSERT( quotient.second <= NumberTraits<Quotient>::ZERO );
389  if ( quotient.second < NumberTraits<Quotient>::ZERO )
390  this->operator=( oneOverZero() );
391  else if ( quotient.first == NumberTraits<Quotient>::ZERO ) // (0,0)
392  this->operator=( zeroOverOne() );
393  else
394  {
395  Fraction f = zeroOverOne();
396  for ( Quotient i = 0; i < quotient.first; ++i )
397  f = f.right();
398  this->operator=( f );
399  }
400  }
401  else if ( NumberTraits<Quotient>::even( quotient.second ) )
402  {
403  Fraction f = left();
404  for ( Quotient i = 1; i < quotient.first; ++i )
405  f = f.right();
406  this->operator=( f );
407  }
408  else
409  {
410  Fraction f = right();
411  for ( Quotient i = 1; i < quotient.first; ++i )
412  f = f.left();
413  this->operator=( f );
414  }
415  // quots.clear();
416  // this->getCFrac( quots );
417  // std::cerr << " => F[";
418  // for ( unsigned int i = 0; i < quots.size(); ++i )
419  // std::cerr << " " << quots[ i ];
420  // std::cerr << "]" << std::endl;
421 }
422 //-----------------------------------------------------------------------------
423 template <typename TInteger, typename TQuotient>
424 inline
425 void
426 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
427 getSplit( Fraction & f1, Fraction & f2 ) const
428 {
429  f1.myNode = myNode->ascendantLeft;
430  f2.myNode = myNode->ascendantRight;
431 }
432 //-----------------------------------------------------------------------------
433 template <typename TInteger, typename TQuotient>
434 inline
435 void
436 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
437 getSplitBerstel( Fraction & f1, Quotient & nb1,
438  Fraction & f2, Quotient & nb2 ) const
439 {
440  if ( odd() )
441  {
442  f1.myNode = myNode->ascendantLeft;
443  nb1 = this->u();
444  f2.myNode = f1.myNode->ascendantRight;
445  nb2 = 1;
446  }
447  else
448  {
449  f2.myNode = myNode->ascendantRight;
450  nb2 = this->u();
451  f1.myNode = f2.myNode->ascendantLeft;
452  nb1 = 1;
453  }
454 }
455 //-----------------------------------------------------------------------------
456 template <typename TInteger, typename TQuotient>
457 inline
458 void
459 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
460 getCFrac( std::vector<Quotient> & quotients ) const
461 {
462  ASSERT( this->k() >= NumberTraits<Quotient>::ZERO );
463  int64_t i = NumberTraits<Quotient>::castToInt64_t( this->k() );
464  quotients.resize( (unsigned int)i + 1 );
465  quotients[ (unsigned int)i-- ] = this->u();
466  Node* n = myNode;
467  bool bleft = odd() ? true : false;
468  while ( i >= 0 )
469  {
470  ASSERT( n->k >= NumberTraits<Quotient>::ZERO );
471  n = bleft ? n->ascendantLeft : n->ascendantRight;
472  quotients[ (unsigned int)i ] =
473  ( i == NumberTraits<Quotient>::castToInt64_t( n->k ) ) ? n->u : 1;
474  --i;
475  bleft = ! bleft;
476  }
477 }
478 //-----------------------------------------------------------------------------
479 template <typename TInteger, typename TQuotient>
480 inline
481 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction::ConstIterator
482 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
483 begin() const
484 {
485  CFracSequence* seq = new CFracSequence;
486  this->getCFrac( *seq );
487  return ConstIterator( seq, seq->begin() );
488 }
489 //-----------------------------------------------------------------------------
490 template <typename TInteger, typename TQuotient>
491 inline
492 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction::ConstIterator
493 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
494 end() const
495 {
496  static CFracSequence dummy;
497  return ConstIterator( 0, dummy.end() );
498 }
499 
500 template <typename TInteger, typename TQuotient>
501 inline
502 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
503 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
504 median(const Fraction & g) const
505 {
506  return Fraction(this->p()+g.p(),this->q()+g.q());
507 }
508 
509 //----------------------------------------------------------------------------
510 template <typename TInteger, typename TQuotient>
511 inline
512 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
513 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
514 simplestFractionInBetween(const Fraction & other) const
515 {
516  Fraction f(*this);
517  Fraction g(other);
518  Fraction res;
519 
520  if(f>g)
521  {
522  res = f; f = g; g = res;
523  }
524  res = Fraction();
525 
526  int i = 0;
527 
528  Value uf, ug;
529  ConstIterator itf=f.begin(), itg=g.begin();
530 
531  DGtal::functors::Abs<TInteger> absComputer;
532 
533  if(absComputer(f.p()*g.q()-f.q()*g.p())==NumberTraits<TInteger>::ONE)
534  return f.median(g);
535 
536  itf = f.begin(); itg = g.begin();
537  uf = *itf; ug = *itg;
538  while(uf.first == ug.first && i != f.k() && i != g.k())
539  {
540  res.push_back(std::make_pair(uf.first,i));
541  i++;
542  itf++;itg++;
543  uf = *itf;
544  ug = *itg;
545  }
546 
547  if(uf.first==ug.first)
548  {
549  if(i == f.k())
550  {
551  res.push_back(std::make_pair(uf.first,i));
552  i++;
553  itg++;
554  ug = *itg;
555  res.push_back(std::make_pair(ug.first+1,i));
556  }
557  else
558  {
559  res.push_back(std::make_pair(uf.first,i));
560  i++;
561  itf++;
562  uf = *itf;
563  res.push_back(std::make_pair(uf.first+1,i));
564  }
565  }
566  else
567  {
568  if(i!=f.k() && i != g.k())
569  (uf.first<ug.first)?res.push_back(std::make_pair(uf.first+1,i)):res.push_back(std::make_pair(ug.first+1,i));
570  else
571  if(i == f.k() && i == g.k())
572  (uf.first<ug.first)?res.push_back(std::make_pair(uf.first+1,i)):res.push_back(std::make_pair(ug.first+1,i));
573  else
574  if(i==f.k())
575  {
576  if(uf.first < ug.first)
577  res.push_back(std::make_pair(uf.first+1,i));
578  else
579  if(uf.first == ug.first + 1)
580  {
581  res.push_back(std::make_pair(ug.first,i));
582  i++;
583  itg++;
584  ug = *itg;
585  if(ug.first==NumberTraits<TInteger>::ONE)
586  {
587  res.push_back(std::make_pair(ug.first,i));
588  i++;
589  itg++;
590  ug = *itg;
591  res.push_back(std::make_pair(ug.first+1,i));
592  }
593  else
594  res.push_back(std::make_pair(2,i));
595  }
596  else
597  res.push_back(std::make_pair(ug.first+1,i));
598  }
599  else
600  {
601  if(ug.first < uf.first)
602  res.push_back(std::make_pair(ug.first+1,i));
603  else
604  if(ug.first == uf.first + 1)
605  {
606  res.push_back(std::make_pair(uf.first,i));
607  i++;
608  itf++;
609  uf = *itf;
610  if(uf.first==NumberTraits<TInteger>::ONE)
611  {
612  res.push_back(std::make_pair(uf.first,i));
613  i++;
614  itf++;
615  uf = *itf;
616  res.push_back(std::make_pair(uf.first+1,i));
617  }
618  else
619  res.push_back(std::make_pair(2,i));
620  }
621  else
622  res.push_back(std::make_pair(uf.first+1,i));
623 
624  }
625  }
626 
627  return res;
628 }
629 
630 
631 //-----------------------------------------------------------------------------
632 template <typename TInteger, typename TQuotient>
633 inline
634 void
635 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
636 selfDisplay( std::ostream & out ) const
637 {
638  if ( this->null() ) out << "[Fraction null]";
639  else
640  {
641  out << "[Fraction f=" << this->p()
642  << "/" << this->q()
643  << " u=" << this->u()
644  << " k=" << this->k()
645  << std::flush;
646  std::vector<Quotient> quotients;
647  if ( this->k() >= 0 )
648  {
649  this->getCFrac( quotients );
650  out << " [" << quotients[ 0 ];
651  for ( unsigned int i = 1; i < quotients.size(); ++i )
652  out << "," << quotients[ i ];
653  out << "]";
654  }
655  out << " ]";
656  }
657 }
658 
659 ///////////////////////////////////////////////////////////////////////////////
660 // DGtal::SternBrocot<TInteger, TQuotient>
661 
662 //-----------------------------------------------------------------------------
663 template <typename TInteger, typename TQuotient>
664 inline
665 DGtal::SternBrocot<TInteger, TQuotient>::~SternBrocot()
666 {
667  if ( myOneOverOne != 0 ) delete myOneOverOne;
668  if ( myOneOverZero != 0 ) delete myOneOverZero;
669  if ( myZeroOverOne != 0 ) delete myZeroOverOne;
670 }
671 //-----------------------------------------------------------------------------
672 template <typename TInteger, typename TQuotient>
673 inline
674 DGtal::SternBrocot<TInteger, TQuotient>::SternBrocot()
675 {
676  myOneOverZero = new Node( NumberTraits<Integer>::ONE,
677  NumberTraits<Integer>::ZERO,
678  NumberTraits<Quotient>::ZERO,
679  -NumberTraits<Quotient>::ONE,
680  myZeroOverOne, 0, myOneOverOne, 0,
681  myZeroOverOne );
682  myZeroOverOne = new Node( NumberTraits<Integer>::ZERO,
683  NumberTraits<Integer>::ONE,
684  NumberTraits<Quotient>::ZERO,
685  NumberTraits<Quotient>::ZERO,
686  myZeroOverOne, myOneOverZero, 0, myOneOverOne,
687  myOneOverZero );
688  myOneOverOne = new Node( NumberTraits<Integer>::ONE,
689  NumberTraits<Integer>::ONE,
690  NumberTraits<Quotient>::ONE,
691  NumberTraits<Quotient>::ZERO,
692  myZeroOverOne, myOneOverZero, 0, 0,
693  myOneOverOne );
694  myOneOverZero->ascendantLeft = myZeroOverOne;
695  myOneOverZero->descendantLeft = myOneOverOne;
696  myOneOverZero->inverse = myZeroOverOne;
697  myZeroOverOne->ascendantLeft = myZeroOverOne;
698  myZeroOverOne->descendantRight = myOneOverOne;
699  myOneOverOne->inverse = myOneOverOne;
700  nbFractions = 3;
701 }
702 //-----------------------------------------------------------------------------
703 template <typename TInteger, typename TQuotient>
704 inline
705 DGtal::SternBrocot<TInteger, TQuotient> &
706 DGtal::SternBrocot<TInteger, TQuotient>::instance()
707 {
708  if ( singleton == 0 )
709  singleton = new SternBrocot;
710  return *singleton;
711 }
712 
713 
714 //-----------------------------------------------------------------------------
715 template <typename TInteger, typename TQuotient>
716 inline
717 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
718 DGtal::SternBrocot<TInteger, TQuotient>::zeroOverOne()
719 {
720  return Fraction( instance().myZeroOverOne );
721 }
722 //-----------------------------------------------------------------------------
723 template <typename TInteger, typename TQuotient>
724 inline
725 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
726 DGtal::SternBrocot<TInteger, TQuotient>::oneOverZero()
727 {
728  return Fraction( instance().myOneOverZero );
729 }
730 
731 ///////////////////////////////////////////////////////////////////////////////
732 // Interface - public :
733 
734 /**
735  * Writes/Displays the object on an output stream.
736  * @param out the output stream where the object is written.
737  */
738 template <typename TInteger, typename TQuotient>
739 inline
740 void
741 DGtal::SternBrocot<TInteger, TQuotient>::display( std::ostream & out,
742  const Fraction & f )
743 {
744  if ( f.null() ) out << "[Fraction null]";
745  else
746  {
747  out << "[Fraction f=" << f.p()
748  << "/" << f.q()
749  << " u=" << f.u()
750  << " k=" << f.k()
751  << std::flush;
752  std::vector<Quotient> quotients;
753  if ( f.k() >= 0 )
754  {
755  f.getCFrac( quotients );
756  out << " [" << quotients[ 0 ];
757  for ( unsigned int i = 1; i < quotients.size(); ++i )
758  out << "," << quotients[ i ];
759  out << "]";
760  }
761  out << " ]";
762  }
763 }
764 
765 /**
766  * Checks the validity/consistency of the object.
767  * @return 'true' if the object is valid, 'false' otherwise.
768  */
769 template <typename TInteger, typename TQuotient>
770 inline
771 bool
772 DGtal::SternBrocot<TInteger, TQuotient>::isValid() const
773 {
774  return true;
775 }
776 
777 ///////////////////////////////////////////////////////////////////////////////
778 // class SternBrocot
779 ///////////////////////////////////////////////////////////////////////////////
780 template <typename TInteger, typename TQuotient>
781 inline
782 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
783 DGtal::SternBrocot<TInteger, TQuotient>::fraction
784 ( Integer p, Integer q,
785  Fraction ancestor )
786 {
787  IntegerComputer<Integer> ic;
788  Integer g = ic.gcd( p, q );
789  if ( g != NumberTraits<Integer>::ZERO )
790  {
791  p /= g;
792  q /= g;
793  }
794  // special case 1/0
795  if ( ( p == NumberTraits<Integer>::ONE )
796  && ( q == NumberTraits<Integer>::ZERO ) ) return oneOverZero();
797  // other positive fractions
798  while ( ! ancestor.equals( p, q ) )
799  {
800  ASSERT( ( p + q ) >= ( ancestor.p() + ancestor.q() )
801  && "[ImaGene::SternBrocot::fraction] bad ancestor." );
802  ancestor = ancestor.lessThan( p, q )
803  ? ancestor.right()
804  : ancestor.left();
805  }
806  return ancestor;
807 }
808 
809 
810 ///////////////////////////////////////////////////////////////////////////////
811 // Implementation of inline functions //
812 
813 // JOL: invalid overloading
814 // template <typename TInteger, typename TQuotient>
815 // inline
816 // std::ostream&
817 // DGtal::operator<< ( std::ostream & out,
818 // const typename SternBrocot<TInteger, TQuotient>::Fraction & object )
819 // {
820 // typedef SternBrocot<TInteger,TQuotient> SB;
821 // SB::display( out, object );
822 // return out;
823 // }
824 
825 // //
826 ///////////////////////////////////////////////////////////////////////////////
827 
828