cfRandomExperiment.hpp

Go to the documentation of this file.

// #include "cfRandomExperiment.h"

// nach  #include "cfRandomExperiment.h"  kann KDE Befehle ergänzen!

using namespace nCenterFocus;
using std::list;

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
        CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::CFRandomExperiment(      
                                                                                                                        CenterFocusExperiment<variant> & cfExperiment,
                                                                                                                        D_CenterfocusParams const*  cfparams,
                                                                                                                        TFrommer1               & frommer1, 
                                                                                                                        TFrommer2               & frommer2,
                                                                                                                        RankStatistic           & fullRankStatistic,
                                                                                                                        RankStatistic           & subRankStatistic,
                                                                                                                        LiftAndQuadricsStatistic        & fullQuadricsStatistic,
                                                                                                                        LiftAndQuadricsStatistic        & subQuadricsStatistic,
                                                                                                                        FailedLiftStatistic             &liftStatistic,
                                                                                                                        TPolynomXY              & randomMinusP_polynom, 
                                                                                                                        TPolynomXY              & randomQ_polynom
                                                                                                                        ):
                                                                                                                                cfparams_m(cfparams),
                                                                                                                                cfExperiment_m(cfExperiment),
                                                                                                                                frommer1_m(frommer1),
                                                                                                                                frommer2_m(frommer2),   
                                                                                                                                randomMinusP_polynom_m(randomMinusP_polynom),
                                                                                                                                randomQ_polynom_m(randomQ_polynom),
                                                                                                                                f1_ring_ref_m( frommer1.getRingRef() ),
                                                                                                                                f2_ring_ref_m( frommer2.getRingRef() ),
                                                                                                                                frommer1_field_ref( frommer1_m.getRing()->getFieldRef() ),
                                                                                                                                polynomialRing_m(frommer1.getRingRef() ),
                                                                                                                                st_m(NULL),
                                                                                                                                fullRankStatistic_m(fullRankStatistic),
                                                                                                                                subRankStatistic_m(subRankStatistic),
                                                                                                                                fullQuadricsStatistic_m(fullQuadricsStatistic),
                                                                                                                                subQuadricsStatistic_m(subQuadricsStatistic),
                                                                                                                                liftStatistic_m(liftStatistic),
                                                                                                                                v1_m ( performFormula23Step_construct_v1(  ) ),
                                                                                                                                v2_m ( performFormula23Step_construct_v2(  ) ),
                                                                                                                                coeffRandomVariablesOrder_m(initRandomVariablesOrder() ),
                                                                                                                                randomCounter_m(0)

                                                                                                                        
                                                                                                                        
{
        #ifdef FORMULA23_TIMER

                formula23_compute_bcd_a1_not_zero_count = 0;
                formula23_compute_bcd_a2_not_zero_count = 0;
                formula23_analyze_all_time_count = 0;
                
                formula23_step_time_m=0;
                formula23_compute_a_time_m=0;
                compute_bcd_a1_not_zero_time_m=0;
                compute_bcd_a2_not_zero_time_m=0;
                computeSolutions_time_m=0;
                analyze_all_time_m = 0;
                random_trials_time_m=0;
                
                formula23_compute_a_getFocalValues_time_m=0;
                compute_bcd_a1_not_zero_getFocalValues_time_m=0;
                compute_bcd_a2_not_zero_getFocalValues_time_m=0;
                tim1_m.clear();
                tim2_m.clear();
                tim3_m.clear();

        #endif

        //assert( randomExperimentWellDefined() );
        psolutions_m = new typename TFrommer1::RingType::ScalarType[frommer1_field_ref.getCharacteristic() ];
        pPointSolutions_m = new pair <TPolynomXY,TPolynomXY>[ frommer1_field_ref.getCharacteristic()*frommer1_field_ref.getCharacteristic() ];

        drittel_m = f1_ring_ref_m.getField()->multInv(f1_ring_ref_m.getField()->Convert(3) ) ;
        zwei_m =  f1_ring_ref_m.getField()->Convert(2)   ;  

        initRandomExperiment();

}



template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void    CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::initRandomExperiment(    )
{
        
        // containsEntry-Parameter sind p_or_q, xExp,yExp.
        bool    bCoeff_Q_30isRandom = cfparams_m->getMonomsWithRandomCoefficients().containsEntry(QCoefficient, 3,0) ;
        bool    bCoeff_P_30isRandom = cfparams_m->getMonomsWithRandomCoefficients().containsEntry(PCoefficient, 3,0) ;


        bUseFormula1_m  = false;
        bUseFormula2_m  = false;
        bUseFormula23_m = false;

        
        if (  cfparams_m->useFormula1() ||   cfparams_m->useFormula2() || cfparams_m->useFormula23()  )
        {
                assert( cfparams_m->getPolynomialDegree()==3 ) ; // Formeln stimmen nur fuer Polynome vom Grad 3
        }

        if (  cfparams_m->useFormula1() )
        {
                assert (bCoeff_Q_30isRandom);
                // @todo was ist der Sinn dieses Tests?
                //assert(cfparams_m->requiredVanishedFocalValuesNum()>0);
                bUseFormula1_m = true;
                assert(cfparams_m->requiredVanishedFocalValuesNum()>0);
                if (cfparams_m->getHamiltonianComponentSwitch()==1)
                        assert(cfparams_m->requiredVanishedFocalValuesNum()>1);// anderenfalls ist korrektes Ermitteln  der Statistik nicht möglich
        }

        if (   cfparams_m->useFormula2() )  
        {
                assert (bCoeff_P_30isRandom);
                assert( cfparams_m->useFormula1() );
                //assert(cfparams_m->requiredVanishedFocalValuesNum()>1);
                
                //assert( cfparams_m->getPolynomialDegree()==3 ) ; // was war nochmal der Grund dafür?
                bUseFormula2_m=true;
                assert(cfparams_m->requiredVanishedFocalValuesNum()>1);
                
                if (cfparams_m->getHamiltonianComponentSwitch()==1)
                        assert(cfparams_m->requiredVanishedFocalValuesNum()>2); // anderenfalls ist korrektes Ermitteln  der Statistik nicht möglich
        }       
        
        if ( cfparams_m->useFormula23() )
        {
        
                assert(cfparams_m->useFormula1() );
                //assert(cfparams_m->requiredVanishedFocalValuesNum()>2);
                
                assert(cfparams_m->getFieldCharacteristic()>=3 );
                // pp_{3,0}=>1,qq_{2,1}=>-3,pp_{1,2}=>-3,qq_{0,3}=>1} 
 
                //qq_{3,0}=>1,pp_{2,1}=>3,qq_{1,2}=>-3,pp_{0,3}=>-1} 

                //p_30, q21,p12,q03,p03,p_21,q12,q30 müsssen variabel sein,
                // und nicht nur das, bei diesen Koeffizienten darf ersmal nur '*' gesetzt sein, und nicht 'A'.
                
                //PCoefficient
               const CoeffList & randomCoeffList =  cfparams_m->getMonomsWithRandomCoefficients();
                
         

                assert( cfparams_m->isCoefficientVariable( 3, 0, PCoefficient) );
                // Konstruktor CoeffListEntry parameter : P_or_q, yexp, degree
                assert(randomCoeffList.containsEntry( CoeffListEntry(PCoefficient, 0,3)  ) );
                assert( cfparams_m->isCoefficientVariable( 1, 2, PCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(PCoefficient, 2,3)  ) );
                assert( cfparams_m->isCoefficientVariable( 2, 1, PCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(PCoefficient, 1,3)  ) );
                assert( cfparams_m->isCoefficientVariable( 0, 3, PCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(PCoefficient, 3,3)  ) );

                assert( cfparams_m->isCoefficientVariable( 3, 0, QCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(QCoefficient, 0,3)  ) );
                assert( cfparams_m->isCoefficientVariable( 1, 2, QCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(QCoefficient, 2,3)  ) );
                assert( cfparams_m->isCoefficientVariable( 2, 1, QCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(QCoefficient, 1,3)  ) );
                assert( cfparams_m->isCoefficientVariable( 0, 3, QCoefficient) );
                assert(randomCoeffList.containsEntry( CoeffListEntry(QCoefficient, 3,3)  ) );
                bUseFormula23_m = true;
                
                assert(cfparams_m->requiredVanishedFocalValuesNum()>2);
                if (cfparams_m->getHamiltonianComponentSwitch()==1)
                        assert(cfparams_m->requiredVanishedFocalValuesNum()>3);// anderenfalls ist korrektes Ermitteln  der Statistik nicht möglich

                //v1_m = performFormula23Step_construct_v1(  );
                //v2_m = performFormula23Step_construct_v2(  );
        }

        if (! bUseFormula1_m && !bUseFormula2_m && !bUseFormula23_m) 
                st_m=new DStatistic( cfparams_m->getMaxFocalValuesToCompute(), 0);
        else if (bUseFormula1_m && !bUseFormula2_m && !bUseFormula23_m)
                st_m=new DStatistic( cfparams_m->getMaxFocalValuesToCompute(), 1);
        else if (bUseFormula1_m &&  bUseFormula2_m  )
                st_m=new DStatistic( cfparams_m->getMaxFocalValuesToCompute(), 2);
        else if (bUseFormula1_m &&  bUseFormula23_m  )
                st_m=new DStatistic( cfparams_m->getMaxFocalValuesToCompute(), 3);

        assert(st_m!=NULL);

}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
list<CoeffListEntry>    CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::initRandomVariablesOrder()
{
        list<CoeffListEntry> coeffRandomVariablesOrder;
        list<CoeffListEntry>::const_iterator coeffIt;
        list<CoeffListEntry> allCoeffVariablesOrder = cfparams_m->getCoeffVariablesOrder();
        coeffIt =allCoeffVariablesOrder.begin();
        while (coeffIt!=allCoeffVariablesOrder.end() )
        {
                if (    cfparams_m->isRandomVariable( *coeffIt ) )
                        coeffRandomVariablesOrder.push_back( *coeffIt   ); 
                coeffIt++;
        }
        return coeffRandomVariablesOrder;
}


