00001
00002
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00062
00063 template <int variant>
00064 CenterFocusExperiment<variant>::CenterFocusExperiment( D_CenterfocusParams const * params,
00065 PrintCenterFocusResults * prn
00066 ) :
00067
00068 cfparams_m(params),
00069 resultPrinter_m(prn),
00070 allCoeffVariablesOrder_m(params->getCoeffVariablesOrder() ),
00071 epsRing_m(new typename D_CenterfocusParams::RingType(params->getFieldCharacteristic(),
00072 max( params->getEpsPrecision(),2 )
00073 )
00074 ),
00075 epsFieldRef_m(epsRing_m->getFieldRef() ),
00076 epsFrommer_m( new D_Eps_Frommer( params->getMaxFocalValuesToCompute(),
00077 epsRing_m
00078 )
00079 ),
00080 epsRegularFrommer_m( new D_Eps_Frommer( params->getMaxFocalValuesToCompute(), params->getRing() )),
00081 ring2_m ( new D_Frommer2::RingType(cfparams_m->getFieldCharacteristic(), 1 )),
00082 ring0_m ( new D_Frommer0::RingType(cfparams_m->getFieldCharacteristic(), 0 )),
00083 polynomialRing_m (epsRegularFrommer_m->getRingRef() ),
00084 polynomialRing_2_m(*ring2_m),
00085 polynomialRing_0_m(*ring0_m)
00086
00087
00088
00089 {
00090 frommer2_m= new D_Frommer( cfparams_m->getMaxFocalValuesToCompute(), ring2_m);
00091 frommer0_m= new D_Frommer0( cfparams_m->getMaxFocalValuesToCompute(), ring0_m);
00092
00093
00094 assert(epsRegularFrommer_m->getRing()->getEpsPrecision()==cfparams_m->getEpsPrecision() );
00095
00096
00097 minusP_polynom_m = new D_Eps_Frommer::TPolynomXY (cfparams_m->getPolynomialDegree() );
00098 Q_polynom_m = new D_Eps_Frommer::TPolynomXY (cfparams_m->getPolynomialDegree() );
00099
00100
00101
00102 st_m=new DStatistic( cfparams_m->getMaxFocalValuesToCompute() , 0);
00103
00104
00105 #ifdef SAFE
00106 assert( typeid( D_Eps_Frommer::RingType ) == typeid( D_CenterfocusParams::RingType) );
00107 assert( typeid( D_Eps_Frommer::TPolynomXY ) == typeid( D_CenterfocusParams::PolynomXY_Type) );
00108 #endif
00109
00110 #ifdef ALGORITHM_TIMER
00111 hamiltonianCheckTime_m=0;
00112 computeVanishedVocalValuesTime_m=0;
00113 computeJacobianMatrixTime_m=0;
00114 computeJacobianRankTime_m=0;
00115 checkSmoothnessTime_m=0;
00116 computeQuadricsTime_m=0;
00117 computeAllFocalValuesTime_m=0;
00118 addQuadricstatistic_m=0;
00119 smoothnessTest_computeJacobiMatrixTime_m=0;
00120 smoothnessTest_solveLGSTime_m=0;
00121 smoothnessTest_getJacobianKernelTime_m=0;
00122 smoothnessTest_addEpsVector_m=0;
00123 smoothnessTest_copyPolynomTime_m=0;
00124 smoothnessTest_computeFocalValuesTime_m=0;
00125 #endif
00126 }
00127
00128
00129
00131 template <int variant>
00132 void CenterFocusExperiment<variant>::internPrintVariableOrder( ostream & os,
00133 const list<CoeffListEntry> & coeffVariablesOrder,
00134 string comment)
00135 {
00136 os << "-- " << comment << " : " << std::endl;
00137 os << "-- " ;
00138 list<CoeffListEntry>::const_iterator coeffIt;
00139
00140 coeffIt= coeffVariablesOrder.begin();
00141
00142 for (unsigned int row = 0; row< coeffVariablesOrder.size(); row++, coeffIt++)
00143 {
00144 if ((*coeffIt).isPMonom() )
00145 os << "pp_";
00146 else os << "qq_";
00147 os << (*coeffIt).x_exp;
00148 os << (*coeffIt).y_exp;
00149 os << ", ";
00150 }
00151 os << std::endl;
00152 return;
00153 }
00154
00155
00156
00158 template <int variant>
00159 void CenterFocusExperiment<variant>::printVariableOrder(ostream & os )
00160 {
00161 internPrintVariableOrder(os, allCoeffVariablesOrder_m, "fullVariableOrder");
00162 internPrintVariableOrder(os, subCoeffVariablesOrder_m, "subVariableOrder");
00163 return;
00164 }
00165
00166 template <int variant>
00167 void purge()
00168 {
00169
00170
00171 }
00172
00175 template <int variant>
00176 list<CoeffListEntry> CenterFocusExperiment<variant>::initSubCoeffVariablesOrder()
00177 {
00178 list<CoeffListEntry> coeffVariablesOrder;
00179 list<CoeffListEntry>::const_iterator coeffIt;
00180
00181 coeffIt =allCoeffVariablesOrder_m.begin();
00182 while (coeffIt!=allCoeffVariablesOrder_m.end() )
00183 {
00184 if ( cfparams_m->isCoefficientVariable( *coeffIt ) )
00185 coeffVariablesOrder.push_back( *coeffIt );
00186 coeffIt++;
00187 }
00188 return coeffVariablesOrder;
00189 }
00190
00191
00194 template <int variant>
00195 void CenterFocusExperiment<variant>::initExperiment( const D_CenterfocusParams * params)
00196 {
00197 cfparams_m = params;
00198
00199 inputPointCounter_m=0;
00200
00201 assert( cfparams_m->getPolynomialDegree()>=2);
00202
00203
00204 subCoeffVariablesOrder_m.clear();
00205
00206 subCoeffVariablesOrder_m = initSubCoeffVariablesOrder( );
00207
00208 if (minusP_polynom_m)
00209 delete minusP_polynom_m;
00210
00211 if (Q_polynom_m)
00212 delete Q_polynom_m;
00213
00214
00215 minusP_polynom_m = new D_Eps_Frommer::TPolynomXY( cfparams_m->getMinusPPolynomRef() ) ;
00216 Q_polynom_m = new D_Eps_Frommer::TPolynomXY( cfparams_m->getQPolynomRef() );
00217 }
00218
00219
00228
00234 template <int variant>
00235 void CenterFocusExperiment<variant>::performExperiment( const D_CenterfocusParams * params,
00236 RankStatistic &fullRankStatistic,
00237 RankStatistic &subRankStatistic,
00238 LiftAndQuadricsStatistic &fullQuadricsStatistic,
00239 LiftAndQuadricsStatistic &subQuadricsStatistic,
00240 FailedLiftStatistic &liftStatistic)
00241 {
00242
00243 initExperiment(params);
00244
00245
00246 if ( !cfparams_m->polynomsHaveCheckAllCoeff() &&
00247 !cfparams_m->polynomsHaveRandomCoeff()
00248 )
00249
00250 {
00251
00252 if (cfparams_m->getRequestedTrialsNum()>0)
00253 {
00254 std::cerr << "Error: it is not possible to perform random trials without random coefficients" <<std::endl;
00255 assert(false);
00256 }
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 if (cfparams_m->getEpsPrecision()==0)
00270 {
00271 D_Frommer0::TPolynomXY * MinusP_polynom = new D_Frommer0::TPolynomXY ( cfparams_m->getPolynomialDegree() );
00272 D_Frommer0::TPolynomXY * Q_polynom = new D_Frommer0::TPolynomXY ( cfparams_m->getPolynomialDegree() );
00273
00274
00275
00276
00277
00278
00279 copyPolynomWithGivenEpsPrecision( *minusP_polynom_m,*MinusP_polynom, cfparams_m->getEpsPrecision() );
00280 copyPolynomWithGivenEpsPrecision( *Q_polynom_m,*Q_polynom , cfparams_m->getEpsPrecision() );
00281
00282 performRegularExperiment(
00283 *frommer0_m,
00284
00285 *frommer2_m,
00286 polynomialRing_0_m,
00287
00288 *MinusP_polynom,
00289 *Q_polynom,
00290 *st_m,
00291 fullRankStatistic,
00292 subRankStatistic,
00293 fullQuadricsStatistic,
00294 subQuadricsStatistic,
00295 liftStatistic );
00296
00297 delete MinusP_polynom; MinusP_polynom=NULL;
00298 delete Q_polynom; Q_polynom=NULL;
00299 }
00300 else if (cfparams_m->getEpsPrecision()==1)
00301 {
00302 D_Frommer::TPolynomXY * MinusP_polynom = new D_Frommer::TPolynomXY ( cfparams_m->getPolynomialDegree() );
00303 D_Frommer::TPolynomXY * Q_polynom = new D_Frommer::TPolynomXY ( cfparams_m->getPolynomialDegree() );
00304
00305
00306
00307
00308
00309
00310 copyPolynomWithGivenEpsPrecision( *minusP_polynom_m,*MinusP_polynom, cfparams_m->getEpsPrecision() );
00311 copyPolynomWithGivenEpsPrecision( *Q_polynom_m,*Q_polynom , cfparams_m->getEpsPrecision() );
00312
00313 performRegularExperiment(
00314 *frommer2_m,
00315
00316 *frommer2_m,
00317 polynomialRing_2_m,
00318
00319 *MinusP_polynom,
00320 *Q_polynom,
00321 *st_m,
00322 fullRankStatistic,
00323 subRankStatistic,
00324 fullQuadricsStatistic,
00325 subQuadricsStatistic,
00326 liftStatistic );
00327
00328 delete MinusP_polynom; MinusP_polynom=NULL;
00329 delete Q_polynom; Q_polynom=NULL;
00330 }
00331 else
00332 {
00333 performRegularExperiment( *epsRegularFrommer_m,
00334 *frommer2_m,
00335 polynomialRing_m,
00336 *minusP_polynom_m,
00337 *Q_polynom_m,
00338 *st_m,
00339 fullRankStatistic,
00340 subRankStatistic,
00341 fullQuadricsStatistic,
00342 subQuadricsStatistic,
00343 liftStatistic );
00344 }
00345
00346
00347
00348
00349 }
00350 else
00351 {
00352
00353
00354
00355
00356
00357 D_Frommer0::TPolynomXY * randomMinusP_polynom = new D_Frommer0::TPolynomXY ( cfparams_m->getPolynomialDegree() );
00358 D_Frommer0::TPolynomXY * randomQ_polynom = new D_Frommer0::TPolynomXY ( cfparams_m->getPolynomialDegree() );
00359
00360
00361 assert ( randomExperimentWellDefined() );
00362
00363
00364
00365 copyPolynomWithGivenEpsPrecision( *minusP_polynom_m,*randomMinusP_polynom, 0 );
00366 copyPolynomWithGivenEpsPrecision( *Q_polynom_m,*randomQ_polynom , 0 );
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 frommer0_m->setName("frommer0_m");
00377 frommer2_m->setName("frommer2_m");
00378
00379
00380 CFRandomExperiment<variant, D_Frommer0::TPolynomXY, D_Frommer0, D_Frommer> cfRandomExperiment(*this,
00381 cfparams_m,
00382 *frommer0_m,
00383 *frommer2_m,
00384 fullRankStatistic,
00385 subRankStatistic,
00386 fullQuadricsStatistic,
00387 subQuadricsStatistic,
00388 liftStatistic,
00389 *randomMinusP_polynom,
00390 *randomQ_polynom);
00391
00392
00393 cfRandomExperiment.performRandomExperiment();
00394 *st_m=cfRandomExperiment.getStatisticRef();
00395
00396
00397
00398 delete randomMinusP_polynom; randomMinusP_polynom=NULL;
00399 delete randomQ_polynom; randomQ_polynom=NULL;
00400
00401 #ifdef FORMULA23_TIMER
00402 cfRandomExperiment.printTimings(cfparams_m->getOsRef() );
00403 #endif
00404
00405
00406 }
00407
00408 }
00409
00410
00411 template <int variant>
00412 void CenterFocusExperiment<variant>::printAlgorithmTimings(std::ostream & os) const
00413 {
00414 #ifdef ALGORITHM_TIMER
00415 os << "-- Algorithm Timings " << std::endl;
00416 os << "--hamiltonianCheckTime : " << hamiltonianCheckTime_m << "sec. " << std::endl;
00417 os << "--computeVanishedVocalValuesTime : " << computeVanishedVocalValuesTime_m << "sec. " << std::endl;
00418 os << "--computeJacobianMatrixTime_m : " << computeJacobianMatrixTime_m << "sec. " << std::endl;
00419 os << "--computeJacobianRankTime_m : " << computeJacobianRankTime_m << "sec. " << std::endl;
00420 os << "--checkSmoothnessTime_m : " << checkSmoothnessTime_m << "sec. " << std::endl;
00421 os << "--computeQuadricsTime_m : " << computeQuadricsTime_m << "sec. " << std::endl;
00422 os << "--addQuadricstatistic_m : " << addQuadricstatistic_m << "sec. " << std::endl;
00423 os << "--computeAllFocalValuesTime_m : " << computeAllFocalValuesTime_m << "sec. " << std::endl;
00424
00425
00426 os << "--smoothnessTest_getJacobianKernelTime : " << smoothnessTest_getJacobianKernelTime_m << "sec. " << std::endl;
00427 os << "--smoothnessTest_copyPolynomTime_m : " << smoothnessTest_copyPolynomTime_m << "sec. " << std::endl;
00428 os << "--smoothnessTest_computeFocalValuesTime_m : " << smoothnessTest_computeFocalValuesTime_m << "sec. " << std::endl;
00429 os << "--smoothnessTest_computeJacobiMatrixTime : " << smoothnessTest_computeJacobiMatrixTime_m << "sec. " << std::endl;
00430 os << "--smoothnessTest_solveLGSTime : " << smoothnessTest_solveLGSTime_m << "sec. " << std::endl;
00431 os << "--smoothnessTest_addEpsVector_m : " << smoothnessTest_addEpsVector_m << "sec. " << std::endl;
00432
00433
00434 #endif
00435 }
00436
00437
00438 template <int variant>
00439 void CenterFocusExperiment<variant>::printStageTimings(std::ostream & os) const
00440 {
00441 os << "\n--FrommerAlgorithm timings:\n" ;
00442 frommer0_m->printStageTimings( os );
00443 frommer2_m->printStageTimings( os );
00444 epsRegularFrommer_m->setName("epsRegularFrommer_m");
00445 epsFrommer_m->setName("epsFrommer_m");
00446 epsRegularFrommer_m->printStageTimings( os );
00447 epsFrommer_m->printStageTimings( os );
00448 }
00449
00450
00451 template <int variant>
00452 DStatistic & CenterFocusExperiment<variant>::getStatisticRef()
00453 {
00454 return *st_m;
00455 }
00456
00457
00458
00459
00467 template <int variant>
00468 template < class PolynomXY_Type, class TFrommer_Type1, class TFrommer_Type2>
00469 void CenterFocusExperiment<variant>::performRegularExperiment( TFrommer_Type1 &frommer1,
00470 TFrommer_Type2 &frommer2,
00471 const PolynomialRing<PolynomXY_Type, typename TFrommer_Type1::RingType> & polynomialRing,
00472 const PolynomXY_Type & minusP_polynom,
00473 const PolynomXY_Type & Q_polynom,
00474 DStatistic & st,
00475 RankStatistic & fullRankStatistic,
00476 RankStatistic & subRankStatistic,
00477 LiftAndQuadricsStatistic & fullQuadricsStatistic,
00478 LiftAndQuadricsStatistic & subQuadricsStatistic,
00479 FailedLiftStatistic & liftStatistic
00480 )
00481 {
00482
00483
00484
00485 inputPointCounter_m++;
00486
00487 #ifdef SAFE
00488 #endif
00489
00490 #ifdef ALGORITHM_TIMER
00491 tim_1_m.clear();
00492 tim_1_m.start();
00493
00494 #endif
00495
00496 bool isHaM= isHamiltonianComponent(minusP_polynom, Q_polynom,frommer1.getRingRef() );
00497
00498 #ifdef ALGORITHM_TIMER
00499 tim_1_m.stop();
00500 hamiltonianCheckTime_m += tim_1_m.usertime();
00501 tim_1_m.clear();
00502 #endif
00503
00504 if (isHaM )
00505 hamiltonianPointCount_g++;
00506
00507
00508
00509 if ( cfparams_m->getHamiltonianComponentSwitch() == 0
00510 || (cfparams_m->getHamiltonianComponentSwitch()== -1 && ! isHaM )
00511 || (cfparams_m->getHamiltonianComponentSwitch()== 1 && isHaM )
00512 )
00513 {
00514
00515 #ifdef ALGORITHM_TIMER
00516 tim_1_m.start();
00517 #endif
00518
00521 frommer1.setPolynoms( &minusP_polynom, &Q_polynom );
00522 frommer1.doit( false );
00523
00524
00526 int successiveVanishedFocalValuesCount = frommer1.getSuccessiveVanishedFocalValuesCount();
00527
00528 #ifdef ALGORITHM_TIMER
00529 tim_1_m.stop();
00530 computeVanishedVocalValuesTime_m += tim_1_m.usertime();
00531 tim_1_m.clear();
00532 #endif
00533
00534 if ( cfparams_m->showProgress() && inputPointCounter_m % cfparams_m->getProgressStepWidth() == 0 )
00535 cerr << ".";
00536
00537 #ifdef DEBUG
00538 std::cerr << "regular frommer finished" << std::endl ;
00539 #endif
00540
00541 st.addStatistic( successiveVanishedFocalValuesCount ) ;
00542
00543 typedef TMatrix< typename D_Eps_Frommer::RingType::FieldType > MatrixTmpType;
00544 typedef Matrix3D<MatrixTmpType> Matrix3DType;
00545 typedef FocalValuesInfo<typename TFrommer_Type1::RingType::ElementType > FocalValuesInfoType;
00546 typedef CFJacobianInfo<MatrixTmpType> CFJacobianInfoType ;
00547 typedef CFQuadricsResult<MatrixTmpType, Matrix3DType> CFQuadricsResultType ;
00548 typedef CPointLiftInfo<D_Eps_Frommer, CFEpsPointType> CPointLiftInfoType;
00549
00550 typedef CFSingleResult< PolynomXY_Type,
00551 FocalValuesInfoType,
00552 CFJacobianInfoType,
00553 CFQuadricsResultType,
00554 CPointLiftInfoType >
00555
00556 CFSingleResultType;
00557
00558
00559 MatrixTmpType* full_jacobi_matrix = NULL;
00560 MatrixTmpType* sub_jacobi_matrix = NULL;
00561
00562
00563
00564 #ifdef SAFE
00565 assert (cfparams_m->requiredVanishedFocalValuesNum()<= cfparams_m->getMaxFocalValuesToCompute());
00566 #endif
00567
00568
00569 if ( successiveVanishedFocalValuesCount >= cfparams_m->requiredVanishedFocalValuesNum() )
00570 {
00571
00572
00573
00574
00575
00576
00577
00578
00579 #ifdef ALGORITHM_TIMER
00580 tim_1_m.start();
00581 #endif
00582 full_jacobi_matrix = createJacobiMatrix( frommer2 ,
00583
00584 minusP_polynom,
00585 Q_polynom,
00586 successiveVanishedFocalValuesCount,
00587 allCoeffVariablesOrder_m
00588 );
00589
00590 #ifdef ALGORITHM_TIMER
00591 tim_1_m.stop();
00592
00593 computeJacobianMatrixTime_m += tim_1_m.usertime();
00594 tim_1_m.clear();
00595 tim_1_m.start();
00596 #endif
00597 unsigned int fullJacobiMatrixRank = full_jacobi_matrix->getRank( );
00598
00599 #ifdef ALGORITHM_TIMER
00600 tim_1_m.stop();
00601
00602 computeJacobianRankTime_m += tim_1_m.usertime();
00603 tim_1_m.clear();
00604 #endif
00605
00606 unsigned int subJacobiMatrixRank ;
00607
00608 if (allCoeffVariablesOrder_m.size()>subCoeffVariablesOrder_m.size() && subCoeffVariablesOrder_m.size()>0 )
00609 {
00610 sub_jacobi_matrix = createJacobiMatrix( frommer2 ,
00611
00612 minusP_polynom, Q_polynom,
00613 successiveVanishedFocalValuesCount,
00614 subCoeffVariablesOrder_m );
00615
00616 subJacobiMatrixRank = sub_jacobi_matrix->getRank( );
00617 }
00618 else if (allCoeffVariablesOrder_m.size()==subCoeffVariablesOrder_m.size())
00619 {
00620 sub_jacobi_matrix = full_jacobi_matrix;
00621 subJacobiMatrixRank = fullJacobiMatrixRank;
00622 }
00623 else
00624 {
00625 sub_jacobi_matrix = NULL;
00626 subJacobiMatrixRank = 0;
00627 }
00628 CFJacobianInfoType jacobianInfo (full_jacobi_matrix, fullJacobiMatrixRank );
00629 CFJacobianInfoType subJacobianInfo(sub_jacobi_matrix, subJacobiMatrixRank );
00630
00631
00632
00633 fullRankStatistic.addRankStatistic( successiveVanishedFocalValuesCount,
00634 fullJacobiMatrixRank
00635 );
00636 if (subCoeffVariablesOrder_m.size()>0)
00637 {
00638 subRankStatistic.addRankStatistic( successiveVanishedFocalValuesCount,
00639 subJacobiMatrixRank
00640 );
00641 }
00642
00643 #ifdef DEBUG
00644 cerr << " full_jacobi_matrix " << fullJacobiMatrixRank << endl;
00645 cerr << " sub_jacobi_matrix " << sub_jacobi_matrix << endl;
00646 #endif
00647
00648
00649 if ( fullJacobiMatrixRank >= cfparams_m->getRequiredFullJacobianMinRank() &&
00650 subJacobiMatrixRank >= cfparams_m->getRequiredSubJacobianMinRank() )
00651 {
00652
00653
00654
00655
00656 #ifdef DEBUG
00657 PolynomXY_Type PE= polynomialRing.addInv(minusP_polynom );
00658 resultPrinter_m->printPolynoms( std::cerr,
00659 PE, Q_polynom ,
00660 string("Debugausgabe: Punkt gefunden 'compute p, q' \n"));
00661 #endif
00662
00663
00664
00665 #ifdef ALGORITHM_TIMER
00666 tim_1_m.start();
00667 #endif
00668
00669
00670 CFQuadricsResultType fullQuadrics(*full_jacobi_matrix, allCoeffVariablesOrder_m);
00671
00672
00673 computeQuadric( minusP_polynom, Q_polynom,
00674 *full_jacobi_matrix,
00675 allCoeffVariablesOrder_m,
00676 fullQuadrics );
00677
00678 #ifdef ALGORITHM_TIMER
00679 tim_1_m.stop();
00680 computeQuadricsTime_m += tim_1_m.usertime();
00681 tim_1_m.clear();
00682 #endif
00683
00684
00685
00686
00687
00688 CPointLiftInfoType pointFullLiftInfo( cfparams_m->getMaxLift(), *epsFrommer_m );
00689 CPointLiftInfoType pointSubLiftInfo ( cfparams_m->getMaxLift(), *epsFrommer_m );
00690 pointSubLiftInfo.setLiftTestPassed(true);
00691
00692 #ifdef ALGORITHM_TIMER
00693 tim_1_m.start();
00694 #endif
00695
00696 static EmptyFailedLiftStatistic emptyLiftStatistic;
00697
00698
00699 #ifdef CF_TEST
00700 if (fullQuadrics.getQuadricsBig()->getZNum()>0 && cfparams_m->getLiftTrials()==0)
00701 {
00702
00705 CPointLiftInfoType pointFullLiftInfoTest( cfparams_m->getMaxLift(), *epsFrommer_m );
00706 performLiftTrials( 2,
00707 1,
00708 frommer2,
00709 minusP_polynom, Q_polynom,
00710 *full_jacobi_matrix,
00711 allCoeffVariablesOrder_m,
00712 pointFullLiftInfoTest,
00713 emptyLiftStatistic );
00714
00715 assert(! pointFullLiftInfoTest.liftTestPassed() && pointFullLiftInfoTest.getMinFailedLiftNr()==2);
00716 }
00717 #endif
00718
00719 performLiftTrials( cfparams_m->getMaxLift(),
00720 cfparams_m->getLiftTrials(),
00721 frommer2,
00722 minusP_polynom, Q_polynom,
00723 *full_jacobi_matrix,
00724 allCoeffVariablesOrder_m,
00725 pointFullLiftInfo,
00726 emptyLiftStatistic );
00727 #ifdef ALGORITHM_TIMER
00728 tim_1_m.stop();
00729 checkSmoothnessTime_m += tim_1_m.usertime();
00730 tim_1_m.clear();
00731 #endif
00732
00733 if (cfparams_m->getRequestedLiftPointNum()>0)
00734 {
00735 liftTest_computeLiftPoints( cfparams_m->getMaxLift(),
00736 cfparams_m->getLiftTrials(),
00737 frommer2,
00738 minusP_polynom, Q_polynom,
00739 *full_jacobi_matrix,
00740 allCoeffVariablesOrder_m,
00741 pointFullLiftInfo.getLiftPointsRef() );
00742 }
00743
00744
00745 #ifdef ALGORITHM_TIMER
00746 tim_1_m.start();
00747 #endif
00748
00751 if ( ! pointFullLiftInfo.liftTestPassed() )
00752 {
00753
00754 FailedLiftStatistic localLiftStatistic( cfparams_m->getExhaustiveMaxLift(), cfparams_m->getExhaustiveLiftTrials() ) ;
00755
00756 pointFullLiftInfo.clear();
00757 pointFullLiftInfo.setMaxLift( cfparams_m->getExhaustiveMaxLift() );
00758
00759
00760
00761
00763 if (cfparams_m->logExtendedLiftStatistic() )
00764 {
00765
00766 FailedLiftStatistic localLiftStatistic( cfparams_m->getExhaustiveMaxLift(), cfparams_m->getExhaustiveLiftTrials() ) ;
00767
00768 performLiftTrials( cfparams_m->getExhaustiveMaxLift(),
00769 cfparams_m->getExhaustiveLiftTrials(),
00770 frommer2,
00771 minusP_polynom, Q_polynom,
00772 *full_jacobi_matrix,
00773 allCoeffVariablesOrder_m,
00774 pointFullLiftInfo,
00775 localLiftStatistic );
00776
00777 localLiftStatistic.setTestFailureCount();
00778 liftStatistic += localLiftStatistic;
00779
00780 g_extLiftStatistic.logFailedLift(localLiftStatistic, fullJacobiMatrixRank);
00781 }
00782 else
00783 {
00784 performLiftTrials( cfparams_m->getExhaustiveMaxLift(),
00785 cfparams_m->getExhaustiveLiftTrials(),
00786 frommer2,
00787 minusP_polynom, Q_polynom,
00788 *full_jacobi_matrix,
00789 allCoeffVariablesOrder_m,
00790 pointFullLiftInfo,
00791 emptyLiftStatistic );
00792 }
00793
00794 }
00795 #ifdef ALGORITHM_TIMER
00796 tim_1_m.stop();
00797 checkSmoothnessTime_m += tim_1_m.usertime();
00798 tim_1_m.clear();
00799 tim_1_m.start();
00800 #endif
00801 fullQuadricsStatistic.addLiftAndQuadricsStatistic( successiveVanishedFocalValuesCount ,
00802 fullJacobiMatrixRank,
00803 fullQuadrics.getQuadricsBig()->getZNum() ,
00804 pointFullLiftInfo.liftTestPassed(),
00805 pointFullLiftInfo.getMinFailedLiftNr()
00806 );
00807
00808
00809
00810 if ( allCoeffVariablesOrder_m.size()==subCoeffVariablesOrder_m.size() )
00811 {
00812 subQuadricsStatistic.addLiftAndQuadricsStatistic( successiveVanishedFocalValuesCount ,
00813 fullJacobiMatrixRank,
00814 fullQuadrics.getQuadricsBig()->getZNum(),
00815 pointSubLiftInfo.liftTestPassed() );
00816 }
00817
00818
00819
00820
00821
00822
00823
00824 #ifdef ALGORITHM_TIMER
00825 tim_1_m.stop();
00826 addQuadricstatistic_m += tim_1_m.usertime();
00827 tim_1_m.clear();
00828 #endif
00829 CFQuadricsResultType subQuadrics(*sub_jacobi_matrix, subCoeffVariablesOrder_m);
00830 if ( allCoeffVariablesOrder_m.size()>subCoeffVariablesOrder_m.size() && subCoeffVariablesOrder_m.size()>0)
00831 {
00832
00833 #ifdef ALGORITHM_TIMER
00834 tim_1_m.start();
00835
00836 #endif
00837 computeQuadric( minusP_polynom,
00838 Q_polynom,
00839 *sub_jacobi_matrix,
00840 subCoeffVariablesOrder_m,
00841 subQuadrics );
00842 #ifdef ALGORITHM_TIMER
00843 tim_1_m.stop();
00844 computeQuadricsTime_m += tim_1_m.usertime();
00845 tim_1_m.clear();
00846 tim_1_m.start();
00847 #endif
00848 subQuadricsStatistic.addLiftAndQuadricsStatistic( successiveVanishedFocalValuesCount ,
00849 subJacobiMatrixRank,
00850 subQuadrics.getQuadricsBig()->getZNum() ,
00851 pointSubLiftInfo.liftTestPassed()
00852 );
00853 #ifdef ALGORITHM_TIMER
00854 tim_1_m.stop();
00855 addQuadricstatistic_m += tim_1_m.usertime();
00856 tim_1_m.clear();
00857 #endif
00858
00859 }
00860
00861 uint16_t fullQuadricsRank=fullQuadrics.getQuadricsBig()->getZNum() ;
00862
00863 if (cfparams_m->pointPassedFilter( successiveVanishedFocalValuesCount,
00864 fullJacobiMatrixRank,
00865 pointFullLiftInfo.liftTestPassed(),
00866 fullQuadricsRank)
00867 )
00868
00870 if ( ( fullJacobiMatrixRank + fullQuadricsRank ) >= cfparams_m->getMinComponentCodim() ||
00871 ( (( fullJacobiMatrixRank + 1) >= cfparams_m->getMinComponentCodim()) && ! pointFullLiftInfo.liftTestPassed() )
00872 )
00873 {
00874
00875 if ( !cfparams_m->nonSmoothPointPassFilter() || ! pointFullLiftInfo.liftTestPassed() )
00876 {
00877
00878 CFSingleResultType singleResult;
00879
00880 PolynomXY_Type pPolynom= polynomialRing.addInv( minusP_polynom );
00881 singleResult.setPolynoms( pPolynom, Q_polynom );
00882
00883 singleResult.setSubJacobianInfo(&subJacobianInfo);
00884 singleResult.setJacobianInfo( &jacobianInfo );
00885
00886 singleResult.setFullQuadrics( &fullQuadrics );
00887 singleResult.setSubQuadrics( &subQuadrics );
00888
00889 singleResult.setFullLiftInfo( & pointFullLiftInfo );
00890
00894
00895 #ifdef SAFE
00896 assert( frommer1.getRingRef().getEpsPrecision() >= cfparams_m->getEpsPrecision() );
00897 #endif
00898 #ifdef ALGORITHM_TIMER
00899 tim_1_m.start();
00900 #endif
00901
00902 frommer1.setPolynoms( &minusP_polynom, &Q_polynom );
00903 frommer1.doit( true );
00904
00905
00906 FocalValuesInfoType fvInfo ;
00907 frommer1.getComputedFocalValues( fvInfo.getFocalValuesRef() );
00908
00909 #ifdef ALGORITHM_TIMER
00910 tim_1_m.stop();
00911 computeAllFocalValuesTime_m += tim_1_m.usertime();
00912 tim_1_m.clear();
00913 #endif
00914
00915 fvInfo.setVanishedFocalValuesNum( frommer1.getSuccessiveVanishedFocalValuesCount() );
00916 singleResult.setFocalValuesInfo( &fvInfo );
00917 singleResult.setPointID(cfparams_m->getPointID() );
00918
00919
00920
00921
00922 resultPrinter_m->beginSingleResultPrint();
00923
00924 if (cfparams_m->pureSmoothnessTest())
00925 singleResult.printSmoothnessTestResult( cfparams_m->getOsRef() );
00926 else
00927 singleResult.print( cfparams_m->getOsRef() );
00928 }
00929 }
00930
00931 }
00932 delete full_jacobi_matrix; full_jacobi_matrix = NULL;
00933 if ( allCoeffVariablesOrder_m.size()>subCoeffVariablesOrder_m.size() )
00934 {
00936 if (sub_jacobi_matrix!=NULL)
00937 {
00938 delete sub_jacobi_matrix; sub_jacobi_matrix=NULL;
00939 }
00940
00941 }
00942 }
00943 }
00944 }
00945
00946
00950 template <int variant>
00951 void CenterFocusExperiment<variant>::liftTest_invertCoeffVectorInPlace (TVector<typename D_Eps_Frommer::RingType::FieldType> & vector)
00952 {
00953
00954 for ( size_t var=0; var < vector.getSize(); var++ )
00955 {
00956
00957 epsFieldRef_m.addInvInPlace(vector.getValRef(var) );
00958 }
00959 }
00960
00961
00965 template <int variant>
00966 typename CenterFocusExperiment<variant>::CFEpsPointType
00967 CenterFocusExperiment<variant>::liftTest_addEpsVector ( const CFEpsPointType & polynomPair,
00968
00969 TVector<typename D_Eps_Frommer::RingType::FieldType> & vector,
00970 const list<CoeffListEntry> & varCoeffOrder,
00971 int epsFactorExp )
00972
00973 {
00974
00975 #ifdef SAFE
00976 assert(varCoeffOrder.size()==vector.getSize() );
00977 #endif
00978 CFEpsPointType res(polynomPair);
00979
00980 int row=0;
00981
00982 for ( std::list<CoeffListEntry>::const_iterator coeffListIt = varCoeffOrder.begin();
00983 coeffListIt != varCoeffOrder.end();
00984 coeffListIt++, row++
00985 )
00986 {
00987 typename D_Eps_Frommer::RingType::ScalarType scalarCoeff = vector.getVal(row);
00988
00989
00990 if ( (*coeffListIt).isPMonom() )
00991 {
00992 epsFieldRef_m.addInvInPlace(scalarCoeff);
00993 res.first.getCoeffRef( (*coeffListIt).x_exp, (*coeffListIt).y_exp)
00994 .setValue( epsFactorExp, scalarCoeff.getX() );
00995 }
00996 else
00997 {
00998 res.second.getCoeffRef( (*coeffListIt).x_exp, (*coeffListIt).y_exp)
00999 .setValue(epsFactorExp, scalarCoeff.getX() );
01000 }
01001 }
01002 return res;
01003 }
01004
01006 template <int variant>
01007 void CenterFocusExperiment<variant>::liftTest_eraseEpsPartInPlace (CFEpsPointType & polynomPair,
01008 unsigned short epsPrecision)
01009 {
01010 for (int deg = polynomPair.first.getDegree(); deg>=0; deg--)
01011 {
01012 for (int x_exp=deg; x_exp>=0; x_exp--)
01013 {
01014 polynomPair.first.getCoeffRef( x_exp, deg-x_exp).setValue(epsPrecision, 0 );
01015 polynomPair.second.getCoeffRef( x_exp, deg-x_exp).setValue(epsPrecision, 0 );
01016 }
01017 }
01018 }
01019
01020
01023 template <int variant>
01024 TVector<typename D_Eps_Frommer::RingType::FieldType>
01025 CenterFocusExperiment<variant>::liftTest_getFocalValuesEpsPart(D_Eps_Frommer & epsFrommer, unsigned int epsPart)
01026 {
01027
01028
01029
01030
01031 TVector<typename D_Eps_Frommer::RingType::FieldType> res( epsFrommer.getSuccessiveVanishedFocalValuesCount(),
01032 &epsFieldRef_m );
01033
01034 for ( int currFocalValue = 1; currFocalValue <= epsFrommer.getSuccessiveVanishedFocalValuesCount() ; currFocalValue++)
01035 {
01036 typename D_Eps_Frommer::TPolynomXY::CoefficientType focalValue = epsFrommer.getFocalValue(currFocalValue);
01037 #ifdef SAFE
01038 typename D_Eps_Frommer::RingType::ScalarType zero(0);
01039 for (int i=0; i<epsPart; i++ )
01040 assert( focalValue.getValue(i) == zero );
01041 #endif
01042
01044 res.setVal( currFocalValue-1, focalValue.getValue(epsPart) );
01045 }
01046 return res;
01047 }
01048
01053
01054 template <int variant>
01055 bool
01056 CenterFocusExperiment<variant>::liftTest_solveLGS( const TMatrix< typename D_Eps_Frommer::RingType::FieldType > & jacobiMatrix,
01057 const TVector<typename D_Eps_Frommer::RingType::FieldType> & rightHandSide,
01058 TVector<typename D_Eps_Frommer::RingType::FieldType> & result )
01059
01060 {
01061
01062 bool transpose =false;
01063
01064 typedef FFpackMatrix< Modular<double> > FFpackMatrixType;
01065 FFpackMatrixType FPJacobi ( jacobiMatrix, jacobiMatrix.getRing()->getCharacteristic(), transpose=false );
01066
01067
01068 FFpackMatrixType::ElementType rightHandSide_FFPack[ rightHandSide.getSize() ];
01069
01070
01071
01072 for (size_t pos=0; pos < rightHandSide.getSize(); pos++)
01073 {
01074 rightHandSide_FFPack[pos] = rightHandSide.getVal(pos);
01075 }
01076 int jacobiRank=-1;
01077
01078
01079 typename FFpackMatrixType::ElementType * ffpackLGSsol = FPJacobi.solveLGS( rightHandSide_FFPack ,jacobiRank);
01080
01081 if (ffpackLGSsol)
01082 {
01083
01084 FFpackMatrixType FPJacobi2 ( jacobiMatrix, jacobiMatrix.getRing()->getCharacteristic(), transpose=false );
01085
01086 int rightKernelDim = -1;
01087 FFpackMatrixType* rightKernel = FPJacobi2.getRightKernel( rightKernelDim ) ;
01088
01089 #ifdef CF_TEST
01090
01091
01092
01093 #endif
01094
01095 for (size_t pos=0; pos < result.getSize(); pos++)
01096 {
01097 result.setVal(pos, ffpackLGSsol[pos] );
01098 }
01099
01100 for (int col=0; col< rightKernel->getColNum(); col++)
01101 {
01102 typename D_Eps_Frommer::RingType::ScalarType factor = random( cfparams_m->getRandomSeedPtr(), cfparams_m->getFieldCharacteristic()-1) ;
01103
01104 for (size_t var=0; var < result.getSize(); var++ )
01105 {
01106 epsFieldRef_m.addInPlace( result.getValRef(var),
01107 epsFieldRef_m.multiply( factor,
01108 rightKernel->getVal( var, col )
01109 )
01110 );
01111 }
01112 }
01113 delete rightKernel;
01114 delete[] ffpackLGSsol;
01115 return true;
01116 }
01117 else
01118 {
01119 return false;
01120 }
01121 }
01122
01123
01124
01132 template <int variant>
01133 template <class PolynomXY_Type, class Matrix_Type, class TFrommer_Type2>
01134 bool
01135 CenterFocusExperiment<variant>::liftTest_performSingleTest( int maxLift,
01136 TFrommer_Type2 & frommer2,
01137 const PolynomXY_Type & minusP_polynom,
01138 const PolynomXY_Type & Q_polynom,
01139 const Matrix_Type & jacobiKernel,
01140 const list<CoeffListEntry> & varCoeffOrder,
01141 CFEpsPointType * & liftResultRef,
01142 int & lastLiftNr)
01143 {
01144
01145 liftResultRef = NULL;
01146 CFEpsPointType liftPolynomVorlage ( D_Eps_Frommer::TPolynomXY(cfparams_m->getPolynomialDegree() ),
01147 D_Eps_Frommer::TPolynomXY(cfparams_m->getPolynomialDegree() ) ) ;
01148
01149 int oldEpsPrecision= epsFrommer_m->getRingRef().getEpsPrecision();
01150
01151 #ifdef ALGORITHM_TIMER
01152 tim_2_m.clear();
01153 tim_2_m.start();
01154 #endif
01155
01156 copyPolynomWithGivenEpsPrecision(minusP_polynom, liftPolynomVorlage.first, 0) ;
01157 copyPolynomWithGivenEpsPrecision(Q_polynom, liftPolynomVorlage.second, 0) ;
01158
01159 polynomSetEpsPrecision(liftPolynomVorlage.first, 2);
01160 polynomSetEpsPrecision(liftPolynomVorlage.second, 2);
01161
01162 #ifdef ALGORITHM_TIMER
01163 tim_2_m.stop();
01164 smoothnessTest_copyPolynomTime_m += tim_2_m.usertime();
01165 tim_2_m.clear();
01166 #endif
01167
01168
01169
01170
01171 TVector<typename D_Eps_Frommer::RingType::FieldType> result( varCoeffOrder.size(), & (epsFieldRef_m ) );
01172 for (unsigned int col=0; col<jacobiKernel.getColNum(); col++)
01173
01174 {
01175 typename D_Eps_Frommer::RingType::ScalarType factor = random( cfparams_m->getRandomSeedPtr(), cfparams_m->getFieldCharacteristic()-1) ;
01176
01177 for (size_t var=0; var < result.getSize(); var++ )
01178 {
01179 epsFieldRef_m.addInPlace( result.getValRef(var),
01180 epsFieldRef_m.multiply( factor,
01181 jacobiKernel.getVal( var, col )
01182 )
01183 );
01184 }
01185 }
01186
01188
01189 bool singleLiftTestPassed = true;
01190 int epsPrecision = 1;
01191 lastLiftNr = epsPrecision ;
01192 for ( epsPrecision=2; epsPrecision <= maxLift; epsPrecision++)
01193 {
01194
01195 lastLiftNr = epsPrecision ;
01196 #ifdef ALGORITHM_TIMER
01197 tim_2_m.clear();
01198 tim_2_m.start();
01199 #endif
01200 liftPolynomVorlage= liftTest_addEpsVector( liftPolynomVorlage,
01201 result, varCoeffOrder,
01202 epsPrecision - 1 );
01203
01204 #ifdef ALGORITHM_TIMER
01205 tim_2_m.stop();
01206 smoothnessTest_addEpsVector_m += tim_2_m.usertime();
01207 tim_2_m.clear();
01208 tim_2_m.start();
01209 #endif
01210
01211 epsFrommer_m->getRingRef().setEpsPrecision( epsPrecision );
01212
01213 polynomSetEpsPrecision(liftPolynomVorlage.first, epsPrecision);
01214 polynomSetEpsPrecision(liftPolynomVorlage.second, epsPrecision);
01215
01216 epsFrommer_m->setPolynoms( &(liftPolynomVorlage.first), &(liftPolynomVorlage.second) );
01217 epsFrommer_m->doit(false);
01218
01219 #ifdef ALGORITHM_TIMER
01220 tim_2_m.stop();
01221 smoothnessTest_computeFocalValuesTime_m += tim_2_m.usertime();
01222 tim_2_m.clear();
01223 #endif
01224
01225 TVector<typename D_Eps_Frommer::RingType::FieldType> rightHandSide = liftTest_getFocalValuesEpsPart(*epsFrommer_m, epsPrecision);
01226
01227
01228 typedef TMatrix< typename D_Eps_Frommer::RingType::FieldType > MatrixTmpType;
01229
01230
01231
01232
01233
01234
01235
01236
01237 #ifdef ALGORITHM_TIMER
01238 tim_2_m.start();
01239 #endif
01240 MatrixTmpType * interJacobiMatrix= createJacobiMatrix( frommer2,
01241 liftPolynomVorlage.first,
01242 liftPolynomVorlage.second,
01243 epsFrommer_m->getSuccessiveVanishedFocalValuesCount(),
01244 varCoeffOrder
01245 );
01246 #ifdef ALGORITHM_TIMER
01247 tim_2_m.stop();
01248 smoothnessTest_computeJacobiMatrixTime_m += tim_2_m.usertime();
01249 tim_2_m.clear();
01250 tim_2_m.start();
01251 #endif
01252
01253
01254 assert( interJacobiMatrix->getRowNum()>0 && interJacobiMatrix->getColNum()>0 );
01255 if ( liftTest_solveLGS (*interJacobiMatrix, rightHandSide, result) )
01256 {
01258 #ifdef SAFE
01259 TVector<typename D_Eps_Frommer::RingType::FieldType>* computedRightHandSide = (*interJacobiMatrix)*result;
01260
01261
01262 assert( (*computedRightHandSide)==rightHandSide );
01263 delete computedRightHandSide; computedRightHandSide=NULL;
01264 #endif
01265 delete interJacobiMatrix;
01266 liftTest_invertCoeffVectorInPlace(result);
01267 }
01268 else
01269 {
01270 #ifdef ALGORITHM_TIMER
01271 tim_2_m.stop();
01272 smoothnessTest_solveLGSTime_m += tim_2_m.usertime();
01273 tim_2_m.clear();
01274
01275 #endif
01276
01277 delete interJacobiMatrix;
01278 singleLiftTestPassed = false;
01279 break;
01280
01281 }
01282 #ifdef ALGORITHM_TIMER
01283 tim_2_m.stop();
01284 smoothnessTest_solveLGSTime_m += tim_2_m.usertime();
01285 tim_2_m.clear();
01286 #endif
01287 }
01288
01289 #ifdef ALGORITHM_TIMER
01290 tim_2_m.clear();
01291 tim_2_m.start();
01292 #endif
01293
01294 if (singleLiftTestPassed)
01295 {
01296 liftPolynomVorlage= liftTest_addEpsVector( liftPolynomVorlage,
01297 result,
01298 varCoeffOrder,
01299 epsPrecision - 1 );
01300 liftResultRef= new CFEpsPointType(liftPolynomVorlage);
01301 }
01302 #ifdef ALGORITHM_TIMER
01303 tim_2_m.stop();
01304 smoothnessTest_addEpsVector_m += tim_2_m.usertime();
01305 tim_2_m.clear();
01306 #endif
01307
01308
01309
01310 epsFrommer_m->getRingRef().setEpsPrecision( oldEpsPrecision );
01311
01312 return singleLiftTestPassed;
01313 }
01314
01315
01316
01320 template <int variant>
01321 template <class PolynomXY_Type, class Matrix_Type, class TFrommer_Type2, class FailedLiftStatisticType>
01322 bool CenterFocusExperiment<variant>::performLiftTrials( int maxLift,
01323 int liftTrials,
01324 TFrommer_Type2 & frommer2,
01325 const PolynomXY_Type & minusP_polynom,
01326 const PolynomXY_Type & Q_polynom,
01327 const Matrix_Type & jacobiMatrix,
01328 const list<CoeffListEntry> & varCoeffOrder,
01329 CPointLiftInfo<D_Eps_Frommer, CFEpsPointType> & pointLiftInfo,
01330 FailedLiftStatisticType & liftStatistic)
01331 {
01332 bool liftTestPassed = true;
01333 if (liftTrials==0)
01334 {
01335 pointLiftInfo.setLiftTestPassed( liftTestPassed );
01336 return liftTestPassed;
01337 }
01338
01339 assert(varCoeffOrder.size()>0 );
01340
01341 #ifdef SAFE
01342
01343 assert(jacobiMatrix.getRowNum()>0 );
01344 assert(jacobiMatrix.getColNum()>0 );
01345 #endif
01346
01347 CFEpsPointType * liftResult;
01348
01349
01350
01351 int minFailedLiftNr = cfparams_m->getMaxLift() + 1;
01352 int maxFailedLiftNr = -1;
01353 int firstFailedTrialNr = -1;
01354 int failedTrialCount = 0;
01355
01356 #ifdef ALGORITHM_TIMER
01357 tim_2_m.clear();
01358 tim_2_m.start();
01359 #endif
01360
01361 Matrix_Type* jacobiKernel = jacobiMatrix.getKernel();
01362 #ifdef ALGORITHM_TIMER
01363 tim_2_m.stop();
01364 smoothnessTest_getJacobianKernelTime_m += tim_2_m.usertime();
01365 tim_2_m.clear();
01366
01367 #endif
01368
01369 for (int currentLiftTrial=1; currentLiftTrial <= liftTrials; currentLiftTrial++)
01370 {
01371 int lastLiftNr;
01372 if (! liftTest_performSingleTest(maxLift, frommer2, minusP_polynom, Q_polynom, *jacobiKernel, varCoeffOrder, liftResult, lastLiftNr) )
01373 {
01374 liftStatistic.logFailedLift(lastLiftNr);
01375 failedTrialCount++;
01376
01377 if (liftTestPassed )
01378 firstFailedTrialNr= currentLiftTrial;
01379 liftTestPassed = false;
01380 if (lastLiftNr<minFailedLiftNr )
01381 {
01382 minFailedLiftNr = lastLiftNr;
01383 }
01384 if (maxFailedLiftNr<lastLiftNr)
01385 {
01386 maxFailedLiftNr=lastLiftNr;
01387 }
01388
01389 }
01390 else
01391 delete liftResult;
01392 }
01393
01394 delete jacobiKernel;
01395
01396 pointLiftInfo.setLiftTestPassed( liftTestPassed );
01397
01398 if (! liftTestPassed)
01399 {
01400
01401 pointLiftInfo.setMinFailedLiftNr( minFailedLiftNr );
01402 pointLiftInfo.setMaxFailedLiftNr( maxFailedLiftNr );
01403 pointLiftInfo.setFirstFailedTrialNr( firstFailedTrialNr );
01404 pointLiftInfo.setFailedTrialCount( failedTrialCount );
01405 }
01406
01407
01408 return liftTestPassed;
01409 }
01410
01411
01412 template <int variant>
01413 template <class PolynomXY_Type, class Matrix_Type, class TFrommer_Type2>
01414 void CenterFocusExperiment<variant>::liftTest_computeLiftPoints( int maxLift,
01415 int liftTrials,
01416 TFrommer_Type2 & frommer2,
01417 const PolynomXY_Type & minusP_polynom,
01418 const PolynomXY_Type & Q_polynom,
01419 const Matrix_Type & jacobiMatrix,
01420 const list<CoeffListEntry> & varCoeffOrder,
01421 vector<CFEpsPointType*> & liftPoints)
01422 {
01423
01424 assert(varCoeffOrder.size()>0 );
01425 #ifdef SAFE
01426 assert(jacobiMatrix.getRowNum()>0 );
01427 assert(jacobiMatrix.getColNum()>0 );
01428 #endif
01429 Matrix_Type* jacobiKernel = jacobiMatrix.getKernel();
01430
01431 for (int liftPointNr= cfparams_m->getRequestedLiftPointNum();liftPointNr>0 ;liftPointNr--)
01432 {
01433 CFEpsPointType * liftResult;
01434 for (int currentLiftTrial = 1; currentLiftTrial <= liftTrials; currentLiftTrial++)
01435 {
01436 int currentLift;
01437 if ( liftTest_performSingleTest(maxLift, frommer2, minusP_polynom, Q_polynom, *jacobiKernel, varCoeffOrder, liftResult, currentLift) )
01438 {
01439 liftPoints.push_back(liftResult);
01440
01441 break;
01442 }
01443 }
01444 }
01445 delete jacobiKernel;
01446 return ;
01447 }
01448
01449
01452 template <int variant>
01453 bool CenterFocusExperiment<variant>::randomExperimentWellDefined()
01454 {
01455 assert(minusP_polynom_m!=0);
01456 assert(Q_polynom_m!=0);
01457
01458 assert( minusP_polynom_m->getDegree() == Q_polynom_m->getDegree() );
01459
01460
01461 int maxEpsDegree = 0;
01462
01463 for (int deg= minusP_polynom_m->getDegree(); deg>=0; deg--)
01464 {
01465 for (int x_exp=deg; x_exp>=0; x_exp--)
01466 {
01467 #ifdef DEBUG
01468 std::cerr << "minusP_polynom_m " << minusP_polynom_m->getCoeffRef(x_exp, deg-x_exp) << std::endl;
01469 std::cerr << "Q_polynom_m " << Q_polynom_m->getCoeffRef(x_exp, deg-x_exp) << std::endl;
01470 #endif
01471 if ( ( minusP_polynom_m->getCoeffRef(x_exp, deg-x_exp) ).getEpsDegree() > maxEpsDegree )
01472
01473 maxEpsDegree = (minusP_polynom_m->getCoeffRef(x_exp, deg-x_exp)).getEpsDegree();
01474
01475 if ( (Q_polynom_m->getCoeffRef(x_exp, deg-x_exp)).getEpsDegree() > maxEpsDegree )
01476
01477 maxEpsDegree = (Q_polynom_m->getCoeffRef(x_exp, deg-x_exp)).getEpsDegree();
01478 }
01479 }
01480
01481
01482 if (maxEpsDegree ==0)
01483 return true;
01484 else
01485 {
01486 std::cerr << "maxEpsDegree " << maxEpsDegree << std::endl;
01487 std::cerr << "randomExperiment is currently not defined if polynom coefficients are in F_p[e] " << std::endl;
01488 return false;
01489 }
01490
01491
01492 }
01493
01494
01495
01519
01520
01521 template <int variant>
01522 template <class PolynomXY_Type, class Ring_Type>
01523 bool CenterFocusExperiment<variant>::isHamiltonianComponent(const PolynomXY_Type & minusP_polynom ,
01524 const PolynomXY_Type & Q_polynom,
01525 const Ring_Type &ring1)
01526 {
01527
01528 static const typename Ring_Type::ScalarType two (2);
01529 static const typename Ring_Type::ScalarType three (3);
01530
01531
01532
01533
01534 if (ring1.addInv(( minusP_polynom.getCoeff(1,1))).nearlyEqual( ring1.scalarMultiply( two,
01535 Q_polynom.getCoeff(2,0) ) ) == 0)
01536 return false;
01537
01538 if ( (Q_polynom.getCoeff(1,1)).nearlyEqual( ring1.addInv( ring1.scalarMultiply(two,
01539 minusP_polynom.getCoeff(0,2) )) ) ==0 )
01540 return false;
01541
01542 if (cfparams_m->getPolynomialDegree()==2)
01543 return true;
01544
01545 if ((ring1.addInv(minusP_polynom.getCoeff(2,1))).nearlyEqual(ring1.scalarMultiply(three,
01546 Q_polynom.getCoeff(3,0))) ==0)
01547 return false;
01548
01549 if ((ring1.addInv(minusP_polynom.getCoeff(1,2))).nearlyEqual( Q_polynom.getCoeff(2,1)) ==0 )
01550 return false;
01551
01552 if ( Q_polynom.getCoeff(1,2).nearlyEqual(ring1.scalarMultiply(three,
01553 ring1.addInv(minusP_polynom.getCoeff(0,3)))) ==0 )
01554 return false;
01555
01556 if (cfparams_m->getPolynomialDegree()==3)
01557 return true;
01558
01559 std::cerr << "HamiltonianComponentCheck is not implemented for polynomial degree greater than 3 " << std::cerr;
01560 assert(false);
01561
01562 return true;
01563 }
01564
01565
01566
01567 template <int variant>
01568 template <class Matrix_Type>
01569 inline Matrix_Type* CenterFocusExperiment<variant>::complementColumns(const Matrix_Type & mat)
01570 {
01571
01572 assert(mat.getRank()==mat.getColNum());
01573
01574 Matrix_Type* result = new Matrix_Type(mat.getRowNum(), mat.getRowNum(), mat.getRing());
01575
01576
01577 for (unsigned int row= 0 ;row< mat.getRowNum(); row++)
01578 {
01579 for (unsigned int col= 0; col< mat.getColNum(); col++)
01580 {
01581 result->setVal(row,col, mat.getVal(row,col) );
01582 }
01583 }
01584
01585 unsigned int colPos = mat.getColNum();
01586 for (unsigned int ev= 0 ;ev< mat.getRowNum() && colPos < result->getColNum(); ev++)
01587 {
01588 result->setVal(ev,colPos, 1 );
01589 if (result->getRank() > colPos )
01590 {
01591 colPos++;
01592 }
01593 else
01594 {
01595 result->setVal(ev, colPos, 0 );
01596 }
01597 }
01598
01599 return result;
01600
01601
01602 }
01603
01605 template <int variant>
01606 template <class TPolynomXY_Type, class Matrix_Type >
01607 inline Matrix3D<Matrix_Type>*
01608 CenterFocusExperiment<variant>::computeQuadric_computeAlpha(const TPolynomXY_Type &minusPpol,
01609 const TPolynomXY_Type &Qpol,
01610 const Matrix_Type & jacobiMat,
01611 const Matrix_Type & jacobiKernel ,
01612 const std::list<CoeffListEntry> & coeffOrder )
01613 {
01614
01615 typename D_Eps_Frommer::TPolynomXY quadricMinusPpol ( minusPpol.getDegree() ) ;
01616 typename D_Eps_Frommer::TPolynomXY quadricQpol ( Qpol.getDegree() );
01617
01618
01619
01620
01621
01622
01623
01624 copyPolynomWithGivenEpsPrecision(minusPpol,quadricMinusPpol,0) ;
01625 copyPolynomWithGivenEpsPrecision(Qpol,quadricQpol,0) ;
01626
01627 polynomSetEpsPrecision(quadricMinusPpol, 2);
01628 polynomSetEpsPrecision(quadricQpol, 2);
01629
01630 int oldEpsPrecision = epsFrommer_m->getRingRef().getEpsPrecision();
01631 epsFrommer_m->getRingRef().setEpsPrecision(2);
01632
01633
01634
01635
01636
01637
01638
01639
01640 int iDim = jacobiMat.getRowNum() ;
01641
01642
01643 int jkDim = jacobiKernel.getColNum();
01644
01645 Matrix3D<Matrix_Type> * alpha = new Matrix3D<Matrix_Type>(iDim,jkDim,jkDim, jacobiMat.getRing() ,"alpha ") ;
01646
01647 typename D_Eps_Frommer::RingType::FieldType const * epsField = epsFrommer_m->getRing()->getField();
01648
01649
01650 for (int j =0; j < jkDim; j++)
01651 for (int k =0; k < jkDim; k++)
01652 {
01653 typename Matrix_Type::TNum * b_j = jacobiKernel.getCol(j);
01654 typename Matrix_Type::TNum * b_k = jacobiKernel.getCol(k);
01655
01656 int row=0;
01657
01658
01659 std::list<CoeffListEntry>::const_iterator coeffListIt;
01660
01661 for ( coeffListIt = coeffOrder.begin(); coeffListIt!=coeffOrder.end(); coeffListIt++, row++)
01662 {
01663
01664 typename D_Eps_Frommer::RingType::ScalarType scalarEpsCoeffPart = epsField->add( b_j[row].getX(), b_k[row].getX() );
01665
01666
01667 if ( (*coeffListIt).isPMonom() )
01668 {
01669 epsField->addInvInPlace(scalarEpsCoeffPart);
01670 quadricMinusPpol.getCoeffRef( (*coeffListIt).x_exp, (*coeffListIt).y_exp)
01671 .setEps( scalarEpsCoeffPart.getX() );
01672 }
01673 else
01674 {
01675 quadricQpol.getCoeffRef( (*coeffListIt).x_exp, (*coeffListIt).y_exp)
01676 .setEps( scalarEpsCoeffPart.getX() );
01677 }
01678 }
01679 delete[] b_j; b_j=NULL;
01680 delete[] b_k; b_k=NULL;
01681
01682 epsFrommer_m->setPolynoms( &quadricMinusPpol, &quadricQpol );
01683 epsFrommer_m->doit(false);
01684
01685 for ( int currFocalValueNumMinusOne = 0; currFocalValueNumMinusOne <iDim; currFocalValueNumMinusOne++)
01686 {
01687 typename D_Eps_Frommer::TPolynomXY::CoefficientType currFocalValue = epsFrommer_m->getFocalValue( currFocalValueNumMinusOne+1 );
01688
01689 assert ( (int)currFocalValue.getX() == 0 );
01690 assert ( (int)currFocalValue.getEps() == 0 );
01691
01692 int epsExponent=2;
01693 alpha->setVal (currFocalValueNumMinusOne, j, k, currFocalValue.getValue( epsExponent ) );
01694 }
01695 }
01696
01697
01698 for ( int i = 0; i < iDim ; i++)
01699 for ( int j = 0; j < jkDim; j++)
01700 for ( int k = 0; k < jkDim ; k++)
01701 assert (alpha->getVal ( i ,j,k) == alpha->getVal ( i,k,j ) );
01702
01703
01704
01705 epsFrommer_m->getRingRef().setEpsPrecision(oldEpsPrecision);
01706
01707 return alpha;
01708
01709 }
01710
01711
01712
01713
01714 template <int variant>
01715 template < class Matrix_Type >
01716
01717 inline Matrix3D<Matrix_Type>*
01718 CenterFocusExperiment<variant>::computeQuadric_computeLambda(const Matrix3D<Matrix_Type> & alpha )
01719 {
01720
01721 assert( &alpha!=0);
01722
01723 int jkDim = alpha.getColNum();
01724
01725 int iDim = alpha.getRowNum();
01726
01727 Matrix3D<Matrix_Type>* lambda = new Matrix3D<Matrix_Type>(iDim, jkDim, jkDim, alpha.getRing(), "lambda " ) ;
01728
01729
01730 typename D_Eps_Frommer::RingType::FieldType const * epsField = epsFrommer_m->getRing()->getField();
01731
01732
01733 for ( int i = 0; i < iDim; i++)
01734 {
01735 for ( int jj = 0; jj < jkDim; jj++)
01736 {
01737 lambda->setVal ( i ,jj ,jj,
01738 epsField ->multiply( epsField->multInv(4) ,
01739 alpha.getVal( i ,jj, jj ) )
01740 );
01741 }
01742 }
01743
01744 for ( int i = 0; i < iDim ; i++)
01745
01746 for ( int j = 0; j < jkDim; j++)
01747
01748 for ( int k = 0; k < jkDim ; k++)
01749
01750 if (k!=j )
01751 {
01752 lambda->setVal(i,j,k, epsField->scalarMultiply( epsField->multInv( 2 ) ,
01753 epsField ->add( alpha.getVal( i, j, k ),
01754 epsField->addInv ( epsField ->add( lambda->getVal( i ,k, k ),
01755 lambda->getVal( i ,j, j )
01756 )
01757 )
01758 )
01759 )
01760
01761 );
01762
01763 }
01764
01765
01766 for ( int i = 0; i < iDim ; i++)
01767 for ( int j = 0; j < jkDim; j++)
01768 for ( int k = 0; k < jkDim ; k++)
01769 assert (lambda->getVal ( i ,j,k)== lambda->getVal ( i,k,j ));
01770
01771 return lambda;
01772 }
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01814 template <int variant>
01815 template < class Matrix_Type >
01816 inline Matrix_Type* CenterFocusExperiment<variant>::computeQuadric_getLeftInverse(const Matrix_Type & jacobiKernel )
01817 {
01818 assert ( jacobiKernel.getRing()->getEpsPrecision()==0 );
01819 assert ( jacobiKernel.getRank() ==jacobiKernel.getColNum() );
01820
01821
01822 Matrix_Type * jacobiKernelTransposed = jacobiKernel.getTransposed();
01823 Matrix_Type * jacobiKernelTransposedKernel = jacobiKernelTransposed->getKernel();
01824
01825 Matrix_Type * composed = complementColumns(jacobiKernel);
01826
01827 assert ( composed->getRank() ==composed->getRowNum() );
01828 assert ( composed->getColNum() ==composed->getRowNum() );
01829
01830 FFpackMatrix< Modular<double> > FPcomposed (*composed,jacobiKernel.getRing()->getCharacteristic(),false );
01831
01832 FFpackMatrix< Modular<double> >* FPcomposedLeftInverse = FPcomposed.getLeftInverse();
01833
01834 assert ( FPcomposedLeftInverse!=NULL );
01835
01836 Matrix_Type reconv(*FPcomposedLeftInverse, jacobiKernel.getRing() );
01837
01838
01839 Matrix_Type * D=NULL;
01840 if ( jacobiKernel.getColNum()>0 )
01841 D = reconv.getSubMatrix(0 ,0,jacobiKernel.getColNum()-1,reconv.getColNum()-1);
01842 else
01843 D = new Matrix_Type(0,jacobiKernel.getRowNum(),reconv.getRing() );
01844
01845 Matrix_Type * testInverse = ((*D)*jacobiKernel);
01846
01847
01848
01849
01850 assert( testInverse->isIdentity() );
01851
01852 delete jacobiKernelTransposed; jacobiKernelTransposed=NULL;
01853 delete jacobiKernelTransposedKernel; jacobiKernelTransposedKernel=NULL;
01854 delete composed; composed=NULL;
01855 delete FPcomposedLeftInverse; FPcomposedLeftInverse=NULL;
01856 delete testInverse; testInverse=NULL;
01857
01858
01859 return D;
01860 }
01861
01862
01863
01869 template <int variant>
01870 template < class Matrix_Type >
01871 inline Matrix3D<Matrix_Type>*
01872 CenterFocusExperiment<variant>::computeQuadric_computeBigQuadric( const Matrix3D<Matrix_Type> & smallQuadric,
01873 const Matrix_Type & jacobiKernel )
01874 {
01875
01876 Matrix_Type * De = computeQuadric_getLeftInverse(jacobiKernel);
01877
01878
01879
01880 Matrix_Type * DeTransposed = De->getTransposed();
01881
01882
01883
01884 Matrix3D<Matrix_Type>* righttmp = smallQuadric.rightMultiply(De);
01885
01886
01887
01888 Matrix3D<Matrix_Type>* bigQuadric = righttmp->leftMultiply(DeTransposed);
01889
01890 delete De; De=NULL;
01891 delete DeTransposed; DeTransposed=NULL;
01892 delete righttmp; righttmp=NULL;
01893
01894 return bigQuadric;
01895 }
01896
01897
01898
01899
01902 template <int variant>
01903 template < class Matrix_Type >
01904 inline Matrix3D<Matrix_Type>*
01905 CenterFocusExperiment<variant>::computeQuadric_computeQuadricSmall( const Matrix_Type & jacobiMat,
01906 const Matrix3D<Matrix_Type> & lambda )
01907 {
01908 assert(&lambda!=0);
01909
01910 int jkDim = lambda.getColNum();
01911
01912 Matrix3D<Matrix_Type>* quadricsSmall = NULL ;
01913
01914 Matrix_Type* jacobiTransposedCopy = jacobiMat.getTransposed();
01915
01916
01917
01918 Matrix_Type* jacobiCoKernel = jacobiTransposedCopy->getKernel();
01919 delete jacobiTransposedCopy; jacobiTransposedCopy=NULL;
01920
01921 #ifdef DEBUG
01922 std::cerr << "jacobiMat rows: " << jacobiMat.getRowNum() << std::endl ;
01923 std::cerr << "jacobiMat Rank: " << jacobiMat.getRank() << std::endl ;
01924 std::cerr << "jacobiMat KernelDim: " << jacobiMat.getKernel()->getColNum() << std::endl ;
01925 if (jacobiCoKernel!=0)
01926 std::cerr << "jacobiCoKernel" << std::endl << *jacobiCoKernel << std::endl;
01927 #endif
01928
01929 if (jacobiCoKernel->getColNum()==0)
01930 {
01931 assert( jacobiMat.getRowNum() == jacobiMat.getRank() );
01932 #ifdef DEBUG
01933 std::cerr << "jacobiCoKernel i s Zero! " << std::endl ;
01934 #endif
01935
01936 quadricsSmall = new Matrix3D<Matrix_Type>(jkDim, jkDim, 0, jacobiMat.getRing(), "quadricsSmall" ) ;
01937 delete jacobiCoKernel; jacobiCoKernel=NULL;
01938 return quadricsSmall;
01939 }
01940
01941 assert( jacobiMat.getRowNum() == jacobiMat.getRank()+jacobiCoKernel->getColNum() );
01942
01943 assert(jacobiCoKernel!=0);
01944
01945
01946 Matrix_Type* jacobiCoKernelTransposed = jacobiCoKernel->getTransposed();
01947 delete jacobiCoKernel; jacobiCoKernel=NULL;
01948
01949 Matrix_Type* prod=(*jacobiCoKernelTransposed)*(jacobiMat);
01950 assert(prod->isZero());
01951 delete prod; prod=NULL;
01952
01953 Matrix3D<Matrix_Type>* CMultLambda= lambda.leftMultiply(jacobiCoKernelTransposed);
01954
01955 CMultLambda->setName("JacobiCoKernel * Lambda");
01956
01957 Matrix3D<Matrix_Type>* preQuadricsSmall = CMultLambda->getTransversalForm();
01958
01959 preQuadricsSmall->setName("preQuadricsSmall");
01960
01961
01962 quadricsSmall = preQuadricsSmall->computeFrontalMatrixBasis();
01963
01964 #ifdef DEBUG
01965
01966
01967 preQuadricsSmall->print3DMatrix();
01968 quadricsSmall->print3DMatrix();
01969 if ( quadricsSmall->getZNum() > 0)
01970 {
01971 quadricsSmall->print3DMatrix(std::cerr);
01972 }
01973
01974 #endif
01975
01976 delete jacobiCoKernelTransposed; jacobiCoKernelTransposed=NULL;
01977 delete CMultLambda; CMultLambda=NULL;
01978 delete preQuadricsSmall; preQuadricsSmall=NULL;
01979
01980 return quadricsSmall;
01981 }
01982
01983
01984
01985
01987 template <int variant>
01988 template <class TPolynomXY_Type, class Matrix_Type >
01989 inline void CenterFocusExperiment<variant>::computeQuadric( const TPolynomXY_Type &minusPpol,
01990 const TPolynomXY_Type &Qpol,
01991 const Matrix_Type & jacobiMat,
01992 const list<CoeffListEntry> & coeffOrder,
01993 CFQuadricsResult<Matrix_Type, Matrix3D<Matrix_Type> > & quadricResult
01994 )
01995 {
01996 Matrix_Type* jacobiKernel = jacobiMat.getKernel();
01997 Matrix3D<Matrix_Type>* quadricSmall = NULL;
01998 Matrix3D<Matrix_Type>* quadricBig = NULL;
01999
02000 if ( jacobiKernel->getColNum()==0)
02001 {
02002 cerr << "jacobiKernel->getRowNum()" << jacobiKernel->getRowNum() << endl;
02003 jacobiKernel= new Matrix_Type(jacobiMat.getColNum(),0,jacobiMat.getRing() );
02004 quadricSmall = new Matrix3D<Matrix_Type>(0, 0, 0, jacobiMat.getRing(), "quadricSmall " ) ;
02005 quadricBig = new Matrix3D<Matrix_Type>(0, 0, 0, jacobiMat.getRing(), "quadricBig " ) ;
02006 quadricResult.setQuadrics(quadricBig, quadricSmall);
02007 return;
02008 }
02009
02010
02011 Matrix_Type* prod=(jacobiMat)*(*jacobiKernel);
02012 assert(prod->isZero());
02013 delete prod; prod=NULL;
02014
02015 assert( jacobiMat.getColNum() == jacobiMat.getRank()+jacobiKernel->getColNum() );
02016
02017
02018
02019
02020
02021 Matrix3D<Matrix_Type>* alpha = computeQuadric_computeAlpha(minusPpol, Qpol, jacobiMat, *jacobiKernel, coeffOrder);
02022 Matrix3D<Matrix_Type>* lambda = computeQuadric_computeLambda(*alpha);
02023
02024 quadricSmall = computeQuadric_computeQuadricSmall(jacobiMat, *lambda );
02025
02026 quadricBig = computeQuadric_computeBigQuadric(*quadricSmall,*jacobiKernel);
02027
02028 quadricBig->setName("big quadric");
02029
02030 #ifdef DEBUG
02031 alpha->printTransversalView();
02032 lambda->printTransversalView();
02033
02034
02035 if (quadricBig->getZNum() > 0 )
02036 {
02037 quadricBig->print3DMatrix(std::cerr);
02038
02040
02041 }
02042 #endif
02043
02044 delete alpha; alpha=NULL;
02045 delete lambda; lambda=NULL;
02046 delete jacobiKernel; jacobiKernel=NULL;
02047
02048 quadricResult.setQuadrics(quadricBig, quadricSmall);
02049 }
02050
02051
02052 template <int variant>
02053 template <class TPolynomXY_SRC_Type>
02054 void CenterFocusExperiment<variant>::polynomSetEpsPrecision( TPolynomXY_SRC_Type & srcDestPol, int epsPrecision)
02055 {
02056 for (int deg = srcDestPol.getDegree(); deg>=0; deg--)
02057 {
02058 for (int x_exp=deg; x_exp>=0; x_exp--)
02059 {
02060 srcDestPol.getCoeffRef( x_exp, deg-x_exp).setEpsPrecision(epsPrecision );
02061 }
02062 }
02063 }
02064
02065
02067 template <int variant>
02068 template <class TPolynomXY_SRC_Type, class TPolynomXY_DEST_Type>
02069
02070 void CenterFocusExperiment<variant>::copyPolynomWithGivenEpsPrecision(const TPolynomXY_SRC_Type & srcPol,
02071 TPolynomXY_DEST_Type & destPol, int epsPrecision)
02072 {
02073 for (int deg = srcPol.getDegree(); deg>=0; deg--)
02074 {
02075 for (int x_exp=deg; x_exp>=0; x_exp--)
02076 {
02077 int MaxEpsPrecision= min(epsPrecision, (int)srcPol.getCoeff(x_exp, deg-x_exp).getEpsPrecision() );
02078 for (int prec=0; prec<= MaxEpsPrecision ; prec++)
02079 {
02080 destPol.getCoeffRef( x_exp, deg-x_exp).setValue(prec,
02081 srcPol.getCoeff(x_exp, deg-x_exp).getValue(prec) );
02082 }
02083 }
02084 }
02085 }
02086
02112
02113
02114
02115
02116 template <int variant>
02117 template <class TFrommer, class TPolynomXY_Type >
02118 TMatrix<typename D_Eps_Frommer::RingType::FieldType >* CenterFocusExperiment<variant>::createJacobiMatrix(
02119 TFrommer &frommer2 ,
02120 const TPolynomXY_Type & minusPpol,
02121 const TPolynomXY_Type & Qpol,
02122 int reqVanishedFocalValuesCount,
02123 const list<CoeffListEntry> & coeffVariablesOrder)
02124 {
02125
02126 #ifdef DEBUG
02127 std::cerr << " createJacobiMatrix" << std::endl;
02128 #endif
02129
02131 typename TFrommer::TPolynomXY minusPTruncatedToPrecision0( minusPpol.getDegree() ) ;
02132 typename TFrommer::TPolynomXY qTruncatedToPrecision0 ( Qpol.getDegree() );
02133
02134
02136 copyPolynomWithGivenEpsPrecision(minusPpol, minusPTruncatedToPrecision0, 0);
02137 copyPolynomWithGivenEpsPrecision(Qpol , qTruncatedToPrecision0 , 0);
02138
02139
02140
02141
02142 int oldEpsPrecision = frommer2.getRingRef().getEpsPrecision();
02143
02144
02145 if (oldEpsPrecision!=1)
02146 {
02147 frommer2.getRingRef().setEpsPrecision(1);
02148 polynomSetEpsPrecision(minusPTruncatedToPrecision0 ,1);
02149 polynomSetEpsPrecision(qTruncatedToPrecision0, 1);
02150 }
02151
02152 TMatrix<typename D_Eps_Frommer::RingType::FieldType >* jacobi_matrix1;
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163 jacobi_matrix1= new TMatrix<typename D_Eps_Frommer::RingType::FieldType >( reqVanishedFocalValuesCount,
02164 coeffVariablesOrder.size() ,
02165 epsFrommer_m->getRing()->getField() );
02166
02167
02168 typename TFrommer::RingType::FieldType::ElementType minusOne = frommer2.getRing()->getField()->addInv(1) ;
02170
02171 int colIdx = 0;
02172
02173 for ( list<CoeffListEntry>::const_iterator coeffIt = coeffVariablesOrder.begin();
02174 coeffIt != coeffVariablesOrder.end();
02175 coeffIt++ ,colIdx++ )
02176 {
02177
02178 if ( (*coeffIt).isPMonom() )
02179 {
02180 minusPTruncatedToPrecision0.getCoeffRef( (*coeffIt).x_exp, (*coeffIt).y_exp).setEps( minusOne.getX() );
02181 }else
02182 {
02183 qTruncatedToPrecision0.getCoeffRef( (*coeffIt).x_exp, (*coeffIt).y_exp).setEps(1);
02184 }
02185
02186 frommer2.setPolynoms( &minusPTruncatedToPrecision0, & qTruncatedToPrecision0 );
02187
02188 bool bComputeMaxFocalValues=false;
02189 bool jacobiVariant=true;
02190 frommer2.doit( bComputeMaxFocalValues=false, jacobiVariant=true );
02191
02192 for (short rowIdxPlus1 = 1; rowIdxPlus1 <= reqVanishedFocalValuesCount; rowIdxPlus1++)
02193 {
02194
02195 typename TFrommer::TPolynomXY::CoefficientType tmpFocalValue = frommer2.getFocalValue(rowIdxPlus1);
02196
02197 if (tmpFocalValue.getX() != TFrommer::TPolynomXY::CoefficientType::Zero.getX() )
02198 {
02199 std::cerr << "createJacobiMatrix, error: focal value not nearly zero !!!" << std::endl;
02200 std::cerr << "tmpFocalValue (" << tmpFocalValue.getX() << "," << tmpFocalValue.getEps() << ")" << std::endl;
02201 exit(0);
02202 }
02203 jacobi_matrix1->setVal(rowIdxPlus1 - 1, colIdx , tmpFocalValue.getEps() );
02204 }
02205
02206 if ( (*coeffIt).isPMonom() )
02207 {
02208 minusPTruncatedToPrecision0.getCoeffRef( (*coeffIt).x_exp, (*coeffIt).y_exp).setEps(0);
02209 }
02210 else
02211 {
02212 qTruncatedToPrecision0.getCoeffRef( (*coeffIt).x_exp, (*coeffIt).y_exp).setEps(0);
02213 }
02214
02215 }
02216 if (oldEpsPrecision!=1)
02217 frommer2.getRingRef().setEpsPrecision(oldEpsPrecision);
02218
02219 return jacobi_matrix1;
02220 }
02221
02222
02223 template <int variant>
02224 CenterFocusExperiment<variant>::~CenterFocusExperiment()
02225 {
02226 if (st_m!=NULL)
02227 {
02228 delete st_m;
02229 st_m=NULL;
02230 }
02231
02232 if (minusP_polynom_m!=NULL)
02233 {
02234 delete minusP_polynom_m;
02235 minusP_polynom_m=NULL;
02236 }
02237
02238 if (Q_polynom_m !=NULL)
02239 {
02240 delete Q_polynom_m;
02241 Q_polynom_m = NULL;
02242 }
02243 delete epsRegularFrommer_m;
02244 delete epsFrommer_m;
02245 delete epsRing_m;
02246 delete ring2_m;
02247
02248 #ifdef ALGORITHM_TIMER
02249 MtcpManager_g.disconnectTimer(tim_1_m);
02250 MtcpManager_g.disconnectTimer(tim_2_m);
02251
02252 #endif
02253
02254 }
02255
