projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023952: Improving thread-safety of intersections, approximations and other modeling...
[occt.git] / src / math / math_FunctionRoots.cxx
CommitLineData
b311480e 1// Copyright (c) 1997-1999 Matra Datavision
2// Copyright (c) 1999-2012 OPEN CASCADE SAS
3//
4// The content of this file is subject to the Open CASCADE Technology Public
5// License Version 6.5 (the "License"). You may not use the content of this file
6// except in compliance with the License. Please obtain a copy of the License
7// at http://www.opencascade.org and read it completely before using this file.
8//
9// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11//
12// The Original Code and all software distributed under the License is
13// distributed on an "AS IS" basis, without warranty of any kind, and the
14// Initial Developer hereby disclaims all such warranties, including without
15// limitation, any warranties of merchantability, fitness for a particular
16// purpose or non-infringement. Please see the License for the specific terms
17// and conditions governing the rights and limitations under the License.
18
7fd59977 19//#ifndef DEB
20#define No_Standard_RangeError
21#define No_Standard_OutOfRange
22#define No_Standard_DimensionError
23//#endif
24
25#include <StdFail_NotDone.hxx>
26#include <Standard_RangeError.hxx>
27#include <math_DirectPolynomialRoots.hxx>
28#include <math_FunctionRoots.ixx>
29#include <math_FunctionWithDerivative.hxx>
232a3086 30#include <TColStd_Array1OfReal.hxx>
7fd59977 31//#ifdef WNT
32#include <stdio.h>
33//#endif
34
35#define ITMAX 100
36#define EPS 1e-14
37#define EPSEPS 2e-14
38#define MAXBIS 100
39
40# ifdef DEB
41static Standard_Boolean myDebug = 0;
42static Standard_Integer nbsolve = 0;
43# endif
44
45static void AppendRoot(TColStd_SequenceOfReal& Sol,
46 TColStd_SequenceOfInteger& NbStateSol,
47 const Standard_Real X,
48 math_FunctionWithDerivative& F,
49// const Standard_Real K,
50 const Standard_Real ,
51 const Standard_Real dX) {
52
53 Standard_Integer n=Sol.Length();
54 Standard_Real t;
55#ifdef DEB
56 if (myDebug) {
57 cout << " Ajout de la solution numero : " << n+1 << endl;
58 cout << " Valeur de la racine :" << X << endl;
59 }
60#endif
61
62 if(n==0) {
63 Sol.Append(X);
64 F.Value(X,t);
65 NbStateSol.Append(F.GetStateNumber());
66 }
67 else {
68 Standard_Integer i=1;
69 Standard_Integer pl=n+1;
70 while(i<=n) {
71 t=Sol.Value(i);
72 if(t >= X) {
73 pl=i;
74 i=n;
75 }
76 if(Abs(X-t) <= dX) {
77 pl=0;
78 i=n;
79 }
80 i++;
81 } //-- while
82 if(pl>n) {
83 Sol.Append(X);
84 F.Value(X,t);
85 NbStateSol.Append(F.GetStateNumber());
86 }
87 else if(pl>0) {
88 Sol.InsertBefore(pl,X);
89 F.Value(X,t);
90 NbStateSol.InsertBefore(pl,F.GetStateNumber());
91 }
92 }
93}
94
95static void Solve(math_FunctionWithDerivative& F,
96 const Standard_Real K,
97 const Standard_Real x1,
98 const Standard_Real y1,
99 const Standard_Real x2,
100 const Standard_Real y2,
101 const Standard_Real tol,
102 const Standard_Real dX,
103 TColStd_SequenceOfReal& Sol,
104 TColStd_SequenceOfInteger& NbStateSol) {
105#ifdef DEB
106 if (myDebug) {
107 cout <<"--> Resolution :" << ++nbsolve << endl;
108 cout <<" x1 =" << x1 << " y1 =" << y1 << endl;
109 cout <<" x2 =" << x2 << " y2 =" << y2 << endl;
110 }
111#endif
112
113 Standard_Integer iter=0;
114 Standard_Real tols2 = 0.5*tol;
115 Standard_Real a,b,c,d=0,e=0,fa,fb,fc,p,q,r,s,tol1,xm,min1,min2;
116 a=x1;b=c=x2;fa=y1; fb=fc=y2;
117 for(iter=1;iter<=ITMAX;iter++) {
118 if((fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0)) {
119 c=a; fc=fa; e=d=b-a;
120 }
121 if(Abs(fc)<Abs(fb)) {
122 a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
123 }
124 tol1 = EPSEPS*Abs(b)+tols2;
125 xm=0.5*(c-b);
126 if(Abs(xm)<tol1 || fb==0) {
127 //-- On tente une iteration de newton
128 Standard_Real Xp,Yp,Dp;
129 Standard_Integer itern=5;
130 Standard_Boolean Ok;
131 Xp=b;
132 do {
133 Ok = F.Values(Xp,Yp,Dp);
134 if(Ok) {
135 Ok=Standard_False;
136 if(Dp>1e-10 || Dp<-1e-10) {
137 Xp = Xp-(Yp-K)/Dp;
138 }
139 if(Xp<=x2 && Xp>=x1) {
140 F.Value(Xp,Yp); Yp-=K;
141 if(Abs(Yp)<Abs(fb)) {
142 b=Xp;
143 fb=Yp;
144 Ok=Standard_True;
145 }
146 }
147 }
148 }
149 while(Ok && --itern>=0);
150
151 AppendRoot(Sol,NbStateSol,b,F,K,dX);
152 return;
153 }
154 if(Abs(e)>=tol1 && Abs(fa)>Abs(fb)) {
155 s=fb/fa;
156 if(a==c) {
157 p=xm*s; p+=p;
158 q=1.0-s;
159 }
160 else {
161 q=fa/fc;
162 r=fb/fc;
163 p=s*((xm+xm)*q*(q-r)-(b-a)*(r-1.0));
164 q=(q-1.0)*(r-1.0)*(s-1.0);
165 }
166 if(p>0.0) {
167 q=-q;
168 }
169 p=Abs(p);
170 min1=3.0*xm*q-Abs(tol1*q);
171 min2=Abs(e*q);
172 if((p+p)<( (min1<min2) ? min1 : min2)) {
173 e=d;
174 d=p/q;
175 }
176 else {
177 d=xm;
178 e=d;
179 }
180 }
181 else {
182 d=xm;
183 e=d;
184 }
185 a=b;
186 fa=fb;
187 if(Abs(d)>tol1) {
188 b+=d;
189 }
190 else {
191 if(xm>=0) b+=Abs(tol1);
192 else b+=-Abs(tol1);
193 }
194 F.Value(b,fb);
195 fb-=K;
196 }
197#ifdef DEB
198 cout<<" Non Convergence dans math_FunctionRoots.cxx "<<endl;
199#endif
200}
201
202#define NEWSEQ 1
203
1ef32e96 204static const Standard_Integer methode = 1; //-- 1:(Nv Traitement) 3:(Nv + Ancien +check) 2:(Ancien)
7fd59977 205
206math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
207 const Standard_Real A,
208 const Standard_Real B,
209 const Standard_Integer NbSample,
210 const Standard_Real _EpsX,
211 const Standard_Real EpsF,
212 const Standard_Real EpsNull,
213 const Standard_Real K )
214{
215#ifdef DEB
216 if (myDebug) {
217 cout << "---- Debut de math_FunctionRoots ----" << endl;
218 nbsolve = 0;
219 }
220#endif
221
1ef32e96
RL		 222#if NEWSEQ
223 TColStd_SequenceOfReal StaticSol;
224#endif
7fd59977 225 Sol.Clear();
226 NbStateSol.Clear();
227 if(methode & 1) {
228 Done = Standard_True;
229 Standard_Real X0=A;
230 Standard_Real XN=B;
231 Standard_Integer N=NbSample;
232 //-- ------------------------------------------------------------
233 //-- Verifications de bas niveau
234 if(B<A) {
235 X0=B;
236 XN=A;
237 }
238 N*=2;
239 if(N < 20) {
240 N=20;
241 }
242 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
243 Standard_Real EpsX = _EpsX;
244 Standard_Real DeltaU = Abs(X0)+Abs(XN);
245 Standard_Real NEpsX = 0.0000000001 * DeltaU;
246 if(EpsX < NEpsX) {
247 EpsX = NEpsX;
248 }
249
250 //-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
251 //-- A .............................................................. B
252 //-- X0 X1 X2 ........................................ Xn-1 Xn
253 Standard_Integer i;
254 Standard_Real X=X0;
255 Standard_Boolean Ok;
256 double dx = (XN-X0)/N;
232a3086 257 TColStd_Array1OfReal ptrval(0, N);
7fd59977 258 Standard_Integer Nvalid = -1;
259 Standard_Real aux = 0;
260 for(i=0; i<=N ; i++,X+=dx) {
261 if( X > XN) X=XN;
262 Ok=F.Value(X,aux);
232a3086
PD		 263 if(Ok) ptrval(++Nvalid) = aux - K;
264// ptrval(i)-=K;
7fd59977 265 }
266 //-- Toute la fonction est nulle ?
267
268 if( Nvalid < N ) {
269 Done = Standard_False;
270 return;
271 }
272
273 AllNull=Standard_True;
274// for(i=0;AllNull && i<=N;i++) {
275 for(i=0;AllNull && i<=N;i++) {
232a3086 276 if(ptrval(i)>EpsNull || ptrval(i)<-EpsNull) {
7fd59977 277 AllNull=Standard_False;
278 }
279 }
280 if(AllNull) {
281 //-- tous les points echantillons sont dans la tolerance
282
283 }
284 else {
285 //-- Il y a des points hors tolerance
286 //-- on detecte les changements de signes STRICTS
287 Standard_Integer ip1;
288// Standard_Boolean chgtsign=Standard_False;
289 Standard_Real tol = EpsX;
290 Standard_Real X2;
291 for(i=0,ip1=1,X=X0;i<N; i++,ip1++,X+=dx) {
292 X2=X+dx;
293 if(X2 > XN) X2 = XN;
232a3086
PD		 294 if(ptrval(i)<0.0) {
295 if(ptrval(ip1)>0.0) {
7fd59977 296 //-- --------------------------------------------------
297 //-- changement de signe dans Xi Xi+1
232a3086 298 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
7fd59977 299 }
300 }
301 else {
232a3086 302 if(ptrval(ip1)<0.0) {
7fd59977 303 //-- --------------------------------------------------
304 //-- changement de signe dans Xi Xi+1
232a3086 305 Solve(F,K,X,ptrval(i),X2,ptrval(ip1),tol,NEpsX,Sol,NbStateSol);
7fd59977 306 }
307 }
308 }
309 //-- On detecte les cas ou la fct s annule sur des Xi et est
310 //-- non nulle au voisinage de Xi
311 //--
312 //-- On prend 2 points u0,u1 au voisinage de Xi
313 //-- Si (F(u0)-K)*(F(u1)-K) <0 on lance une recherche
314 //-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
315 for(i=0; i<=N; i++) {
232a3086 316 if(ptrval(i)==0) {
7fd59977 317// Standard_Real Val,Deriv;
318 X=X0+i*dx;
319 if (X>XN) X=XN;
320 Standard_Real u0,u1;
321 u0=dx*0.5; u1=X+u0; u0+=X;
322 if(u0<X0) u0=X0;
323 if(u0>XN) u0=XN;
324 if(u1<X0) u1=X0;
325 if(u1>XN) u1=XN;
326
327 Standard_Real y0,y1;
328 F.Value(u0,y0); y0-=K;
329 F.Value(u1,y1); y1-=K;
330 if(y0*y1 < 0.0) {
331 Solve(F,K,u0,y0,u1,y1,tol,NEpsX,Sol,NbStateSol);
332 }
333 else {
334 if(y0!=0.0 || y1!=0.0) {
335 AppendRoot(Sol,NbStateSol,X,F,K,NEpsX);
336 }
337 }
338 }
339 }
340 //-- --------------------------------------------------------------------------------
341 //-- Il faut traiter differement le cas des points en bout :
232a3086 342 if(ptrval(0)<=EpsF && ptrval(0)>=-EpsF) {
7fd59977 343 AppendRoot(Sol,NbStateSol,X0,F,K,NEpsX);
344 }
232a3086 345 if(ptrval(N)<=EpsF && ptrval(N)>=-EpsF) {
7fd59977 346 AppendRoot(Sol,NbStateSol,XN,F,K,NEpsX);
347 }
348
349 //-- --------------------------------------------------------------------------------
350 //-- --------------------------------------------------------------------------------
351 //-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
352 //-- un maximum avec f(x)<0
353 //-- On reprend une discretisation plus fine au voisinage de ces extremums
354 //--
355 //-- Recherche d un minima positif
356 Standard_Real xm,ym,dym,xm1,xp1;
357 Standard_Real majdx = 5.0*dx;
358 Standard_Boolean Rediscr;
359// Standard_Real ptrvalbis[MAXBIS];
360 Standard_Integer im1=0;
361 ip1=2;
362 for(i=1,xm=X0+dx; i<N; xm+=dx,i++,im1++,ip1++) {
363 Rediscr = Standard_False;
364 if (xm > XN) xm=XN;
232a3086
PD		 365 if(ptrval(i)>0.0) {
366 if((ptrval(im1)>ptrval(i)) && (ptrval(ip1)>ptrval(i))) {
7fd59977 367 //-- Peut on traverser l axe Ox
368 //-- -------------- Estimation a partir de Xim1
369 xm1=xm-dx;
370 if(xm1 < X0) xm1=X0;
371 F.Values(xm1,ym,dym); ym-=K;
372 if(dym<-1e-10 || dym>1e-10) { // normalement dym < 0
373 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
374 if(t<majdx && t > -majdx) {
375 Rediscr = Standard_True;
376 }
377 }
378 //-- -------------- Estimation a partir de Xip1
379 if(Rediscr==Standard_False) {
380 xp1=xm+dx;
381 if(xp1 > XN ) xp1=XN;
382 F.Values(xp1,ym,dym); ym-=K;
383 if(dym<-1e-10 || dym>1e-10) { // normalement dym > 0
384 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
385 if(t<majdx && t > -majdx) {
386 Rediscr = Standard_True;
387 }
388 }
389 }
390 }
391 }
232a3086
PD		 392 else if(ptrval(i)<0.0) {
393 if((ptrval(im1)<ptrval(i)) && (ptrval(ip1)<ptrval(i))) {
7fd59977 394 //-- Peut on traverser l axe Ox
395 //-- -------------- Estimation a partir de Xim1
396 xm1=xm-dx;
397 if(xm1 < X0) xm1=X0;
398 F.Values(xm1,ym,dym); ym-=K;
399 if(dym>1e-10 || dym<-1e-10) { // normalement dym > 0
400 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
401 if(t<majdx && t > -majdx) {
402 Rediscr = Standard_True;
403 }
404 }
405 //-- -------------- Estimation a partir de Xim1
406 if(Rediscr==Standard_False) {
407 xm1=xm-dx;
408 if(xm1 < X0) xm1=X0;
409 F.Values(xm1,ym,dym); ym-=K;
410 if(dym>1e-10 || dym<-1e-10) { // normalement dym < 0
411 Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
412 if(t<majdx && t > -majdx) {
413 Rediscr = Standard_True;
414 }
415 }
416 }
417 }
418 }
419 if(Rediscr) {
420 //-- --------------------------------------------------
421 //-- On recherche un extrema entre x0 et x3
422 //-- x1 et x2 sont tels que x0<x1<x2<x3
423 //-- et |f(x0)| > |f(x1)| et |f(x3)| > |f(x2)|
424 //--
425 //-- En entree : a=xm-dx b=xm c=xm+dx
426 Standard_Real x0,x1,x2,x3,f0,f3;
427 Standard_Real R=0.61803399;
428 Standard_Real C=1.0-R;
429 Standard_Real tolCR=NEpsX*10.0;
232a3086
PD		 430 f0=ptrval(im1);
431 f3=ptrval(ip1);
7fd59977 432 x0=xm-dx;
433 x3=xm+dx;
434 if(x0 < X0) x0=X0;
435 if(x3 > XN) x3=XN;
436 Standard_Boolean recherche_minimum = (f0>0.0);
437
438 if(Abs(x3-xm) > Abs(x0-xm)) { x1=xm; x2=xm+C*(x3-xm); }
439 else { x2=xm; x1=xm-C*(xm-x0); }
440 Standard_Real f1,f2;
441 F.Value(x1,f1); f1-=K;
442 F.Value(x2,f2); f2-=K;
443 //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
444 while(Abs(x3-x0) > tolCR*(Abs(x1)+Abs(x2)) && (Abs(x1 -x2) > 0)) {
445 //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
446 //-- x0,f0,x1,f1,x2,f2,x3,f3);
447 if(recherche_minimum) {
448 if(f2<f1) {
449 x0=x1; x1=x2; x2=R*x1+C*x3;
450 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
451 }
452 else {
453 x3=x2; x2=x1; x1=R*x2+C*x0;
454 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
455 }
456 }
457 else {
458 if(f2>f1) {
459 x0=x1; x1=x2; x2=R*x1+C*x3;
460 f0=f1; f1=f2; F.Value(x2,f2); f2-=K;
461 }
462 else {
463 x3=x2; x2=x1; x1=R*x2+C*x0;
464 f3=f2; f2=f1; F.Value(x1,f1); f1-=K;
465 }
466 }
467 //-- On ne fait pas que chercher des extremas. Il faut verifier
468 //-- si on ne tombe pas sur une racine
469 if(f1*f0 <0.0) {
470 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
471 Solve(F,K,x0,f0,x1,f1,tol,NEpsX,Sol,NbStateSol);
472 }
473 if(f2*f3 <0.0) {
474 //-- printf("\n Recherche entre (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
475 Solve(F,K,x2,f2,x3,f3,tol,NEpsX,Sol,NbStateSol);
476 }
477 }
478 if(f1<f2) {
479 //-- x1,f(x1) minimum
480 if(Abs(f1) < EpsF) {
481 AppendRoot(Sol,NbStateSol,x1,F,K,NEpsX);
482 }
483 }
484 else {
485 //-- x2.f(x2) minimum
486 if(Abs(f2) < EpsF) {
487 AppendRoot(Sol,NbStateSol,x2,F,K,NEpsX);
488 }
489 }
490 } //-- Recherche d un extrema
491 } //-- for
492 }
493
7fd59977 494#if NEWSEQ
495 if(methode==3) {
496 StaticSol.Clear();
497 Standard_Integer n=Sol.Length();
498 for(Standard_Integer ii=1;ii<=n;ii++) {
499 StaticSol.Append(Sol.Value(ii));
500 }
501 Sol.Clear();
502 NbStateSol.Clear();
503 }
504#endif
505 }
506 if(methode & 2) {
507 //-- ********************************************************************************
508 //-- ANCIEN TRAITEMENT
509 //-- ********************************************************************************
510
511
512 // calculate all the real roots of a function within the range
513 // A..B. whitout condition on A and B
514 // a solution X is found when
515 // abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
516 // The function is considered as null between A and B if
517 // abs(F-K) <= EpsNull within this range.
518 Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
519 //-- Il ne faut pas EpsX = 0.000...001 car dans ce cas
520 //-- U + Nn*EpsX == U
521 Standard_Real Lowr,Upp;
522 Standard_Real Increment;
523 Standard_Real Null2;
524 Standard_Real FLowr,FUpp,DFLowr,DFUpp;
525 Standard_Real U,Xu;
526 Standard_Real Fxu,DFxu,FFxu,DFFxu;
527 Standard_Real Fyu,DFyu,FFyu,DFFyu;
528 Standard_Boolean Finish;
529 Standard_Real FFi;
530 Standard_Integer Nbiter = 30;
531 Standard_Integer Iter;
532 Standard_Real Ambda,T;
533 Standard_Real AA,BB,CC;
534 Standard_Integer Nn;
535 Standard_Real Alfa1=0,Alfa2;
536 Standard_Real OldDF = RealLast();
537 Standard_Real Standard_Underflow = 1e-32 ; //-- RealSmall();
538 Standard_Boolean Ok;
539
540 Done = Standard_False;
541
542 StdFail_NotDone_Raise_if(NbSample <= 0, " ");
543
544 // initialisation
545
546 if (A > B) {
547 Lowr=B;
548 Upp=A;
549 }
550 else {
551 Lowr=A;
552 Upp=B;
553 }
554
555 Increment = (Upp-Lowr)/NbSample;
556 StdFail_NotDone_Raise_if(Increment < EpsX, " ");
557 Done = Standard_True;
558 //-- On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
559 Standard_Real DeltaU = Abs(Upp)+Abs(Lowr);
560 Standard_Real NEpsX = 0.0000000001 * DeltaU;
561 if(EpsX < NEpsX) {
562 EpsX = NEpsX;
563 //-- cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" ) EpsX = "<<EpsX<<endl;
564 }
565 //--
566 Null2 = EpsNull*EpsNull;
567
568 Ok = F.Values(Lowr,FLowr,DFLowr);
569
570 if(!Ok) {
571 Done = Standard_False;
572 return;
573 }
574
575 FLowr = FLowr - K;
576
577 Ok = F.Values(Upp,FUpp,DFUpp);
578
579 if(!Ok) {
580 Done = Standard_False;
581 return;
582 }
583
584 FUpp = FUpp - K;
585
586 // Calcul sur U
587
588 U = Lowr-EpsX;
589 Fyu = FLowr-EpsX*DFLowr; // extrapolation lineaire
590 DFyu = DFLowr;
591 FFyu = Fyu*Fyu;
592 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;
593 AllNull = ( FFyu <= Null2 );
594
595 while ( U < Upp) {
596
597 Xu = U;
598 Fxu = Fyu;
599 DFxu = DFyu;
600 FFxu = FFyu;
601 DFFxu = DFFyu;
602
603 U = Xu + Increment;
604 if (U <= Lowr) {
605 Fyu = FLowr + (U-Lowr)*DFLowr;
606 DFyu = DFLowr;
607 }
608 else if (U >= Upp) {
609 Fyu = FUpp + (U-Upp)*DFUpp;
610 DFyu = DFUpp;
611 }
612 else {
613 Ok = F.Values(U,Fyu,DFyu);
614
615 if(!Ok) {
616 Done = Standard_False;
617 return;
618 }
619
620 Fyu = Fyu - K;
621 }
622 FFyu = Fyu*Fyu;
623 DFFyu = Fyu*DFyu; DFFyu+=DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
624
625 if ( !AllNull || ( FFyu > Null2 && U <= Upp) ){
626
627 if (AllNull) { //rechercher les vraix zeros depuis le debut
628
629 AllNull = Standard_False;
630 Xu = Lowr-EpsX;
631 Fxu = FLowr-EpsX*DFLowr;
632 DFxu = DFLowr;
633 FFxu = Fxu*Fxu;
634 DFFxu = Fxu*DFxu; DFFxu+=DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
635 U = Xu + Increment;
636 Ok = F.Values(U,Fyu,DFyu);
637
638 if(!Ok) {
639 Done = Standard_False;
640 return;
641 }
642
643 Fyu = Fyu - K;
644 FFyu = Fyu*Fyu;
645 DFFyu = Fyu*DFyu; DFFyu+=DFFyu;//-- DFFyu = 2.*Fyu*DFyu;
646 }
647 Standard_Real FxuFyu=Fxu*Fyu;
648
649 if ( (DFFyu > 0. && DFFxu <= 0.)
650 || (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
651 || (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.)
652 || (FxuFyu <= 0.) ) {
653 // recherche d 1 minimun possible
654 Finish = Standard_False;
655 Ambda = Increment;
656 T = 0.;
657 Iter=0;
658 FFi=0.;
659
660 if (FxuFyu >0.) {
661 // chercher si f peut s annuler pour eviter
662 // des iterations inutiles
663 if ( Fxu*(Fxu + 2.*DFxu*Increment) > 0. &&
664 Fyu*(Fyu - 2.*DFyu*Increment) > 0. ) {
665
666 Finish = Standard_True;
667 FFi = Min ( FFxu , FFyu); //pour ne pas recalculer yu
668 }
669 else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow) ||
670 (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)) {
671
672 Finish = Standard_True;
673 FFxu = 0.0;
674 FFi = FFyu; // pour recalculer yu
675 }
676 else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow) ||
677 (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)) {
678
679 Finish = Standard_True;
680 FFyu =0.0;
681 FFi = FFxu; // pour recalculer U
682 }
683 }
684 else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow) {
685
686 Finish = Standard_True;
687 FFxu = 0.0;
688 FFi = FFyu;
689 }
690 else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow) {
691
692 Finish = Standard_True;
693 FFyu =0.0;
694 FFi = FFxu;
695 }
696 while (!Finish) {
697
698 // calcul des 2 solutions annulant la derivee de l interpolation cubique
699 // Ambda*t=(U-Xu) F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
700 // df=aa*t*t+2*bb*t+cc
701
702 AA = 3.*(Ambda*(DFFxu+DFFyu)+2.*(FFxu-FFyu));
703 BB = -2*(Ambda*(DFFyu+2.*DFFxu)+3.*(FFxu-FFyu));
704 CC = Ambda*DFFxu;
705
706 if(Abs(AA)<1e-14 && Abs(BB)<1e-14 && Abs(CC)<1e-14) {
707 AA=BB=CC=0;
708 }
709
710
711 math_DirectPolynomialRoots Solv(AA,BB,CC);
712 if ( !Solv.InfiniteRoots() ) {
713 Nn=Solv.NbSolutions();
714 if (Nn <= 0) {
715 if (Fxu*Fyu < 0.) { Alfa1 = 0.5;}
716 else Finish = Standard_True;
717 }
718 else {
719 Alfa1 = Solv.Value(1);
720 if (Nn == 2) {
721 Alfa2 = Solv.Value(2);
722 if (Alfa1 > Alfa2){
723 AA = Alfa1;
724 Alfa1 = Alfa2;
725 Alfa2 = AA;
726 }
727 if (Alfa1 > 1. || Alfa2 < 0.){
728 // resolution par dichotomie
729 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
730 else Finish = Standard_True;
731 }
732 else if ( Alfa1 < 0. ||
733 ( DFFxu > 0. && DFFyu >= 0.) ) {
734 // si 2 derivees >0
735 //(cas changement de signe de la distance signee sans
736 // changement de signe de la derivee:
737 //cas de 'presque'tangence avec 2
738 // solutions proches) ,on prend la plus grane racine
739 if (Alfa2 > 1.) {
740 if (Fxu*Fyu < 0.) Alfa1 = 0.5;
741 else Finish = Standard_True;
742 }
743 else Alfa1 = Alfa2;
744 }
745 }
746 }
747 }
748 else if (Fxu*Fyu < -1e-14) Alfa1 = 0.5;
749 //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
750 else Finish = Standard_True;
751
752 if (!Finish) {
753 // petits tests pour diminuer le nombre d iterations
754 if (Alfa1 <= EpsX ) {
755 Alfa1+=Alfa1;
756 }
757 else if (Alfa1 >= (1.-EpsX) ) {
758 Alfa1 = Alfa1+Alfa1-1.;
759 }
760 Alfa1 = Ambda*(1. - Alfa1);
761 U = U + T - Alfa1;
762 if ( U <= Lowr ) {
763 AA = FLowr + (U - Lowr)*DFLowr;
764 BB = DFLowr;
765 }
766 else if ( U >= Upp ) {
767 AA = FUpp + (U - Upp)*DFUpp;
768 BB = DFUpp;
769 }
770 else {
771 Ok = F.Values(U,AA,BB);
772
773 if(!Ok) {
774 Done = Standard_False;
775 return;
776 }
777
778 AA = AA - K;
779 }
780 FFi = AA*AA;
781 CC = AA*BB; CC+=CC;
782 if ( ( ( CC < 0. && FFi < FFxu ) || DFFxu > 0.)
783 && AA*Fxu > 0. ) {
784 FFxu = FFi;
785 DFFxu = CC;
786 Fxu = AA;
787 DFxu = BB;
788 T = Alfa1;
789 if (Alfa1 > Ambda*0.5) {
790 // remarque (1)
791 // determination d 1 autre borne pour diviser
792 //le nouvel intervalle par 2 au -
793 Xu = U + Alfa1*0.5;
794 if ( Xu <= Lowr ) {
795 AA = FLowr + (Xu - Lowr)*DFLowr;
796 BB = DFLowr;
797 }
798 else if ( Xu >= Upp ) {
799 AA = FUpp + (Xu - Upp)*DFUpp;
800 BB = DFUpp;
801 }
802 else {
803 Ok = F.Values(Xu,AA,BB);
804
805 if(!Ok) {
806 Done = Standard_False;
807 return;
808 }
809
810 AA = AA -K;
811 }
812 FFi = AA*AA;
813 CC = AA*BB; CC+=CC;
814 if ( (( CC >= 0. || FFi >= FFxu ) && DFFxu <0.)
815 || Fxu * AA < 0. ) {
816 Fyu = AA;
817 DFyu = BB;
818 FFyu = FFi;
819 DFFyu = CC;
820 T = Alfa1*0.5;
821 Ambda = Alfa1*0.5;
822 }
823 else if ( AA*Fyu < 0. && AA*Fxu >0.) {
824 // changement de signe sur l intervalle u,U
825 Fxu = AA;
826 DFxu = BB;
827 FFxu = FFi;
828 DFFxu = CC;
829 FFi = Min(FFxu,FFyu);
830 T = Alfa1*0.5;
831 Ambda = Alfa1*0.5;
832 U = Xu;
833 }
834 else { Ambda = Alfa1;}
835 }
836 else { Ambda = Alfa1;}
837 }
838 else {
839 Fyu = AA;
840 DFyu = BB;
841 FFyu = FFi;
842 DFFyu = CC;
843 if ( (Ambda-Alfa1) > Ambda*0.5 ) {
844 // meme remarque (1)
845 Xu = U - (Ambda-Alfa1)*0.5;
846 if ( Xu <= Lowr ) {
847 AA = FLowr + (Xu - Lowr)*DFLowr;
848 BB = DFLowr;
849 }
850 else if ( Xu >= Upp ) {
851 AA = FUpp + (Xu - Upp)*DFUpp;
852 BB = DFUpp;
853 }
854 else {
855 Ok = F.Values(Xu,AA,BB);
856
857 if(!Ok) {
858 Done = Standard_False;
859 return;
860 }
861
862 AA = AA - K;
863 }
864 FFi = AA*AA;
865 CC = AA*BB; CC+=CC;
866 if ( AA*Fyu <=0. && AA*Fxu > 0.) {
867 FFxu = FFi;
868 DFFxu = CC;
869 Fxu = AA;
870 DFxu = BB;
871 Ambda = ( Ambda - Alfa1 )*0.5;
872 T = 0.;
873 }
874 else if ( AA*Fxu < 0. && AA*Fyu > 0.) {
875 FFyu = FFi;
876 DFFyu = CC;
877 Fyu = AA;
878 DFyu = BB;
879 Ambda = ( Ambda - Alfa1 )*0.5;
880 T = 0.;
881 FFi = Min(FFxu , FFyu);
882 U = Xu;
883 }
884 else {
885 T =0.;
886 Ambda = Ambda - Alfa1;
887 }
888 }
889 else {
890 T =0.;
891 Ambda = Ambda - Alfa1;
892 }
893 }
894 // tests d arrets
895 if (Abs(FFxu) <= Standard_Underflow ||
896 (Abs(DFFxu) <= Standard_Underflow && Fxu*Fyu > 0.)
897 ){
898 Finish = Standard_True;
899 if (Abs(FFxu) <= Standard_Underflow ) {FFxu =0.0;}
900 FFi = FFyu;
901 }
902 else if (Abs(FFyu) <= Standard_Underflow ||
903 (Abs(DFFyu) <= Standard_Underflow && Fxu*Fyu > 0.)
904 ){
905 Finish = Standard_True;
906 if (Abs(FFyu) <= Standard_Underflow ) {FFyu =0.0;}
907 FFi = FFxu;
908 }
909 else {
910 Iter = Iter + 1;
911 Finish = Iter >= Nbiter || (Ambda <= EpsX &&
912 ( Fxu*Fyu >= 0. || FFi <= EpsF*EpsF));
913 }
914 }
915 } // fin interpolation cubique
916
917 // restitution du meilleur resultat
918
919 if ( FFxu < FFi ) { U = U + T -Ambda;}
920 else if ( FFyu < FFi ) { U = U + T;}
921
922 if ( U >= (Lowr -EpsX) && U <= (Upp + EpsX) ) {
923 U = Max( Lowr , Min( U , Upp));
924 Ok = F.Value(U,FFi);
925 FFi = FFi - K;
926 if ( Abs(FFi) < EpsF ) {
927 //coherence
928 if (Abs(Fxu) <= Standard_Underflow) { AA = DFxu;}
929 else if (Abs(Fyu) <= Standard_Underflow) { AA = DFyu;}
930 else if (Fxu*Fyu > 0.) { AA = 0.;}
931 else { AA = Fyu-Fxu;}
932 if (!Sol.IsEmpty()) {
933 if (Abs(Sol.Last() - U) > 5.*EpsX
934 || OldDF*AA < 0. ) {
935 Sol.Append(U);
936 NbStateSol.Append(F.GetStateNumber());
937 }
938 }
939 else {
940 Sol.Append(U);
941 NbStateSol.Append(F.GetStateNumber());
942 }
943 OldDF = AA ;
944 }
945 }
946 DFFyu = 0.;
947 Fyu = 0.;
948 Nn = 1;
949 while ( Nn < 1000000 && DFFyu <= 0.) {
950 // on repart d 1 pouyem plus loin
951 U = U + Nn*EpsX;
952 if ( U <= Lowr ) {
953 AA = FLowr + (U - Lowr)*DFLowr;
954 BB = DFLowr;
955 }
956 else if ( U >= Upp ) {
957 AA = FUpp + (U - Upp)*DFUpp;
958 BB = DFUpp;
959 }
960 else {
961 Ok = F.Values(U,AA,BB);
962 AA = AA - K;
963 }
964 if ( AA*Fyu < 0.) {
965 U = U - Nn*EpsX;
966 Nn = 1000001;
967 }
968 else {
969 Fyu = AA;
970 DFyu = BB;
971 FFyu = AA*AA;
972 DFFyu = 2.*DFyu*Fyu;
973 Nn = 2*Nn;
974 }
975 }
976 }
977 }
978 }
979#if NEWSEQ
980 if(methode == 3) {
981 Standard_Integer n1=StaticSol.Length();
982 Standard_Integer n2=Sol.Length();
983 if(n1!=n2) {
984 printf("\n mathFunctionRoots : n1=%d n2=%d EpsF=%g EpsX=%g\n",n1,n2,EpsF,NEpsX);
985 for(Standard_Integer x1=1; x1<=n1;x1++) {
986 Standard_Real v; F.Value(StaticSol(x1),v); v-=K;
987 printf(" (%+13.8g:%+13.8g) ",StaticSol(x1),v);
988 }
989 printf("\n");
990 for(Standard_Integer x2=1; x2<=n2; x2++) {
991 Standard_Real v; F.Value(Sol(x2),v); v-=K;
992 printf(" (%+13.8g:%+13.8g) ",Sol(x2),v);
993 }
994 printf("\n");
995 }
996 Standard_Integer n=n1;
997 if(n1>n2) n=n2;
998 for(Standard_Integer i=1;i<=n;i++) {
999 Standard_Real t = Sol(i)-StaticSol(i);
7fd59977 1000 if(Abs(t)>NEpsX) {
1001 printf("\n mathFunctionRoots : i:%d/%d delta: %g",i,n,t);
1002 }
1003 }
1004 }
1005#endif
1006}
1007}
1008
1009
1010void math_FunctionRoots::Dump(Standard_OStream& o) const
1011{
1012
1013 o << "math_FunctionRoots ";
1014 if(Done) {
1015 o << " Status = Done \n";
1016 o << " Number of solutions = " << Sol.Length() << endl;
1017 for (Standard_Integer i = 1; i <= Sol.Length(); i++) {
1018 o << " Solution Number " << i << "= " << Sol.Value(i) << endl;
1019 }
1020 }
1021 else {
1022 o<< " Status = not Done \n";
1023 }
1024}
Open CASCADE Technology repository