template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
DStatistic&             CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::getStatisticRef()
{
        return *st_m;
}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
pair<TPolynomXY,TPolynomXY>             CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::performFormula23Step_construct_v1(    )
{

        // TODO: nicht vergessen, dass es sich bei v1.first um das MINUSPPolynom handelt!
        pair<TPolynomXY,TPolynomXY>     v1 ( TPolynomXY ( cfparams_m->getPolynomialDegree() ), TPolynomXY ( cfparams_m->getPolynomialDegree() ) );

        if (cfparams_m->getPolynomialDegree()<3)
                //tuhe nichts
                return v1;

        typename TPolynomXY::CoefficientType coeff;

        coeff.setX(1)  ;
                v1.first.setCoeff(3,0, f1_ring_ref_m.addInv( coeff) );
        //coeff.setX(1)  ;      
                v1.second.setCoeff( 0,3, coeff ) ;
        
        coeff.setX(3) ;
                v1.first.setCoeff(1,2,   coeff  );
        //coeff.setX(3) ;
                v1.second.setCoeff( 2,1,  f1_ring_ref_m.addInv( coeff  ) );

        return v1;
}



template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
pair<TPolynomXY,TPolynomXY>             CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::performFormula23Step_construct_v2(    )
{
        pair<TPolynomXY,TPolynomXY> v2 ( TPolynomXY ( cfparams_m->getPolynomialDegree() ), TPolynomXY ( cfparams_m->getPolynomialDegree() ) );

        if (cfparams_m->getPolynomialDegree()<3)
                //tuhe nichts
                return v2;

        typename TPolynomXY::CoefficientType coeff;


        // TODO: nicht vergessen, dass es sich bei v2.first um das MINUSPPolynom handelt!

        coeff.setX(1);  
                v2.first.setCoeff( 0,3,  coeff );
        //coeff.setX(1)  ;
                v2.second.setCoeff( 3,0,  coeff ) ;
        coeff.setX(3);  
                v2.first.setCoeff( 2,1,  f1_ring_ref_m.addInv(coeff) );
        //coeff.setX(3) ;       
                v2.second.setCoeff( 1,2, f1_ring_ref_m.addInv( coeff )  );

        return v2;
}



template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void            CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::initRandomCoefficients(  TPolynomXY  & p, 
                                                                                                                                        TPolynomXY  & q
                                                                                                                                )
{
        typename TPolynomXY::CoefficientType  value(0);


        list<CoeffListEntry>::const_iterator    coeffIt = coeffRandomVariablesOrder_m.begin();
        while (coeffIt!=coeffRandomVariablesOrder_m.end() )
        {
                value.setX( f1_ring_ref_m.getField()->ConvertScalar(   random(  cfparams_m->getRandomSeedPtr(), 
                                                                                                        f1_ring_ref_m.getCharacteristic() -1 
                                                                                                )   
                                                                                ) 
                        ); 

                if ((*coeffIt).porq == PCoefficient)
                {
                        //p.setCoeff((*coeffIt)->x_exp, (*coeffIt)->y_exp, value );
                        p.setCoeff((*coeffIt).x_exp, (*coeffIt).y_exp, f1_ring_ref_m.addInv(value) );
                }
                else
                {
                        q.setCoeff((*coeffIt).x_exp,(*coeffIt).y_exp, value);
                }
                coeffIt++;
        }
}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
 void           CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::performSingleALLExperiment(TPolynomXY  randomMinusP_polynom,  TPolynomXY randomQ_polynom )
{

        bool finished = false;
        assert( cfparams_m->polynomsHaveCheckAllCoeff() );
        assert( ! bUseFormula23_m );

        while (!finished)
        {
                #ifdef DEBUG
                        TPolynomXY  PE= polynomialRing_m.addInv(randomMinusP_polynom  );
                        cfExperiment_m.getResultPrinterConstRef().printPolynoms(std::cerr,PE, randomMinusP_polynom ,string("Testausgabe: performSingleALLExperiment \n "));
                        getchar();
                #endif

                        
                #ifdef SAFE
                for (int deg = randomMinusP_polynom.getDegree(); deg>=0; deg--)
                {
                        for (int x_exp=deg; x_exp>=0; x_exp--)
                        {
                                assert( (int)randomQ_polynom.getCoeff    (x_exp, deg-x_exp).getEps()==0 );
                                assert( (int)randomMinusP_polynom.getCoeff(x_exp, deg-x_exp).getEps()==0 );
                        }
                }
                #endif
                #ifdef MTCP
                        mtcp_no();
                #endif
                correctStatistic();
                cfExperiment_m.performRegularExperiment(        frommer1_m ,
                                                                        frommer2_m,
                                                                        polynomialRing_m,
                                                                        randomMinusP_polynom, 
                                                                        randomQ_polynom,
                                                                        *st_m, 
                                                                        fullRankStatistic_m, 
                                                                        subRankStatistic_m, 
                                                                        fullQuadricsStatistic_m, 
                                                                        subQuadricsStatistic_m,
                                                                        liftStatistic_m );
//Alle Zahlen an geforderten Koeffizientenstellen durchgehen. (x^ay^b='A')
                
                finished = setNextAllCoefficient( randomMinusP_polynom, randomQ_polynom );
                #ifdef MTCP
                        MtcpManager_g.checkpointingTimeFrame();
                #endif
        }
        return;
}


