frommer.cpp

Go to the documentation of this file.

// Frommer.cpp: Implementierung der Klasse Frommer.
//

 
using namespace nCenterFocus;

template < class TFrommerDefs, int variant>
const typename TFrommerDefs::RingType* Frommer< TFrommerDefs, variant >::getRing() const
{
        return ring;
};

template < class TFrommerDefs, int variant>
const typename TFrommerDefs::RingType& Frommer< TFrommerDefs, variant >::getRingRef() const
{
        return *ring;
};
// Konstruktion/Destruktion
template < class TFrommerDefs, int variant >
Frommer< TFrommerDefs, variant >::Frommer(short m,const TRing *_ring):  
                                                        name_m(""),
                                                        ring(_ring), 
                                                        field_ref_m( _ring->getFieldRef() ),
                                                        maxFocalValuesToCompute(m),
                                                        maxDegree( m *2 + 3),
                                                        A  ( new TMatrixXYType ("A"  ,maxDegree)),
                                                        B  ( new TMatrixXYType ("B"  ,maxDegree)),

                                                        dxA( new TMatrixXYType ("dxA",maxDegree)),

                                                        dyA( new TMatrixXYType ("dyA",maxDegree)),
                                                        a_table(createLowATable() ),

                                                        currFocalValuePos(0)
                                                        #ifdef FROMMERSTAGE_TIMER
                                                                , frommerTimer_m(m)
                                                        #endif

{
        #ifdef FROMMERSTAGE_TIMER
                tim_m.clear();
                tim2_m.clear();
        #endif

        
        #ifdef SAFE
                initialized=false;
        #endif

        minusp = NULL;
        q      = NULL;
        lastDegree = maxDegree;

        assert(minusp == NULL);

        //computedFocalValues = new typename TRing::ElementType[ maxFocalValuesToCompute+1 ];
        computedFocalValues = new typename TRing::ElementType[ maxFocalValuesToCompute+1 ];
        

        assert(minusp == NULL);

        A->setCoeff (2,0,  TRing::ElementType::One);
        A->setCoeff (1,1,  TRing::ElementType::Zero);
        A->setCoeff (0,2,  TRing::ElementType::One);


}


template < class TFrommerDefs, int variant>
Frommer< TFrommerDefs, variant >::~Frommer()
{
        delete[] computedFocalValues;
        delete   B;
        delete   A;
        delete   dxA;
        delete   dyA;
        delete[] a_table;
        #ifdef FROMMERSTAGE_TIMER
                MtcpManager_g.disconnectTimer(tim_m);
                MtcpManager_g.disconnectTimer(tim2_m);
        #endif
}


template < class TFrommerDefs, int variant>
void Frommer< TFrommerDefs, variant >::printLowATable(typename TFrommerDefs::RingType::ScalarType* aTable)
{

        for (int row=0;row<maxDegree+1;row++)
        {
                for (int col=0;col<maxDegree+1;col++)
                {
                        cerr << aTable[(maxDegree+1)*row+col] << " " ;
                }
                cerr << "\n" << endl;
        }
}

template < class TFrommerDefs, int variant>
typename TFrommerDefs::RingType::ScalarType* Frommer< TFrommerDefs, variant >::createLowATable()
{

        int tableSize= (maxDegree+1)*(maxDegree+1);

        typename TRing::ScalarType* t_a_table = new typename TRing::ScalarType[ tableSize ];

        short   i,j, h;

        for ( h=0; h<tableSize; h++)
        {
                t_a_table[h ] = TRing::ScalarType::Zero;
        }

        for ( h=0; h<maxDegree+1; h++)
        {
                t_a_table[ h*(maxDegree+1) ] =  TRing::ScalarType::One;
                        
                for ( i=1; i<= h/2; i++)
                {
                        typename TRing::ScalarType  Produkt =  TRing::ScalarType::One;
                        for ( j=1; j<=i; j++)
                        {
                                field_ref_m.multiplyInPlace(    Produkt,
                                                                                field_ref_m.multiply(   h - (2*j-1), 
                                                                                                                        field_ref_m.multInv(2*j-1)
                                                                                                ) 
                                                                );
                        }
                        t_a_table[ h*(maxDegree+1)+i ]= Produkt;
                }
        }
        //printLowATable(t_a_table); getchar();
        return t_a_table;
}


