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.
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.
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/>.
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
26 * Implementation of inline methods defined in SternBrocot.h
28 * This file is part of the DGtal library.
32 //////////////////////////////////////////////////////////////////////////////
34 #include "DGtal/arithmetic/IntegerComputer.h"
35 //////////////////////////////////////////////////////////////////////////////
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
44 ///////////////////////////////////////////////////////////////////////////////
45 // DGtal::SternBrocot<TInteger, TQuotient>::Node
46 //-----------------------------------------------------------------------------
47 template <typename TInteger, typename TQuotient>
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,
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 ),
61 // std::cerr << "(" << p1 << "/" << q1 << "," << u1 << "," << k1 << ")";
64 //////////////////////////////////////////////////////////////////////////////
65 // DGtal::SternBrocot<TInteger, TQuotient>::Fraction
66 //-----------------------------------------------------------------------------
67 template <typename TInteger, typename TQuotient>
69 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
70 Fraction( Integer aP, Integer aQ, Fraction ancestor )
72 this->operator=( SternBrocotTree::fraction( aP, aQ, ancestor ) );
74 //-----------------------------------------------------------------------------
75 template <typename TInteger, typename TQuotient>
77 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
78 Fraction( Node* sb_node )
82 //-----------------------------------------------------------------------------
83 template <typename TInteger, typename TQuotient>
85 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
86 Fraction( const Self & other )
87 : myNode( other.myNode )
90 //-----------------------------------------------------------------------------
91 template <typename TInteger, typename TQuotient>
93 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction &
94 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
95 operator=( const Self & other )
99 myNode = other.myNode;
103 //-----------------------------------------------------------------------------
104 template <typename TInteger, typename TQuotient>
107 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
112 //-----------------------------------------------------------------------------
113 template <typename TInteger, typename TQuotient>
115 typename DGtal::SternBrocot<TInteger, TQuotient>::Integer
116 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
119 return myNode ? myNode->p : 0;
121 //-----------------------------------------------------------------------------
122 template <typename TInteger, typename TQuotient>
124 typename DGtal::SternBrocot<TInteger, TQuotient>::Integer
125 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
128 return myNode ? myNode->q : 0;
130 //-----------------------------------------------------------------------------
131 template <typename TInteger, typename TQuotient>
133 typename DGtal::SternBrocot<TInteger, TQuotient>::Quotient
134 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
137 ASSERT( myNode != 0 );
140 //-----------------------------------------------------------------------------
141 template <typename TInteger, typename TQuotient>
143 typename DGtal::SternBrocot<TInteger, TQuotient>::Quotient
144 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
147 ASSERT( myNode != 0 );
150 //-----------------------------------------------------------------------------
151 template <typename TInteger, typename TQuotient>
154 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
155 equals( Integer p1, Integer q1 ) const
157 return ( this->p() == p1 ) && ( this->q() == q1 );
159 //-----------------------------------------------------------------------------
160 template <typename TInteger, typename TQuotient>
163 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
164 lessThan( Integer p1, Integer q1 ) const
166 Integer d = p() * q1 - q() * p1;
167 return d < NumberTraits<Integer>::ZERO;
169 //-----------------------------------------------------------------------------
170 template <typename TInteger, typename TQuotient>
173 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
174 moreThan( Integer p1, Integer q1 ) const
176 Integer d = p() * q1 - q() * p1;
177 return d > NumberTraits<Integer>::ZERO;
179 //-----------------------------------------------------------------------------
180 template <typename TInteger, typename TQuotient>
183 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
184 operator==( const Fraction & other ) const
186 return myNode == other.myNode;
188 //-----------------------------------------------------------------------------
189 template <typename TInteger, typename TQuotient>
192 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
193 operator!=( const Fraction & other ) const
195 return myNode != other.myNode;
197 //-----------------------------------------------------------------------------
198 template <typename TInteger, typename TQuotient>
201 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
202 operator<( const Fraction & other ) const
204 return this->lessThan( other.p(), other.q() );
206 //-----------------------------------------------------------------------------
207 template <typename TInteger, typename TQuotient>
210 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
211 operator>( const Fraction & other ) const
213 return this->moreThan( other.p(), other.q() );
215 //-----------------------------------------------------------------------------
216 template <typename TInteger, typename TQuotient>
218 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
219 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
222 if ( myNode->descendantLeft == 0 )
224 Node* pleft = myNode->ascendantLeft;
225 Node* n = new Node( p() + pleft->p,
227 odd() ? u() + 1 : (Quotient) 2,
228 odd() ? k() : k() + 1,
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,
240 myNode->inverse->descendantRight = invn;
241 myNode->descendantLeft = n;
242 instance().nbFractions += 2;
244 return Fraction( myNode->descendantLeft );
246 //-----------------------------------------------------------------------------
247 template <typename TInteger, typename TQuotient>
249 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
250 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
253 if ( myNode->descendantRight == 0 )
255 Fraction inv( myNode->inverse );
257 ASSERT( myNode->descendantRight != 0 );
259 return Fraction( myNode->descendantRight );
261 //-----------------------------------------------------------------------------
262 template <typename TInteger, typename TQuotient>
265 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
268 return ( k() & 1 ) == 0;
270 //-----------------------------------------------------------------------------
271 template <typename TInteger, typename TQuotient>
274 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
277 return ( k() & 1 ) != 0;
279 //-----------------------------------------------------------------------------
280 template <typename TInteger, typename TQuotient>
282 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
283 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
286 return Fraction( odd() ? myNode->ascendantRight : myNode->ascendantLeft );
288 //-----------------------------------------------------------------------------
289 template <typename TInteger, typename TQuotient>
291 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
292 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
293 father( Quotient m ) const
295 if ( m > NumberTraits<Quotient>::ONE ) // > 1
299 n = odd() ? n->ascendantRight : n->ascendantLeft;
300 return Fraction( n );
302 else if ( m != NumberTraits<Quotient>::ZERO ) // == 1
304 return odd() ? previousPartial().right() : previousPartial().left();
307 return reduced( 2 ); //previousPartial().previousPartial();
309 //-----------------------------------------------------------------------------
310 template <typename TInteger, typename TQuotient>
312 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
313 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
314 previousPartial() const
316 return Fraction( odd() ? myNode->ascendantLeft : myNode->ascendantRight );
317 // return Fraction( odd()
318 // ? myNode->ascendantLeft->descendantRight
319 // : myNode->ascendantRight->descendantLeft );
321 //-----------------------------------------------------------------------------
322 template <typename TInteger, typename TQuotient>
324 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
325 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
328 return Fraction( myNode->inverse );
330 //-----------------------------------------------------------------------------
331 template <typename TInteger, typename TQuotient>
333 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
334 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
335 partial( Quotient kp ) const
337 ASSERT( ( ((Quotient)-2) <= kp ) && ( kp <= k() ) );
338 return reduced( k() - kp );
340 //-----------------------------------------------------------------------------
341 template <typename TInteger, typename TQuotient>
343 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
344 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
345 reduced( Quotient i ) const
347 ASSERT( ( ((Quotient)0) <= i ) && ( i <= ( k()+((Quotient)2) ) ) );
348 Node* n = this->myNode;
350 bool bleft = ( n->k & NumberTraits<Quotient>::ONE )
351 != NumberTraits<Quotient>::ZERO;
352 while ( i-- > NumberTraits<Quotient>::ZERO )
354 n = bleft ? n->ascendantLeft : n->ascendantRight;
357 return Fraction( n );
359 //-----------------------------------------------------------------------------
360 template <typename TInteger, typename TQuotient>
363 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
364 push_back( const std::pair<Quotient, Quotient> & quotient )
366 pushBack( quotient );
368 //-----------------------------------------------------------------------------
369 template <typename TInteger, typename TQuotient>
372 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
373 pushBack( const std::pair<Quotient, Quotient> & quotient )
375 // std::vector<Quotient> quots;
378 // this->getCFrac( quots );
379 // std::cerr << "F[";
380 // for ( unsigned int i = 0; i < quots.size(); ++i )
381 // std::cerr << " " << quots[ i ];
383 // else std::cerr << "[";
384 // std::cerr << "] + " << "(" << quotient.first
385 // << "," << quotient.second << ")";
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() );
395 Fraction f = zeroOverOne();
396 for ( Quotient i = 0; i < quotient.first; ++i )
398 this->operator=( f );
401 else if ( NumberTraits<Quotient>::even( quotient.second ) )
404 for ( Quotient i = 1; i < quotient.first; ++i )
406 this->operator=( f );
410 Fraction f = right();
411 for ( Quotient i = 1; i < quotient.first; ++i )
413 this->operator=( f );
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;
422 //-----------------------------------------------------------------------------
423 template <typename TInteger, typename TQuotient>
426 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
427 getSplit( Fraction & f1, Fraction & f2 ) const
429 f1.myNode = myNode->ascendantLeft;
430 f2.myNode = myNode->ascendantRight;
432 //-----------------------------------------------------------------------------
433 template <typename TInteger, typename TQuotient>
436 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
437 getSplitBerstel( Fraction & f1, Quotient & nb1,
438 Fraction & f2, Quotient & nb2 ) const
442 f1.myNode = myNode->ascendantLeft;
444 f2.myNode = f1.myNode->ascendantRight;
449 f2.myNode = myNode->ascendantRight;
451 f1.myNode = f2.myNode->ascendantLeft;
455 //-----------------------------------------------------------------------------
456 template <typename TInteger, typename TQuotient>
459 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
460 getCFrac( std::vector<Quotient> & quotients ) const
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();
467 bool bleft = odd() ? true : false;
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;
478 //-----------------------------------------------------------------------------
479 template <typename TInteger, typename TQuotient>
481 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction::ConstIterator
482 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
485 CFracSequence* seq = new CFracSequence;
486 this->getCFrac( *seq );
487 return ConstIterator( seq, seq->begin() );
489 //-----------------------------------------------------------------------------
490 template <typename TInteger, typename TQuotient>
492 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction::ConstIterator
493 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
496 static CFracSequence dummy;
497 return ConstIterator( 0, dummy.end() );
500 template <typename TInteger, typename TQuotient>
502 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
503 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
504 median(const Fraction & g) const
506 return Fraction(this->p()+g.p(),this->q()+g.q());
509 //----------------------------------------------------------------------------
510 template <typename TInteger, typename TQuotient>
512 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
513 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
514 simplestFractionInBetween(const Fraction & other) const
522 res = f; f = g; g = res;
529 ConstIterator itf=f.begin(), itg=g.begin();
531 DGtal::functors::Abs<TInteger> absComputer;
533 if(absComputer(f.p()*g.q()-f.q()*g.p())==NumberTraits<TInteger>::ONE)
536 itf = f.begin(); itg = g.begin();
537 uf = *itf; ug = *itg;
538 while(uf.first == ug.first && i != f.k() && i != g.k())
540 res.push_back(std::make_pair(uf.first,i));
547 if(uf.first==ug.first)
551 res.push_back(std::make_pair(uf.first,i));
555 res.push_back(std::make_pair(ug.first+1,i));
559 res.push_back(std::make_pair(uf.first,i));
563 res.push_back(std::make_pair(uf.first+1,i));
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));
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));
576 if(uf.first < ug.first)
577 res.push_back(std::make_pair(uf.first+1,i));
579 if(uf.first == ug.first + 1)
581 res.push_back(std::make_pair(ug.first,i));
585 if(ug.first==NumberTraits<TInteger>::ONE)
587 res.push_back(std::make_pair(ug.first,i));
591 res.push_back(std::make_pair(ug.first+1,i));
594 res.push_back(std::make_pair(2,i));
597 res.push_back(std::make_pair(ug.first+1,i));
601 if(ug.first < uf.first)
602 res.push_back(std::make_pair(ug.first+1,i));
604 if(ug.first == uf.first + 1)
606 res.push_back(std::make_pair(uf.first,i));
610 if(uf.first==NumberTraits<TInteger>::ONE)
612 res.push_back(std::make_pair(uf.first,i));
616 res.push_back(std::make_pair(uf.first+1,i));
619 res.push_back(std::make_pair(2,i));
622 res.push_back(std::make_pair(uf.first+1,i));
631 //-----------------------------------------------------------------------------
632 template <typename TInteger, typename TQuotient>
635 DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
636 selfDisplay( std::ostream & out ) const
638 if ( this->null() ) out << "[Fraction null]";
641 out << "[Fraction f=" << this->p()
643 << " u=" << this->u()
644 << " k=" << this->k()
646 std::vector<Quotient> quotients;
647 if ( this->k() >= 0 )
649 this->getCFrac( quotients );
650 out << " [" << quotients[ 0 ];
651 for ( unsigned int i = 1; i < quotients.size(); ++i )
652 out << "," << quotients[ i ];
659 ///////////////////////////////////////////////////////////////////////////////
660 // DGtal::SternBrocot<TInteger, TQuotient>
662 //-----------------------------------------------------------------------------
663 template <typename TInteger, typename TQuotient>
665 DGtal::SternBrocot<TInteger, TQuotient>::~SternBrocot()
667 if ( myOneOverOne != 0 ) delete myOneOverOne;
668 if ( myOneOverZero != 0 ) delete myOneOverZero;
669 if ( myZeroOverOne != 0 ) delete myZeroOverOne;
671 //-----------------------------------------------------------------------------
672 template <typename TInteger, typename TQuotient>
674 DGtal::SternBrocot<TInteger, TQuotient>::SternBrocot()
676 myOneOverZero = new Node( NumberTraits<Integer>::ONE,
677 NumberTraits<Integer>::ZERO,
678 NumberTraits<Quotient>::ZERO,
679 -NumberTraits<Quotient>::ONE,
680 myZeroOverOne, 0, myOneOverOne, 0,
682 myZeroOverOne = new Node( NumberTraits<Integer>::ZERO,
683 NumberTraits<Integer>::ONE,
684 NumberTraits<Quotient>::ZERO,
685 NumberTraits<Quotient>::ZERO,
686 myZeroOverOne, myOneOverZero, 0, myOneOverOne,
688 myOneOverOne = new Node( NumberTraits<Integer>::ONE,
689 NumberTraits<Integer>::ONE,
690 NumberTraits<Quotient>::ONE,
691 NumberTraits<Quotient>::ZERO,
692 myZeroOverOne, myOneOverZero, 0, 0,
694 myOneOverZero->ascendantLeft = myZeroOverOne;
695 myOneOverZero->descendantLeft = myOneOverOne;
696 myOneOverZero->inverse = myZeroOverOne;
697 myZeroOverOne->ascendantLeft = myZeroOverOne;
698 myZeroOverOne->descendantRight = myOneOverOne;
699 myOneOverOne->inverse = myOneOverOne;
702 //-----------------------------------------------------------------------------
703 template <typename TInteger, typename TQuotient>
705 DGtal::SternBrocot<TInteger, TQuotient> &
706 DGtal::SternBrocot<TInteger, TQuotient>::instance()
708 if ( singleton == 0 )
709 singleton = new SternBrocot;
714 //-----------------------------------------------------------------------------
715 template <typename TInteger, typename TQuotient>
717 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
718 DGtal::SternBrocot<TInteger, TQuotient>::zeroOverOne()
720 return Fraction( instance().myZeroOverOne );
722 //-----------------------------------------------------------------------------
723 template <typename TInteger, typename TQuotient>
725 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
726 DGtal::SternBrocot<TInteger, TQuotient>::oneOverZero()
728 return Fraction( instance().myOneOverZero );
731 ///////////////////////////////////////////////////////////////////////////////
732 // Interface - public :
735 * Writes/Displays the object on an output stream.
736 * @param out the output stream where the object is written.
738 template <typename TInteger, typename TQuotient>
741 DGtal::SternBrocot<TInteger, TQuotient>::display( std::ostream & out,
744 if ( f.null() ) out << "[Fraction null]";
747 out << "[Fraction f=" << f.p()
752 std::vector<Quotient> quotients;
755 f.getCFrac( quotients );
756 out << " [" << quotients[ 0 ];
757 for ( unsigned int i = 1; i < quotients.size(); ++i )
758 out << "," << quotients[ i ];
766 * Checks the validity/consistency of the object.
767 * @return 'true' if the object is valid, 'false' otherwise.
769 template <typename TInteger, typename TQuotient>
772 DGtal::SternBrocot<TInteger, TQuotient>::isValid() const
777 ///////////////////////////////////////////////////////////////////////////////
779 ///////////////////////////////////////////////////////////////////////////////
780 template <typename TInteger, typename TQuotient>
782 typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
783 DGtal::SternBrocot<TInteger, TQuotient>::fraction
784 ( Integer p, Integer q,
787 IntegerComputer<Integer> ic;
788 Integer g = ic.gcd( p, q );
789 if ( g != NumberTraits<Integer>::ZERO )
795 if ( ( p == NumberTraits<Integer>::ONE )
796 && ( q == NumberTraits<Integer>::ZERO ) ) return oneOverZero();
797 // other positive fractions
798 while ( ! ancestor.equals( p, q ) )
800 ASSERT( ( p + q ) >= ( ancestor.p() + ancestor.q() )
801 && "[ImaGene::SternBrocot::fraction] bad ancestor." );
802 ancestor = ancestor.lessThan( p, q )
810 ///////////////////////////////////////////////////////////////////////////////
811 // Implementation of inline functions //
813 // JOL: invalid overloading
814 // template <typename TInteger, typename TQuotient>
817 // DGtal::operator<< ( std::ostream & out,
818 // const typename SternBrocot<TInteger, TQuotient>::Fraction & object )
820 // typedef SternBrocot<TInteger,TQuotient> SB;
821 // SB::display( out, object );
826 ///////////////////////////////////////////////////////////////////////////////