template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void            CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::performSingleRandomTrial(const TPolynomXY & randomMinusP_polynom,const   TPolynomXY &randomQ_polynom )
{

        if ( cfparams_m->polynomsHaveCheckAllCoeff() )
        {
                performSingleALLExperiment( randomMinusP_polynom, randomQ_polynom );
        }
        else
        {
                #ifdef DEBUG
                        TPolynomXY  PE= polynomialRing_m.addInv(randomMinusP_polynom  );
                        cfExperiment_m.getResultPrinterConstRef().printPolynoms(std::cerr,PE, randomMinusP_polynom ,string("Testausgabe: performSingleRandomExperiment \n "));
                        getchar();
                #endif
                        
                #ifdef SAFE
                for (int deg = randomMinusP_polynom.getDegree(); deg>=0; deg--)
                {
                        for (int x_exp=deg; x_exp>=0; x_exp--)
                        {
                                assert( (int)randomQ_polynom.getCoeff    (x_exp, deg-x_exp).getEps()==0 );
                                assert( (int)randomMinusP_polynom.getCoeff(x_exp, deg-x_exp).getEps()==0 );
                        }
                }
                #endif
                
                // Hier darf kein CorrectStatistic aufgerufen werden!
                cfExperiment_m.performRegularExperiment(        frommer1_m ,
                                                                        frommer2_m,
                                                                        polynomialRing_m,
                                                                        randomMinusP_polynom, 
                                                                        randomQ_polynom,
                                                                        *st_m, 
                                                                        fullRankStatistic_m, 
                                                                        subRankStatistic_m, 
                                                                        fullQuadricsStatistic_m, 
                                                                        subQuadricsStatistic_m,
                                                                        liftStatistic_m );
                
         
        }  
}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
inline void     CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>:: correctStatistic()
{
        #ifdef SAFE
                assert( (bUseFormula2_m &&bUseFormula1_m && !bUseFormula23_m )|| !bUseFormula2_m );
                assert( (bUseFormula23_m &&bUseFormula1_m && !bUseFormula2_m)|| !bUseFormula23_m );
        #endif

        if (bUseFormula1_m  )
        {
                if ( bUseFormula23_m) // 
                {
                        st_m->addFormula1And23Try( cfparams_m->getFieldCharacteristic() );
                                
                }
                else if ( bUseFormula2_m && bUseFormula1_m )
                {
                                
                                st_m->addFormula1And2Try( cfparams_m->getFieldCharacteristic() );
                                        
                }
                else if (       ! bUseFormula2_m && ! bUseFormula23_m)
                        st_m->addFormula1Try( cfparams_m->getFieldCharacteristic() );
                                
        }
}



template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void            CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::performRandomExperiment(  )
{
        // experimentCounter_m=0;

        long64 randomExperimentsMaxNum = 0;     

        if ( cfparams_m->polynomsHaveRandomCoeff() )
        {
                randomExperimentsMaxNum = cfparams_m->getRequestedTrialsNum()-1;
                assert(randomExperimentsMaxNum>=0);
                //cerr << "randomExperimentsMaxNum" << randomExperimentsMaxNum << endl; 
        }
        else
        {
                cerr << "error: no randomCoeff detected" << endl;
                exit(0);
        }
        
        // für Debug-Zwecke kann man damit bewirken, dass die ersten 'startRandomcounter' Versuche nicht durchgeführt werden.
        long64 startRandomcounter = cfparams_m->getRandomCounterStartValue();
        // fuehrt einen Zufaellsversuch durch. Ist auch eventuell  mit einem All-Versuch gekoppelt:
        // fuer jeden Zufallsversuch werden alle All-Koeffizienten ausprobiert.
        for ( randomCounter_m = 0; randomCounter_m <= randomExperimentsMaxNum; randomCounter_m++ )
        {       
        
                if (cfparams_m->polynomsHaveRandomCoeff() ) /* Zufallszahlen setzen */
                {
                        initRandomCoefficients(  randomMinusP_polynom_m, randomQ_polynom_m);
                } 
                
                // 
                if (randomCounter_m >= startRandomcounter )
                {
                        #ifdef MTCP
                                mtcp_no();
                        #endif
                        #ifdef DEBUG
                                TPolynomXY  PE= polynomialRing_m.addInv(randomMinusP_polynom_m  );
                                cfExperiment_m.getResultPrinterConstRef().printPolynoms(std::cerr,PE, randomQ_polynom_m ,string("Testausgabe: vor (compute p, q) \n "));
                                
                                getchar();
                        #endif
                        correctStatistic();
                        //Formel 1 berechnet q(3,0) so, dass die erste Strudelgroesse 0 isttch          
                        if (bUseFormula1_m && ! bUseFormula2_m && ! bUseFormula23_m)
                        {
                                compute_q_coeff_x_3_y_0 (randomMinusP_polynom_m, randomQ_polynom_m  );
                                performSingleRandomTrial(       randomMinusP_polynom_m,
                                                                                randomQ_polynom_m       );
                        }
                        
                        //Formel 2 berechnet p(3,0) so, dass die zweite Strudelgroesse 0 ibtst  
                        if ( bUseFormula2_m && bUseFormula1_m )
                        {
                                        compute_q_coeff_x_3_y_0 (randomMinusP_polynom_m, randomQ_polynom_m  );
                                        // compute_p_coeff_x_3_y_0 calls performSingleRandomTrial.
                                        compute_p_coeff_x_3_y_0(randomMinusP_polynom_m, randomQ_polynom_m );
                        }

                        if ( bUseFormula23_m) // 
                        {
                                //cerr << "bUseFormula23_m step" << endl;
                                compute_q_coeff_x_3_y_0 (randomMinusP_polynom_m, randomQ_polynom_m  );
                                // performFormula23Step has to set number of vanished second focal values.
                                performFormula23Step(randomMinusP_polynom_m, randomQ_polynom_m );
                        }
                        
                        if (! bUseFormula1_m && ! bUseFormula2_m && ! bUseFormula23_m)
                        {

                                //cerr << "noformulastep " << endl;
                                performSingleRandomTrial(       randomMinusP_polynom_m,
                                                                                randomQ_polynom_m       );
                        }
                        #ifdef MTCP
                                MtcpManager_g.checkpointingTimeFrame();
                        #endif
                        // Fortschrittsanzeige!
                        //#ifdef SHOWPROGRESS
                        //if (experimentCounter%schrittweite == 0 && experimentCounter!=0)
                        //      std::cerr << experimentCounter << " . Berechnung durchgefuehrt" << std::endl;
                        //#endif
                }

        } // end for Zufall
}





template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
inline bool     CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::setNextAllCoefficient( TPolynomXY        & allMinusP_polynom, 
                                                                                                                                          TPolynomXY    & allQ_polynom )
{
        bool finished = false;
        
        const CoeffListEntry *          currentMonomWithAllCoefficients = cfparams_m->getMonomsWithAllCoefficients().getFirstConst();

        typename TFrommer1::RingType::ElementType value = TFrommer1::RingType::ElementType::Zero;

        /* zaehle ersten Koeffizienten hoch, wenn der 0 ist, dann zweiten, wenn der 0, dann dritten , ... */
        while ((  value.getX()==0)&&( value.getEps()==0)&& (!finished) ) 
        {
                if (currentMonomWithAllCoefficients->porq == PCoefficient)
                {
                        value = allMinusP_polynom.getCoeff(currentMonomWithAllCoefficients->x_exp, currentMonomWithAllCoefficients->y_exp);
                        
                        f1_ring_ref_m.addInPlaceRef(value, TFrommer1::RingType::ElementType::One); 

                        allMinusP_polynom.setCoeff(currentMonomWithAllCoefficients->x_exp, currentMonomWithAllCoefficients->y_exp, value);
                }
                else
                {
                        value = allQ_polynom.getCoeff(currentMonomWithAllCoefficients->x_exp, currentMonomWithAllCoefficients->y_exp);
                        
                        f1_ring_ref_m.addInPlaceRef(value, TFrommer1::RingType::ElementType::One);

                        allQ_polynom.setCoeff(currentMonomWithAllCoefficients->x_exp, currentMonomWithAllCoefficients->y_exp, value);
                }
                // ähnlich zu einem Übertrag:
                if (( value.getX() == 0) && (  value.getEps()==0)) 
                {//nimm die naechste Zahl und add sie hoch. Wenn du bei der letzten bist, sind wir fertig (ende=1)
                        currentMonomWithAllCoefficients = currentMonomWithAllCoefficients->getNextConst();
                        if ( currentMonomWithAllCoefficients == NULL )
                                finished=true;
                }
        } //end while
        return finished;

}