template < class TFrommerDefs,int variant>
inline void Frommer< TFrommerDefs, variant >::init()
{

        #ifdef SAFE
                assert(A->getCoeff(2,0) == TRing::ElementType::One);
                assert(A->getCoeff(1,1) == TRing::ElementType::Zero);
                assert(A->getCoeff(0,2) == TRing::ElementType::One);
        
                A->clear();
                B->clear();
                dxA->clear();
                dyA->clear();

                A->setCoeff (2,0,  TRing::ElementType::One);
                A->setCoeff (1,1,  TRing::ElementType::Zero);
                A->setCoeff (0,2,  TRing::ElementType::One);
                
        #else
                // A->clear(); // nicht mehr notwendig !
                // dyA->clear(lastDegree);//  nicht mehr notwendig!
                B->clear(lastDegree); //notwendig bei Wiederverwendung!
                // dxA->clear(lastDegree); // nicht mehr notwendig, zumindest nicht für Grad 3, weil
                //  vor dem C-Schritt dxa und dya vom aktuellenGrad-2(zeile+spalte) initialisiert  werden 
        #endif
}



template < class TFrommerDefs, int variant>
inline void Frommer< TFrommerDefs, variant >::setPolynoms(const TPolynomXY * const minuspp, const TPolynomXY *const qq) 
{

        // setze  berechneten focal Values = nix ! (computedFocalValus = currFocalValuePos-1 )
        currFocalValuePos=1;
        #ifdef SAFE
                initialized=true;
        #endif
        #ifdef SSAFE
                assert( minuspp->getDegree() == qq->getDegree() );
                if (minusp!=NULL)
                        assert( minusp->getDegree() == qq->getDegree() );
                initialized=true;
        #endif

        q       =       qq;
        minusp  =       minuspp;

        

};


template < class TFrommerDefs, int variant >
inline  void Frommer< TFrommerDefs, variant >::compute_dyA_and_dxA()
{
        register const short    currentDegree_minus_2           = currentDegree -2; //fest
        register short          currentDegree_minus_2_minus_i   = currentDegree_minus_2; //Laufvariable

        for (register short i=0; i<=currentDegree_minus_2; ++i,--currentDegree_minus_2_minus_i)
        {

                //dyA[i, currentDegree-2-i]:=(currentDegree-1-i)*A[i, currentDegree-1-i];
                dyA->setCoeff(  i, currentDegree_minus_2_minus_i, 
                                        ring->scalarMultiply (  currentDegree_minus_2_minus_i + 1 ,
                                                                        A->getCoeffConstRef( i   , currentDegree_minus_2_minus_i + 1)
                                                        
                                                                )
                                );

                //dxA[i, aktueller_grad-2-i]:=(i+1)*A[i+1, aktueller_grad-2-i];
                dxA->setCoeff(  i, currentDegree_minus_2_minus_i, 
                                        ring->scalarMultiply (  i+1 ,
                                                                        A->getCoeffConstRef( i+1 , currentDegree_minus_2_minus_i  )
                                                        
                                                                )
                                );
        }
}


