1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
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.
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.
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.
20 #define No_Standard_RangeError
21 #define No_Standard_OutOfRange
22 #define No_Standard_DimensionError
25 #include <math_DirectPolynomialRoots.ixx>
27 #include <Standard_RangeError.hxx>
28 #include <StdFail_InfiniteSolutions.hxx>
32 // Reference pour solution equation 3ieme degre et 2ieme degre :
33 // ALGORITHMES NUMERIQUES ANALYSE ET MISE EN OEUVRE, tome 2
34 // (equations et systemes non lineaires)
35 // J. VIGNES editions TECHNIP.
37 const Standard_Real ZERO = 1.0e-30;
38 const Standard_Real EPSILON = RealEpsilon();
39 const Standard_Real RADIX = 2;
40 const Standard_Real Un_Sur_Log_RADIX = 1.0/log(2.0);
42 static Standard_Real Value(const Standard_Integer N, Standard_Real *Poly, const Standard_Real X) {
44 Standard_Real Result = Poly[0];
45 for(Standard_Integer Index = 1; Index < N; Index++) {
46 Result = Result * X + Poly[Index];
52 static void Values(const Standard_Integer N, Standard_Real *Poly, const Standard_Real X,
53 Standard_Real& Val, Standard_Real& Der) {
55 Val = Poly[0] * X + Poly[1];
57 for(Standard_Integer Index = 2; Index < N; Index++) {
59 Val = Val * X + Poly[Index];
63 static Standard_Real Improve(const Standard_Integer N, Standard_Real *Poly, const Standard_Real IniSol) {
65 Standard_Real Val = 0., Der, Delta;
66 Standard_Real Sol = IniSol;
67 Standard_Real IniVal = Value(N, Poly, IniSol);
68 Standard_Integer Index;
70 // cout << "Improve\n";
71 for(Index = 1; Index < 10; Index++) {
72 Values(N, Poly, Sol, Val, Der);
73 if(Abs(Der) <= ZERO) break;
75 if(Abs(Delta) <= EPSILON * Abs(Sol)) break;
77 // cout << " Iter = " << Index << " Delta = " << Delta
78 // << " Val = " << Val << " Der = " << Der << "\n";
80 if(Abs(Val) <= Abs(IniVal)) {
88 Standard_Real Improve(const Standard_Real A, const Standard_Real B, const Standard_Real C,
89 const Standard_Real D, const Standard_Real E, const Standard_Real IniSol) {
91 Standard_Real Poly[5];
97 return Improve(5, Poly, IniSol);
100 Standard_Real Improve(const Standard_Real A, const Standard_Real B,
101 const Standard_Real C, const Standard_Real D, const Standard_Real IniSol) {
103 Standard_Real Poly[4];
108 return Improve(4, Poly, IniSol);
111 Standard_Real Improve(const Standard_Real A, const Standard_Real B,
112 const Standard_Real C, const Standard_Real IniSol) {
114 Standard_Real Poly[3];
118 return Improve(3, Poly, IniSol);
121 Standard_Integer BaseExponent(const Standard_Real X) {
124 return (Standard_Integer)(log(X) * Un_Sur_Log_RADIX);
127 return (Standard_Integer)(-log(-X) * Un_Sur_Log_RADIX);
135 math_DirectPolynomialRoots::math_DirectPolynomialRoots(const Standard_Real A,
136 const Standard_Real B,
137 const Standard_Real C,
138 const Standard_Real D,
139 const Standard_Real E) {
140 InfiniteStatus = Standard_False;
141 Done = Standard_True;
142 Solve(A, B, C, D, E);
145 math_DirectPolynomialRoots::math_DirectPolynomialRoots(const Standard_Real A,
146 const Standard_Real B,
147 const Standard_Real C,
148 const Standard_Real D) {
149 Done = Standard_True;
150 InfiniteStatus = Standard_False;
154 math_DirectPolynomialRoots::math_DirectPolynomialRoots(const Standard_Real A,
155 const Standard_Real B,
156 const Standard_Real C) {
157 Done = Standard_True;
158 InfiniteStatus = Standard_False;
162 math_DirectPolynomialRoots::math_DirectPolynomialRoots(const Standard_Real A,
163 const Standard_Real B) {
164 Done = Standard_True;
165 InfiniteStatus = Standard_False;
170 void math_DirectPolynomialRoots::Solve(const Standard_Real a,
171 const Standard_Real b,
172 const Standard_Real c,
173 const Standard_Real d,
174 const Standard_Real e) {
180 //// modified by jgv, 22.01.09 ////
181 Standard_Real aZero = ZERO;
182 Standard_Real Abs_b = Abs(b), Abs_c = Abs(c), Abs_d = Abs(d), Abs_e = Abs(e);
193 aZero = Epsilon(100.*aZero);
195 if(Abs(a) <= aZero) {
196 Standard_Real aZero1000 = 1000.*aZero;
197 Standard_Boolean with_a = Standard_False;
198 if (Abs_b > ZERO && Abs_b <= aZero1000)
199 with_a = Standard_True;
200 if (Abs_c > ZERO && Abs_c <= aZero1000)
201 with_a = Standard_True;
202 if (Abs_d > ZERO && Abs_d <= aZero1000)
203 with_a = Standard_True;
204 if (Abs_e > ZERO && Abs_e <= aZero1000)
205 with_a = Standard_True;
213 ///////////////////////////////////
215 Standard_Real A, B, C, D, R3, S3, T3, Q3, Y0, P0, Q0, P, Q, P1, Q1;
216 Standard_Real Discr, Sdiscr;
217 Standard_Integer Index;
218 Standard_Integer Exp;
219 Standard_Real PowRadix1,PowRadix2;
225 Exp = BaseExponent(D) / 4;
227 //-- A = A / pow(RADIX, Exp);
228 //-- B = B / pow(RADIX, 2 * Exp);
229 //-- C = C / pow(RADIX, 3 * Exp);
230 //-- D = D / pow(RADIX, 4 * Exp);
231 PowRadix1 = pow(RADIX,Exp);
232 A/= PowRadix1; PowRadix2 = PowRadix1 * PowRadix1;
234 C/= PowRadix2 * PowRadix1;
235 D/= PowRadix2 * PowRadix2;
238 S3 = A * C - 4.0 * D;
239 T3 = D * (4.0 * B - A * A) - C * C;
241 math_DirectPolynomialRoots Sol3(Q3, R3, S3, T3);
242 //-- ################################################################################
243 if(Sol3.IsDone() == Standard_False) { Done = Standard_False; return; }
244 //-- ################################################################################
248 for(Index = 2; Index <= Sol3.NbSolutions(); Index++) {
249 if(Sol3.Value(Index) > Y0) Y0 = Sol3.Value(Index);
251 Discr = A * Y0 * 0.5 - C;
258 P0 = A * A * 0.25 - B + Y0;
259 if(P0 < 0.0) P0 = 0.0;
261 Q0 = Y0 * Y0 * 0.25 - D;
262 if(Q0 < 0.0) Q0 = 0.0;
265 Standard_Real Ademi = A * 0.5;
266 Standard_Real Ydemi = Y0 * 0.5;
267 Standard_Real SdiscrQ0 = Sdiscr * Q0;
270 Q = Ydemi + SdiscrQ0;
272 Q1 = Ydemi - SdiscrQ0;
273 // Modified by skv - Wed Apr 14 16:05:24 2004 IDEM(Airbus) Begin
276 eps = Epsilon(100.*Max(Ademi, P0));
282 eps = Epsilon(100.*Max(Ydemi, SdiscrQ0));
287 // Modified by skv - Wed Apr 14 16:05:24 2004 IDEM(Airbus) End
290 math_DirectPolynomialRoots ASol2(Ademi, P, Q);
291 //-- ################################################################################
292 if(ASol2.IsDone() == Standard_False) { Done = Standard_False; return; }
293 //-- ################################################################################
294 math_DirectPolynomialRoots BSol2(Ademi, P1, Q1);
295 //-- ################################################################################
296 if(BSol2.IsDone() == Standard_False) { Done = Standard_False; return; }
297 //-- ################################################################################
299 NbSol = ASol2.NbSolutions() + BSol2.NbSolutions();
300 for(Index = 0; Index < ASol2.NbSolutions(); Index++) {
301 TheRoots[Index] = ASol2.TheRoots[Index];
303 for(Index = 0; Index < BSol2.NbSolutions(); Index++) {
304 TheRoots[ASol2.NbSolutions() + Index] = BSol2.TheRoots[Index];
306 for(Index = 0; Index < NbSol; Index++) {
307 TheRoots[Index] = TheRoots[Index] * PowRadix1;
308 TheRoots[Index] = Improve(a, b, c, d, e, TheRoots[Index]);
312 void math_DirectPolynomialRoots::Solve(const Standard_Real A,
313 const Standard_Real B,
314 const Standard_Real C,
315 const Standard_Real D) {
322 Standard_Real Beta, Gamma, Del, P1, P2, P, Ep, Q1, Q2, Q3, Q, Eq, A1, A2, Discr;
323 Standard_Real Sigma, Psi, D1, D2, Sb, Omega, Sp3, Y1, Dbg, Sdbg, Den1, Den2;
324 Standard_Real U, H, Sq;
325 Standard_Integer Exp;
331 Exp = BaseExponent(Del) / 3;
333 Standard_Real PowRadix1 = pow(RADIX,Exp);
334 Standard_Real PowRadix2 = PowRadix1*PowRadix1;
337 Del/= PowRadix2 * PowRadix1;
338 //-- Beta = Beta / pow(RADIX, Exp);
339 //-- Gamma = Gamma / pow(RADIX, 2 * Exp);
340 //-- Del = Del / pow(RADIX, 3 * Exp);
343 P2 = - (Beta * Beta) / 3.0;
345 Ep = 5.0 * EPSILON * (Abs(P1) + Abs(P2));
346 if(Abs(P) <= Ep) P = 0.0;
348 Q2 = - Beta * Gamma / 3.0;
349 Q3 = 2.0 * (Beta * Beta * Beta) / 27.0;
351 Eq = 10.0 * EPSILON * (Abs(Q1) + Abs(Q2) + Abs(Q3));
352 if(Abs(Q) <= Eq) Q = 0.0;
353 //-- ############################################################
354 Standard_Real AbsP = P; if(P<0.0) AbsP = -P;
355 if(AbsP>1e+80) { Done = Standard_False; return; }
356 //-- ############################################################
357 A1 = (P * P * P) / 27.0;
362 Psi = Gamma * Gamma * (4.0 * Gamma - Beta * Beta) / 27.0;
364 D1 = Sigma + 2.0 * sqrt(-A1);
367 D1 = Sigma - 2.0 * sqrt(-A1);
371 if(Abs(Del - D1) >= 18.0 * EPSILON * (Abs(Del) + Abs(D1)) &&
372 Abs(Del - D2) >= 24.0 * EPSILON * (Abs(Del) + Abs(D2))) {
373 Discr = (Del - D1) * (Del - D2) / 4.0;
384 if(Beta == 0.0 && Q == 0.0) {
385 TheRoots[0] = sqrt(-P);
386 TheRoots[1] = -TheRoots[0];
390 Omega = atan(0.5 * Q / sqrt(- Discr));
391 Sp3 = sqrt(-P / 3.0);
392 Y1 = -2.0 * Sb * Sp3 * cos(M_PI / 6.0 - Sb * Omega / 3.0);
393 TheRoots[0] = - Beta / 3.0 + Y1;
394 if(Beta * Q <= 0.0) {
395 TheRoots[1] = - Beta / 3.0 + 2.0 * Sp3 * sin(Omega / 3.0);
398 Dbg = Del - Beta * Gamma;
405 Den1 = 8.0 * Beta * Beta / 9.0 - 4.0 * Beta * Y1 / 3.0
407 Den2 = 2.0 * Y1 * Y1 - Q / Y1;
408 TheRoots[1] = Dbg / Den1 + Sdbg * sqrt(-27.0 * Discr) / Den2;
410 TheRoots[2] = - Del / (TheRoots[0] * TheRoots[1]);
413 else if(Discr > 0.0) {
415 U = sqrt(Discr) + Abs(Q / 2.0);
417 U = pow(U, 1.0 / 3.0);
420 U = - pow(Abs(U), 1.0 / 3.0);
423 H = U * U + P / 3.0 + (P / U) * (P / U) / 9.0;
426 H = U * Abs(Q) / (U * U - P / 3.0);
428 if(Beta * Q >= 0.0) {
429 if(Abs(H) <= RealSmall() && Abs(Q) <= RealSmall()) {
430 TheRoots[0] = - Beta / 3.0 - U + P / (3.0 * U);
433 TheRoots[0] = - Beta / 3.0 - Q / H;
437 TheRoots[0] = - Del / (Beta * Beta / 9.0 + H - Beta * Q / (3.0 * H));
448 Sp3 = sqrt(-P / 3.0);
449 if(Beta * Q <= 0.0) {
450 TheRoots[0] = -Beta / 3.0 + Sq * Sp3;
451 TheRoots[1] = TheRoots[0];
452 if(Beta * Q == 0.0) {
453 TheRoots[2] = -Beta / 3.0 - 2.0 * Sq * Sp3;
456 TheRoots[2] = - Del / (TheRoots[0] * TheRoots[1]);
460 TheRoots[0] = -Gamma / (Beta + 3.0 * Sq * Sp3);
461 TheRoots[1] = TheRoots[0];
462 TheRoots[2] = -Beta / 3.0 - 2.0 * Sq * Sp3;
465 for(Standard_Integer Index = 0; Index < NbSol; Index++) {
466 TheRoots[Index] = TheRoots[Index] * pow(RADIX, Exp);
467 TheRoots[Index] = Improve(A, B, C, D, TheRoots[Index]);
471 void math_DirectPolynomialRoots::Solve(const Standard_Real A,
472 const Standard_Real B,
473 const Standard_Real C) {
480 Standard_Real EpsD = 3.0 * EPSILON * (B * B + Abs(4.0 * A * C));
481 Standard_Real Discrim = B * B - 4.0 * A * C;
483 if(Abs(Discrim) <= EpsD) Discrim = 0.0;
487 else if(Discrim == 0.0) {
489 TheRoots[0] = -0.5 * B / A;
490 TheRoots[0] = Improve(A, B, C, TheRoots[0]);
491 TheRoots[1] = TheRoots[0];
496 TheRoots[0] = - (B + sqrt(Discrim)) / (2.0 * A);
499 TheRoots[0] = - (B - sqrt(Discrim)) / (2.0 * A);
501 TheRoots[0] = Improve(A, B, C, TheRoots[0]);
502 TheRoots[1] = C / (A * TheRoots[0]);
503 TheRoots[1] = Improve(A, B, C, TheRoots[1]);
507 void math_DirectPolynomialRoots::Solve(const Standard_Real A,
508 const Standard_Real B) {
511 if (Abs(B) <= ZERO) {
512 InfiniteStatus = Standard_True;
519 TheRoots[0] = -B / A;
522 void math_DirectPolynomialRoots::Dump(Standard_OStream& o) const {
523 o << "math_DirectPolynomialRoots ";
527 else if(InfiniteStatus) {
528 o << " Status = Infinity Roots \n";
530 else if (!InfiniteStatus) {
531 o << " Status = Not Infinity Roots \n";
532 o << " Number of solutions = " << NbSol <<"\n";
533 for (Standard_Integer i = 1; i <= NbSol; i++) {
534 o << " Solution number " << i << " = " << TheRoots[i-1] <<"\n";