//OriginalFormel 1
                                /*q(3,0) = -1/3*p(2,0)*p(1,1)           .
                                                 + -1/3*p(1,1)*p(0,2)           .
                                                 + 2/3*p(2,0)*q(2,0)            .
                                                 + 1/3*q(2,0)*q(1,1)            .               
                                                 + -2/3*p(0,2)*q(0,2)           .
                                                 + 1/3*q(1,1)*q(0,2)            .
                                                 + 1/3*p(2,1)                           .
                                                 + p(0,3)                                       .
                                                 + -1/3*q(1,2)                          .*/
template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
inline void     
CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::compute_q_coeff_x_3_y_0(         const   TPolynomXY      & minusP_polynom, 
                                                                                                                        TPolynomXY                      & Q_polynom 
                                                                                                                )
{       
        //typename TRing::ElementType   formel1erg= ring1->Convert( 0 );
        //typename TRing::ElementType   formel1zw = ring1->Convert( 0 );

        
        typename TFrommer1::RingType::ElementType       formel1zw ( minusP_polynom.getCoeff(2,0) );

        f1_ring_ref_m.addInPlaceRef(formel1zw, minusP_polynom.getCoeffConstRef(0,2));
        f1_ring_ref_m.multiplyInPlaceRef(formel1zw, minusP_polynom.getCoeffConstRef(1,1));
        
        f1_ring_ref_m.addInPlaceRef(formel1zw,  minusP_polynom.getCoeffConstRef(2,1) );
        f1_ring_ref_m.addInPlaceRef(formel1zw,  Q_polynom.getCoeffConstRef(1,2) );
        
        f1_ring_ref_m.addInvInPlace(formel1zw);
        // here formel1zw is = -(  mp(2,0)*mp(1,1) +   mp(1,1)*mp(0,2) + mp(2,1)+  q(1,2) )

        f1_ring_ref_m.accMultRef(formel1zw,  Q_polynom.getCoeffConstRef(2,0),  Q_polynom.getCoeffConstRef(1,1) );
        f1_ring_ref_m.accMultRef(formel1zw,  Q_polynom.getCoeffConstRef(1,1),  Q_polynom.getCoeffConstRef(0,2) );
        // here formel1zw is = -(  mp(2,0)*mp(1,1) +   mp(1,1)*mp(0,2) +  q(1,2) ) + q(2,0)*q(1,1) + q(1,1)*q(0,2)

        typename TFrommer1::RingType::ElementType       formel1erg( minusP_polynom.getCoeff(2,0) );
        f1_ring_ref_m.multiplyInPlaceRef(formel1erg,  Q_polynom.getCoeffConstRef(2,0) );
        f1_ring_ref_m.addInvInPlace( formel1erg );
        f1_ring_ref_m.accMultRef(formel1erg,  minusP_polynom.getCoeffConstRef(0,2),  Q_polynom.getCoeffConstRef(0,2) );
        // here formel1erg is =  mp(2,0)*q(2,0) + mp(0,2)*q(0,2)

        f1_ring_ref_m.scalarMultiplyInPlaceRef(zwei_m, formel1erg);
        // here formel1erg is = 2* ( mp(2,0)*q(2,0) + mp(0,2)*q(0,2) )
        f1_ring_ref_m.addInPlaceRef(formel1erg,formel1zw);
        // here formel1erg is = -(  mp(2,0)*mp(1,1) +   mp(1,1)*mp(0,2) +  q(1,2) ) + q(2,0)*q(1,1) + q(1,1)*q(0,2) + 2*( mp(2,0)*q(2,0) + mp(0,2)*q(0,2) )

        //formel1erg = 1/3*formel1erg
        f1_ring_ref_m.scalarMultiplyInPlaceRef(drittel_m, formel1erg);

        //formel1erg = formel1erg  - mp(0,3)
        f1_ring_ref_m.addInPlace(formel1erg, f1_ring_ref_m.addInv( minusP_polynom.getCoeffConstRef(0,3) ) );
        #ifdef SAFE
                assert(formel1erg.getEps()==0);
        #endif

        #ifdef CF_TEST 
                
                for (int i=0; i<f1_ring_ref_m.getCharacteristic(); i++)
                {
                        if ( i != formel1erg.getX() )
                        {
                                //cerr << "test Formula 1" << endl;
                                Q_polynom.setCoeff(3,0, i);
                                frommer1_m.setPolynoms(&(minusP_polynom), &(Q_polynom) )        ;
                                //frommer1_m.doit2();
                                frommer1_m.template  doitx<1>(1);
        
                                assert(frommer1_m.getFocalValue(1)!=TFrommer1::RingType::ElementType::Zero);
        
                        }
                        else
                        {
                                Q_polynom.setCoeff(3,0, i);
                                frommer1_m.setPolynoms(&(minusP_polynom), &(Q_polynom) )        ;
                                //frommer1_m.doit2();
                                frommer1_m.template doitx<1>(1);
        
                                assert(frommer1_m.getFocalValue(1)==TFrommer1::RingType::ElementType::Zero);
                        }
                        
                }
        #endif 

        Q_polynom.setCoeff(3,0, formel1erg);

        return; 
}