template < class TFrommerDefs, int variant >
inline void Frommer< TFrommerDefs, variant >::perform_inner_C_Step(register const short j, 
                                                                                         register const short l,
                                                                                         register  short n       )
{
        // m läuft, j,l hier fest, n =(currentDegree-m-j-l) folgt aus Bedingung (m+n+j+l==currentDegree).
        // Berechnungen finden nur statt,
        //  wenn (m>=0), (n>=0) und  (m+n)>0 und 
        //(m+n)<=currentDegree-2  (letztere Bedingung ist nach Opel 
        //zumindest für  Polynome dritten Grades berechtigt, siehe :
        // implizite Annahmen: bei Funktionseintritt ist n>0; (n<=currentDegree-2), (j+l)>1
        // im ersten Schritt sind alle Bedingungen erfüllt, erst ab dem zweiten (m>=1) muss 
        //lediglich geprüft werden, ob (n>=0) ist.
        
        //-------------------
        // Starte also bei m=0:
        #ifdef SAFE
                assert(n<=currentDegree-2);
                assert(n>0);
        #endif 

        const   typename TRing::ElementType &   ppp = minusp->getCoeffConstRef(j, l); 

        

        register        typename TRing::ElementType & b   = B->getCoeffRef(j, l+n);     

        //ACCUMULATE MULTIPLICATION: b = b + p(j,l) * dya(m,n)
        ring->accMultRef(b, ppp, dyA->getCoeffConstRef(0, n));
        

        const   typename TRing::ElementType &  qqq = q->getCoeffConstRef(j, l); 
        //ACCUMULATE MULTIPLICATION:  b = b + q(j,l) * dxa(m,n)
        ring->accMultRef(b, qqq, dxA->getCoeffConstRef(0, n));

        //----------m>0
        --n; // das dazugehoerige m++ steht in der Scheifeninitialisierung (m=1)

        // bei erfullten Vorbedingungen reicht es zu prüfen ob n>=0 .
        int tmpn = n;
        assert( n>=0);
        
        for (register unsigned short m=1; n>=0; ++m,--n)
        {
                
                //register typename TRing::ElementType & b2 = ; 
                //ACCUMULATE mULTIPLICATION:   b2 = b2 + q(j, l) * dxa(m,n)
                ring->accMultRef(B->getCoeffRef(j+m, l+n), qqq, dxA->getCoeffConstRef(m, n));
        
                // Formel moeglicherweiwe so implementieren,
                // dass erst am Schluss eines C-schritts, also beim Plusgleich 
                //das Ergebnis in das des vorgegebenen Koerpers umgerechnet wird.
                // dazu ist 
                //1 .eine Analyse notwendig, welche Operation wie lange dauert.
                //2. Ein entsprechend geraeumiger Datentyp fuer ein B-Eintrag, welcher nicht gleich ueberlaeuft.
        }
        n=tmpn;
        for (register unsigned short m=1; n>=0; ++m,--n)
        
        {
        //      register typename TRing::ElementType & b2 = B->getCoeffRef(j+m, l+n);   

                //ACCUMULATE mULTIPLICATION:   b2 = b2 + p(j, l) * dya(m,n)
                ring->accMultRef(B->getCoeffRef(j+m, l+n), ppp, dyA->getCoeffConstRef(m, n));
        }
}




template < class TFrommerDefs, int variant >
template <int depth,int stage>
inline void Frommer< TFrommerDefs, variant >::perform_inner_C_Step_differ_stages(register const short j, 
                                                                                         register const short l,
                                                                                         register  short n       )
{
        // m läuft, j,l hier fest, n =(currentDegree-m-j-l) folgt aus Bedingung (m+n+j+l==currentDegree).
        // Berechnungen finden nur statt,
        //  wenn (m>=0), (n>=0) und  (m+n)>0 und 
        //(m+n)<=currentDegree-2  (letztere Bedingung ist nach Opel 
        //zumindest für  Polynome dritten Grades berechtigt, siehe :
        // implizite Annahmen: bei Funktionseintritt ist n>0; (n<=currentDegree-2), (j+l)>1
        // im ersten Schritt sind alle Bedingungen erfüllt, erst ab dem zweiten (m>=1) muss 
        //lediglich geprüft werden, ob (n>=0) ist.
        
        //-------------------
        // Starte also bei m=0:
        #ifdef SAFE
                assert(n<=currentDegree-2);
                assert(n>0);
        #endif 

        const   typename TRing::ElementType &   ppp = minusp->getCoeffConstRef(j, l); 

        

        register        typename TRing::ElementType & b   = B->getCoeffRef(j, l+n);     

        //ACCUMULATE MULTIPLICATION: b = b + p(j,l) * dya(m,n)
        ring->accMultRef(b, ppp, dyA->getCoeffConstRef(0, n));
        

        const   typename TRing::ElementType &  qqq = q->getCoeffConstRef(j, l); 
        //ACCUMULATE MULTIPLICATION:  b = b + q(j,l) * dxa(m,n)
        ring->accMultRef(b, qqq, dxA->getCoeffConstRef(0, n));

        //----------m>0
        --n; // das dazugehoerige m++ steht in der Scheifeninitialisierung (m=1)

        // bei erfullten Vorbedingungen reicht es zu prüfen ob n>=0 .

        cerr << "ppp" << ppp << endl;
        cerr << " qqq" << qqq << endl;
        cerr << " bo" << B->getCoeffRef(j, l+n) << endl;

        for (register unsigned short m=1; n>=0; ++m,--n)
        {
                
                //register typename TRing::ElementType & b2 = ; 
                //ACCUMULATE mULTIPLICATION:   b2 = b2 + q(j, l) * dxa(m,n)
                ring->accMultRef(B->getCoeffRef(j+m, l+n), qqq, dxA->getCoeffConstRef(m, n));
        
                ring->accMultRef(B->getCoeffRef(j+m, l+n), ppp, dyA->getCoeffConstRef(m, n));

                // Formel moeglicherweiwe so implementieren,
                // dass erst am Schluss eines C-schritts, also beim Plusgleich 
                //das Ergebnis in das des vorgegebenen Koerpers umgerechnet wird.
                // dazu ist 
                //1 .eine Analyse notwendig, welche Operation wie lange dauert.
                //2. Ein entsprechend geraeumiger Datentyp fuer ein B-Eintrag, welcher nicht gleich ueberlaeuft.
        }
}


