File indexing completed on 2024-04-06 12:26:30
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0015
0016 #include "RecoLocalTracker/SiStripRecHitConverter/interface/VVIObj.h"
0017 #else
0018 #include "VVIObj.h"
0019 #endif
0020
0021 #include <cmath>
0022 #include <algorithm>
0023 #include <functional>
0024
0025 namespace sistripvvi {
0026
0027 namespace VVIObjDetails {
0028 void sincosint(double x, double& sint, double& cint);
0029 double cosint(double x);
0030 double sinint(double x);
0031 double expint(double x);
0032
0033 inline double f1(double x, double const* h_) { return h_[0] + h_[1] * std::log(h_[2] * x) - h_[3] * x; }
0034 inline double f2(double x, double const* h_) {
0035 return h_[4] - x + h_[5] * (std::log(std::abs(x)) + expint(x)) - h_[6] * std::exp(-x);
0036 }
0037 template <typename F>
0038 int dzero(double a, double b, double& x0, double& rv, double eps, int mxf, F func);
0039 }
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049 VVIObj::VVIObj(double kappa, double beta2, int mode) : mode_(mode) {
0050 const double xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
0051 const double xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
0052 double h_[7];
0053 double q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
0054 int lp, lq, k, l, n;
0055
0056
0057
0058 if (kappa < 0.01)
0059 kappa = 0.01;
0060 if (kappa > 10.)
0061 kappa = 10.;
0062 if (beta2 < 0.)
0063 beta2 = 0.;
0064 if (beta2 > 1.)
0065 beta2 = 1.;
0066
0067 h_[4] = 1. - beta2 * 0.42278433999999998 + 7.6 / kappa;
0068 h_[5] = beta2;
0069 h_[6] = 1. - beta2;
0070 h4 = -7.6 / kappa - (beta2 * .57721566 + 1);
0071 h5 = log(kappa);
0072 h6 = 1. / kappa;
0073 t0_ = (h4 - h_[4] * h5 - (h_[4] + beta2) * (log(h_[4]) + VVIObjDetails::expint(h_[4])) + exp(-h_[4])) / h_[4];
0074
0075
0076
0077 for (lp = 0; lp < 9; ++lp) {
0078 if (kappa >= xp[lp])
0079 break;
0080 }
0081 ll = -lp - 1.5;
0082 for (lq = 0; lq < 7; ++lq) {
0083 if (kappa <= xq[lq])
0084 break;
0085 }
0086 ul = lq - 6.5;
0087
0088 VVIObjDetails::dzero(ll, ul, u, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f2, std::placeholders::_1, h_));
0089 q = 1. / u;
0090 t1_ = h4 * q - h5 - (beta2 * q + 1) * (log((fabs(u))) + VVIObjDetails::expint(u)) + exp(-u) * q;
0091 t_ = t1_ - t0_;
0092 omega_ = 6.2831853000000004 / t_;
0093 h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
0094 if (kappa >= .07) {
0095 h_[0] += 6.90775527;
0096 }
0097 h_[1] = beta2 * kappa;
0098 h_[2] = h6 * omega_;
0099 h_[3] = omega_ * 1.5707963250000001;
0100
0101 VVIObjDetails::dzero(5., 155., x0_, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f1, std::placeholders::_1, h_));
0102 n = x0_ + 1.;
0103 d = exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
0104 a_[n - 1] = 0.;
0105 if (mode_ == 0) {
0106 a_[n - 1] = omega_ * .31830988654751274;
0107 }
0108 q = -1.;
0109 q2 = 2.;
0110 for (k = 1; k < n; ++k) {
0111 l = n - k;
0112 x = omega_ * k;
0113 x1 = h6 * x;
0114 VVIObjDetails::sincosint(x1, c2, c1);
0115 c1 = log(x) - c1;
0116 c3 = sin(x1);
0117 c4 = cos(x1);
0118 xf1 = kappa * (beta2 * c1 - c4) - x * c2;
0119 xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
0120 if (mode_ == 0) {
0121 d1 = q * d * omega_ * exp(xf1);
0122 a_[l - 1] = d1 * cos(xf2);
0123 b_[l - 1] = -d1 * sin(xf2);
0124 } else {
0125 d1 = q * d * exp(xf1) / k;
0126 a_[l - 1] = d1 * sin(xf2);
0127 b_[l - 1] = d1 * cos(xf2);
0128 a_[n - 1] += q2 * a_[l - 1];
0129 }
0130 q = -q;
0131 q2 = -q2;
0132 }
0133
0134 }
0135
0136
0137
0138
0139
0140
0141
0142 double VVIObj::fcn(double x) const {
0143
0144
0145 double f, u, y, a0, a1;
0146 double a2 = 0.;
0147 double b1, b0, b2, cof;
0148 int k, n, n1;
0149
0150 n = x0_;
0151 if (x < t0_) {
0152 f = 0.;
0153 } else if (x <= t1_) {
0154 y = x - t0_;
0155 u = omega_ * y - 3.141592653589793;
0156 cof = cos(u) * 2.;
0157 a1 = 0.;
0158 a0 = a_[0];
0159 n1 = n + 1;
0160 for (k = 2; k <= n1; ++k) {
0161 a2 = a1;
0162 a1 = a0;
0163 a0 = a_[k - 1] + cof * a1 - a2;
0164 }
0165 b1 = 0.;
0166 b0 = b_[0];
0167 for (k = 2; k <= n; ++k) {
0168 b2 = b1;
0169 b1 = b0;
0170 b0 = b_[k - 1] + cof * b1 - b2;
0171 }
0172 f = (a0 - a2) * .5 + b0 * sin(u);
0173 if (mode_ != 0) {
0174 f += y / t_;
0175 }
0176 } else {
0177 f = 0.;
0178 if (mode_ != 0) {
0179 f = 1.;
0180 }
0181 }
0182 return f;
0183 }
0184
0185
0186
0187
0188
0189
0190
0191 void VVIObj::limits(double& xl, double& xu) const {
0192 xl = t0_;
0193 xu = t1_;
0194 return;
0195 }
0196
0197 namespace VVIObjDetails {
0198 double cosint(double x) {
0199
0200
0201 const double zero = 0.;
0202 const double one = 1.;
0203 const double two = 2.;
0204 const double eight = 8.;
0205 const double ce = .57721566490153;
0206 const double c__[14] = {1.9405491464836,
0207 .9413409132865,
0208 -.579845034293,
0209 .3091572011159,
0210 -.0916101792208,
0211 .0164437407515,
0212 -.0019713091952,
0213 1.692538851e-4,
0214 -1.09393296e-5,
0215 5.522386e-7,
0216 -2.23995e-8,
0217 7.465e-10,
0218 -2.08e-11,
0219 5e-13};
0220 const double p[23] = {
0221 .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
0222 6.0069212e-7, -1.2586794e-7, 2.932563e-8, -7.45696e-9, 2.04105e-9, -5.9502e-10,
0223 1.8323e-10, -5.921e-11, 1.997e-11, -7e-12, 2.54e-12, -9.5e-13,
0224 3.7e-13, -1.4e-13, 6e-14, -2e-14, 1e-14};
0225 const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
0226 -4.2110196e-7, 6.907245e-8, -1.318321e-8, 2.83697e-9, -6.7329e-10,
0227 1.734e-10, -4.787e-11, 1.403e-11, -4.33e-12, 1.4e-12,
0228 -4.7e-13, 1.7e-13, -6e-14, 2e-14, -1e-14};
0229
0230
0231 double d__1;
0232
0233
0234 double h__;
0235 int i__;
0236 double r__, y, b0, b1, b2, pp, qq, alfa;
0237
0238
0239
0240 if (x == zero) {
0241 return zero;
0242 }
0243 if (fabs(x) <= eight) {
0244 y = x / eight;
0245
0246 d__1 = y;
0247 h__ = two * (d__1 * d__1) - one;
0248 alfa = -two * h__;
0249 b1 = zero;
0250 b2 = zero;
0251 for (i__ = 13; i__ >= 0; --i__) {
0252 b0 = c__[i__] - alfa * b1 - b2;
0253 b2 = b1;
0254 b1 = b0;
0255 }
0256 b1 = ce + log((fabs(x))) - b0 + h__ * b2;
0257 } else {
0258 r__ = one / x;
0259 y = eight * r__;
0260
0261 d__1 = y;
0262 h__ = two * (d__1 * d__1) - one;
0263 alfa = -two * h__;
0264 b1 = zero;
0265 b2 = zero;
0266 for (i__ = 22; i__ >= 0; --i__) {
0267 b0 = p[i__] - alfa * b1 - b2;
0268 b2 = b1;
0269 b1 = b0;
0270 }
0271 pp = b0 - h__ * b2;
0272 b1 = zero;
0273 b2 = zero;
0274 for (i__ = 19; i__ >= 0; --i__) {
0275 b0 = q[i__] - alfa * b1 - b2;
0276 b2 = b1;
0277 b1 = b0;
0278 }
0279 qq = b0 - h__ * b2;
0280 b1 = r__ * (qq * sin(x) - r__ * pp * cos(x));
0281 }
0282 return b1;
0283 }
0284
0285 double sinint(double x) {
0286
0287
0288 const double zero = 0.;
0289 const double one = 1.;
0290 const double two = 2.;
0291 const double eight = 8.;
0292 const double pih = 1.5707963267949;
0293 const double s[14] = {1.9522209759531,
0294 -.6884042321257,
0295 .4551855132256,
0296 -.1804571236838,
0297 .0410422133759,
0298 -.0059586169556,
0299 6.001427414e-4,
0300 -4.44708329e-5,
0301 2.5300782e-6,
0302 -1.141308e-7,
0303 4.1858e-9,
0304 -1.273e-10,
0305 3.3e-12,
0306 -1e-13};
0307 const double p[23] = {
0308 .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
0309 6.0069212e-7, -1.2586794e-7, 2.932563e-8, -7.45696e-9, 2.04105e-9, -5.9502e-10,
0310 1.8323e-10, -5.921e-11, 1.997e-11, -7e-12, 2.54e-12, -9.5e-13,
0311 3.7e-13, -1.4e-13, 6e-14, -2e-14, 1e-14};
0312 const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
0313 -4.2110196e-7, 6.907245e-8, -1.318321e-8, 2.83697e-9, -6.7329e-10,
0314 1.734e-10, -4.787e-11, 1.403e-11, -4.33e-12, 1.4e-12,
0315 -4.7e-13, 1.7e-13, -6e-14, 2e-14, -1e-14};
0316
0317
0318 double d__1;
0319
0320
0321 double h__;
0322 int i__;
0323 double r__, y, b0, b1, b2, pp, qq, alfa;
0324
0325 if (fabs(x) <= eight) {
0326 y = x / eight;
0327 d__1 = y;
0328 h__ = two * (d__1 * d__1) - one;
0329 alfa = -two * h__;
0330 b1 = zero;
0331 b2 = zero;
0332 for (i__ = 13; i__ >= 0; --i__) {
0333 b0 = s[i__] - alfa * b1 - b2;
0334 b2 = b1;
0335 b1 = b0;
0336 }
0337 b1 = y * (b0 - b2);
0338 } else {
0339 r__ = one / x;
0340 y = eight * r__;
0341 d__1 = y;
0342 h__ = two * (d__1 * d__1) - one;
0343 alfa = -two * h__;
0344 b1 = zero;
0345 b2 = zero;
0346 for (i__ = 22; i__ >= 0; --i__) {
0347 b0 = p[i__] - alfa * b1 - b2;
0348 b2 = b1;
0349 b1 = b0;
0350 }
0351 pp = b0 - h__ * b2;
0352 b1 = zero;
0353 b2 = zero;
0354 for (i__ = 19; i__ >= 0; --i__) {
0355 b0 = q[i__] - alfa * b1 - b2;
0356 b2 = b1;
0357 b1 = b0;
0358 }
0359 qq = b0 - h__ * b2;
0360 d__1 = fabs(pih);
0361 if (x < 0.)
0362 d__1 = -d__1;
0363 b1 = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
0364 }
0365
0366 return b1;
0367 }
0368
0369 void sincosint(double x, double& sint, double& cint) {
0370
0371
0372 const double zero = 0.;
0373 const double one = 1.;
0374 const double two = 2.;
0375 const double eight = 8.;
0376 const double ce = .57721566490153;
0377 const double pih = 1.5707963267949;
0378 const double s__[14] = {1.9522209759531,
0379 -.6884042321257,
0380 .4551855132256,
0381 -.1804571236838,
0382 .0410422133759,
0383 -.0059586169556,
0384 6.001427414e-4,
0385 -4.44708329e-5,
0386 2.5300782e-6,
0387 -1.141308e-7,
0388 4.1858e-9,
0389 -1.273e-10,
0390 3.3e-12,
0391 -1e-13};
0392
0393 const double c__[14] = {1.9405491464836,
0394 .9413409132865,
0395 -.579845034293,
0396 .3091572011159,
0397 -.0916101792208,
0398 .0164437407515,
0399 -.0019713091952,
0400 1.692538851e-4,
0401 -1.09393296e-5,
0402 5.522386e-7,
0403 -2.23995e-8,
0404 7.465e-10,
0405 -2.08e-11,
0406 5e-13};
0407
0408 const double p[23] = {
0409 .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
0410 6.0069212e-7, -1.2586794e-7, 2.932563e-8, -7.45696e-9, 2.04105e-9, -5.9502e-10,
0411 1.8323e-10, -5.921e-11, 1.997e-11, -7e-12, 2.54e-12, -9.5e-13,
0412 3.7e-13, -1.4e-13, 6e-14, -2e-14, 1e-14};
0413 const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
0414 -4.2110196e-7, 6.907245e-8, -1.318321e-8, 2.83697e-9, -6.7329e-10,
0415 1.734e-10, -4.787e-11, 1.403e-11, -4.33e-12, 1.4e-12,
0416 -4.7e-13, 1.7e-13, -6e-14, 2e-14, -1e-14};
0417
0418
0419 double d__1;
0420
0421
0422 double h__;
0423 int i__;
0424 double r__, y, b0, b1, b2, pp, qq, alfa;
0425
0426 sint = 0;
0427 cint = 0;
0428
0429 if (fabs(x) <= eight) {
0430 y = x / eight;
0431
0432 d__1 = y;
0433 h__ = two * (d__1 * d__1) - one;
0434 alfa = -two * h__;
0435
0436
0437 if (x != 0) {
0438 b1 = zero;
0439 b2 = zero;
0440 for (i__ = 13; i__ >= 0; --i__) {
0441 b0 = c__[i__] - alfa * b1 - b2;
0442 b2 = b1;
0443 b1 = b0;
0444 }
0445 cint = ce + log((fabs(x))) - b0 + h__ * b2;
0446 }
0447
0448 b1 = zero;
0449 b2 = zero;
0450 for (i__ = 13; i__ >= 0; --i__) {
0451 b0 = s__[i__] - alfa * b1 - b2;
0452 b2 = b1;
0453 b1 = b0;
0454 }
0455 sint = y * (b0 - b2);
0456
0457 } else {
0458 r__ = one / x;
0459 y = eight * r__;
0460
0461 d__1 = y;
0462 h__ = two * (d__1 * d__1) - one;
0463 alfa = -two * h__;
0464 b1 = zero;
0465 b2 = zero;
0466 for (i__ = 22; i__ >= 0; --i__) {
0467 b0 = p[i__] - alfa * b1 - b2;
0468 b2 = b1;
0469 b1 = b0;
0470 }
0471 pp = b0 - h__ * b2;
0472 b1 = zero;
0473 b2 = zero;
0474 for (i__ = 19; i__ >= 0; --i__) {
0475 b0 = q[i__] - alfa * b1 - b2;
0476 b2 = b1;
0477 b1 = b0;
0478 }
0479 qq = b0 - h__ * b2;
0480
0481 cint = r__ * (qq * sin(x) - r__ * pp * cos(x));
0482
0483 d__1 = pih;
0484 if (x < 0.)
0485 d__1 = -d__1;
0486 sint = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
0487 }
0488 }
0489
0490 double expint(double x) {
0491
0492
0493 const double zero = 0.;
0494 const double q2[7] = {
0495 .10340013040487, 3.319092135933, 20.449478501379, 41.280784189142, 32.426421069514, 10.041164382905, 1.};
0496 const double p3[6] = {
0497 -2.3909964453136, -147.98219500504, -254.3763397689, -119.55761038372, -19.630408535939, -.9999999999036};
0498 const double q3[6] = {177.60070940351, 530.68509610812, 462.23027156148, 156.81843364539, 21.630408494238, 1.};
0499 const double p4[8] = {-8.6693733995107,
0500 -549.14226552109,
0501 -4210.0161535707,
0502 -249301.39345865,
0503 -119623.66934925,
0504 -22174462.775885,
0505 3892804.213112,
0506 -391546073.8091};
0507 const double q4[8] = {34.171875,
0508 -1607.0892658722,
0509 35730.029805851,
0510 -483547.43616216,
0511 4285596.2461175,
0512 -24903337.574054,
0513 89192576.757561,
0514 -165254299.72521};
0515 const double a1[8] = {-2.1808638152072,
0516 -21.901023385488,
0517 9.3081638566217,
0518 25.076281129356,
0519 -33.184253199722,
0520 60.121799083008,
0521 -43.253113287813,
0522 1.0044310922808};
0523 const double b1[8] = {0.,
0524 3.9370770185272,
0525 300.89264837292,
0526 -6.2504116167188,
0527 1003.6743951673,
0528 14.325673812194,
0529 2736.2411988933,
0530 .52746885196291};
0531 const double a2[8] = {-3.4833465360285,
0532 -18.65454548834,
0533 -8.2856199414064,
0534 -32.34673303054,
0535 17.960168876925,
0536 1.7565631546961,
0537 -1.9502232128966,
0538 .99999429607471};
0539 const double b2[8] = {0.,
0540 69.500065588743,
0541 57.283719383732,
0542 25.777638423844,
0543 760.76114800773,
0544 28.951672792514,
0545 -3.4394226689987,
0546 1.0008386740264};
0547 const double a3[6] = {
0548 -27.780928934438, -10.10479081576, -9.1483008216736, -5.0223317461851, -3.0000077799358, 1.0000000000704};
0549 const double one = 1.;
0550 const double b3[6] = {0., 122.39993926823, 2.7276100778779, -7.1897518395045, -2.9990118065262, 1.999999942826};
0551 const double two = 2.;
0552 const double three = 3.;
0553 const double x0 = .37250741078137;
0554 const double xl[6] = {-24., -12., -6., 0., 1., 4.};
0555 const double p1[5] = {4.293125234321, 39.894153870321, 292.52518866921, 425.69682638592, -434.98143832952};
0556 const double q1[5] = {1., 18.899288395003, 150.95038744251, 568.05252718987, 753.58564359843};
0557 const double p2[7] = {.43096783946939,
0558 6.9052252278444,
0559 23.019255939133,
0560 24.378408879132,
0561 9.0416155694633,
0562 .99997957705159,
0563 4.656271079751e-7};
0564
0565
0566 double v, y, ap, bp, aq, dp, bq, dq;
0567
0568 if (x <= xl[0]) {
0569 ap = a3[0] - x;
0570 for (int i__ = 2; i__ <= 5; ++i__) {
0571
0572 ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
0573 }
0574 y = exp(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
0575 } else if (x <= xl[1]) {
0576 ap = a2[0] - x;
0577 for (int i__ = 2; i__ <= 7; ++i__) {
0578 ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
0579 }
0580 y = exp(-x) / x * (a2[7] + b2[7] / ap);
0581 } else if (x <= xl[2]) {
0582 ap = a1[0] - x;
0583 for (int i__ = 2; i__ <= 7; ++i__) {
0584 ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
0585 }
0586 y = exp(-x) / x * (a1[7] + b1[7] / ap);
0587 } else if (x < xl[3]) {
0588 v = -two * (x / three + one);
0589 bp = zero;
0590 dp = p4[0];
0591 for (int i__ = 2; i__ <= 8; ++i__) {
0592 ap = bp;
0593 bp = dp;
0594 dp = p4[i__ - 1] - ap + v * bp;
0595 }
0596 bq = zero;
0597 dq = q4[0];
0598 for (int i__ = 2; i__ <= 8; ++i__) {
0599 aq = bq;
0600 bq = dq;
0601 dq = q4[i__ - 1] - aq + v * bq;
0602 }
0603 y = -log(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
0604 } else if (x == xl[3]) {
0605 return zero;
0606 } else if (x < xl[4]) {
0607 ap = p1[0];
0608 aq = q1[0];
0609 for (int i__ = 2; i__ <= 5; ++i__) {
0610 ap = p1[i__ - 1] + x * ap;
0611 aq = q1[i__ - 1] + x * aq;
0612 }
0613 y = -log(x) + ap / aq;
0614 } else if (x <= xl[5]) {
0615 y = one / x;
0616 ap = p2[0];
0617 aq = q2[0];
0618 for (int i__ = 2; i__ <= 7; ++i__) {
0619 ap = p2[i__ - 1] + y * ap;
0620 aq = q2[i__ - 1] + y * aq;
0621 }
0622 y = exp(-x) * ap / aq;
0623 } else {
0624 y = one / x;
0625 ap = p3[0];
0626 aq = q3[0];
0627 for (int i__ = 2; i__ <= 6; ++i__) {
0628 ap = p3[i__ - 1] + y * ap;
0629 aq = q3[i__ - 1] + y * aq;
0630 }
0631 y = exp(-x) * y * (one + y * ap / aq);
0632 }
0633 return y;
0634 }
0635
0636 template <typename F>
0637 int dzero(double a, double b, double& x0, double& rv, double eps, int mxf, F func) {
0638
0639 double d__1, d__2, d__3, d__4;
0640
0641
0642 double f1, f2, f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc, fa, fb, ee, ff;
0643 int mc;
0644 double xa, xb, fx, xx, su4;
0645
0646 xa = std::min(a, b);
0647 xb = std::max(a, b);
0648 fa = func(xa);
0649 fb = func(xb);
0650 if (fa * fb > 0.) {
0651 rv = (xb - xa) * -2;
0652 x0 = 0.;
0653 return 1;
0654 }
0655 mc = 0;
0656 L1:
0657 x0 = (xa + xb) * .5;
0658 rv = x0 - xa;
0659 ee = eps * (fabs(x0) + 1);
0660 if (rv <= ee) {
0661 rv = ee;
0662 ff = func(x0);
0663 return 0;
0664 }
0665 f1 = fa;
0666 x1 = xa;
0667 f2 = fb;
0668 x2 = xb;
0669 L2:
0670 fx = func(x0);
0671 ++mc;
0672 if (mc > mxf) {
0673 rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
0674 x0 = 0.;
0675 return 0;
0676 }
0677 if (fx * fa > 0.) {
0678 xa = x0;
0679 fa = fx;
0680 } else {
0681 xb = x0;
0682 fb = fx;
0683 }
0684 L3:
0685 u1 = f1 - f2;
0686 u2 = x1 - x2;
0687 u3 = f2 - fx;
0688 u4 = x2 - x0;
0689 if (u2 == 0. || u4 == 0.) {
0690 goto L1;
0691 }
0692 f3 = fx;
0693 x3 = x0;
0694 u1 /= u2;
0695 u2 = u3 / u4;
0696 ca = u1 - u2;
0697 cb = (x1 + x2) * u2 - (x2 + x0) * u1;
0698 cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
0699 if (ca == 0.) {
0700 if (cb == 0.) {
0701 goto L1;
0702 }
0703 x0 = -cc / cb;
0704 } else {
0705 u3 = cb / (ca * 2);
0706 u4 = u3 * u3 - cc / ca;
0707 if (u4 < 0.) {
0708 goto L1;
0709 }
0710 su4 = fabs(u4);
0711 if (x0 + u3 < 0.f) {
0712 su4 = -su4;
0713 }
0714 x0 = -u3 + su4;
0715 }
0716 if (x0 < xa || x0 > xb) {
0717 goto L1;
0718 }
0719
0720 d__3 = (d__1 = x0 - x3, fabs(d__1)), d__4 = (d__2 = x0 - x2, fabs(d__2));
0721 rv = std::min(d__3, d__4);
0722 ee = eps * (fabs(x0) + 1);
0723 if (rv > ee) {
0724 f1 = f2;
0725 x1 = x2;
0726 f2 = f3;
0727 x2 = x3;
0728 goto L2;
0729 }
0730 fx = func(x0);
0731 if (fx == 0.) {
0732 rv = ee;
0733 ff = func(x0);
0734 return 0;
0735 }
0736 if (fx * fa < 0.) {
0737 xx = x0 - ee;
0738 if (xx <= xa) {
0739 rv = ee;
0740 ff = func(x0);
0741 return 0;
0742 }
0743 ff = func(xx);
0744 fb = ff;
0745 xb = xx;
0746 } else {
0747 xx = x0 + ee;
0748 if (xx >= xb) {
0749 rv = ee;
0750 ff = func(x0);
0751 return 0;
0752 }
0753 ff = func(xx);
0754 fa = ff;
0755 xa = xx;
0756 }
0757 if (fx * ff > 0.) {
0758 mc += 2;
0759 if (mc > mxf) {
0760 rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
0761 x0 = 0.;
0762 return 0;
0763 }
0764 f1 = f3;
0765 x1 = x3;
0766 f2 = fx;
0767 x2 = x0;
0768 x0 = xx;
0769 fx = ff;
0770 goto L3;
0771 }
0772
0773 rv = ee;
0774 ff = func(x0);
0775 return 0;
0776 }
0777
0778 }
0779 }