// interessant waere der Anteil von Frommer f2 zu Frommer f1, oder?
//ist momentan nur fuer epsPrecision 0 und 1 zulässig! Abfragen!
// Frommer fuer Formel 1 und Frommer fuer Formel2 brauchen etwa gleich lange; warum weiss ich nicht
template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
inline void     CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::
					compute_p_coeff_x_3_y_0(   TPolynomXY              & randomMinusP_polynom, 
                                                                                const TPolynomXY        & randomQ_polynom
                                                                        )
{
        typename TFrommer2::RingType::ElementType       epsilon1(0,1); // = 0+1*eps, ist ein erlaubter Aufruf.

        typename TFrommer2::TPolynomXY randomMinusP_polynom_epsPrecision_1( randomMinusP_polynom.getDegree() ) ; 
        typename TFrommer2::TPolynomXY randomQ_polynom_epsPrecision_1  (      randomQ_polynom.getDegree()  );

        // Kopie erstellen mit epsPrecision 0. (zweiter Parameter)
        cfExperiment_m.copyPolynomWithGivenEpsPrecision(randomMinusP_polynom, randomMinusP_polynom_epsPrecision_1, 0);
        cfExperiment_m.copyPolynomWithGivenEpsPrecision(randomQ_polynom     , randomQ_polynom_epsPrecision_1     , 0);
        

        //#ifdef SAFE
                // gleichgleich1 , damit nicht unnötig genau gerechnet wird.
        assert(f2_ring_ref_m.getEpsPrecision() == 1 );
        //#endif 

        randomMinusP_polynom_epsPrecision_1.setCoeff(3,0, f2_ring_ref_m.addInv( epsilon1 ) ); //muss  nochmal überdacht werden

        frommer2_m.setPolynoms (& randomMinusP_polynom_epsPrecision_1, & randomQ_polynom_epsPrecision_1);
        // Fuehlre den Frommer-Algorithmus  bis Strudelgroesse 2 aus !!!!! 
        //frommer2_m.doit2() ;
        frommer2_m.template doitx<2>(2) ;

        typename TFrommer2::RingType::ElementType focalValue_2 = frommer2_m.getFocalValue(2); // kein Convert - weil die Zahl schon Element des Rings ist!

        //der Koeeffizient von p(3,0) soll -a/b sein, wenn die zweite Strudelgröße = a+b*eps ist
        //das geht natürlich nur, wenn b <> 0 ; wenn es nicht geht, setzen wir es auf 0
        
        typename TFrommer1::RingType::ElementType resCoeff(0);

        if (focalValue_2.getEps() != 0)
        {
                typename TFrommer2::RingType::FieldType const   & f2field_ref = f2_ring_ref_m.getFieldRef();
                
                resCoeff.setX(  (f2field_ref.addInv(    f2field_ref.multiply( focalValue_2.getX(),
                                                                                        f2field_ref.multInv(focalValue_2.getEps() ) 
                                                                                )        
                                                        ) 
                                ).getX() 
                        );
                f1_ring_ref_m.addInvInPlace(resCoeff);
                
                // Statistik: p-1 Werte sind auf Jeden Fall !=0.
        }
        else
        {
                // Statistik: p-1 Werte sind auf Jeden Fall !=0., wen .getX()!=0., Wenn x==0 musst du alle pruefen!
        }


        #ifdef CF_TEST 
                for (int i=0; i<f1_ring_ref_m.getCharacteristic(); i++)
                {
                        if ( typename TFrommer1::RingType::ElementType(i) != resCoeff )
                        {
                                //cerr << "test Formula 2" << endl;
                                randomMinusP_polynom_m.setCoeff(3,0, i);
                                frommer1_m.setPolynoms(&(randomMinusP_polynom), &(randomQ_polynom) )    ;
                                //frommer1_m.doit2();
                                frommer1_m.template doitx<2>(2);
                                if (frommer1_m.getFocalValue(2)==TFrommer1::RingType::ElementType::Zero)
                                {
                                //      cerr << "Formel 2 Sonderfall" << endl;
                                //      cerr << "focalValue_2.getX()" << (int)focalValue_2.getX() << endl;
                                //      cerr << "focalValue_2.getEps()" << (int)focalValue_2.getEps() << endl;
                                }
                                if (focalValue_2 != 0)
                                        assert(frommer1_m.getFocalValue(2)!=TFrommer1::RingType::ElementType::Zero);
                        }
                        else
                        {
                                randomMinusP_polynom_m.setCoeff(3,0, i);
                                frommer1_m.setPolynoms(&(randomMinusP_polynom), &(randomQ_polynom) )    ;
                                //frommer1_m.doit2();
                                frommer1_m.template doitx<2>(2);
                                if (focalValue_2.getEps() != 0)
                                {
                                        assert(frommer1_m.getFocalValue(2)==TFrommer1::RingType::ElementType::Zero);
                                }
                                else
                                {
                                        if (focalValue_2.getX() != 0) // no solutions
                                                assert(frommer1_m.getFocalValue(2)!=TFrommer1::RingType::ElementType::Zero);
                                        // if focalValue_2 ==0 there is no Information and all cases must be tested.
                                        
                                }
                        }
                }

        #endif 
        if (focalValue_2 != 0)
        {
                randomMinusP_polynom_m.setCoeff(3,0,  resCoeff);
                // In die Statistik (p-1)-Werte Eintragen bei Zweite Strudelgroesse =0 und die Anzahl der Versuche mit p multiplizieren?
                // 
                assert(randomMinusP_polynom_m.getCoeff(3,0).getEps()==0);
                performSingleRandomTrial(       randomMinusP_polynom_m, randomQ_polynom_m       );
        }
        else
        {
                int characteristic = f1_ring_ref_m.getCharacteristic();
                for (int i=0; i< characteristic; i++)
                {
                        randomMinusP_polynom_m.setCoeff(3,0, i);
                        frommer1_m.setPolynoms(&(randomMinusP_polynom), &(randomQ_polynom) )    ;
                                //frommer1_m.doit2();
                        frommer1_m.template doitx<2>(2);
                        int solutionCounter=0;
                        if (frommer1_m.getFocalValue(2)== typename TFrommer1::RingType::ElementType(0) )
                        {

                                solutionCounter++;
                                performSingleRandomTrial(       randomMinusP_polynom_m,
                                                                                randomQ_polynom_m       );
                        }
                }
        }

        
        #ifdef DEBUG
                std::cerr << "compute_p_koeff_x_3_y_0 finished" << std::endl ;
                getchar();
        #endif

}


template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
 typename TFrommer1::RingType::ScalarType 
        CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::
		performFormula23Step_getThirdFocalValue       (       const pair< TPolynomXY, TPolynomXY >    &       polynompair  ) const
{
        frommer1_m.setPolynoms(&(polynompair.first), &(polynompair.second) )    ;

        //frommer1_m.doit3();
        frommer1_m.template doitx < 3> ( 3 );
        
        #ifdef SAFE
                typename  TPolynomXY::CoefficientType firstFocalValue = frommer1_m.getFocalValue(1);
                
                assert( firstFocalValue.isZero() );
        #endif

        //typename  TPolynomXY::CoefficientType secondFocalValue = frommer1_m.getFocalValue(2);
        //assert( secondFocalValue.isZero() );
 
        return frommer1_m.getFocalValue(3);
}


template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
 typename TFrommer1::RingType::ScalarType 
        CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::
		performFormula23Step_getSecondFocalValue      (       const pair< TPolynomXY , TPolynomXY > & polynompair ) const
{
        frommer1_m.setPolynoms(&(polynompair.first), &(polynompair.second) )    ;
        //frommer1_m.doit2();
        frommer1_m.template  doitx<2>(2);
        //frommer1_m.doit2();
        #ifdef SAFE
                typename  TPolynomXY::CoefficientType firstFocalValue = frommer1_m.getFocalValue(1);
                assert( firstFocalValue.isZero() );
        #endif
        //typename  TPolynomXY::CoefficientType secondFocalValue = frommer1_m.getFocalValue(2);

        //return secondFocalValue.getX();
        return  frommer1_m.getFocalValue(2);
}


template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
inline  pair<TPolynomXY,TPolynomXY> 
                CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::
							performFormula23Step_computePoint        (typename   TFrommer1::RingType::ScalarType     v1coeff,
                                                                                                         typename    TFrommer1::RingType::ScalarType    v2coeff,
                                                                                                         const pair<const TPolynomXY ,const TPolynomXY   > & PQPolynomPair
                                                                                                        ) const
{
                pair<TPolynomXY,TPolynomXY> ret         (     polynomialRing_m.scalarMultiply( v1coeff, v1_m.first ), 
                                                                        polynomialRing_m.scalarMultiply( v1coeff, v1_m.second) );
        
                polynomialRing_m.addInPlace( ret.first,  PQPolynomPair.first  );
                polynomialRing_m.addInPlace( ret.second, PQPolynomPair.second  );
                
                polynomialRing_m.addInPlace(ret.first,   polynomialRing_m.scalarMultiply( v2coeff, v2_m.first )  );
                polynomialRing_m.addInPlace(ret.second,  polynomialRing_m.scalarMultiply( v2coeff, v2_m.second) );
                 

                return ret;     
}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
inline  void
                CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::
							performFormula23Step_computePointNew(typename   TFrommer1::RingType::ScalarType  v1coeff,
                                                                                                         typename    TFrommer1::RingType::ScalarType    v2coeff,
                                                                                                         const pair<const TPolynomXY ,const TPolynomXY   > & PQPolynomPair,
                                                                                                        pair<TPolynomXY,TPolynomXY> & ret
                                                                                                        ) const
{
                 polynomialRing_m.scalarMultiplyInPlace( v1coeff, ret.first );
                polynomialRing_m.scalarMultiplyInPlace( v1coeff, ret.second) ;
        
                polynomialRing_m.addInPlace( ret.first,  PQPolynomPair.first  );
                polynomialRing_m.addInPlace( ret.second, PQPolynomPair.second  );
                
                polynomialRing_m.addInPlace(ret.first,   polynomialRing_m.scalarMultiply( v2coeff, v2_m.first )  );
                polynomialRing_m.addInPlace(ret.second,  polynomialRing_m.scalarMultiply( v2coeff, v2_m.second) );

}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void    CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::
		performFormula23Step_solveQ(unsigned  int & solutions,typename TFrommer1::RingType::ScalarType  * b  )
{
        
         
        typename TFrommer1::RingType::ScalarType  frommer1_Null(0);

        typename TFrommer1::RingType::ScalarType bp,bq;

        solutions=0;

        if ( b[2] != frommer1_Null)
        {
                //cerr << " b[2] != 0" << endl;
                // -- case 1.1: bb2 <> 0, i.e. really a quadratic equation
                //-- divide by bb2

                // bq = bb0/bb2;
                bq = frommer1_field_ref.multiply(b[0], frommer1_field_ref.multInv( b[2] ) );
                // bp = bb1/bb2/2;
                bp = frommer1_field_ref.multiply(       b[1], 
                                                                frommer1_field_ref.multiply( frommer1_field_ref.multInv( b[2] ),
                                                                                                        frommer1_field_ref.multInv(2)
                                                                                                ) 
                                                        );
                // -> x^2 + 2*bp*x + bq 
                //  -- decide on number of solution
                //        discr = bp^2-bq;
        
                typename TFrommer1::RingType::ScalarType discriminant = frommer1_field_ref.multiply(bp,bp);
                frommer1_field_ref.addInPlace(discriminant, frommer1_field_ref.addInv( bq) );
        
                
                //if discr == 0 then solQ = { -bp };    
                if ( discriminant== frommer1_Null )
                {
                        psolutions_m[0]=frommer1_field_ref.addInv( bp);
                        solutions=1;
                        //cerr << " discriminant==0 " << endl;
                        //res.push_back  (frommer1_field_ref.addInv( bp) );
                        return ;
                }

                typename TFrommer1::RingType::sqrtInf_t sqrt=f1_ring_ref_m.sqrt( discriminant );
        
 
                //if sqrtFp#discr#0 == 2 then (
              // x1 = -bp + sqrtFp#discr#1;
              // x2 = -bp - sqrtFp#discr#1;
              // solQ = {x1,x2} )
               
                if ( sqrt.solutions ==2 )
                {
                        psolutions_m[0]=frommer1_field_ref.add ( frommer1_field_ref.addInv( bp ),
                                                                          sqrt.sqrt 
                                                                        ) ;

                        psolutions_m[1]=frommer1_field_ref.add ( frommer1_field_ref.addInv( bp),
                                                                                frommer1_field_ref.addInv( sqrt.sqrt )
                                                                        ) ;
                        solutions=2;
                        return ;
                }
                //if sqrtFp#discr#0 == 0 then solQ = {};
                #ifdef SAFE
                        assert ( sqrt.solutions ==0 ); // nur der obige Fall sollte noch uebrig sein.
                #endif
                return ;

        }
        else if ( b[1] != frommer1_Null)
        {
                //cerr << " b[1] != 0" << endl;
                //-- case 1.2: bb2=0,bb1<>0
                //-- ie bb0+bb1*x = 0
                // solQ = { -bb0/bb1};
                solutions=1;
                psolutions_m[0]=frommer1_field_ref.addInv ( frommer1_field_ref.multiply( b[0],
                                                                                                                frommer1_field_ref.multInv(b[1] )
                                                                                                        )
                                                                );
                return ;
        }
        else if ( b[0] == frommer1_Null )
        {
                solutions = frommer1_field_ref.getCharacteristic();
                //cerr << " b[0] == 0" << endl;
                //-- case 1.4: bb2=bb1=bb0=0
                //-- ie all x are solutions
                int characteristic=frommer1_field_ref.getCharacteristic();
                for (int i=0;i< characteristic; i++)
                        psolutions_m[i] = i;
                return ;
        }
        
        //if bb2 == 0 and bb1 == 0 and bb0 !=0 then (
        // -- case 1.3: bb2=bb1=0,b0<>0
        //-- ie no solution
        //solQ = {} 
        return ;
}