template < class TFrommerDefs, int variant >
template <int depth, int stage>
inline void Frommer< TFrommerDefs, variant >::perform_C_Step_differ_stages()
{

        register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.

        switch ( q->getDegree() )
        {
                                
        case 3: 
        
                #ifdef SAFE
                        assert(n>0);//gilt normalerweise immer, da currentDegree bei 3 startet!
                #endif
                //if ( n>0) gilt immer!
                {       
                        //(0,2)
                                 perform_inner_C_Step_differ_stages<depth, stage>(0,2,n);
                        //(1,1)
                                 perform_inner_C_Step_differ_stages<depth,stage>(1,1,n);
                
                        //(2,0)
                                 perform_inner_C_Step_differ_stages<depth,stage>(2,0,n);
                        --n;
                        //if (n>0) gilt immer fuer currentDegree >=4
                        if (n>0)
                        {
                                //(0,3)
                                         perform_inner_C_Step_differ_stages<depth,stage>(0,3,n);
                                
                                //(1,2)
                                         perform_inner_C_Step_differ_stages<depth,stage>(1,2,n);
                                        
                                //(2,1)
                                         perform_inner_C_Step_differ_stages<depth,stage>(2,1,n);
                                
                                //(3,0)
                                         perform_inner_C_Step_differ_stages<depth,stage>(3,0,n);
                        }
                }
        
                break;
        
        case 2:
                if (n>0) // gilt immer !
                {
                
                        //(0,2)
                                 perform_inner_C_Step_differ_stages<depth,stage>(0,2,n);
        #ifdef DEBUG2
                                                this->outputMatrix(B);
                                        #endif
                        
                        //(1,1)
                                 perform_inner_C_Step_differ_stages<depth,stage>(1,1,n);
        #ifdef DEBUG2
                                                this->outputMatrix(B);
                                        #endif
                        
                        //(2,0)
                                 perform_inner_C_Step_differ_stages<depth,stage>(2,0,n);        
        #ifdef DEBUG2
                                                this->outputMatrix(B);
                                        #endif
                }
                
                break;
        
        case 1:
        
                break;
        
        
        default:
                perform_generic_C_Step();
        }
}



template <class TFrommerDefs, int variant>
inline void Frommer< TFrommerDefs, variant >::perform_generic_C_Step()
{
         
        short d = q->getDegree(); // p->getgrad kann nicht genommen werden!! sonst maxDegree evtl falsch!!
        
                if (currentDegree<d)
                        d=currentDegree;
        
                for (register short  j=0; j<=d; j++)
                        for (register short l=0; l<=d-j; l++)
                        {
                                short nn=currentDegree-j-l;

                        

                                if (nn<=currentDegree-2 && nn>0)
                                {
                                         perform_inner_C_Step(j, l, nn);

                                        #ifdef DEBUG2
                                                this->outputMatrix(B);
                                        #endif

                                }
                        }
}


template < class TFrommerDefs, int variant >
inline void Frommer< TFrommerDefs, variant >::perform_C_Step()
{

        register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.

        switch ( q->getDegree() )
        {
                                
        case 3: 
        
                #ifdef SAFE
                        assert(n>0);//gilt normalerweise immer, da currentDegree bei 3 startet!
                #endif
                //if ( n>0) gilt immer!
                {       
                        //(0,2)
                                perform_inner_C_Step(0,2,n);
                        //(1,1)
                                perform_inner_C_Step(1,1,n);
                
                        //(2,0)
                                perform_inner_C_Step(2,0,n);
                        --n;
                        //if (n>0)
                        if (n>0)
                        {
                                //(0,3)
                                        perform_inner_C_Step(0,3,n);
                                
                                //(1,2)
                                        perform_inner_C_Step(1,2,n);
                                        
                                //(2,1)
                                        perform_inner_C_Step(2,1,n);
                                
                                //(3,0)
                                        perform_inner_C_Step(3,0,n);
                        }
                }
        
                break;
        
        case 2:
                if (n>0)
                {
                
                        //(0,2)
                                perform_inner_C_Step(0,2,n);
                                #ifdef DEBUG2
                                                this->outputMatrix(B);
                                #endif
                        
                        //(1,1)
                                perform_inner_C_Step(1,1,n);
                                        #ifdef DEBUG2
                                                this->outputMatrix(B);
                                        #endif
                        
                        //(2,0)
                                perform_inner_C_Step(2,0,n);    
                                #ifdef DEBUG2
                                                this->outputMatrix(B);
                                #endif
                }
                
                break;
        
        case 1:
        
                break;
        
        
        default:
                perform_generic_C_Step();
        }
}


