00001
00002
00003
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 #define POLEPSILON 1e-12
00035 #define POLINFINITY 9.0E99
00036 #define ITERATION 12 // 20 But 8 is really enough for Newton with a good guess
00037
00038 template <class T, class F>
00039 G4PolynomialSolver<T,F>::G4PolynomialSolver (T* typeF, F func, F deriv,
00040 G4double precision)
00041 {
00042 Precision = precision ;
00043 FunctionClass = typeF ;
00044 Function = func ;
00045 Derivative = deriv ;
00046 }
00047
00048 template <class T, class F>
00049 G4PolynomialSolver<T,F>::~G4PolynomialSolver ()
00050 {
00051 }
00052
00053 template <class T, class F>
00054 G4double G4PolynomialSolver<T,F>::solve(G4double IntervalMin,
00055 G4double IntervalMax)
00056 {
00057 return Newton(IntervalMin,IntervalMax);
00058 }
00059
00060
00061
00062
00063
00064
00065 #define NBBEZIER 5
00066
00067 template <class T, class F>
00068 G4int
00069 G4PolynomialSolver<T,F>::BezierClipping(
00070 G4double *IntervalMin,
00071 G4double *IntervalMax)
00072 {
00076 G4double P[NBBEZIER][2],D[2];
00077 G4double NewMin,NewMax;
00078
00079 G4int IntervalIsVoid = 1;
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 D[0] = (FunctionClass->*Derivative)(*IntervalMin);
00100
00101 P[0][0] = (*IntervalMin);
00102 P[0][1] = (FunctionClass->*Function)(*IntervalMin);
00103
00104
00105 if (std::fabs(P[0][1]) < Precision) {
00106 return 1;
00107 }
00108
00109 if (((*IntervalMax) - (*IntervalMin)) < POLEPSILON) {
00110 return 1;
00111 }
00112
00113 P[1][0] = (*IntervalMin) + ((*IntervalMax) - (*IntervalMin))/4;
00114 P[1][1] = P[0][1] + (((*IntervalMax) - (*IntervalMin))/4.0) * D[0];
00115
00116 D[1] = (FunctionClass->*Derivative)(*IntervalMax);
00117
00118 P[4][0] = (*IntervalMax);
00119 P[4][1] = (FunctionClass->*Function)(*IntervalMax);
00120
00121 P[3][0] = (*IntervalMax) - ((*IntervalMax) - (*IntervalMin))/4;
00122 P[3][1] = P[4][1] - ((*IntervalMax) - (*IntervalMin))/4 * D[1];
00123
00124 P[2][0] = ((*IntervalMax) + (*IntervalMin))/2;
00125 P[2][1] = (16*(FunctionClass->*Function)(((*IntervalMax)+(*IntervalMin))/2)
00126 - (P[0][1] + 4*P[1][1] + 4*P[3][1] + P[4][1]))/6 ;
00127
00128 {
00129 G4double Intersection ;
00130 G4int i,j;
00131
00132 NewMin = (*IntervalMax) ;
00133 NewMax = (*IntervalMin) ;
00134
00135 for (i=0;i<5;i++)
00136 for (j=i+1;j<5;j++)
00137 {
00138
00139 if (((P[j][1] > -Precision) && (P[i][1] < Precision)) ||
00140 ((P[j][1] < Precision) && (P[i][1] > -Precision))) {
00141 IntervalIsVoid = 0;
00142 Intersection = P[j][0] - P[j][1]*((P[i][0] - P[j][0])/
00143 (P[i][1] - P[j][1]));
00144 if (Intersection < NewMin) {
00145 NewMin = Intersection;
00146 }
00147 if (Intersection > NewMax) {
00148 NewMax = Intersection;
00149 }
00150 }
00151 }
00152
00153
00154 if (IntervalIsVoid != 1) {
00155 (*IntervalMax) = NewMax;
00156 (*IntervalMin) = NewMin;
00157 }
00158 }
00159
00160 if (IntervalIsVoid == 1) {
00161 return -1;
00162 }
00163
00164 return 0;
00165 }
00166
00167 template <class T, class F>
00168 G4double G4PolynomialSolver<T,F>::Newton (G4double IntervalMin,
00169 G4double IntervalMax)
00170 {
00171
00172
00173
00174 G4double Value = 0;
00175 G4double Gradient = 0;
00176 G4double Lambda ;
00177
00178 G4int i=0;
00179 G4int j=0;
00180
00181
00182
00183 {
00184 G4int NewtonIsSafe ;
00185
00186 while ((NewtonIsSafe = BezierClipping(&IntervalMin,&IntervalMax)) == 0) ;
00187
00188 if (NewtonIsSafe == -1) {
00189 return POLINFINITY;
00190 }
00191 }
00192
00193 Lambda = IntervalMin;
00194 Value = (FunctionClass->*Function)(Lambda);
00195
00196
00197
00198 while (j != -1) {
00199
00200 Value = (FunctionClass->*Function)(Lambda);
00201
00202 Gradient = (FunctionClass->*Derivative)(Lambda);
00203
00204 Lambda = Lambda - Value/Gradient ;
00205
00206 if (std::fabs(Value) <= Precision) {
00207 j ++;
00208 if (j == 2) {
00209 j = -1;
00210 }
00211 } else {
00212 i ++;
00213
00214 if (i > ITERATION)
00215 return POLINFINITY;
00216 }
00217 }
00218 return Lambda ;
00219 }