template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void    CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::performFormula23Step(    TPolynomXY              & minusP_polynom, 
                                                                                        TPolynomXY              & Q_polynom )
{
        #ifdef MTCP
                mtcp_no();
        #endif

        #ifdef FORMULA23_TIMER
                tim1_m.start();
                tim2_m.start();
        #endif

//      cerr <<         "performFormula23Step" << endl;
//      cerr <<         "Q_polynom" << Q_polynom;

        // TODO: wieso Frommer2, und funktioniert das auch immer?

        pair <TPolynomXY,TPolynomXY>    PQpolynom_pair_refs(minusP_polynom, Q_polynom);

        pair <TPolynomXY ,TPolynomXY >  tmppolynoms  (minusP_polynom, Q_polynom);


        typename TFrommer1::RingType::FieldType const & field_ref = frommer1_m.getRing()->getFieldRef();

        typename TFrommer1::RingType::ScalarType a[3];
        
 

        #ifdef FORMULA23_TIMER
                tim3_m.start();
        #endif

        a[0]    = performFormula23Step_getSecondFocalValue(  tmppolynoms );

        #ifdef FORMULA23_TIMER
                tim3_m.stop();
                formula23_compute_a_getFocalValues_time_m += tim3_m.realtime();
                tim3_m.clear();
        #endif

        polynomialRing_m.addInPlace(tmppolynoms.first,          v1_m.first);    
        polynomialRing_m.addInPlace(tmppolynoms.second ,      v1_m.second);


        // a1 = (focalValuesCF(L#(i*3+1)))_1_0-a0;
        #ifdef FORMULA23_TIMER
                tim3_m.start();
        #endif

        a[1]    = field_ref.add ( performFormula23Step_getSecondFocalValue( tmppolynoms ) , field_ref.addInv( a[0] ));
        #ifdef FORMULA23_TIMER
                tim3_m.stop();
                formula23_compute_a_getFocalValues_time_m += tim3_m.realtime();
                tim3_m.clear();
        #endif
        tmppolynoms.first       = polynomialRing_m.add(minusP_polynom, v2_m.first);
        tmppolynoms.second      = polynomialRing_m.add(Q_polynom,      v2_m.second);

        #ifdef FORMULA23_TIMER
                tim3_m.start();
        #endif
        // a2 = (focalValuesCF(L#(i*3+2)))_1_0-a0;

        a[2]    = field_ref.add (performFormula23Step_getSecondFocalValue( tmppolynoms ) , field_ref.addInv( a[0] ));

        #ifdef FORMULA23_TIMER
                tim3_m.stop();
                formula23_compute_a_getFocalValues_time_m += tim3_m.realtime();
                tim3_m.clear();
        #endif
        //cerr <<       "a[0]" << a[0]<<endl;
        //cerr <<       "a[1]" << a[1]<<endl;
        //cerr <<       "a[2]" << a[2]<<endl;
        //cerr <<       "v1_m.mp" << polynomialRing_m.addInv(v1_m.first);
        //cerr <<       "v1_m.q" <<  v1_m.second;
        //cerr <<       "v2_m" << v2_m;

        #ifdef FORMULA23_TIMER
                tim2_m.stop();
                formula23_compute_a_time_m +=  tim2_m.realtime();
                tim2_m.clear();
        #endif


        typename  TFrommer1::RingType::ScalarType b[3], c, d;

        if (a[1] != TFrommer1::RingType::ScalarType::Zero ) 
        {
                
                
                #ifdef FORMULA23_TIMER
                        formula23_compute_bcd_a1_not_zero_count ++;
                        tim2_m.start();
                #endif
        //      cerr << "a[1]!=0   " << endl;
                //then focal3 = b0 + b1*v2 + b2*v2^2

                typename  TFrommer1::RingType::ScalarType v1coeff, v2coeff , ma0_div_a1;

                //vv1 := (vv2) -> -a2/a1*vv2 - a0/a1;

                //vv1(0) 
                ma0_div_a1 = field_ref.multiply(        field_ref.addInv(a[0]), field_ref.multInv(a[1]) );
                v1coeff = ma0_div_a1;
                v2coeff = TFrommer1::RingType::ScalarType::Zero;
                

                pair<TPolynomXY,TPolynomXY> tmppoint (v1_m );
                performFormula23Step_computePointNew(v1coeff, v2coeff , PQpolynom_pair_refs, tmppoint);
        
        //      cerr << "tmppoint.first =  " << tmppoint.first << endl;
        //      cerr << "tmppoint.second =  " << tmppoint.second << endl;
        
                #ifdef FORMULA23_TIMER
                        tim3_m.start();
                #endif
                b[0]= performFormula23Step_getThirdFocalValue(tmppoint);
        //      cerr << "b[0] =  " << b[0] << endl;
                #ifdef FORMULA23_TIMER
                        tim3_m.stop();
                        compute_bcd_a1_not_zero_getFocalValues_time_m += tim3_m.realtime();
                        tim3_m.clear();
                #endif

                //vv1(1)
                v1coeff = field_ref.add(        ma0_div_a1,     
                                        field_ref.multiply(     field_ref.addInv(a[2]),
                                                                        field_ref.multInv(a[1])
                                         ));
                
                v2coeff= typename TFrommer1::RingType::ScalarType(1);

                #ifdef FORMULA23_TIMER
                        tim3_m.start();
                #endif

                tmppoint = v1_m;
                performFormula23Step_computePointNew(v1coeff, v2coeff , PQpolynom_pair_refs, tmppoint);

                //-- c = focal3(1) = b0+b1+b2
                c= performFormula23Step_getThirdFocalValue(tmppoint);
                #ifdef FORMULA23_TIMER
                        tim3_m.stop();
                        compute_bcd_a1_not_zero_getFocalValues_time_m += tim3_m.realtime();
                        tim3_m.clear();
                #endif
                //vv1(-1)
                v1coeff = field_ref.add(        ma0_div_a1,     
                                                        field_ref.multiply(      a[2],
                                                                                        field_ref.multInv(a[1])
                                                                        )
                                        );

                v2coeff= field_ref.addInv(v2coeff);

                tmppoint = v1_m;
                performFormula23Step_computePointNew(v1coeff, v2coeff , PQpolynom_pair_refs, tmppoint);

                #ifdef FORMULA23_TIMER
                        tim3_m.start();
                #endif
                //-- d = focal3(-1) = b0-b1+b2
                d= performFormula23Step_getThirdFocalValue(tmppoint);

                #ifdef FORMULA23_TIMER
                        tim3_m.stop();
                        compute_bcd_a1_not_zero_getFocalValues_time_m += tim3_m.realtime();
                        tim3_m.clear();
                 
                        tim2_m.stop();
                        compute_bcd_a1_not_zero_time_m +=  tim2_m.realtime();
                        tim2_m.clear();
                #endif
 
        }else if  ( a[2] != TFrommer1::RingType::ScalarType::Zero) 
        {

                
                #ifdef FORMULA23_TIMER
                        formula23_compute_bcd_a2_not_zero_count++;
                        tim2_m.start();
                #endif
                //cerr << "a[1]==0, a[2]!=0  " << endl;
                //vv2 = -a0/a2;
                //    focal3(point,0,vv2),
                //-- c = focal3(1) = b0+b1+b2
                //focal3(point,1,vv2),
                //-- d = focal3(-1) = b0-b1+b2
                //focal3(point,-1,vv2)

                typename  TFrommer1::RingType::ScalarType v1coeff, v2coeff ;
                
                //vv2 = -a0/a2;
                v2coeff = field_ref.multiply(   field_ref.addInv(a[0]),
                                                        field_ref.multInv(a[2]) 
                                        );

                v1coeff = TFrommer1::RingType::ScalarType::Zero;


                pair<TPolynomXY,TPolynomXY>     tmppoint ( v1_m );
                performFormula23Step_computePointNew(v1coeff, v2coeff , PQpolynom_pair_refs, tmppoint);
                #ifdef FORMULA23_TIMER
                        tim3_m.start();
                #endif

                b[0] = performFormula23Step_getThirdFocalValue(tmppoint);
                #ifdef FORMULA23_TIMER
                        tim3_m.stop();
                        compute_bcd_a2_not_zero_getFocalValues_time_m += tim3_m.realtime();
                        tim3_m.clear();
                #endif
                v1coeff= typename TFrommer1::RingType::ScalarType(1);


                tmppoint = v1_m;
                performFormula23Step_computePointNew(v1coeff, v2coeff , PQpolynom_pair_refs, tmppoint);

                #ifdef FORMULA23_TIMER
                        tim3_m.start();
                #endif
                //-- c = focal3(1) = b0+b1+b2
                c = performFormula23Step_getThirdFocalValue(tmppoint);
                #ifdef FORMULA23_TIMER
                        tim3_m.stop();
                        compute_bcd_a2_not_zero_getFocalValues_time_m += tim3_m.realtime();
                        tim3_m.clear();
                #endif
                v1coeff= field_ref.addInv(v1coeff);

                tmppoint = v1_m;
                performFormula23Step_computePointNew(v1coeff, v2coeff , PQpolynom_pair_refs, tmppoint);


                #ifdef FORMULA23_TIMER
                        tim3_m.start();
                #endif
                //-- d = focal3(-1) = b0-b1+b2
                d = performFormula23Step_getThirdFocalValue(tmppoint);
                #ifdef FORMULA23_TIMER
                        tim3_m.stop();
                        compute_bcd_a2_not_zero_getFocalValues_time_m += tim3_m.realtime();
                        tim3_m.clear();
                #endif

                #ifdef FORMULA23_TIMER
                        tim2_m.stop();
                        compute_bcd_a2_not_zero_time_m +=   tim2_m.realtime();
                        tim2_m.clear();
                #endif
        }
        else if  ( a[0] != TFrommer1::RingType::ScalarType::Zero) 
        {

                //cerr << "a[0]!=0  " << endl;
                // a1=a2=0 a0<>0. No solutions
                // also weder f2 noch f3 =0. 
        }
        unsigned int pPointSolutions=0;

        if  ( (a[1] != TFrommer1::RingType::ScalarType::Zero ) || (a[2] != TFrommer1::RingType::ScalarType::Zero )  )
        {
                #ifdef FORMULA23_TIMER
                        tim2_m.start();
                #endif

                st_m->correctSecondFocalValueStatistic( field_ref.getCharacteristic() );
                //cerr << "a[1]!=0 || a[2] !=0" << endl;
                //b1 = (c-d)/2;
                //  b2 = (c+d)/2-b0;
                b[1] = field_ref.multiply( field_ref.multInv(2),
                                        field_ref.add(  c,
                                                        field_ref.addInv(d)
                                                )
                                );

                b[2] = field_ref.add (field_ref.multiply(       field_ref.multInv(2),
                                                                        field_ref.add(  c,
                                                                                                d
                                                                                        )
                                                                ),
                                field_ref.addInv(b[0])
                        );

                //      cerr <<         "b[0]" << b[0]<<endl;
                //      cerr <<         "b[1]" << b[1]<<endl;
                //      cerr <<         "b[2]" << b[2]<<endl;
                unsigned int solutions=0;
                // Antworten stehen in psolutions_m
                  performFormula23Step_solveQ( solutions, & (b[0]) );
                pPointSolutions = solutions;
                //cerr << "rr.size" << rr.size()<< endl;
                if  ( a[1] != TFrommer1::RingType::ScalarType::Zero )  
                {

                        //cerr << "a[1]!=0  " << endl;
                        /*if a1 != 0 then (
                        -- solve for v1
                        vv1 := (vv2) -> -a2/a1*vv2 - a0/a1;
                        -- eliminate v1 
                        -- then focal3 = b0 + b1*v2 + b2*v2^2
                        solutions = apply(solveQ(b0,b1,b2),i->(vv1(i),i));
                        //apply(solutions,i->point+i#0*v1+i#1*v2)
                -- end for a1 <> 0 
                        );*/
                        typename TFrommer1::RingType::ScalarType ma_2div_a1, ma0_div_a1;
                        ma0_div_a1=field_ref.addInv(field_ref.multiply(a[0],
                                                                                                 field_ref.multInv( a[1] )
                                                                                )
                                                        );
                        ma_2div_a1=field_ref.addInv(field_ref.multiply(a[2],
                                                                                                 field_ref.multInv( a[1] )
                                                                                )
                                                        );

                        for (size_t i=0; i<solutions; i++ )
                        {
                                        typename  TFrommer1::RingType::ScalarType v1coeff_= field_ref.add (     ma0_div_a1, 
                                                                                                                                        field_ref.multiply(     ma_2div_a1, 
                                                                                                                                                                        psolutions_m[i]
                                                                                                                                                                )
                                                                                                                                );
                                
                                //cerr << "Solutions: " << endl;
                                //cerr << "v1coeff_ =  " << v1coeff_ << endl;
                                //cerr << " rr[i] =  " << rr[i] << endl;
                                
                                pPointSolutions_m[i] = v1_m;
                                performFormula23Step_computePointNew(v1coeff_,  psolutions_m[i] , PQpolynom_pair_refs, pPointSolutions_m[i] );

                                //cerr << "solution.first =  " << tmppoint.first << endl;
                                //cerr << "solution.second =  " << tmppoint.second << endl;
                        }
                        
                }
                else if  ( a[2] != TFrommer1::RingType::ScalarType::Zero )  
                {
                        //cerr << "a[2]!=0  " << endl;
                        /*
                        -- case 2: a1==0, a2<>0
                        if (a1==0) and (a2!=0) then (
                                vv2 = -a0/a2;
                        -- eliminate v2 
                        -- then focal3 = b0 + b1*v1 + b2*v1^2 
                        solutions = apply(solveQ(b0,b1,b2),i->(i,vv2));
                        );*/


                        typename  TFrommer1::RingType::ScalarType v2coeff_ = field_ref.addInv(field_ref.multiply(       a[0], 
                                                                                                                                                        field_ref.multInv(a[2])
                                                                                                                                        )
                                                                                                                );
                        for (size_t i=0;i<solutions;i++)
                        {
                                        pPointSolutions_m[i]=  v1_m;
                                        performFormula23Step_computePointNew( psolutions_m[i], v2coeff_, PQpolynom_pair_refs, pPointSolutions_m[i]) ;   
                        }
                
                
                }
                #ifdef FORMULA23_TIMER
                        tim2_m.stop();
                        computeSolutions_time_m         +=  tim2_m.realtime();
                        tim2_m.clear();
                #endif
        }
        else if (a[0] == TFrommer1::RingType::ScalarType::Zero )        //(a2==0) and (a1==0) and (a0==0)
        {

                #ifdef FORMULA23_TIMER
                        formula23_analyze_all_time_count ++;
                        tim2_m.start();
                #endif

                //cerr << "a[0]==0" << endl;
                //cerr << "minusP_polynom" << polynomialRing_m.addInv(minusP_polynom) << endl;
                //cerr << "Q" << Q_polynom << endl;
 
                //-- in this case the solutions are not really solutions
                //-- they have to be analysed even for focal3=0. Macaulay-code:
                //solutions = flatten apply(p,i->apply(p,j->(sub(i,Fp),sub(j,Fp))));
                //apply(solutions,i->point+i#0*v1+i#1*v2)
                
                int pos=0;
                int characteristic = field_ref.getCharacteristic();
                for (int i=0; i< characteristic; i++)
                for (int j=0; j< characteristic; j++)
                        {
                                pair <TPolynomXY ,TPolynomXY > tmppoint  (v1_m );
                                 performFormula23Step_computePointNew( i,j, PQpolynom_pair_refs, tmppoint) ;    

                                frommer1_m.setPolynoms(&(tmppoint.first), &(tmppoint.second) )  ;
                        
                                frommer1_m.template doitx<3>(3);
                                typename  TFrommer1::RingType::ScalarType s2,s3;
                                        
                                s2=frommer1_m.getFocalValue(2);
                                s3=frommer1_m.getFocalValue(3);

                                if (s2==TFrommer1::RingType::ScalarType::Zero && s3==TFrommer1::RingType::ScalarType::Zero)
                                {
                                        pPointSolutions_m[pos] = tmppoint ;
                                        pos++;
                                        pPointSolutions++;
                                }
                                if (s2==TFrommer1::RingType::ScalarType::Zero) 
                                        st_m->correctSecondFocalValueStatistic(1);
                        }
                #ifdef FORMULA23_TIMER
                        tim2_m.stop();
                        analyze_all_time_m              +=  tim2_m.realtime();
                        tim2_m.clear();
                #endif
        }
        //cerr << "solutions.size()"<<solutions.size()<< endl;
        #ifdef CF_TEST
                int countSecondFocalValues=0;
                int characteristic = field_ref.getCharacteristic();

                for (int i=0; i< characteristic; i++)
                {       for (int j=0; j< characteristic; j++)
                        {
                                pair <TPolynomXY ,TPolynomXY > tmppoint  (v1_m );
                                performFormula23Step_computePointNew( i,j, PQpolynom_pair_refs, tmppoint) ;     

                                bool isSolution=false;
                                for (size_t currSolution=0; currSolution<pPointSolutions; currSolution++)
                                {
                                        if ( tmppoint==pPointSolutions_m[currSolution]  )
                                                isSolution=true;
                                }
                                typename  TFrommer1::RingType::ScalarType s2,s3;
                                s2 = performFormula23Step_getSecondFocalValue(tmppoint);
                                s3 = performFormula23Step_getThirdFocalValue(tmppoint);
                                if (s2==TFrommer1::RingType::ElementType::Zero) countSecondFocalValues++;
                                if (!isSolution)
                                {
                                        //solutionscomplement.append(tmppoint);
                                        assert( s2!=TFrommer1::RingType::ElementType::Zero || s3 !=TFrommer1::RingType::ElementType::Zero );
                                }
                                else
                                {
                                        //cerr << "solution test" << endl;
                                        assert( s2==TFrommer1::RingType::ElementType::Zero && s3 == TFrommer1::RingType::ElementType::Zero );
                                }
                        }
                }
                if (a[1] != TFrommer1::RingType::ElementType::Zero || a[2] != TFrommer1::RingType::ElementType::Zero )
                {
                        assert(countSecondFocalValues == field_ref.getCharacteristic() );
                }
                
        
        #endif

        #ifdef FORMULA23_TIMER
                tim1_m.stop();
                formula23_step_time_m +=  tim1_m.realtime() ;
                tim1_m.clear();
                tim1_m.start();
        #endif

        for (size_t i=0; i<pPointSolutions; i++)
        {
                performSingleRandomTrial(       pPointSolutions_m[i].first,
                                                                                pPointSolutions_m[i].second     );
        }

        #ifdef FORMULA23_TIMER
                tim1_m.stop();
                random_trials_time_m    +=  tim1_m.realtime();
                tim1_m.clear();
        #endif

        #ifdef MTCP
                MtcpManager_g.checkpointingTimeFrame();
        #endif

}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
void    CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::printTimings(std::ostream & os) const
{
        #ifdef FORMULA23_TIMER
                os << "---timer_formula23_step: "                       <<formula23_step_time_m << " sec." << std::endl;
                os << "---timer_formula23_compute_a: "          << formula23_compute_a_time_m << " sec." << std::endl;
                os << "---timer_compute_bcd_a1_not_zero: "      << compute_bcd_a1_not_zero_time_m << " sec." << std::endl;
                os << "---timer_compute_bcd_a2_not_zero: "      << compute_bcd_a2_not_zero_time_m << " sec." << std::endl;
                os << "---timer_computeSolutions: "             << computeSolutions_time_m << " sec." << std::endl;
                os << "---timer_analyze_all: "                  << analyze_all_time_m << " sec." << std::endl;
                os << "---timer_random_trials: "                        << random_trials_time_m << " sec." << std::endl;

                os << "---timer_formula23_compute_a_getFocalValues: "           << formula23_compute_a_getFocalValues_time_m << "sec." << std::endl;
                os << "---timer_compute_bcd_a1_not_zero_getFocalValues_time: "  << compute_bcd_a1_not_zero_getFocalValues_time_m << "sec." << std::endl;
                os << "---timer_compute_bcd_a2_not_zero_getFocalValues_time: "  << compute_bcd_a2_not_zero_getFocalValues_time_m << "sec." << std::endl;

                
                os << "---compute_bcd_a1_not_zero_count: \t"    << formula23_compute_bcd_a1_not_zero_count << " " << std::endl;
                os << "---compute_bcd_a2_not_zero_count: \t"    << formula23_compute_bcd_a2_not_zero_count << "" << std::endl;
                os << "---analyze_all_count: \t\t"                      << formula23_analyze_all_time_count << "" << std::endl;


        #endif

}

template <int variant,  class TPolynomXY, class  TFrommer1, class  TFrommer2>
        CFRandomExperiment<variant, TPolynomXY, TFrommer1, TFrommer2>::~CFRandomExperiment()
{
        if (st_m!=NULL) delete st_m;
        delete[] psolutions_m;
        delete[] pPointSolutions_m;
        
        #ifdef FORMULA23_TIMER
                MtcpManager_g.disconnectTimer(tim1_m);
                MtcpManager_g.disconnectTimer(tim2_m);
                MtcpManager_g.disconnectTimer(tim3_m);
        #endif
}