template < class TFrommerDefs, int variant >
inline  void Frommer< TFrommerDefs, variant >::perform_the_A_Step()
{
        //ungerader A-Schritt, siehe http://www.uni-bayreuth.de/departments/math/org/mathe6/publ/da/hoehn/node13.html
        if (currentDegree % 2 == 1)  
        {
                perform_the_odd_A_Step();
        }
        else //Gerader A-Schritt, siehe z.B. http://www.uni-bayreuth.de/departments/math/org/mathe6/publ/da/hoehn/node14.html
        {
                perform_the_even_A_Step();
                
        }

}
template < class TFrommerDefs, int variant >
inline  void Frommer< TFrommerDefs, variant >::perform_the_odd_A_Step()
{
        //A(n-1,1) = B(n,0)
                A->setCoeff(currentDegree-1, 1, B->getCoeffRef(currentDegree,0));
        //      typename TRing::ScalarType* test = new typename TRing::ScalarType(5);

        //      assert (typename TRing::ScalarType(5) == *test );
        //      typename TRing::ScalarType* test2 =&(typename TRing::ScalarType(5));
        //      assert (*test2 == *test);
                //j=3,5,7,..,n:
                //A(n-j,j) = ((n-j+2)A(n-j+2,j-2)+B(n-j+1,j-1))/j
                for(register short j=3;j<=currentDegree;j+=2)
                {
                        short j_minus_zwei=j-2;
                        A->setCoeff(    currentDegree-j,j,
                                        ring->scalarMultiplyRef(        field_ref_m.multInv( typename TRing::ScalarType(j) ),
                                                                        ring->addRef(   ring->scalarMultiplyRef(        currentDegree-j_minus_zwei, 
                                                                                                                        A->getCoeffConstRef(currentDegree-j_minus_zwei,j_minus_zwei)
                                                                                                                ),
                                                                                        B->getCoeffConstRef(currentDegree-j+1,j-1)
                                                                                )
                                                                )
                                );
                }

                //A(1,n-1) = -B(0,n)
                A->setCoeff     (       1, currentDegree-1, 
                                        ring->addInvRef( B->getCoeffConstRef( 0,currentDegree) ) 
                                );
                
                //j=3,5,7,...n:

                //A(j,n-j) = ((n-j+2)A(j-2,n-j+2)-B(j-1,n-j+1))/j
                for(register short j=3;j<=currentDegree;j+=2)
                { 
                        short j_minus_zwei=j-2;
                        A->setCoeff(j,currentDegree-j,
                                        ring->scalarMultiplyRef(        field_ref_m.multInv(typename TRing::ScalarType(j)),
                                                                        ring->addRef(   ring->scalarMultiplyRef( currentDegree - j_minus_zwei,
                                                                                                                                A->getCoeffConstRef( j_minus_zwei, currentDegree-j_minus_zwei )
                                                                                                                ),
                                                                                        ring->addInvRef( B->getCoeffConstRef(j-1,currentDegree-j+1) )
                                                                                )
                                                                )
                                        );
                }
}



