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// DEFINITION of static data members
39///////////////////////////////////////////////////////////////////////////////
41template <typename TInteger, typename TQuotient>
42DGtal::SternBrocot<TInteger, TQuotient>*
43DGtal::SternBrocot<TInteger, TQuotient>::singleton = 0;
45///////////////////////////////////////////////////////////////////////////////
46// IMPLEMENTATION of inline methods.
47///////////////////////////////////////////////////////////////////////////////
49///////////////////////////////////////////////////////////////////////////////
50// ----------------------- Standard services ------------------------------
52///////////////////////////////////////////////////////////////////////////////
53// DGtal::SternBrocot<TInteger, TQuotient>::Node
54//-----------------------------------------------------------------------------
55template <typename TInteger, typename TQuotient>
57DGtal::SternBrocot<TInteger, TQuotient>::Node::
58Node( Integer p1, Integer q1, Quotient u1, Quotient k1,
59 Node* ascendant_left1, Node* ascendant_right1,
60 Node* descendant_left1, Node* descendant_right1,
62 : p( p1 ), q( q1 ), u( u1 ), k( k1 ),
63 ascendantLeft( ascendant_left1 ),
64 ascendantRight( ascendant_right1 ),
65 descendantLeft( descendant_left1 ),
66 descendantRight( descendant_right1 ),
69 // std::cerr << "(" << p1 << "/" << q1 << "," << u1 << "," << k1 << ")";
72//////////////////////////////////////////////////////////////////////////////
73// DGtal::SternBrocot<TInteger, TQuotient>::Fraction
74//-----------------------------------------------------------------------------
75template <typename TInteger, typename TQuotient>
77DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
78Fraction( Integer aP, Integer aQ, Fraction ancestor )
80 this->operator=( SternBrocotTree::fraction( aP, aQ, ancestor ) );
82//-----------------------------------------------------------------------------
83template <typename TInteger, typename TQuotient>
85DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
86Fraction( Node* sb_node )
90//-----------------------------------------------------------------------------
91template <typename TInteger, typename TQuotient>
93DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
94Fraction( const Self & other )
95 : myNode( other.myNode )
98//-----------------------------------------------------------------------------
99template <typename TInteger, typename TQuotient>
101typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction &
102DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
103operator=( const Self & other )
105 if ( this != &other )
107 myNode = other.myNode;
111//-----------------------------------------------------------------------------
112template <typename TInteger, typename TQuotient>
115DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
120//-----------------------------------------------------------------------------
121template <typename TInteger, typename TQuotient>
123typename DGtal::SternBrocot<TInteger, TQuotient>::Integer
124DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
127 return myNode ? myNode->p : 0;
129//-----------------------------------------------------------------------------
130template <typename TInteger, typename TQuotient>
132typename DGtal::SternBrocot<TInteger, TQuotient>::Integer
133DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
136 return myNode ? myNode->q : 0;
138//-----------------------------------------------------------------------------
139template <typename TInteger, typename TQuotient>
141typename DGtal::SternBrocot<TInteger, TQuotient>::Quotient
142DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
145 ASSERT( myNode != 0 );
148//-----------------------------------------------------------------------------
149template <typename TInteger, typename TQuotient>
151typename DGtal::SternBrocot<TInteger, TQuotient>::Quotient
152DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
155 ASSERT( myNode != 0 );
158//-----------------------------------------------------------------------------
159template <typename TInteger, typename TQuotient>
162DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
163equals( Integer p1, Integer q1 ) const
165 return ( this->p() == p1 ) && ( this->q() == q1 );
167//-----------------------------------------------------------------------------
168template <typename TInteger, typename TQuotient>
171DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
172lessThan( Integer p1, Integer q1 ) const
174 Integer d = p() * q1 - q() * p1;
175 return d < NumberTraits<Integer>::ZERO;
177//-----------------------------------------------------------------------------
178template <typename TInteger, typename TQuotient>
181DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
182moreThan( Integer p1, Integer q1 ) const
184 Integer d = p() * q1 - q() * p1;
185 return d > NumberTraits<Integer>::ZERO;
187//-----------------------------------------------------------------------------
188template <typename TInteger, typename TQuotient>
191DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
192operator==( const Fraction & other ) const
194 return myNode == other.myNode;
196//-----------------------------------------------------------------------------
197template <typename TInteger, typename TQuotient>
200DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
201operator!=( const Fraction & other ) const
203 return myNode != other.myNode;
205//-----------------------------------------------------------------------------
206template <typename TInteger, typename TQuotient>
209DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
210operator<( const Fraction & other ) const
212 return this->lessThan( other.p(), other.q() );
214//-----------------------------------------------------------------------------
215template <typename TInteger, typename TQuotient>
218DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
219operator>( const Fraction & other ) const
221 return this->moreThan( other.p(), other.q() );
223//-----------------------------------------------------------------------------
224template <typename TInteger, typename TQuotient>
226typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
227DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
230 if ( myNode->descendantLeft == 0 )
232 Node* pleft = myNode->ascendantLeft;
233 Node* n = new Node( p() + pleft->p,
235 odd() ? u() + 1 : (Quotient) 2,
236 odd() ? k() : k() + 1,
239 Fraction inv = Fraction( myNode->inverse );
240 Node* invpright = inv.myNode->ascendantRight;
241 Node* invn = new Node( inv.p() + invpright->p,
242 inv.q() + invpright->q,
243 inv.even() ? inv.u() + 1 : (Quotient) 2,
244 inv.even() ? inv.k() : inv.k() + 1,
245 myNode->inverse, invpright,
248 myNode->inverse->descendantRight = invn;
249 myNode->descendantLeft = n;
250 instance().nbFractions += 2;
252 return Fraction( myNode->descendantLeft );
254//-----------------------------------------------------------------------------
255template <typename TInteger, typename TQuotient>
257typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
258DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
261 if ( myNode->descendantRight == 0 )
263 Fraction inv( myNode->inverse );
265 ASSERT( myNode->descendantRight != 0 );
267 return Fraction( myNode->descendantRight );
269//-----------------------------------------------------------------------------
270template <typename TInteger, typename TQuotient>
273DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
276 return ( k() & 1 ) == 0;
278//-----------------------------------------------------------------------------
279template <typename TInteger, typename TQuotient>
282DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
285 return ( k() & 1 ) != 0;
287//-----------------------------------------------------------------------------
288template <typename TInteger, typename TQuotient>
290typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
291DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
294 return Fraction( odd() ? myNode->ascendantRight : myNode->ascendantLeft );
296//-----------------------------------------------------------------------------
297template <typename TInteger, typename TQuotient>
299typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
300DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
301father( Quotient m ) const
303 if ( m > NumberTraits<Quotient>::ONE ) // > 1
307 n = odd() ? n->ascendantRight : n->ascendantLeft;
308 return Fraction( n );
310 else if ( m != NumberTraits<Quotient>::ZERO ) // == 1
312 return odd() ? previousPartial().right() : previousPartial().left();
315 return reduced( 2 ); //previousPartial().previousPartial();
317//-----------------------------------------------------------------------------
318template <typename TInteger, typename TQuotient>
320typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
321DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
322previousPartial() const
324 return Fraction( odd() ? myNode->ascendantLeft : myNode->ascendantRight );
325 // return Fraction( odd()
326 // ? myNode->ascendantLeft->descendantRight
327 // : myNode->ascendantRight->descendantLeft );
329//-----------------------------------------------------------------------------
330template <typename TInteger, typename TQuotient>
332typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
333DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
336 return Fraction( myNode->inverse );
338//-----------------------------------------------------------------------------
339template <typename TInteger, typename TQuotient>
341typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
342DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
343partial( Quotient kp ) const
345 ASSERT( ( ((Quotient)-2) <= kp ) && ( kp <= k() ) );
346 return reduced( k() - kp );
348//-----------------------------------------------------------------------------
349template <typename TInteger, typename TQuotient>
351typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
352DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
353reduced( Quotient i ) const
355 ASSERT( ( ((Quotient)0) <= i ) && ( i <= ( k()+((Quotient)2) ) ) );
356 Node* n = this->myNode;
358 bool bleft = ( n->k & NumberTraits<Quotient>::ONE )
359 != NumberTraits<Quotient>::ZERO;
360 while ( i-- > NumberTraits<Quotient>::ZERO )
362 n = bleft ? n->ascendantLeft : n->ascendantRight;
365 return Fraction( n );
367//-----------------------------------------------------------------------------
368template <typename TInteger, typename TQuotient>
371DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
372push_back( const std::pair<Quotient, Quotient> & quotient )
374 pushBack( quotient );
376//-----------------------------------------------------------------------------
377template <typename TInteger, typename TQuotient>
380DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
381pushBack( const std::pair<Quotient, Quotient> & quotient )
383 // std::vector<Quotient> quots;
386 // this->getCFrac( quots );
387 // std::cerr << "F[";
388 // for ( unsigned int i = 0; i < quots.size(); ++i )
389 // std::cerr << " " << quots[ i ];
391 // else std::cerr << "[";
392 // std::cerr << "] + " << "(" << quotient.first
393 // << "," << quotient.second << ")";
396 ASSERT( quotient.second <= NumberTraits<Quotient>::ZERO );
397 if ( quotient.second < NumberTraits<Quotient>::ZERO )
398 this->operator=( oneOverZero() );
399 else if ( quotient.first == NumberTraits<Quotient>::ZERO ) // (0,0)
400 this->operator=( zeroOverOne() );
403 Fraction f = zeroOverOne();
404 for ( Quotient i = 0; i < quotient.first; ++i )
406 this->operator=( f );
409 else if ( NumberTraits<Quotient>::even( quotient.second ) )
412 for ( Quotient i = 1; i < quotient.first; ++i )
414 this->operator=( f );
418 Fraction f = right();
419 for ( Quotient i = 1; i < quotient.first; ++i )
421 this->operator=( f );
424 // this->getCFrac( quots );
425 // std::cerr << " => F[";
426 // for ( unsigned int i = 0; i < quots.size(); ++i )
427 // std::cerr << " " << quots[ i ];
428 // std::cerr << "]" << std::endl;
430//-----------------------------------------------------------------------------
431template <typename TInteger, typename TQuotient>
434DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
435getSplit( Fraction & f1, Fraction & f2 ) const
437 f1.myNode = myNode->ascendantLeft;
438 f2.myNode = myNode->ascendantRight;
440//-----------------------------------------------------------------------------
441template <typename TInteger, typename TQuotient>
444DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
445getSplitBerstel( Fraction & f1, Quotient & nb1,
446 Fraction & f2, Quotient & nb2 ) const
450 f1.myNode = myNode->ascendantLeft;
452 f2.myNode = f1.myNode->ascendantRight;
457 f2.myNode = myNode->ascendantRight;
459 f1.myNode = f2.myNode->ascendantLeft;
463//-----------------------------------------------------------------------------
464template <typename TInteger, typename TQuotient>
467DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
468getCFrac( std::vector<Quotient> & quotients ) const
470 ASSERT( this->k() >= NumberTraits<Quotient>::ZERO );
471 int64_t i = NumberTraits<Quotient>::castToInt64_t( this->k() );
472 quotients.resize( (unsigned int)i + 1 );
473 quotients[ (unsigned int)i-- ] = this->u();
475 bool bleft = odd() ? true : false;
478 ASSERT( n->k >= NumberTraits<Quotient>::ZERO );
479 n = bleft ? n->ascendantLeft : n->ascendantRight;
480 quotients[ (unsigned int)i ] =
481 ( i == NumberTraits<Quotient>::castToInt64_t( n->k ) ) ? n->u : 1;
486//-----------------------------------------------------------------------------
487template <typename TInteger, typename TQuotient>
489typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction::ConstIterator
490DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
493 CFracSequence* seq = new CFracSequence;
494 this->getCFrac( *seq );
495 return ConstIterator( seq, seq->begin() );
497//-----------------------------------------------------------------------------
498template <typename TInteger, typename TQuotient>
500typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction::ConstIterator
501DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
504 static CFracSequence dummy;
505 return ConstIterator( 0, dummy.end() );
508template <typename TInteger, typename TQuotient>
510typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
511DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
512median(const Fraction & g) const
514 return Fraction(this->p()+g.p(),this->q()+g.q());
517//----------------------------------------------------------------------------
518template <typename TInteger, typename TQuotient>
520typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
521DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
522simplestFractionInBetween(const Fraction & other) const
530 res = f; f = g; g = res;
537 ConstIterator itf=f.begin(), itg=g.begin();
539 DGtal::functors::Abs<TInteger> absComputer;
541 if(absComputer(f.p()*g.q()-f.q()*g.p())==NumberTraits<TInteger>::ONE)
544 itf = f.begin(); itg = g.begin();
545 uf = *itf; ug = *itg;
546 while(uf.first == ug.first && i != f.k() && i != g.k())
548 res.push_back(std::make_pair(uf.first,i));
555 if(uf.first==ug.first)
559 res.push_back(std::make_pair(uf.first,i));
563 res.push_back(std::make_pair(ug.first+1,i));
567 res.push_back(std::make_pair(uf.first,i));
571 res.push_back(std::make_pair(uf.first+1,i));
576 if(i!=f.k() && i != g.k())
577 (uf.first<ug.first)?res.push_back(std::make_pair(uf.first+1,i)):res.push_back(std::make_pair(ug.first+1,i));
579 if(i == f.k() && i == g.k())
580 (uf.first<ug.first)?res.push_back(std::make_pair(uf.first+1,i)):res.push_back(std::make_pair(ug.first+1,i));
584 if(uf.first < ug.first)
585 res.push_back(std::make_pair(uf.first+1,i));
587 if(uf.first == ug.first + 1)
589 res.push_back(std::make_pair(ug.first,i));
593 if(ug.first==NumberTraits<TInteger>::ONE)
595 res.push_back(std::make_pair(ug.first,i));
599 res.push_back(std::make_pair(ug.first+1,i));
602 res.push_back(std::make_pair(2,i));
605 res.push_back(std::make_pair(ug.first+1,i));
609 if(ug.first < uf.first)
610 res.push_back(std::make_pair(ug.first+1,i));
612 if(ug.first == uf.first + 1)
614 res.push_back(std::make_pair(uf.first,i));
618 if(uf.first==NumberTraits<TInteger>::ONE)
620 res.push_back(std::make_pair(uf.first,i));
624 res.push_back(std::make_pair(uf.first+1,i));
627 res.push_back(std::make_pair(2,i));
630 res.push_back(std::make_pair(uf.first+1,i));
639//-----------------------------------------------------------------------------
640template <typename TInteger, typename TQuotient>
643DGtal::SternBrocot<TInteger, TQuotient>::Fraction::
644selfDisplay( std::ostream & out ) const
646 if ( this->null() ) out << "[Fraction null]";
649 out << "[Fraction f=" << this->p()
651 << " u=" << this->u()
652 << " k=" << this->k()
654 std::vector<Quotient> quotients;
655 if ( this->k() >= 0 )
657 this->getCFrac( quotients );
658 out << " [" << quotients[ 0 ];
659 for ( unsigned int i = 1; i < quotients.size(); ++i )
660 out << "," << quotients[ i ];
667///////////////////////////////////////////////////////////////////////////////
668// DGtal::SternBrocot<TInteger, TQuotient>
670//-----------------------------------------------------------------------------
671template <typename TInteger, typename TQuotient>
673DGtal::SternBrocot<TInteger, TQuotient>::~SternBrocot()
675 if ( myOneOverOne != 0 ) delete myOneOverOne;
676 if ( myOneOverZero != 0 ) delete myOneOverZero;
677 if ( myZeroOverOne != 0 ) delete myZeroOverOne;
679//-----------------------------------------------------------------------------
680template <typename TInteger, typename TQuotient>
682DGtal::SternBrocot<TInteger, TQuotient>::SternBrocot()
684 myOneOverZero = new Node( NumberTraits<Integer>::ONE,
685 NumberTraits<Integer>::ZERO,
686 NumberTraits<Quotient>::ZERO,
687 -NumberTraits<Quotient>::ONE,
688 myZeroOverOne, 0, myOneOverOne, 0,
690 myZeroOverOne = new Node( NumberTraits<Integer>::ZERO,
691 NumberTraits<Integer>::ONE,
692 NumberTraits<Quotient>::ZERO,
693 NumberTraits<Quotient>::ZERO,
694 myZeroOverOne, myOneOverZero, 0, myOneOverOne,
696 myOneOverOne = new Node( NumberTraits<Integer>::ONE,
697 NumberTraits<Integer>::ONE,
698 NumberTraits<Quotient>::ONE,
699 NumberTraits<Quotient>::ZERO,
700 myZeroOverOne, myOneOverZero, 0, 0,
702 myOneOverZero->ascendantLeft = myZeroOverOne;
703 myOneOverZero->descendantLeft = myOneOverOne;
704 myOneOverZero->inverse = myZeroOverOne;
705 myZeroOverOne->ascendantLeft = myZeroOverOne;
706 myZeroOverOne->descendantRight = myOneOverOne;
707 myOneOverOne->inverse = myOneOverOne;
710//-----------------------------------------------------------------------------
711template <typename TInteger, typename TQuotient>
713DGtal::SternBrocot<TInteger, TQuotient> &
714DGtal::SternBrocot<TInteger, TQuotient>::instance()
716 if ( singleton == 0 )
717 singleton = new SternBrocot;
722//-----------------------------------------------------------------------------
723template <typename TInteger, typename TQuotient>
725typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
726DGtal::SternBrocot<TInteger, TQuotient>::zeroOverOne()
728 return Fraction( instance().myZeroOverOne );
730//-----------------------------------------------------------------------------
731template <typename TInteger, typename TQuotient>
733typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
734DGtal::SternBrocot<TInteger, TQuotient>::oneOverZero()
736 return Fraction( instance().myOneOverZero );
739///////////////////////////////////////////////////////////////////////////////
740// Interface - public :
743 * Writes/Displays the object on an output stream.
744 * @param out the output stream where the object is written.
746template <typename TInteger, typename TQuotient>
749DGtal::SternBrocot<TInteger, TQuotient>::display( std::ostream & out,
752 if ( f.null() ) out << "[Fraction null]";
755 out << "[Fraction f=" << f.p()
760 std::vector<Quotient> quotients;
763 f.getCFrac( quotients );
764 out << " [" << quotients[ 0 ];
765 for ( unsigned int i = 1; i < quotients.size(); ++i )
766 out << "," << quotients[ i ];
774 * Checks the validity/consistency of the object.
775 * @return 'true' if the object is valid, 'false' otherwise.
777template <typename TInteger, typename TQuotient>
780DGtal::SternBrocot<TInteger, TQuotient>::isValid() const
785///////////////////////////////////////////////////////////////////////////////
787///////////////////////////////////////////////////////////////////////////////
788template <typename TInteger, typename TQuotient>
790typename DGtal::SternBrocot<TInteger, TQuotient>::Fraction
791DGtal::SternBrocot<TInteger, TQuotient>::fraction
792( Integer p, Integer q,
795 IntegerComputer<Integer> ic;
796 Integer g = ic.gcd( p, q );
797 if ( g != NumberTraits<Integer>::ZERO )
803 if ( ( p == NumberTraits<Integer>::ONE )
804 && ( q == NumberTraits<Integer>::ZERO ) ) return oneOverZero();
805 // other positive fractions
806 while ( ! ancestor.equals( p, q ) )
808 ASSERT( ( p + q ) >= ( ancestor.p() + ancestor.q() )
809 && "[ImaGene::SternBrocot::fraction] bad ancestor." );
810 ancestor = ancestor.lessThan( p, q )
818///////////////////////////////////////////////////////////////////////////////
819// Implementation of inline functions //
821// JOL: invalid overloading
822// template <typename TInteger, typename TQuotient>
825// DGtal::operator<< ( std::ostream & out,
826// const typename SternBrocot<TInteger, TQuotient>::Fraction & object )
828// typedef SternBrocot<TInteger,TQuotient> SB;
829// SB::display( out, object );
834///////////////////////////////////////////////////////////////////////////////