fast_frommer.cpp

Go to the documentation of this file.
00001 
00018 #include <assert.h>
00019 
00020 
00021 template <class TFrommerDefs, int variant>
00022 const  typename TFrommerDefs::RingType * fast_Frommer< TFrommerDefs, variant>::getRing() const
00023 {
00024         return ring1;
00025 };
00026 
00027 template <class TFrommerDefs, int variant>
00028 const  typename TFrommerDefs::RingType & fast_Frommer< TFrommerDefs, variant>::getRingRef() const
00029 {
00030         return *ring1;
00031 };
00032 
00034 // Konstruktion/Destruktion
00036 template <class TFrommerDefs, int variant>
00037 fast_Frommer< TFrommerDefs, variant>::
00038 fast_Frommer(short _maxFocalValuesToCompute,     const TRing *_ring                     ) :
00039                                                                         name_m(""),
00040                                                                         ring1                           ( _ring),
00041                                                                         field_ref_m(_ring->getFieldRef()),
00042                                                                         maxFocalValuesToCompute ( _maxFocalValuesToCompute ),
00043                                                                         maxDegree                       ( maxFocalValuesToCompute *2 + 3 ),
00044                                                                         pq                              ( NULL ), 
00045                                                                         A                               ( new TMatrixXYType("A",maxDegree)),
00046                                                                         dA                              ( new TMatrixXYPairType ("dA",maxDegree) ),
00047                                                                         B                               ( new TMatrixXYType ("B",maxDegree) ),
00048                                                                         B0                              ( new  TMatrixXY_0Type ("B",maxDegree) ),
00049                                                                         a_table                 ( createLowATable(maxDegree, _ring) ),
00050                                                                         currFocalValuePos(0)
00051 
00052                 #ifdef FROMMERSTAGE_TIMER
00053                         ,frommerTimer_m(_maxFocalValuesToCompute)
00054                 #endif
00055 {
00056         #ifdef SAFE
00057                  initialized=false;
00058         #endif
00059 
00060         #ifdef FROMMERSTAGE_TIMER
00061                 tim_m.clear();
00062                 tim2_m.clear();
00063         #endif
00064 
00065         minusp = NULL;
00066         q = NULL;
00067         lastDegree = maxDegree;
00068         computedFocalValues = new typename TRing::ElementType[ maxFocalValuesToCompute+1 ];
00069 
00070         A->setCoeff (2,0, TRing::ElementType::One);
00071         A->setCoeff (1,1, TRing::ElementType::Zero);
00072         A->setCoeff (0,2, TRing::ElementType::One);
00073 
00074         wellDefined();
00075 }
00076 
00077 template <class TFrommerDefs, int variant>
00078 void fast_Frommer< TFrommerDefs, variant>::wellDefined()
00079 {
00080 
00081         assert ( ring1->getLookupModuloTableSize()>14 );
00082 
00084         TPolynomXY * testp = new TPolynomXY(3);
00085 
00086         register typename TPolynomXY::CoefficientType * p_ElementAddr = testp->getCoeffAddr(0,2);
00087         
00088         //p_ElementAddr += 2; 
00089         //assert (testp->getCoeffAddr(3,0)==p_ElementAddr);
00090 
00091 
00092         for (int degree = 0; degree<=3; degree++)
00093         {
00094                 p_ElementAddr = testp->getCoeffAddr(degree,0);
00095                 for (int yExp = 0; yExp<=degree; yExp++)
00096                 {
00097                         if ( testp->getCoeffAddr(degree-yExp, yExp) !=p_ElementAddr )
00098                         {
00099                                 cerr << "degree" << degree << endl;
00100                                 cerr << "degree-yExp" << degree-yExp << endl;
00101                                 cerr << " yExp" << yExp << endl;
00102                                 cerr << "testp->getCoeffAddr(degree-yExp, yExp)" << testp->getCoeffAddr(degree-yExp, yExp) << endl;
00103                                 cerr << "p_ElementAddr" << p_ElementAddr << endl;
00104                         }
00105                         assert( testp->getCoeffAddr(degree-yExp, yExp) ==p_ElementAddr );
00106                         p_ElementAddr++;
00107                 }
00108         }
00109 
00110 
00112         for (int degree = 0; degree<=maxDegree; degree++)
00113         {
00114                 register typename TRing::ElementType * p_A = A->getCoeffAddr(degree,0);
00115                 for (int yExp = 0; yExp<=degree; yExp++)
00116                 {
00117                         assert( A->getCoeffAddr(degree-yExp, yExp) ==p_A );
00118                         p_A++;
00119                 }
00120         }
00121 
00123         for (int degree = 0; degree<=maxDegree; degree++)
00124         {
00125                 register typename TRing::ElementType * p_dA = dA->getCoeffAddr(degree,0);
00126                 for (int yExp = 0; yExp<=degree; yExp++)
00127                 {
00128                         assert( dA->getCoeffAddr(degree-yExp, yExp) ==p_dA );
00129                         p_dA++;
00130                         p_dA++;
00131                 }
00132         }
00133 
00134         delete testp;
00135 }
00136 
00137 template <class TFrommerDefs, int variant>
00138 fast_Frommer< TFrommerDefs, variant>::~fast_Frommer()
00139 {
00140         delete[] computedFocalValues;
00141         delete  B;
00142         delete  A;
00143         delete  dA;
00144         delete[] a_table;
00145 
00146         #ifdef FROMMERSTAGE_TIMER
00147                 MtcpManager_g.disconnectTimer(tim_m);
00148                 MtcpManager_g.disconnectTimer(tim2_m);
00149         #endif
00150 
00151  
00152 }
00153 
00154 
00155 
00166 template <class TFrommerDefs, int variant>
00167 const typename TFrommerDefs::RingType::ScalarType * const  fast_Frommer< TFrommerDefs, variant>::createLowATable(int maxDegree,const TRing *ring1 )
00168 {
00169         short i,j,h;
00170 
00171         int tableSize=(maxDegree+1)*(maxDegree+1) ;
00172 
00173          typename TRing::ScalarType * local_a_table = new  typename TRing::ScalarType[ tableSize ];
00174 
00175         for ( h=0; h<tableSize; h++)
00176         {
00177                 local_a_table[h ] = TRing::ScalarType::Zero;
00178         }
00179 
00180         //for h from 0 to (alphaTableRowNum-1) do
00181         for ( h=0; h<=maxDegree; h++)
00182         {
00183                 //alphaTable_(h,0)= 1_cfRing; 
00184                 local_a_table[h*(maxDegree+1)]= TRing::ScalarType::One;
00185                         
00186                 for ( i=1; i<= h/2; i++)
00187                 {
00188                          typename TRing::ScalarType  Produkt=TRing::ScalarType::One;
00189                         for ( j=1; j<=i; j++)
00190                         {
00191                                 ring1->getField()->multiplyInPlace(Produkt, 
00192                                                                                 ring1->getField()->multiply( h - (2*j-1), 
00193                                                                                                                         ring1->getField()->multInv(2*j-1)
00194                                                                                                                 ) 
00195                                                                         );
00196                         }
00197                         local_a_table[ h*(maxDegree+1)+i ] = Produkt;
00198                 }
00199                 
00200         }
00201         return local_a_table;
00202 }
00203 
00204 
00205 
00212 template <class TFrommerDefs, int variant>
00213 inline void fast_Frommer< TFrommerDefs, variant>::init()
00214 {
00215 
00216         #ifdef SAFE
00217                 assert(A->getCoeff(2,0) == TRing::ElementType::One);
00218                 assert(A->getCoeff(1,1) == TRing::ElementType::Zero);
00219                 assert(A->getCoeff(0,2) == TRing::ElementType::One);
00220                  
00221                 A->clear();
00222                 //B->clear();
00223                 B->clear(lastDegree); //notwendig bei Wiederverwendung!
00224                 B0->clear(lastDegree);
00225                 dA->clear();
00226 
00227                 A->setCoeff (2,0, TRing::ElementType::One);
00228                 A->setCoeff (1,1, TRing::ElementType::Zero);
00229                 A->setCoeff (0,2, TRing::ElementType::One);
00230         #else
00231                 //A->clear();             // nicht mehr notwendig !
00232                 //dxA->clear(lastDegree); // nicht mehr notwendig!
00233                 //dyA->clear(lastDegree); //  nicht mehr notwendig!
00234                 B->clear(lastDegree); //notwendig bei Wiederverwendung! // wuerde bei linearer Speicherung auch schneller gehen!
00235                 B0->clear(lastDegree);
00236         #endif
00237         // A is intitialized in constructor.
00238 }
00239 
00240 
00241 
00242 
00254 template <class TFrommerDefs, int variant>
00255 void fast_Frommer< TFrommerDefs, variant>::setPolynoms( const TPolynomXY * const minuspp,
00256                                          const TPolynomXY * const qq) 
00257 {
00258         //1. delete computedFocalValues: (not nessecary, they are overwritten later)!
00259         
00260         #ifdef SSAFE
00261 
00262                 assert( minuspp->getDegree() == qq->getDegree() );
00263 
00264                 if (minusp!=NULL)
00265                         assert( minusp->getDegree() == qq->getDegree() );
00266                 initialized=true;
00267         #endif
00268         #ifdef SAFE
00269         initialized=true;
00270         #endif
00271         q      = qq;
00272         minusp = minuspp;
00273 
00274         //      init();
00275 };
00276 
00277 
00291 template <class TFrommerDefs, int variant>
00292  inline void fast_Frommer< TFrommerDefs, variant>::compute_dyA_and_dxA()
00293 {
00294         register short currentDegree_minus_2_minus_i=0;
00295 
00296         register typename TRing::ElementType * p_dA = dA->getCoeffAddr(currentDegree-2, currentDegree_minus_2_minus_i);
00297 
00298         register typename TRing::ElementType * pA = A->getCoeffAddr(currentDegree-1, currentDegree_minus_2_minus_i);
00299 
00300 
00301         for (register short i = currentDegree-2; i>=0;  i--, currentDegree_minus_2_minus_i++ )
00302         {
00303                  
00304                 *p_dA = ring1->scalarMultiply( i+1,*pA );
00305 
00306                 p_dA++;
00307                 pA++;
00308                  
00309                 *p_dA = ring1->scalarMultiply ( currentDegree_minus_2_minus_i+1 ,*pA );
00310                 p_dA++;
00311         }
00312 }
00313 
00314 
00315 
00337 template <class TFrommerDefs, int variant>
00338 inline void fast_Frommer< TFrommerDefs, variant>::perform_inner_C_Step( register typename TRing::ElementType * const b0,
00339                                                         register const typename TRing::ElementType * const minusp0,
00340                                                         register const typename TRing::ElementType * const q0,
00341                                                         register const typename TRing::ElementType * const da0,
00342                                                          short _n )
00343 {
00344         // m läuft, j,l fest, n gleich (currentDegree-m-j-l) folgt aus Bedingung (m+n+j+l==currentDegree).
00345         // Berechnungen finden nur statt,
00346         //  wenn (m>=0), (n>=0) und  (m+n)>0 und (m+n)<=currentDegree-2
00347 
00348         // Wenn eine Schleife bei m=1 statt m=0 startet , ist Summe(m,n) >0,
00349         // also kann man sich den Vergleich "(m+n)>0" (und viel Rechenzeit ) sparen.
00350 
00351         //-------------------
00352         // Starte bei m=0:
00353         //register short n=currentDegree-j-l;  //hiermit ist  m+n+j+l==currentDegree für alle nachfolgenden Schritte {m++;n--} erfüllt
00354         // erforderliche Tests (m+n)>0 UND (m+n)<=currentDegree-2 werden mit n<=(currentDegree-2), m=0 erfüllt
00355         //if (n<=currentDegree-2) // Da die nächsten Schritte immer (m++;n--) lauten, reicht diese einmalige Prüfung. 
00356                                 //Es muss in jedem nachfolgenden Schritt nur noch (n>=0) getestet werden!
00357         #ifdef SAFE
00358                 assert(_n<=currentDegree-2);
00359                 assert(_n>0);
00360         #endif 
00361         register short          n = _n;
00362         register  typename TRing::ElementType const *    datmp   = da0;
00363         register  typename TRing::ElementType   polCoeff = *minusp0; // hole p, pointer++, hole q, pointer++
00364         register  typename TRing::ElementType *   b1tmp  = b0;
00365 
00366         /* SSE-Konforme Variante: so wird die Schleife 'countable':
00367         for (register int i=_n;i>=0;i--)
00368         {
00369                 datmp++;
00370                 ring1->accMultSpec(b1tmp, polCoeff, datmp++);   //ACCUMULATE MULTIPLICATION: b = b + p(j,l) * dya(m,n)
00371                 ++b1tmp;
00372         }*/
00373         do 
00374         {
00375                 datmp++;
00376                 //ring1->accMultSpec(b1tmp, polCoeff, datmp++); //ACCUMULATE MULTIPLICATION: b = b + p(j,l) * dya(m,n)
00377                 ring1->accMultRef(*b1tmp,polCoeff, *datmp);
00378                 datmp++;
00379                 ++b1tmp;
00380                 --n;
00381         }
00382         while (n>=0);
00383 
00384         // separated into two loops, because it is faster. (less local variables for the CPU or more cache hits?)
00385         n        = _n;
00386         b1tmp    = b0;
00387         polCoeff = *q0; 
00388         datmp    = da0;
00389 
00390         /* so wird die schleife countable:
00391         for (register int i=_n;i>=0;i--)
00392         {
00393                 ring1->accMultSpec(b1tmp,polCoeff, datmp++);    //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00394                 datmp++;
00395                 ++b1tmp;
00396         }*/
00397 
00398         do 
00399         {
00400                 //ring1->accMultSpec(b1tmp,polCoeff, datmp++);  //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00401                 ring1->accMultRef(*b1tmp,polCoeff, *datmp);
00402                 datmp++;
00403                 datmp++;
00404                 ++b1tmp;
00405                 --n;
00406         }
00407         while (n>=0);
00408         // that's all, but it is not fast enough... 
00409 }
00410 
00412 template <class TFrommerDefs, int variant>
00413 inline void fast_Frommer< TFrommerDefs, variant>::perform_fast_inner_C_Step(    register typename TMatrixXY_0Type::CoefficientType &  b0,
00414                                                         register const typename TRing::ElementType &  minusp0,
00415                                                         register const typename TRing::ElementType &  q0,
00416                                                         register const typename TRing::ElementType &  da0,
00417                                                          short _n )
00418 {
00419         #ifdef SAFE
00420                 assert(_n<=currentDegree-2);
00421                 assert(_n>0);
00422         #endif 
00423 
00424         register  typename TRing::ElementType const *     datmp          = &da0;
00425         register  typename TMatrixXY_0Type::CoefficientType *     b1tmp  = &b0;
00426 
00427         for (register short i=_n; i>=0; i--)
00428         {
00429                 *b1tmp += q0.getX()*datmp->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00430                 datmp++;
00431                 *b1tmp += minusp0.getX()*datmp->getX() ;        //ACCUMULATE MULTIPLICATION: b = b + p(j,l) * dya(m,n)
00432                 datmp++;
00433                  
00434                 ++b1tmp;
00435         }
00436 
00437  
00438         /*
00439         register short coeff       = q0.getX();
00440 
00441         for (register short i=_n;i>=0;i--)
00442         {
00443                 *b1tmp+= coeff*datmp->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00444                 datmp++;
00445                 datmp++;
00446                  
00447                 ++b1tmp;
00448         }
00449 
00450         coeff       = minusp0.getX();
00451         datmp    = &da0;
00452         b1tmp    = &b0;
00453 
00454         for (register short i=_n;i>=0;i--)
00455         {
00456                 // *b1tmp+= qcoeff*datmp->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00457                 datmp++;
00458                 *b1tmp+= coeff*datmp->getX() ;  //ACCUMULATE MULTIPLICATION: b = b + p(j,l) * dya(m,n)
00459                 datmp++;
00460                  
00461                 ++b1tmp;
00462         }*/
00463 
00464 
00465 }
00466 
00467 
00477 template <class TFrommerDefs, int variant>
00478 inline void fast_Frommer< TFrommerDefs, variant>::perform_C_Step(   )
00479 {
00480   
00481         // Es muss nicht geprueft werden, ob p->getGrad<=currentDegree ist, da im perform_inner_C_Step 
00482         //  sowieso geprft wird, ob n>=0 ist, n=currentDegree-j-l-m;, m \in {0,currentDegree-2}
00483 
00484         // Verbesserung? statt p und q ein pq-Polynom erstellen, weil
00485         // immer auf den gleichen Index zugegriffen wird!!!
00486         
00487 
00488 
00489         #ifdef SAFE
00490                 assert(currentDegree-2>0);
00491         #endif
00492 switch ( q->getDegree() )
00493         {
00494                         
00495 case 3: 
00496 
00497         
00498         //if ( n>0) -  gilt immer , da currentDegree bei 3 startet!
00499         {
00500                 register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.     
00501                         //setzte p0,q0 pointer
00502 
00503                 register  typename TRing::ElementType const*    po = minusp->getCoeffConstAddr( 2,0 );// setzte b1 auf b0
00504                 register  typename TRing::ElementType const*    qo = q->getCoeffConstAddr( 2,0 );
00505                 register  typename TRing::ElementType *         bo = B->getCoeffAddr( currentDegree,0 );
00506                 register  typename TRing::ElementType const *   da0 = dA->getCoeffGroupAddr(n);
00507 
00508                 //(2,0)
00509                         perform_inner_C_Step(bo,po,qo,da0,n);
00510                         
00511                 //(1,1)
00512                         po++;qo++;bo++;
00513                         #ifdef SAFE
00514                                 assert( *po==minusp->getCoeff(1,1) );
00515                                 assert( *qo==q->getCoeff(1,1) );
00516                                 assert( *bo==B->getCoeff(currentDegree-1 , 1) );
00517                         #endif
00518 
00519                         perform_inner_C_Step(bo,po,qo,da0,n);
00520                 
00521         
00522                 //(0,2)
00523                         po++;qo++;bo++;
00524                         #ifdef SAFE
00525                                 assert(*po==minusp->getCoeff(0,2));
00526                                 assert(*qo==q->getCoeff(0,2));
00527                         #endif
00528                         perform_inner_C_Step(bo,po,qo,da0,n);
00529                 
00530                 --n;
00531                 if (n>0)
00532                 {
00533                         // setzte b0 pointer.
00534                         //bo=bo-2; 
00535                         bo = B->getCoeffAddr( currentDegree,0 );
00536                         //setzte p0,q0 pointer
00537                         po=minusp->getCoeffConstAddr(3,0);
00538                         qo=q->getCoeffConstAddr(3,0);
00539                         //po+=2; qo+=2;
00540                 
00541                         da0 = dA->getCoeffGroupAddr(n);
00542                  
00543                         //(3,0)
00544                                 #ifdef SAFE
00545                                         assert(bo == B->getCoeffAddr( 3+n, 0 ));
00546                                         assert(bo == B->getCoeffAddr( currentDegree, 0 ));
00547         
00548                                         assert(*po==minusp->getCoeff(3,0));
00549                                         assert(*qo==q->getCoeff(3,0));
00550                                 #endif
00551 
00552                                 perform_inner_C_Step(bo,po,qo,da0,n);
00553                                 
00554 
00555                         //(2,1)
00556                                 po++;qo++;bo++;
00557                                 perform_inner_C_Step(bo,po,qo,da0,n);
00558 
00559                         //(1,2)
00560                                 po++;qo++;bo++;
00561                                 #ifdef SAFE
00562                                         assert(*po==minusp->getCoeff(1,2));
00563                                         assert(*qo==q->getCoeff(1,2));
00564                                 #endif
00565 
00566                                 perform_inner_C_Step(bo,po,qo,da0,n);
00567                                 
00568                         
00569                         //(0,3)
00570                                 po++;qo++;bo++;
00571                                 #ifdef SAFE
00572                                         assert(*po==minusp->getCoeff(0,3));
00573                                         assert(*qo==q->getCoeff(0,3));
00574                                 #endif
00575 
00576                                 perform_inner_C_Step(bo,po,qo,da0,n);   
00577         
00578                 }
00579         }
00580 
00581 
00582 
00583         break;
00584 /*
00585 case 2:
00586         {
00587 
00588         register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.
00589         if (n>0)
00590         {
00591                 
00592                 register  typename TRing::ElementType const*    po = minusp->getCoeffConstAddr(2,0);// setzte b1 auf b0
00593                 register  typename TRing::ElementType const*    qo = q->getCoeffConstAddr(2,0);
00594                 register  typename TRing::ElementType *         bo = B->getCoeffAddr(currentDegree,0);
00595                 register  typename TRing::ElementType const *   da0 = dA->getCoeffGroupAddr(n);
00596                 //(2,0)
00597                         perform_inner_C_Step(bo,po,qo,da0,n);   
00598                         po++;qo++;bo++;
00599                 //(1,1)
00600                         #ifdef SAFE
00601                                         assert(*po==minusp->getCoeff(1,1));
00602                                         assert(*qo==q->getCoeff(1,1));
00603                         #endif
00604                         perform_inner_C_Step(bo,po,qo,da0,n);   
00605                         po++;qo++;bo++;
00606                 //(0,2)
00607                         #ifdef SAFE
00608                                 assert(*po==minusp->getCoeff(0,2));
00609                                 assert(*qo==q->getCoeff(0,2));
00610                         #endif
00611                         perform_inner_C_Step(bo,po,qo,da0,n);   
00612         }
00613         
00614         break;
00615         }*/
00616 
00617 case 1:
00618 
00619         break;
00620 
00621 
00622 default:
00623         perform_generic_C_Step();
00624         }
00625 }
00626 
00627 /*for monomDegree from 2 to pointDegree do  
00628         (
00629                 if ( (currentDegree - monomDegree)>0  ) then
00630                 (
00631                         for yExp from 0 to monomDegree do
00632                         (
00633                                 n := currentDegree - monomDegree;
00634                                 xExp := monomDegree - yExp;
00635 
00636 
00637 */
00638 
00640 template <class TFrommerDefs, int variant>
00641  void fast_Frommer< TFrommerDefs, variant>::perform_generic_C_Step()
00642 {
00643         short d = q->getDegree();  
00644 
00645         cerr << "q->getDegree()  " << d << endl;
00646         
00647         //if (currentDegree<d)
00648         //      d=currentDegree;
00649 
00650         for (register short  monomDegree=2; monomDegree<=d; monomDegree++)
00651 
00652                 if ( (currentDegree - monomDegree)>0  )
00653                 for (register short yExp=0; yExp<=monomDegree; yExp++)
00654                 {
00655                         short nn = currentDegree-monomDegree ;
00656                         short xExp =  monomDegree - yExp;
00657                         //std::cerr << "xExp: " << xExp << endl;
00658                         //std::cerr << "yExp: " << yExp << endl;
00659                         assert( nn<=currentDegree-2  &&  nn>0 );
00660                         if ( nn<=currentDegree-2  &&  nn>0 )
00661                         {
00662                                 register typename TRing::ElementType const*     po  = minusp->getCoeffConstAddr(xExp,yExp);// setzte b1 auf b0
00663                                 register typename TRing::ElementType const*     qo  = q->getCoeffConstAddr(xExp,yExp);
00664                                 
00665                                 register typename TRing::ElementType *          bo  = B->getCoeffAddr(currentDegree - yExp, yExp);
00666                                 cerr << "currentDegree - yEx " << currentDegree - yExp << endl;
00667                                 cerr << "xExp  " << xExp << endl;
00668                                 register typename TRing::ElementType const *    da0 = dA->getCoeffGroupAddr(nn);
00669                                 #ifdef SAFE
00670                                         
00671                                         assert( *po==minusp->getCoeff(xExp,yExp) );
00672                                         assert( *qo==q->getCoeff(xExp,yExp) );
00673                                         assert( *bo==B->getCoeff(currentDegree - yExp, yExp ) );
00674 
00675                                         cerr << " *po" << *po << endl;
00676                                         cerr << " *qo" << *qo << endl;
00677                                         cerr << " *bo" << *bo << endl;
00678                                 #endif
00679                                 perform_inner_C_Step(bo,po,qo,da0,nn);  
00680                                 #ifdef DEBUG2
00681                                         this->outputMatrix(B);
00682                                         //getchar();
00683                                 #endif
00684                         }
00685                 }
00686 }
00687 
00688 
00689 template <class TFrommerDefs, int variant>
00690 void fast_Frommer< TFrommerDefs, variant>::perform_first_C_Step()
00691 {
00692 
00693         switch ( q->getDegree() )
00694         {
00695                         
00696                 case 3: 
00697                 {
00698 
00699                         //      register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.
00700 
00701                         register  typename TRing::ElementType const*    coeff = q->getCoeffConstAddr( 2,0 );// setzte b1 auf b0
00702                         
00703                         register  typename TMatrixXY_0Type::CoefficientType *   bo = B0->getCoeffAddr( currentDegree,0 );
00704 
00705                         register  typename TRing::ElementType const *   da0 = dA->getCoeffGroupAddr(currentDegree-2);
00706 
00707                         register  typename TMatrixXY_0Type::CoefficientType *     tmpbo = bo;
00708 
00709                         //assert(da0->getX()!=0);
00710                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00711                         tmpbo++;coeff++;
00712                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00713                         tmpbo++;coeff++;
00714                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00715                         
00716                         
00717                         da0++;
00718                         //assert(da0->getX()==0);
00719                         da0++;
00720                         //assert(da0->getX()==0);
00721                         da0++;
00722                         //assert(da0->getX()!=0);
00723                         tmpbo = bo; tmpbo++;
00724                         coeff = minusp->getCoeffConstAddr( 2,0 );
00725                         *tmpbo += coeff->getX()*da0->getX() ; 
00726                         tmpbo++;coeff++;
00727                         *tmpbo += coeff->getX()*da0->getX() ; 
00728                         tmpbo++;coeff++;
00729                         *tmpbo = coeff->getX()*da0->getX() ; 
00730 
00731                         //da0++;
00732                         //assert(da0->getX()==0);
00733 
00734                         tmpbo = bo;
00735                         register  typename TRing::ElementType *         boorig = B->getCoeffAddr(currentDegree,0);
00736                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00737                         tmpbo++;boorig++;
00738                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00739                         tmpbo++;boorig++;
00740                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00741                         tmpbo++;boorig++;
00742                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00743                 }
00744         
00745         case 1:
00746         
00747                 break;
00748         
00749         
00750         default:
00751                 perform_generic_C_Step();
00752         }
00753 }
00754 
00755 
00756 /*
00757 template <class TFrommerDefs, int variant>
00758  void fast_Frommer< TFrommerDefs, variant>::perform_second_C_Step()
00759 {
00760 
00761         #ifdef SAFE
00762                 assert (currentDegree==4);
00763         #endif
00764 
00765         switch ( q->getDegree() )
00766         {
00767                         
00768                 case 3: 
00769                 {
00770 
00771                         //      register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.
00772 
00773                         register  typename TRing::ElementType const*    coeff = q->getCoeffConstAddr( 2,0 );// setzte b1 auf b0
00774                         
00775                         register  typename TMatrixXY_0Type::CoefficientType *   bo = B0->getCoeffAddr( 4,0 );
00776 
00777                         //register  typename TRing::ElementType const *         da0 = dA->getCoeffGroupAddr(currentDegree-2);
00778                         register  typename TRing::ElementType const *   da0 = dA->getCoeffGroupAddr(2);
00779 
00780                         register  typename TMatrixXY_0Type::CoefficientType *   tmpbo = bo;
00781 
00782                         //assert(da0->getX()!=0);
00783                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00784                         tmpbo++;coeff++;
00785                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00786                         tmpbo++;coeff++;
00787                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00788                         
00789                         
00790                         da0++;
00791                         tmpbo = bo;
00792                         coeff = minusp->getCoeffConstAddr( 2,0 );       
00793                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00794                         tmpbo++;coeff++;
00795                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00796                         tmpbo++;coeff++;
00797                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00798 
00799                         //assert(da0->getX()==0);
00800                         da0++;
00801                         bo++;
00802                         tmpbo = bo;
00803                         coeff = q->getCoeffConstAddr( 2,0 );    
00804                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00805                         tmpbo++;coeff++;
00806                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00807                         tmpbo++;coeff++;
00808                         *tmpbo = coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00809 
00810                         //assert(da0->getX()==0);
00811                         da0++;
00812                         tmpbo = bo;
00813                         coeff = minusp->getCoeffConstAddr( 2,0 );       
00814                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00815                         tmpbo++;coeff++;
00816                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00817                         tmpbo++;coeff++;
00818                         *tmpbo += coeff->getX()*da0->getX() ; //ACCUMULATE MULTIPLICATION: b = b + q(j,l) * dxa(m,n)
00819 
00820                         //assert(da0->getX()!=0);
00821                         da0++;
00822                         bo++;
00823                         tmpbo = bo;
00824                         coeff = q->getCoeffConstAddr( 2,0 );
00825                         *tmpbo += coeff->getX()*da0->getX() ; 
00826                         tmpbo++;coeff++;
00827                         *tmpbo += coeff->getX()*da0->getX() ; 
00828                         tmpbo++;coeff++;
00829                         *tmpbo = coeff->getX()*da0->getX() ; 
00830 
00831                         da0++;
00832 
00833                         tmpbo = bo;
00834                         coeff = minusp->getCoeffConstAddr( 2,0 );
00835                         *tmpbo += coeff->getX()*da0->getX() ; 
00836                         tmpbo++;coeff++;
00837                         *tmpbo += coeff->getX()*da0->getX() ; 
00838                         tmpbo++;coeff++;
00839                         *tmpbo += coeff->getX()*da0->getX() ;
00840 
00842                         coeff = q->getCoeffConstAddr( 3,0 ); 
00843                         //da0 = dA->getCoeffGroupAddr(currentDegree-3);
00844                         da0 = dA->getCoeffGroupAddr(1);
00845                         bo = B0->getCoeffAddr( 4,0 );
00846                         tmpbo = bo;
00847 
00848                 
00849                         *tmpbo += coeff->getX()*da0->getX() ; 
00850                         tmpbo++;coeff++;
00851                         *tmpbo += coeff->getX()*da0->getX() ; 
00852                         tmpbo++;coeff++;
00853                         *tmpbo += coeff->getX()*da0->getX() ;
00854                         tmpbo++;coeff++;
00855                         *tmpbo += coeff->getX()*da0->getX() ;
00856 
00857 
00858                         coeff = minusp->getCoeffConstAddr( 3,0 ); 
00859                         //assert(da0->getX()!=0);
00860                         da0++;  
00861 
00862                         //assert(da0->getX()==0);
00863                         da0++;
00864                         //assert(da0->getX()==0);
00865                         da0++;
00866                         //assert(da0->getX()!=0);
00867                         tmpbo = bo; tmpbo ++;
00868                         *tmpbo += coeff->getX()*da0->getX() ; 
00869                         tmpbo++;coeff++;
00870                         *tmpbo += coeff->getX()*da0->getX() ; 
00871                         tmpbo++;coeff++;
00872                         *tmpbo += coeff->getX()*da0->getX() ;
00873                         tmpbo++;coeff++;
00874                         *tmpbo += coeff->getX()*da0->getX() ;
00875 
00876                         //da0++;
00877                         //assert(da0->getX()==0);
00878 
00879                         tmpbo =  bo;
00880                         register  typename TRing::ElementType *         boorig = B->getCoeffAddr(4,0);
00881                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00882                         tmpbo++;boorig++;
00883                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00884                         tmpbo++;boorig++;
00885                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00886                         tmpbo++;boorig++;
00887                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00888                         tmpbo++;boorig++;
00889                         *boorig = field_ref_m.lookupModuloTable( *tmpbo     ) ;
00890                 }
00891         
00892         case 1:
00893         
00894                 break;
00895         
00896         
00897         default:
00898                 perform_generic_C_Step();
00899         }
00900 }*/
00901 
00902 
00904 template <class TFrommerDefs, int variant>
00905  void fast_Frommer< TFrommerDefs, variant>::perform_fast_C_Step()
00906 {
00907   
00908         // Es muss nicht geprueft werden, ob p->getGrad<=currentDegree ist, da im perform_inner_C_Step 
00909         //  sowieso geprft wird, ob n>=0 ist, n=currentDegree-j-l-m;, m \in {0,currentDegree-2}
00910 
00911         // Verbesserung? statt p und q ein pq-Polynom erstellen, weil
00912         // immer auf den gleichen Index zugegriffen wird!!!
00913         
00914         
00915 
00916         #ifdef SAFE
00917                 assert( currentDegree-2>0);
00918         #endif
00919 switch ( q->getDegree() )
00920         {
00921                         
00922 case 3: 
00923 
00924                 
00925          
00926         {
00927                 register short n = currentDegree-2; // n muss n=currentDegree-l-j erfüllen ! una, ach vergesst es.     
00928                 //if ( n>0) -  gilt immer , da currentDegree bei 3 startet!
00929 
00930                         //setzte p0,q0 pointer
00931 
00932                 register  typename TRing::ElementType const*    po = minusp->getCoeffConstAddr( 2,0 );// setzte b1 auf b0
00933                 register  typename TRing::ElementType const*    qo = q->getCoeffConstAddr( 2,0 );
00934                 register  typename TMatrixXY_0Type::CoefficientType *   bo = B0->getCoeffAddr( currentDegree,0 );
00935                 register  typename TRing::ElementType const *   da0 = dA->getCoeffGroupAddr(n);
00936         
00937                 /*      manuell ausgerollte Schleife ist schneller.
00938                         for (int i=0;i<3; i++)
00939                         {
00940                                 perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
00941                                 po++;qo++;bo++;
00942                         }*/
00943 
00944                 //(2,0)
00945                         perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
00946                         
00947                 //(1,1)
00948                         po++;qo++;bo++;
00949                         #ifdef SAFE
00950                                 assert( *po==minusp->getCoeff(1,1) );
00951                                 assert( *qo==q->getCoeff(1,1) );
00952                                 assert( *bo==B0->getCoeff(currentDegree-1 , 1) );
00953                         #endif
00954 
00955                         perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
00956                 
00957         
00958                 //(0,2)
00959                         po++;qo++;bo++;
00960                         #ifdef SAFE
00961                                 assert(*po==minusp->getCoeff(0,2));
00962                                 assert(*qo==q->getCoeff(0,2));
00963                         #endif
00964                         perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
00965                 
00966                 --n;
00967                 if (n>0)
00968                 {
00969                         // setzte b0 pointer.
00970                         //bo=bo-2; 
00971                         bo = B0->getCoeffAddr( currentDegree,0 );
00972                         //setzte p0,q0 pointer
00973                         po=minusp->getCoeffConstAddr(3,0);
00974                         qo=q->getCoeffConstAddr(3,0);
00975                         //po+=2; qo+=2;
00976                 
00977                         da0 = dA->getCoeffGroupAddr(n);
00978 
00979                         #ifdef SAFE
00980                                         assert(bo == B0->getCoeffAddr( 3+n, 0 ));
00981                                         assert(bo == B0->getCoeffAddr( currentDegree, 0 ));
00982         
00983                                         assert(*po==minusp->getCoeff(3,0));
00984                                         assert(*qo==q->getCoeff(3,0));
00985                                 #endif
00986 
00987                         /*      manuell ausgerollte Schleife ist schneller. 
00988                                 for (int i=0;i<4; i++)
00989                         {
00990                                 perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
00991                                 po++;qo++;bo++;
00992                         }*/
00993                 
00994                  
00995                         //(3,0)
00996                                 #ifdef SAFE
00997                                         assert(bo == B0->getCoeffAddr( 3+n, 0 ));
00998                                         assert(bo == B0->getCoeffAddr( currentDegree, 0 ));
00999         
01000                                         assert(*po==minusp->getCoeff(3,0));
01001                                         assert(*qo==q->getCoeff(3,0));
01002                                 #endif
01003 
01004                                 perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
01005                                 
01006 
01007                         //(2,1)
01008                                 po++;qo++;bo++;
01009                                 perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
01010 
01011                         //(1,2)
01012                                 po++;qo++;bo++;
01013                                 #ifdef SAFE
01014                                         assert(*po==minusp->getCoeff(1,2));
01015                                         assert(*qo==q->getCoeff(1,2));
01016                                 #endif
01017 
01018                                 perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
01019                                 
01020                         
01021                         //(0,3)
01022                                 po++;qo++;bo++;
01023                                 #ifdef SAFE
01024                                         assert(*po==minusp->getCoeff(0,3));
01025                                         assert(*qo==q->getCoeff(0,3));
01026                                 #endif
01027 
01028                                 perform_fast_inner_C_Step(*bo,*po,*qo,*da0,n);
01029         
01030                 }
01031                 register  typename TRing::ElementType *         boorig = B->getCoeffAddr(currentDegree,0);
01032                  bo = B0->getCoeffAddr(currentDegree,0);
01033         
01034                 for (register short i=0;i<=currentDegree;i++ )
01035                 {
01036                         *boorig = field_ref_m.lookupModuloTable( *bo     ) ;
01037                         boorig++;
01038                         bo++;
01039                 }
01040         
01041                 
01042 
01043         }
01044 
01045         break;
01046 
01047 
01048 case 1:
01049 
01050         break;
01051 
01052 
01053 default:
01054         perform_generic_C_Step();
01055         }
01056 }
01057 
01058 
01067 template <class TFrommerDefs, int variant>
01068 inline void fast_Frommer< TFrommerDefs, variant>::perform_the_A_Step()
01069 {
01070         if (currentDegree % 2 == 1)  
01071         {
01072                 perform_the_odd_A_Step();
01073         }
01074         else 
01075         {
01076                 perform_the_even_A_Step();
01077         }
01078 }
01079 
01080 
01081 
01082 
01091 
01092 template <class TFrommerDefs, int variant>
01093 inline  void fast_Frommer< TFrommerDefs, variant>::perform_the_odd_A_Step()
01094 {
01095         
01096         //pA = &(A(n-1,1))
01097         typename TRing::ElementType * pA = A->getCoeffAddr(currentDegree-1, 1);
01098         typename TRing::ElementType const * pB = B->getCoeffConstAddr(currentDegree,0);
01099         //A(n-1,1) = B(n,0)     
01100         typename TRing::ElementType prevAVal=* pB;
01101         * pA= prevAVal;
01102         
01103         // die auskommentierte Anweisungen sind richtig, falls im For-Schleifen-Kopf pA+=2 und pB+=2 gerechnet wird, statt im Schlefeninneren.
01104         //pA+=2; pB+=2;
01105         short tmp=currentDegree-1;
01106         //j=3,5,7,..,n:
01107         //A(n-j,j) = ((n-j+2)A(n-j+2,j-2)+B(n-j+1,j-1))/j
01108         for(register short j=3;j<=currentDegree;j+=2,tmp-=2 )
01109         {       
01110                 pA+=2;
01111                 pB+=2;
01112                 #ifdef SAFE     
01113                         assert (pA == A->getCoeffAddr(currentDegree-j, j) );
01114                         assert (pB == B->getCoeffAddr(currentDegree-j+1,j-1) );
01115                         assert (prevAVal==A->getCoeffConstRef(currentDegree-j+2,j-2) );
01116                         assert (tmp ==  currentDegree-j+2 );
01117                 #endif
01118                 prevAVal=       ring1->scalarMultiply( field_ref_m.multInv((typename TRing::ScalarType)(j) ),
01119                                                                 ring1->add( ring1->scalarMultiply(      tmp, 
01120                                                                                                                 prevAVal
01121                                                                                                         ),
01122                                                                                 *pB
01123                                                                         )
01124                                                 );
01125                 *pA=prevAVal;
01126         
01127         }
01128         // die auskommentierte Anweisungen sind richtig, falls im (vorherigen) for-Schleifen-Kopf pA+=2 und pB+=2 gerechnet wird.
01129         //pA-=3;        pB-=1;  
01130         pA-=1;  pB+=1;
01131         #ifdef SAFE     
01132                 assert (pA == A->getCoeffAddr(1,currentDegree-1 ) );
01133                 assert (pB == B->getCoeffAddr(0,currentDegree));
01134         #endif
01135         //A(1,n-1) = -B(0,n)
01136         prevAVal= ring1->addInv(*pB);
01137         *pA=    prevAVal;
01138 
01139         // die auskommentierte Anweisungen sind richtig, falls im (folgenden) For-Schleifen-Kopf pA-=2 und pB-=2 gerechnet wird, statt im Schlefeninneren.
01140         //pA-=2; pB-=2;
01141         tmp=currentDegree-1;// = ((n-j+2)
01142         //j=3,5,7,...n:
01143         for(register short j=3;j<=currentDegree; j+=2, tmp-=2 )
01144         { 
01145                 //A(j,n-j) = ((n-j+2)A(j-2,n-j+2)-B(j-1,n-j+1))/j
01146                 pA-=2;
01147                 pB-=2;
01148                 #ifdef SAFE                     
01149                         assert (pA == A->getCoeffAddr(j, currentDegree-j) );
01150                         assert (pB == B->getCoeffAddr(j-1,currentDegree-j+1) );
01151                         assert(prevAVal==A->getCoeffConstRef(j-2,currentDegree-j+2) );
01152                         assert (tmp ==  currentDegree-j+2 );
01153                 #endif
01154                 prevAVal=               ring1->scalarMultiply( field_ref_m.multInv(     (typename TRing::ScalarType)(j) ),
01155                                                                 ring1->add( ring1->scalarMultiply( tmp,
01156                                                                                                                 prevAVal
01157                                                                                                                 ),
01158                                                                                 ring1->addInv(*pB)
01159                                                                         )
01160                                                         );
01161                 *pA=prevAVal;
01162         }
01163 }
01164 
01165 
01174 
01175 
01176 
01177 template <class TFrommerDefs, int variant>
01178 inline  void fast_Frommer< TFrommerDefs, variant>::perform_the_even_A_Step()
01179 {
01180         // 1. A[n-l,l] berechnen; l=1, siehe Paper
01181         typename  TRing::ScalarType const *  p_a_table_loc = &( a_table[(maxDegree+1)*currentDegree] );
01182 
01183         typename TRing::ElementType Summe = B->getCoeff(currentDegree,0);
01184 
01185         typename TRing::ElementType const * pB = B->getCoeffConstAddr(currentDegree,0);
01186         // jetzt in der Schleife alle uebrigen negativen Werte aufaddieren:
01187         #ifdef SAFE
01188                 assert(currentDegree%2==0);
01189         #endif
01190         for(register short j=1; j<=currentDegree/2; j++)
01191         {       
01192                 pB+=2;
01193 
01194                 #ifdef SAFE     
01195                         assert(pB==B->getCoeffAddr(currentDegree-j*2, j*2) );
01196                 #endif
01197                 ring1->addInPlace( Summe, 
01198                                         ring1->addInv (ring1->scalarMultiply( field_ref_m.multInv( p_a_table_loc[j] ),
01199                                                                                         * pB
01200                                                                                 
01201                                                                                 )
01202                                                 )
01203                                 );
01204         }
01205         ring1->scalarMultiplyInPlace(   field_ref_m.multInv( (typename TRing::ScalarType )(2) ) ,
01206                                                 Summe
01207                                         );
01208 
01209         
01210         typename TRing::ElementType  * pA = A->getCoeffAddr(currentDegree-1,1);
01211 
01212         *pA=Summe;
01213         //A->setCoeff(currentDegree-1, 1, Summe);
01214         // so. jetzt is A[currentDegree-1,1] berechnet.
01215 
01216         pB = B->getCoeffConstAddr(currentDegree,0);
01217         // vielleicht Tabelle anlegen C[x..,..]/a[i]
01218         for(register int j=3; j<=currentDegree-1; j+=2)
01219         {
01220                 short mini = j-1;
01221                 
01222                 pB+=2;
01223                 #ifdef SAFE
01224                         assert(pB==B->getCoeffAddr(currentDegree-mini,mini));
01225                 #endif
01226                 ring1->addInPlace( Summe, 
01227                                         ring1->scalarMultiply( field_ref_m.multInv( p_a_table_loc[mini/2]),
01228                                                                         *pB
01229                                                                 )
01230                                 );
01231 
01232                 mini = min( j, currentDegree-j );
01233                 
01234                 typename TRing::ElementType Summe2 = ring1->scalarMultiply( field_ref_m.multInv( mini), Summe );        
01235                 //  Summe2 kann man sich nicht sparen , da Summe oben  in der Schleife verwendet wird und nicht ueberschrieben werden kann!
01236 
01237                 mini--;  mini /= 2;
01238 
01239                 pA+=2;
01240                 #ifdef SAFE
01241                         assert(pA==A->getCoeffAddr(currentDegree-j,j));
01242                 #endif
01243 
01244                 *pA= ring1->scalarMultiply(  p_a_table_loc[mini], Summe2 ) ;
01245         }
01246         //puhhh, geschafft.
01247 
01248         //------------
01249         //A(0,n) = 0;
01250         pA++;
01251         #ifdef SAFE
01252                 assert(pA==A->getCoeffAddr(0,currentDegree));
01253         #endif
01254         *pA=TRing::ElementType::Zero;
01255         //A->setCoeff(0,currentDegree, TRing::ElementType::Zero);
01256         typename TRing::ElementType prevAVal=* pA;
01257 
01258         pB = B->getCoeffConstAddr(1,currentDegree-1);
01259         
01260         //j=2,4,6,...n:
01261         //A(j,n-j) = ((n-j+2)A(j-2,n-j+2)-B(j-1,n-j+1))/j
01262         for(register short j=2;j<=currentDegree;j+=2)
01263         { 
01264                 pA-=2;
01265                 #ifdef SAFE
01266                         short j_minus_zwei=j-2;         
01267                         assert(pA==A->getCoeffAddr(j, currentDegree-j));
01268                         assert(pB==B->getCoeffAddr(j-1, currentDegree-j+1));
01269                         assert(prevAVal==A->getCoeffConstRef(j_minus_zwei,currentDegree-j_minus_zwei));
01270                 #endif
01271                 prevAVal=       ring1->scalarMultiply ( field_ref_m.multInv( (typename TRing::ScalarType)(j) ),
01272                                                         ring1->add( ring1->scalarMultiply(      currentDegree-j+2,
01273                                                                                                         prevAVal
01274                                                                                                         ),
01275                                                                                 ring1->addInv(*pB )
01276                                                                         )
01277                                                         
01278                         );
01279                 *pA=prevAVal;
01280                 pB-=2;
01281         }
01282 }
01283 
01284 
01286 template <class TFrommerDefs, int variant>
01287 template <class TemplateMatrix>
01288 inline void fast_Frommer< TFrommerDefs, variant>::outputMatrix(TemplateMatrix *mat) const
01289 {
01290         std::cerr << mat->getName() << std::endl;
01291         mat->outputMatrix();
01292         return;
01293 }
01294 
01296 template <class TFrommerDefs, int variant>
01297 inline void fast_Frommer< TFrommerDefs, variant>::outputMatrices() const
01298 {
01299         
01300         std::cout << "A: " << std::endl;
01301                 A->outputMatrix();
01302                 //getchar();
01303         
01304         std::cout << "dA: " << std::endl;
01305                 dA->outputMatrix();
01306         
01307         
01308         std::cout << "B: " << std::endl;
01309                 B->outputMatrix();
01310                 //getchar();
01311         return;
01312 }
01313 
01314 
01315 template <class TFrommerDefs, int variant>
01316 template<int tHowManyTimes>
01317   void fast_Frommer< TFrommerDefs, variant>::doitx(short howManyFocalValues)
01318 {
01319         #ifdef CF_TEST
01320                 //std::cerr << "testDoitX" << std::endl;
01321                 doitxOld<tHowManyTimes>(howManyFocalValues);
01322                 TMatrixXYType B=getB();
01323                 //TMatrixXYType A=getA();
01324 
01325                 doitxNew<tHowManyTimes>(howManyFocalValues);
01326                 assert(B==getB());
01327                 //assert(A==getA());
01328 
01329         #else
01330                 doitxNew<tHowManyTimes>(howManyFocalValues);
01331         #endif
01332 }
01333 
01338 
01339 template <class TFrommerDefs, int variant>
01340 template<int tHowManyTimes>
01341   void fast_Frommer< TFrommerDefs, variant>::doitxNew(short howManyFocalValues)
01342 {
01343         init();
01344         #ifdef SAFE
01345                 assert(initialized);
01346                 initialized=false;
01347         #endif
01348 
01349         currFocalValuePos = 1;
01350         
01351         currentDegree = 2; 
01352 
01353 
01354         if  (ring1->getEpsPrecision()==0)
01355         {
01356                 #ifdef FROMMERSTAGE_TIMER
01357                         tim_m.start();
01358                         #endif
01359 
01360                 currentDegree++; 
01361                                 
01362                         
01363                 compute_dyA_and_dxA();
01364 
01365                 #ifdef FROMMERSTAGE_TIMER
01366                                         tim2_m.start();
01367                 #endif
01368 
01369                 perform_first_C_Step();
01370                 //perform_C_Step();
01371 
01372                 #ifdef DEBUG2
01373                         this->outputMatrix(B);
01374                 #endif
01375                         #ifdef FROMMERSTAGE_TIMER
01376                                         tim2_m.stop();
01377                                         frommerTimer_m.logCStepTime(tHowManyTimes,currentDegree-2 , tim2_m );
01378                                         tim2_m.clear();
01379                                         tim2_m.start();
01380                                 #endif
01381 
01382 
01383                 perform_the_A_Step();
01384 
01385                         #ifdef FROMMERSTAGE_TIMER
01386                                 tim2_m.stop();
01387                                 tim_m.stop();   
01388                                 frommerTimer_m.logAStepTime(tHowManyTimes,currentDegree-2, tim2_m );
01389                                 frommerTimer_m.logStepTime(tHowManyTimes, currentDegree-2, tim_m );
01390                                 tim2_m.clear();
01391                                 tim_m.clear();
01392                         #endif
01393 
01394                 while(true)
01395                 {
01396                                 #ifdef FROMMERSTAGE_TIMER
01397                                         tim_m.start();
01398                                 #endif
01399         
01400                         currentDegree++; 
01401                                 
01402                         compute_dyA_and_dxA();
01403                                 #ifdef DEBUG2
01404                                         this->outputMatrices();
01405                                 #endif
01406                                 #ifdef FROMMERSTAGE_TIMER
01407                                         tim2_m.start();
01408                                 #endif
01409                          
01410                                 perform_fast_C_Step();
01411                                 //perform_C_Step();
01412                                 
01413                         
01414                                 #ifdef DEBUG2
01415                                         this->outputMatrix(B);
01416                                 #endif
01417                                 #ifdef FROMMERSTAGE_TIMER
01418                                         tim2_m.stop();
01419                                         frommerTimer_m.logCStepTime(tHowManyTimes,currentDegree-2 , tim2_m );
01420                                         tim2_m.clear();
01421                                         tim2_m.start();
01422                                 #endif
01423                         perform_the_A_Step();
01424                                 #ifdef DEBUG2
01425                                         this->outputMatrix(A);
01426                                 #endif
01427                         #ifdef FROMMERSTAGE_TIMER
01428                                 tim2_m.stop();
01429                                 tim_m.stop();   
01430                                 frommerTimer_m.logAStepTime(tHowManyTimes,currentDegree-2, tim2_m );
01431                                 frommerTimer_m.logStepTime(tHowManyTimes, currentDegree-2, tim_m );
01432                                 tim2_m.clear();
01433                                 tim_m.clear();
01434                         #endif
01435                         if ( (currentDegree %2 )==0)//gerade --> Strudelgroesse berechnen
01436                         {
01437                                 computedFocalValues[currFocalValuePos] = ring1->add(    A->getCoeffConstRef(1, currentDegree-1),
01438                                                                                                         B->getCoeffConstRef(0, currentDegree));
01439         
01440                                 ++currFocalValuePos;
01441         
01442                                 if (  currFocalValuePos > howManyFocalValues )  break;
01443                         }
01444                 }
01445         }
01446         else
01447                 while(true)
01448                 {
01449                                 #ifdef FROMMERSTAGE_TIMER
01450                                         tim_m.start();
01451                                 #endif
01452         
01453                         currentDegree++; 
01454                                 
01455                         compute_dyA_and_dxA();
01456                                 #ifdef DEBUG2
01457                                         this->outputMatrices();
01458                                 #endif
01459                                 #ifdef FROMMERSTAGE_TIMER
01460                                         tim2_m.start();
01461                                 #endif
01462                                 perform_C_Step();
01463                                 
01464 
01465                                 #ifdef DEBUG2
01466                                         this->outputMatrix(B);
01467                                 #endif
01468                                 #ifdef FROMMERSTAGE_TIMER
01469                                         tim2_m.stop();
01470                                         frommerTimer_m.logCStepTime(tHowManyTimes,currentDegree-2 , tim2_m );
01471                                         tim2_m.clear();
01472                                         tim2_m.start();
01473                                 #endif
01474                         perform_the_A_Step();
01475                                 #ifdef DEBUG2
01476                                         this->outputMatrix(A);
01477                                 #endif
01478                         #ifdef FROMMERSTAGE_TIMER
01479                                 tim2_m.stop();
01480                                 tim_m.stop();   
01481                                 frommerTimer_m.logAStepTime(tHowManyTimes,currentDegree-2, tim2_m );
01482                                 frommerTimer_m.logStepTime(tHowManyTimes, currentDegree-2, tim_m );
01483                                 tim2_m.clear();
01484                                 tim_m.clear();
01485                         #endif
01486                         if ( (currentDegree %2 )==0)//gerade --> Strudelgroesse berechnen
01487                         {
01488                                 computedFocalValues[currFocalValuePos] = ring1->add(    A->getCoeffConstRef(1, currentDegree-1),
01489                                                                                                         B->getCoeffConstRef(0, currentDegree));
01490         
01491                                 ++currFocalValuePos;
01492         
01493                                 if (  currFocalValuePos > howManyFocalValues )  break;
01494                         }
01495                 }
01496 
01497         // Merke mir aktuellen Grad, um bei einem Neudurchlauf nur verwendete Datenbereiche 
01498         // zu löschen, da dies eine extrem teuere Operation ist!
01499         lastDegree = currentDegree;
01500         #ifdef SAFE
01501                 assert( (currFocalValuePos-1) <= maxFocalValuesToCompute );
01502         #endif 
01503 }
01504 
01505 
01506 
01507 template <class TFrommerDefs, int variant>
01508 template<int tHowManyTimes>
01509   void fast_Frommer< TFrommerDefs, variant>::doitxOld(short howManyFocalValues)
01510 {
01511         init();
01512         #ifdef SAFE
01513                 assert(initialized);
01514                 initialized=false;
01515         #endif
01516 
01517         currFocalValuePos = 1;
01518         
01519         currentDegree = 2; 
01520 
01521 
01522         while(true)
01523         {
01524                         #ifdef FROMMERSTAGE_TIMER
01525                                 tim_m.start();
01526                         #endif
01527 
01528                 currentDegree++; 
01529                         
01530                 compute_dyA_and_dxA();
01531                         #ifdef DEBUG2
01532                                 this->outputMatrices();
01533                         #endif
01534 
01535                         #ifdef FROMMERSTAGE_TIMER
01536                                 tim2_m.start();
01537                         #endif
01538                         perform_C_Step();
01539                         
01540 
01541                         #ifdef DEBUG2
01542                                 this->outputMatrix(B);
01543                         #endif
01544 
01545                         #ifdef FROMMERSTAGE_TIMER
01546                                 tim2_m.stop();
01547                                 frommerTimer_m.logCStepTime(tHowManyTimes,currentDegree-2 , tim2_m );
01548                                 tim2_m.clear();
01549                                 tim2_m.start();
01550                         #endif
01551 
01552                 perform_the_A_Step();
01553 
01554                         #ifdef DEBUG2
01555                                 this->outputMatrix(A);
01556                         #endif
01557 
01558                 #ifdef FROMMERSTAGE_TIMER
01559                         tim2_m.stop();
01560                         tim_m.stop();   
01561                         frommerTimer_m.logAStepTime(tHowManyTimes,currentDegree-2, tim2_m );
01562                         frommerTimer_m.logStepTime(tHowManyTimes, currentDegree-2, tim_m );
01563                         tim2_m.clear();
01564                         tim_m.clear();
01565                 #endif
01566                 if ( (currentDegree %2 )==0)//gerade --> Strudelgroesse berechnen
01567                 {
01568                         computedFocalValues[currFocalValuePos] = ring1->add(    A->getCoeffConstRef(1, currentDegree-1),
01569                                                                                                 B->getCoeffConstRef(0, currentDegree));
01570 
01571                         ++currFocalValuePos;
01572 
01573                         if (  currFocalValuePos > howManyFocalValues )  break;
01574                 }
01575 
01576         }
01577         // Merke mir aktuellen Grad, um bei einem Neudurchlauf nur verwendete Datenbereiche 
01578         // zu löschen, da dies eine extrem teuere Operation ist!
01579         lastDegree = currentDegree;
01580         #ifdef SAFE
01581                 assert( (currFocalValuePos-1) <= maxFocalValuesToCompute );
01582         #endif 
01583 }
01584 
01585 
01586 
01587 
01588 template <class TFrommerDefs, int variant>
01589   void fast_Frommer< TFrommerDefs, variant>::doit(short bComputeAllFocalValues, int jacobian )
01590 {
01591         #ifdef CF_TEST
01592 
01593                 doitOld(bComputeAllFocalValues,jacobian);
01594                 TMatrixXYType B=getB();
01595                 //TMatrixXYType A=getA();
01596 
01597                 #ifdef SAFE
01598                         initialized=true;
01599                 #endif
01600                 doitNew(bComputeAllFocalValues,jacobian);
01601                 assert(B==getB());
01602                 //assert(A==getA());
01603 
01604         #else
01605                 doitNew(bComputeAllFocalValues, jacobian);
01606         #endif
01607 }
01608 
01619 template <class TFrommerDefs, int variant>
01620    void fast_Frommer< TFrommerDefs, variant>::doitNew(short bComputeAllFocalValues, int jacobian )
01621 {
01622         // init gehoert hier rein, da SetPolynoms nicht immer aufgerufen werden muss!
01623         
01624         #ifdef FROMMERSTAGE_TIMER
01625                  
01626                 int hm=0;
01627                 if (jacobian)
01628                         hm=-2 ;
01629                 else if (bComputeAllFocalValues)
01630                         hm=maxFocalValuesToCompute ;
01631                 else
01632                         hm = 0 ;
01633         #endif
01634 
01635 
01636         #ifdef SAFE
01637                 assert(initialized); 
01638                 initialized = false;
01639         #endif
01640         init();
01641          currFocalValuePos = 1;
01642         
01643         currentDegree = 2;
01644 
01645         typename TRing::ElementType res;
01646 
01647         try
01648         {
01649 
01650                 #ifdef FROMMERSTAGE_TIMER
01651                         tim_m.start();
01652                 #endif
01653 
01654                 if  (ring1->getEpsPrecision()==0)
01655                 {
01656                         currentDegree++;        
01657                         compute_dyA_and_dxA();
01658 
01659                         perform_first_C_Step();
01660                         //perform_C_Step();
01661 
01662                         #ifdef FROMMERSTAGE_TIMER
01663                                 tim2_m.stop();
01664                                 frommerTimer_m.logCStepTime(hm, currentDegree-2, tim2_m );
01665                                 tim2_m.clear();
01666                                 tim2_m.start();
01667                         #endif
01668 
01669                         #ifdef DEBUG2
01670                                         this->outputMatrix(B);
01671                         #endif
01672 
01673                         perform_the_A_Step();
01674                         #ifdef DEBUG2
01675                                         this->outputMatrix(A);
01676                         #endif
01677 
01678                         #ifdef FROMMERSTAGE_TIMER
01679                                 tim2_m.stop();
01680                                 tim_m.stop();   
01681                                 frommerTimer_m.logAStepTime(hm, currentDegree-2, tim2_m );
01682                                 frommerTimer_m.logStepTime( hm, currentDegree-2, tim_m );
01683                                 tim2_m.clear();
01684                                 tim_m.clear();
01685                         #endif
01686 
01687                 }
01688                 else
01689                 {
01690 
01691                 
01692                         currentDegree++;        
01693                         compute_dyA_and_dxA();
01694 
01695                 
01696 
01697                         perform_C_Step();
01698 
01699 
01700                         #ifdef DEBUG2
01701                                 this->outputMatrix(B);
01702                         #endif
01703 
01704                         #ifdef FROMMERSTAGE_TIMER
01705                                 tim2_m.stop();
01706                                 frommerTimer_m.logCStepTime(hm, currentDegree-2, tim2_m );
01707                                 tim2_m.clear();
01708                                 tim2_m.start();
01709                         #endif
01710 
01711 
01712                         perform_the_A_Step();
01713                         #ifdef DEBUG2
01714                                         this->outputMatrix(A);
01715                         #endif
01716 
01717                         #ifdef FROMMERSTAGE_TIMER
01718                                 tim2_m.stop();
01719                                 tim_m.stop();   
01720                                 frommerTimer_m.logAStepTime(hm, currentDegree-2, tim2_m );
01721                                 frommerTimer_m.logStepTime( hm, currentDegree-2, tim_m );
01722                                 tim2_m.clear();
01723                                 tim_m.clear();
01724                         #endif
01725 
01726                 }
01727 
01728                 while(true)
01729                 {
01730                         #ifdef FROMMERSTAGE_TIMER
01731                                 tim_m.start();
01732                         #endif
01733 
01734                         currentDegree++; 
01735                         
01736                         compute_dyA_and_dxA();
01737 
01738                                 #ifdef DEBUG2
01739                                         this->outputMatrices();
01740                                 #endif
01741                         #ifdef FROMMERSTAGE_TIMER
01742                                 tim2_m.start();
01743                         #endif
01744 
01745                         if  (ring1->getEpsPrecision()!=0)
01746                         {
01747 
01748                                         perform_C_Step();
01749                         }
01750                                 else
01751                         {
01752                          
01753                                         perform_fast_C_Step();
01754                                         //perform_C_Step();
01755                         }
01756 
01757                                 #ifdef DEBUG2
01758                                         this->outputMatrix(B);
01759                                 #endif
01760 
01761                         #ifdef FROMMERSTAGE_TIMER
01762                                 tim2_m.stop();
01763                                 frommerTimer_m.logCStepTime(hm, currentDegree-2, tim2_m );
01764                                 tim2_m.clear();
01765                                 tim2_m.start();
01766                         #endif
01767                         perform_the_A_Step();
01768 
01769                                 #ifdef DEBUG2
01770                                         this->outputMatrix(A);
01771                                 #endif
01772                         #ifdef FROMMERSTAGE_TIMER
01773                                 tim2_m.stop();
01774                                 tim_m.stop();   
01775                                 frommerTimer_m.logAStepTime(hm, currentDegree-2, tim2_m );
01776                                 frommerTimer_m.logStepTime( hm, currentDegree-2, tim_m );
01777                                 tim2_m.clear();
01778                                 tim_m.clear();
01779                         #endif
01780 
01781                         if ((currentDegree %2) ==0)
01782                         {
01783 
01784                                 //currentDegree gerade --> Strudelgroesse berechnen
01785 
01786                                 res = ring1->add( A->getCoeffConstRef(1, currentDegree-1), 
01787                                                                                    B->getCoeffConstRef(0, currentDegree));
01788 
01789                                 computedFocalValues[currFocalValuePos] =res;
01790                                 #ifdef DEBUG
01791                                         std::cerr << currFocalValuePos << ". Strudelgroesse: ";
01792                                         std::cerr << (computedFocalValues[currFocalValuePos].getX()) << std::endl;
01793                                 #endif
01794 
01795 
01796                                 if (++currFocalValuePos > maxFocalValuesToCompute  )
01797                                 {                                                               
01798                                         break;
01799                                 }
01800 
01801                                 //++currFocalValuePos gehoert ganz ans ende der if-Abfrage !, 
01802                                 if ( (bComputeAllFocalValues==0) && (res.getX() != 0)    )
01803                                 {                                                               
01804                                         break;
01805                                 }
01806                                 // else: bComputeAllFocalValues !=0 or computedFocalValues[currFocalValuePos].getX()==0 : compute next focal Value!
01807                         }
01808                 }
01809                 // Merke mir aktuellen Grad, um bei einem Neudurchlauf nur verwendete Datenbereiche 
01810                 // zu löschen, da dies eine extrem teuere Operation ist!
01811                 lastDegree = currentDegree;
01812 
01813         }
01814         catch (const char* str)
01815         {       
01816                 #ifdef DEBUG
01817                         std::cerr <<  typeid(this).name() << std::endl;
01818                 #endif 
01819                 std::cerr <<"fast_Frommer< TFrommerDefs, variant> :: doit error during computation of " << currFocalValuePos << "-th focal value: " << std::endl << str << std::endl;
01820                 exit(0);
01821         }
01822 }
01823 
01824 
01825 
01826 template <class TFrommerDefs, int variant>
01827    void fast_Frommer< TFrommerDefs, variant>::doitOld(short bComputeAllFocalValues, int jacobian )
01828 {
01829         // init gehoert hier rein, da SetPolynoms nicht immer aufgerufen werden muss!
01830         
01831         #ifdef FROMMERSTAGE_TIMER
01832                  
01833                 int hm=0;
01834                 if (jacobian)
01835                         hm=-2 ;
01836                 else if (bComputeAllFocalValues)
01837                         hm=maxFocalValuesToCompute ;
01838                 else
01839                         hm = 0 ;
01840         #endif
01841 
01842 
01843         #ifdef SAFE
01844                 assert(initialized);
01845                 initialized = false;
01846         #endif
01847         init();
01848          currFocalValuePos = 1;
01849         
01850         currentDegree = 2;
01851 
01852         typename TRing::ElementType res;
01853 
01854         try
01855         {
01856 
01857                 while(true)
01858                 {
01859                         #ifdef FROMMERSTAGE_TIMER
01860                                 tim_m.start();
01861                         #endif
01862 
01863                         currentDegree++; 
01864                         
01865                         compute_dyA_and_dxA();
01866 
01867                                 #ifdef DEBUG2
01868                                         this->outputMatrices();
01869                                 #endif
01870                         #ifdef FROMMERSTAGE_TIMER
01871                                 tim2_m.start();
01872                         #endif
01873 
01874                                 perform_C_Step();
01875 
01876                                 #ifdef DEBUG2
01877                                         this->outputMatrix(B);
01878                                 #endif
01879 
01880                         #ifdef FROMMERSTAGE_TIMER
01881                                 tim2_m.stop();
01882                                 frommerTimer_m.logCStepTime(hm, currentDegree-2, tim2_m );
01883                                 tim2_m.clear();
01884                                 tim2_m.start();
01885                         #endif
01886 
01887                         perform_the_A_Step();
01888 
01889                                 #ifdef DEBUG2
01890                                         this->outputMatrix(A);
01891                                 #endif
01892 
01893                         #ifdef FROMMERSTAGE_TIMER
01894                                 tim2_m.stop();
01895                                 tim_m.stop();   
01896                                 frommerTimer_m.logAStepTime(hm, currentDegree-2, tim2_m );
01897                                 frommerTimer_m.logStepTime( hm, currentDegree-2, tim_m );
01898                                 tim2_m.clear();
01899                                 tim_m.clear();
01900                         #endif
01901 
01902                         if ((currentDegree %2) ==0)
01903                         {
01904 
01905                                 //currentDegree gerade --> Strudelgroesse berechnen
01906 
01907                                 res = ring1->add( A->getCoeffConstRef(1, currentDegree-1), 
01908                                                                                   B->getCoeffConstRef(0, currentDegree));
01909 
01910                                 computedFocalValues[currFocalValuePos] =res;
01911                                 #ifdef DEBUG
01912                                         std::cerr << currFocalValuePos << ". Strudelgroesse: ";
01913                                         std::cerr << (computedFocalValues[currFocalValuePos].getX()) << std::endl;
01914                                 #endif
01915 
01916 
01917                                 if (++currFocalValuePos > maxFocalValuesToCompute  )
01918                                 {                                                               
01919                                         break;
01920                                 }
01921 
01922                                 //++currFocalValuePos gehoert ganz ans ende der if-Abfrage !, 
01923                                 if ( (bComputeAllFocalValues==0) && (res.getX() != 0)    )
01924                                 {                                                               
01925                                         break;
01926                                 }
01927                                 // else: bComputeAllFocalValues !=0 or computedFocalValues[currFocalValuePos].getX()==0 : compute next focal Value!
01928                         }
01929                 }
01930                 // Merke mir aktuellen Grad, um bei einem Neudurchlauf nur verwendete Datenbereiche 
01931                 // zu löschen, da dies eine extrem teuere Operation ist!
01932                 lastDegree = currentDegree;
01933 
01934         }
01935         catch (const char* str)
01936         {       
01937                 #ifdef DEBUG
01938                         std::cerr <<  typeid(this).name() << std::endl;
01939                 #endif 
01940                 std::cerr <<"fast_Frommer< TFrommerDefs, variant> :: doit error during computation of " << currFocalValuePos << "-th focal value: " << std::endl << str << std::endl;
01941                 exit(0);
01942         }
01943 }
01944 
01945 
01955 template <class TFrommerDefs, int variant>
01956 inline typename TFrommerDefs::RingType::ElementType & fast_Frommer< TFrommerDefs, variant>::getFocalValue( int pos) const
01957 {
01958 
01959         if (currFocalValuePos>pos && pos>0)
01960                 return computedFocalValues[pos];
01961         else 
01962                 throw " requested focal value was not computed!" ;
01963 
01964 }
01965 
01966 template < class TFrommerDefs, int variant >
01967 void            fast_Frommer< TFrommerDefs, variant>::getComputedFocalValues( vector<typename TFrommerDefs::RingType::ElementType> &computedFV ) const
01968 {
01969         computedFV.clear();
01970         computedFV.resize( currFocalValuePos-1 );
01971         for (int pos=1; pos<currFocalValuePos; pos++)
01972         {
01973                 computedFV[ pos - 1 ] = computedFocalValues[ pos ] ;
01974         }
01975 }
01976 
01977 
01982 template <class TFrommerDefs, int variant>
01983 inline short fast_Frommer< TFrommerDefs, variant>::getSuccessiveVanishedFocalValuesCount() const
01984 {
01985         short i=0;
01986         while ( (i<maxFocalValuesToCompute) && (i+1<currFocalValuePos) && (computedFocalValues[i+1].getX()==0 )/*&&(computedFocalValues[i+1].eps==0)*/)
01987         {
01988                 
01989                 i++;
01990         }
01991         assert(currFocalValuePos>=i);
01992         return i;
01993 }
01994 
01995 
01996 template < class TFrommerDefs, int variant >
01997 void fast_Frommer< TFrommerDefs, variant>::printStageTimings(std::ostream & os) const
01998 {
01999         #ifdef FROMMERSTAGE_TIMER
02000                 os << "-- " << name_m << ": " << endl;
02001                 frommerTimer_m.print(os);
02002         #endif
02003 }