template < class TFrommerDefs, int variant >
inline  void Frommer< TFrommerDefs, variant >::perform_the_even_A_Step()
{
                // 1. compute A[n-l,l];  l=1, see Paper
                register typename TRing::ScalarType const * p_a_table_loc       =& (a_table[(maxDegree+1)*currentDegree]);

                typename TRing::ElementType             Summe           = B->getCoeff(currentDegree,0);

                // jetzt in der Schleife alle uebrigen negativen  (was?) aufaddieren:
                for( register short j=2; j<=currentDegree; j+=2 )
                {       
                        
                        ring->addInPlace( Summe, 
                                          ring->addInv  (  ring->scalarMultiplyRef( field_ref_m.multInv(p_a_table_loc[j/2]) ,
                                                                                        B->getCoeffConstRef(currentDegree-j, j)
                                                                                        
                                                                                )   
                                                                )
                                        );
                }
                // Summe = summe *2. - kommentar stimmt nicht. in Wirklichkeit wird berechet: Summe=Summe/2.
                ring->scalarMultiplyInPlace( field_ref_m.multInv(typename TRing::ScalarType(2) ), Summe  );

                A->setCoeff(currentDegree-1, 1, Summe );
                
                // so. jetzt is A[currentDegree-1,1] berechnet.

                //beim naechsten A muss ich ein positives B zweimal dazuaddieren und das Ergebnis mit p_a_table[ ...?
                // vielleicht Tabelle  C[x..,..]/a[i] anlegen
                for(register int j=3; j<=currentDegree-1; j+=2)
                {
                        register short mini=j-1;
                        

                        ring->addInPlace( Summe, 
                                          ring->scalarMultiplyRef(       field_ref_m.multInv(p_a_table_loc[ mini/2 ]),
                                                                 B->getCoeffConstRef(currentDegree-mini, mini)
                                                        
                                                      )
                                        );

                        mini = min( j, currentDegree-j );

                        typename TRing::ElementType     Summe2 = ring->scalarMultiplyRef(field_ref_m.multInv(mini), Summe  );
                        
                        mini--; mini/=2;

                        A->setCoeff(    currentDegree-j, j, 
                                        ring->scalarMultiplyRef( p_a_table_loc[mini], Summe2    ) 
                                   );

                }
                //puhhh, geschafft.

                //------------
                //A(0,n) = 0;
                A->setCoeff( 0, currentDegree,  TRing::ElementType::Zero );

                //j=2,4,6,...n:
                for(register short j=2;j<=currentDegree;j+=2)
                { 
                        register short j_minus_zwei = j-2;

                        //A(j,n-j) = ((n-j+2)A(j-2,n-j+2)-B(j-1,n-j+1))/j
                        A->setCoeff(j,currentDegree-j, 
                                        ring->scalarMultiplyRef(   field_ref_m.multInv(typename TRing::ScalarType(j) ) ,
                                                                        ring->addRef(   ring->scalarMultiplyRef( currentDegree-j_minus_zwei,
                                                                                                                         A->getCoeffConstRef(j_minus_zwei, currentDegree-j_minus_zwei) 
                                                                                                        ),
                                                                                        ring->addInvRef( B->getCoeffConstRef(j-1,currentDegree-j+1) ) 
                                                                                ) 
                                                                )
                                  );
                }
}






template < class TFrommerDefs, int variant >
template <class TemplateMatrix>
inline void Frommer< TFrommerDefs, variant >::outputMatrix(TemplateMatrix *mat) const
{
        std::cerr << mat->getName() << std::endl;
        mat->outputMatrix();
        return;
}

template < class TFrommerDefs, int variant >
inline void Frommer< TFrommerDefs, variant >::outputMatrices() const
{
        std::cerr << "A: " << std::endl;
        A->outputMatrix();
                //getchar();
        
        std::cerr << "dxA: " << std::endl;
                dxA->outputMatrix();
                //getchar();
        
        std::cerr << "dyA: " << std::endl;
                dyA->outputMatrix();
                //getchar();
        
        std::cerr << "B: " << std::endl;
                B->outputMatrix();
                //getchar();
}




 
template < class TFrommerDefs, int variant >
template < int hm >
  void Frommer< TFrommerDefs, variant >::doitx(short howMany)
{
        init(); //- clears b and initializes some Coefficients of A

        #ifdef SAFE
                assert(initialized);
                initialized=false;
                assert(howMany<=maxFocalValuesToCompute && howMany>0);
        #endif


        currFocalValuePos=1;
        
        // Idee: merke mir aktuellen Grad und beim Neustart Daten nur bis currentDegree loeschen!
        currentDegree = 2; 

        while(true)
        {
                #ifdef FROMMERSTAGE_TIMER
                        tim_m.start();
                #endif

                currentDegree++; 
                        
                compute_dyA_and_dxA();

                #ifdef FROMMERSTAGE_TIMER
                        tim2_m.start();
                #endif

                //perform_C_Step_wrapper< hm >(currentDegree);
                perform_C_Step();
                #ifdef FROMMERSTAGE_TIMER
                        tim2_m.stop();
                        frommerTimer_m.logCStepTime(hm,currentDegree-2,tim2_m  );
                        tim2_m.clear();
                        tim2_m.start();
                #endif
                
                

                perform_the_A_Step();

                #ifdef FROMMERSTAGE_TIMER
                        tim2_m.stop();
                        tim_m.stop();   
                        frommerTimer_m.logAStepTime(hm,currentDegree-2, tim2_m );
                        frommerTimer_m.logStepTime(hm, currentDegree-2, tim_m  );
                        tim2_m.clear();
                        tim_m.clear();
                #endif

                
                if ( (currentDegree %2 )==0)//gerade --> Strudelgroesse berechnen
                {
                        computedFocalValues[currFocalValuePos] = ring->add(     A->getCoeffConstRef( 1, currentDegree-1 ), 
                                                                                                B->getCoeffConstRef( 0, currentDegree   )
                                                                                        );

                        ++currFocalValuePos;

                        if (  currFocalValuePos > hm )  break;
                }
        }
        lastDegree=currentDegree;
}


template < class TFrommerDefs, int variant >
template <int depth>
inline void 
Frommer< TFrommerDefs, variant >::perform_C_Step_wrapper(const int stage)
{
        switch (stage)
        {
                case 3:
                         perform_C_Step_differ_stages<depth, 4>();
                        break;  
                case 4:
                         perform_C_Step_differ_stages<depth,4>();
                        break;  
                case 5:
                         perform_C_Step_differ_stages<depth,6>();
                        break;  
                case 6:
                         perform_C_Step_differ_stages<depth,6>();
                        break;  
                case 7:
                         perform_C_Step_differ_stages<depth,8>();
                        break;  
                case 8:
                         perform_C_Step_differ_stages<depth,8>();
                        break;  
                case 9:
                         perform_C_Step_differ_stages<depth,10>();
                        break;  
                case 10:
                         perform_C_Step_differ_stages<depth,10>();
                        break;  
                case 11:
                         perform_C_Step_differ_stages<depth,12>();
                        break;  
                case 12:
                         perform_C_Step_differ_stages<depth,12>();
                        break;  
                case 13:
                         perform_C_Step_differ_stages<depth,14>();
                        break;  
                case 14:
                         perform_C_Step_differ_stages<depth,14>();
                        break;  
                case 15:
                         perform_C_Step_differ_stages<depth,16>();
                        break;  
                case 16:
                         perform_C_Step_differ_stages<depth,16>();
                        break;  
                case 17:
                         perform_C_Step_differ_stages<depth,18>();
                        break;  
                case 18:
                         perform_C_Step_differ_stages<depth,18>();
                        break;  
                case 19:
                         perform_C_Step_differ_stages<depth,20>();
                        break;  
                case 20:
                         perform_C_Step_differ_stages<depth,20>();
                        break;  
                case 21:
                         perform_C_Step_differ_stages<depth,22>();
                        break;  
                case 22:
                         perform_C_Step_differ_stages<depth,22>();
                        break;
                case 23:
                         perform_C_Step_differ_stages<depth,24>();
                        break;  
                case 24:
                         perform_C_Step_differ_stages<depth,24>();
                        break;  
                case 25:
                         perform_C_Step_differ_stages<depth,26>();
                        break;          
                case 26:
                         perform_C_Step_differ_stages<depth,26>();
                        break;
                default:
                        perform_C_Step();
        }
}

template < class TFrommerDefs, int variant >
 
inline void Frommer< TFrommerDefs, variant >::doit(int bComputeAllFocalValues, int jacobi)
//inline void Frommer< TFrommerDefs, variant >::doit(int currNumValuesToCompute, bool bComputeAllFocalValues)
{

        
        init();//- clears b and initializes some Coefficients of A

        

        currFocalValuePos=1;
        
        currentDegree = 2; 
        
        try     // evtl ist try-catch teuer!
        {
                while(true)
                {
                        #ifdef FROMMERSTAGE_TIMER
                                tim_m.start();
                        #endif

                        currentDegree++; 
                

        
                        compute_dyA_and_dxA();

                                #ifdef DEBUG2
                                        this->outputMatrices();
                                #endif
                
                        #ifdef FROMMERSTAGE_TIMER
                                tim2_m.start();
                        #endif

                                perform_C_Step();
                        
                                /*if (jacobi)
                                        perform_C_Step_wrapper<-2>(currentDegree);
                                else if (bComputeAllFocalValues)
                                        perform_C_Step_wrapper<-1>(currentDegree);
                                else
                                        perform_C_Step_wrapper<0>(currentDegree);
                                */

                                #ifdef DEBUG2
                                        outputMatrix(B);
                                #endif

                        #ifdef FROMMERSTAGE_TIMER
                                tim2_m.stop();
                                int hm=-1000;
                                if (jacobi)
                                        hm=-2 ;
                                else if (bComputeAllFocalValues)
                                        hm=maxFocalValuesToCompute ;
                                else
                                        hm = 0 ;
                                frommerTimer_m.logCStepTime(hm,currentDegree-2, tim2_m  );

                                tim2_m.clear();
                                tim2_m.start();
                        #endif

                        perform_the_A_Step();

                                #ifdef DEBUG2
                                        outputMatrix(A);
                                        //getchar();
                                #endif
                
                        // Problem: um A(1, currentDegree-1) zu berechnen braucht man alle B-Werte...
                        // Man kann lediglich die letzte for-Schleife der A-Stufe weglassen...
                        // Aber: man kann die Fälle focalValues=2 ,=3 Spezialisieren 

                        //A(1,n-1) = -B(0,n)
                        #ifdef FROMMERSTAGE_TIMER
                                tim2_m.stop();
                                tim_m.stop();   
                                frommerTimer_m.logAStepTime(hm,currentDegree-2 , tim2_m  );
                                frommerTimer_m.logStepTime( hm, currentDegree-2, tim_m   );
                                tim2_m.clear();
                                tim_m.clear();
                        #endif
                 


                        if (currentDegree % 2 == 0)
                        {
                                //gerade --> Strudelgroesse berechnen

                                typename TRing::ElementType res = ring->add(    A->getCoeffConstRef(1, currentDegree-1), 
                                                                                                B->getCoeffConstRef(0, currentDegree)
                                                                                );

                                computedFocalValues[currFocalValuePos] = res;

                                if (++currFocalValuePos > maxFocalValuesToCompute  )
                                {                                                               
                                        break;
                                }
                                //++currFocalValuePos gehoert ganz ans ende der if-Abfrage !, 
                                if (      (     ( bComputeAllFocalValues==0 )   &&      
                                                ( res.getX() != typename TRing::ElementType(0).getX() ) )
                                )
                                {                                                               
                                        break;
                                        
                                        #ifdef DEBUG
                                                std::cerr << "finished" << std::endl;
                                        #endif
                                }
                                
                        }
                }
                //  merke mir aktuellen Grad um bei einem Neustart möglichst wenig neu zu initialisieren
                lastDegree = currentDegree;
                #ifdef DEBUG2
                        std::cerr << "doit::getChar" << std::endl;
                        //getchar();
                #endif
                /*
                if 
                */
        
        }
        catch (const char* str)
        {
                std::cerr << " error during computation of " << currFocalValuePos << "-th focal value: ";
                std::cerr << std::endl << str << std::endl;
                exit(0);
        }
        
}


template < class TFrommerDefs, int variant >
typename TFrommerDefs::RingType::ElementType
Frommer< TFrommerDefs, variant >::getFocalValue( int pos) const
{
        assert(pos<=maxFocalValuesToCompute );

        if (currFocalValuePos>pos && pos>0)
                return computedFocalValues[pos];
        else 
        {
                std::cerr << "currFocalValuePos" << currFocalValuePos << std::endl;
                throw " requested focal value was not computed!" ;
        }
}


template < class TFrommerDefs, int variant >
void            Frommer< TFrommerDefs, variant >::getComputedFocalValues( vector<typename TRing::ElementType> &computedFV ) const
{
        computedFV.clear();
        computedFV.resize( currFocalValuePos-1 );
        for (int pos=1; pos<currFocalValuePos; pos++)
        {
                computedFV[ pos - 1 ] = computedFocalValues[ pos ] ;
        }
}



template < class TFrommerDefs, int variant >
short Frommer< TFrommerDefs, variant >::getSuccessiveVanishedFocalValuesCount() const
{
        short i=0;
        while (         (i<maxFocalValuesToCompute)
                        &&      (i+1<currFocalValuePos)
                        &&      ((int) computedFocalValues[i+1].getX()==0)
                  )
        {
                i++;
        }

        return i;
}

template < class TFrommerDefs, int variant >
void Frommer< TFrommerDefs, variant >::printStageTimings(std::ostream & os) const
{
        #ifdef FROMMERSTAGE_TIMER
                os  << "--" << name_m << ": " << endl;
                frommerTimer_m.print(os);
        #endif
}