1 // Created on: 1996-01-22
2 // Created by: Philippe MANGIN
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #ifndef OCCT_DEBUG
18 #define No_Standard_RangeError
19 #define No_Standard_OutOfRange
20 #endif
21
22
23 #include <FairCurve_Energy.hxx>
24 #include <gp_Pnt2d.hxx>
25 #include <gp_Vec2d.hxx>
26 #include <math_IntegerVector.hxx>
27 #include <math_Matrix.hxx>
28
29 //=======================================================================
FairCurve_Energy(const Handle (TColgp_HArray1OfPnt2d)& Poles,const Standard_Integer ContrOrder1,const Standard_Integer ContrOrder2,const Standard_Boolean WithAuxValue,const Standard_Real Angle1,const Standard_Real Angle2,const Standard_Integer Degree,const Standard_Real Curvature1,const Standard_Real Curvature2)30 FairCurve_Energy::FairCurve_Energy(const Handle(TColgp_HArray1OfPnt2d)& Poles,
31 const Standard_Integer ContrOrder1,
32 const Standard_Integer ContrOrder2,
33 const Standard_Boolean WithAuxValue,
34 const Standard_Real Angle1,
35 const Standard_Real Angle2,
36 const Standard_Integer Degree,
37 const Standard_Real Curvature1,
38 const Standard_Real Curvature2 )
39 //=======================================================================
40 : MyPoles (Poles),
41 MyContrOrder1(ContrOrder1),
42 MyContrOrder2(ContrOrder2),
43 MyWithAuxValue(WithAuxValue),
44 MyNbVar(2*(Poles->Length()-2) - ContrOrder2 - ContrOrder1 + WithAuxValue),
45 MyNbValues(2*Poles->Length() + WithAuxValue),
46 MyLinearForm(0, 1),
47 MyQuadForm(0, 1),
48 MyGradient( 0, MyNbValues),
49 MyHessian( 0, MyNbValues + MyNbValues*(MyNbValues+1)/2 )
50 {
51 // chesk angles in reference (Ox,Oy)
52 gp_XY L0 (Cos(Angle1), Sin(Angle1)), L1 (-Cos(Angle2), Sin(Angle2));
53 MyLinearForm.SetValue(0, L0);
54 MyLinearForm.SetValue(1, L1);
55 gp_XY Q0(-Sin(Angle1), Cos(Angle1)), Q1 (Sin(Angle2), Cos(Angle2));
56 MyQuadForm.SetValue(0, ((double)Degree) / (Degree-1) * Curvature1 * Q0);
57 MyQuadForm.SetValue(1, ((double)Degree) / (Degree-1) * Curvature2 * Q1);
58 }
59
60 //=======================================================================
Value(const math_Vector & X,Standard_Real & E)61 Standard_Boolean FairCurve_Energy::Value(const math_Vector& X,
62 Standard_Real& E)
63 //=======================================================================
64 {
65 Standard_Boolean IsDone;
66 math_Vector Energie(0,0);
67 ComputePoles(X);
68 IsDone = Compute(0, Energie);
69 E = Energie(0);
70 return IsDone;
71 }
72
73 //=======================================================================
Gradient(const math_Vector & X,math_Vector & G)74 Standard_Boolean FairCurve_Energy::Gradient(const math_Vector& X,
75 math_Vector& G)
76 //=======================================================================
77 {
78 Standard_Boolean IsDone;
79 Standard_Real E;
80
81 IsDone = Values(X, E, G);
82 return IsDone;
83 }
84
85 //=======================================================================
Gradient1(const math_Vector & Vect,math_Vector & Grad)86 void FairCurve_Energy::Gradient1(const math_Vector& Vect,
87 math_Vector& Grad)
88 //=======================================================================
89 {
90 Standard_Integer ii,
91 DebG = Grad.Lower(), FinG = Grad.Upper();
92 Standard_Integer Vdeb = 3,
93 Vfin = 2*MyPoles->Length()-2;
94
95 // .... by calculation
96 if (MyContrOrder1 >= 1) {
97 gp_XY DPole (Vect(Vdeb), Vect(Vdeb+1));
98 Grad(DebG) = MyLinearForm(0) * DPole;
99 Vdeb += 2;
100 DebG += 1;
101 }
102 if(MyContrOrder1 == 2) {
103 Standard_Real Lambda0 = MyPoles->Value(MyPoles->Lower())
104 .Distance( MyPoles->Value(MyPoles->Lower()+1) );
105 gp_XY DPole (Vect(Vdeb), Vect(Vdeb+1));
106 Grad(DebG-1) += (MyLinearForm(0) + 2*Lambda0*MyQuadForm(0)) * DPole;
107 Grad(DebG) = MyLinearForm(0) * DPole;
108 Vdeb += 2;
109 DebG += 1;
110 }
111 if (MyWithAuxValue) {
112 Grad(FinG) = Vect( 2*MyPoles->Length()+1 );
113 FinG -= 1;
114 }
115 if (MyContrOrder2 >= 1) {
116 gp_XY DPole (Vect(Vfin-1), Vect(Vfin));
117 Grad(FinG) = MyLinearForm(1) * DPole;
118 FinG -= 1;
119 }
120 if(MyContrOrder2 == 2) {
121 Standard_Real Lambda1 = MyPoles->Value(MyPoles->Upper())
122 .Distance(MyPoles->Value(MyPoles->Upper()-1) );
123 gp_XY DPole (Vect(Vfin-3), Vect(Vfin-2));
124 Grad(FinG) = Grad(FinG+1) + (MyLinearForm(1) + 2*Lambda1*MyQuadForm(1)) * DPole;
125 Grad(FinG+1) = MyLinearForm(1) * DPole;
126 FinG -= 1;
127 }
128 // ... by recopy
129 for (ii=DebG; ii<=FinG; ii++) {
130 Grad(ii) = Vect(Vdeb);
131 Vdeb += 1;
132 }
133 }
134
135 //=======================================================================
Values(const math_Vector & X,Standard_Real & E,math_Vector & G)136 Standard_Boolean FairCurve_Energy::Values(const math_Vector& X,
137 Standard_Real& E,
138 math_Vector& G)
139 //=======================================================================
140 {
141 Standard_Boolean IsDone;
142
143 ComputePoles(X);
144 IsDone = Compute(1, MyGradient);
145 if (IsDone) {
146 E = MyGradient(0);
147 Gradient1(MyGradient, G);
148 }
149 return IsDone;
150 }
151
152 //=======================================================================
Values(const math_Vector & X,Standard_Real & E,math_Vector & G,math_Matrix & H)153 Standard_Boolean FairCurve_Energy::Values(const math_Vector& X,
154 Standard_Real& E,
155 math_Vector& G,
156 math_Matrix& H)
157 //=======================================================================
158 {
159 Standard_Boolean IsDone;
160
161 ComputePoles(X);
162 IsDone = Compute(2, MyHessian);
163 if (IsDone) {
164 E = MyHessian(0);
165 Gradient1(MyHessian, G);
166 Hessian1 (MyHessian, H);
167 }
168 return IsDone;
169 }
170
171
172 //=======================================================================
Hessian1(const math_Vector & Vect,math_Matrix & H)173 void FairCurve_Energy::Hessian1(const math_Vector& Vect,
174 math_Matrix& H)
175 //=======================================================================
176 {
177
178 Standard_Integer ii, jj, kk, Vk;
179 Standard_Integer Vdeb = 3 + 2*MyContrOrder1,
180 Vfin = 2*MyPoles->Length() - 2*(MyContrOrder2+1),
181 Vup = 2*MyPoles->Length()+MyWithAuxValue;
182 Standard_Integer DebH = 1+MyContrOrder1,
183 FinH = MyNbVar - MyWithAuxValue - MyContrOrder2 ;
184 Standard_Real Cos0 = pow(MyLinearForm(0).X(),2),
185 Sin0 = pow(MyLinearForm(0).Y(),2),
186 CosSin0 = 2 * MyLinearForm(0).X() * MyLinearForm(0).Y(),
187 Cos1 = pow(MyLinearForm(1).X(),2),
188 Sin1 = pow(MyLinearForm(1).Y(),2),
189 CosSin1 = 2 * MyLinearForm(1).X() * MyLinearForm(1).Y() ;
190 Standard_Real Lambda0=0, Lambda1=0;
191
192 if (MyContrOrder1 >= 1) {
193 Lambda0 = MyPoles->Value(MyPoles->Lower())
194 .Distance( MyPoles->Value(MyPoles->Lower()+1) );}
195
196 if (MyContrOrder2 >= 1) {
197 Lambda1 = MyPoles->Value(MyPoles->Upper())
198 .Distance(MyPoles->Value(MyPoles->Upper()-1) );}
199
200
201 if (MyContrOrder1 >= 1) {
202
203 // calculate the left lambda column --------------------------------
204
205 jj = Vdeb-2*MyContrOrder1;
206 kk=Indice(jj, jj); // X2X2
207 ii=Indice(jj+1, jj); // X2Y2
208 H(1, 1) = Cos0 * Vect(kk) + CosSin0*Vect(ii) + Sin0 * Vect(ii+1);
209
210 if (MyContrOrder1 >= 2) {
211 gp_XY Laux = (MyLinearForm(0) + 2*Lambda0*MyQuadForm(0));
212 jj = Vdeb-2*(MyContrOrder1-1);
213 kk=Indice(jj, jj-2); // X1X2
214 ii=Indice(jj+1, jj-2); //X1Y2
215 gp_XY Aux(Vect(kk+2), Vect(ii+3));
216
217 H (1, 1) += 2 * (
218 ( MyQuadForm(0).X() * Vect(5) + MyQuadForm(0).Y() * Vect(6) )
219 + ( Laux.X() * ( MyLinearForm(0).X()*Vect(kk) + MyLinearForm(0).Y()*Vect(kk+1))
220 + Laux.Y() * ( MyLinearForm(0).X()*Vect(ii) + MyLinearForm(0).Y()*Vect(ii+1)) )
221 + Laux.X() * Laux.Y() * Vect(ii+2) )
222 + ( Pow(Laux.X(),2) * Vect(kk+2) + Pow(Laux.Y(),2) * Vect(ii+3));
223
224 H(2,1) = (Cos0 * Vect(kk) + CosSin0*(Vect(ii)+Vect(kk+1))/2 + Sin0 * Vect(ii+1))
225 + Laux * MyLinearForm(0).Multiplied(Aux)
226 + (Laux.X()*MyLinearForm(0).Y() + Laux.Y()*MyLinearForm(0).X()) * Vect(ii+2);
227 }
228
229
230 if (MyWithAuxValue) {
231 kk =Indice(Vup, Vdeb-2*MyContrOrder1);
232 H(MyNbVar, 1) = MyLinearForm(0).X() * Vect(kk)
233 + MyLinearForm(0).Y() * Vect(kk+1);
234 }
235
236 if (MyContrOrder2 >= 1) {
237 H(MyNbVar-MyWithAuxValue, 1) = 0; // correct if there are less than 3 nodes
238 if (MyContrOrder2 == 2) {H(MyNbVar-MyWithAuxValue-1, 1) = 0;}
239 }
240
241
242 Vk = Vdeb;
243 kk = Indice(Vk, Vdeb-2*MyContrOrder1);
244 for (ii=DebH; ii<=FinH; ii++) {
245 H(ii, 1) = MyLinearForm(0).X() * Vect(kk)
246 + MyLinearForm(0).Y() * Vect(kk+1);
247 kk += Vk;
248 Vk++;
249 }
250 }
251
252 // calculate the left line mu ----------------------
253 if (MyContrOrder1 >= 2) {
254 jj = Vdeb-2*(MyContrOrder1-1);
255 kk=Indice(jj, jj); // X3X3
256 ii=Indice(jj+1, jj); // X3Y3
257 H(2, 2) = Cos0 * Vect(kk) + CosSin0*Vect(ii) + Sin0 * Vect(ii+1);
258
259 if (MyWithAuxValue) {
260 kk =Indice(Vup, Vdeb-2*(MyContrOrder1-1));
261 gp_XY Pole (Vect(kk), Vect(kk+1));
262 H(MyNbVar, 1) += (MyLinearForm(0) + 2*Lambda0*MyQuadForm(0)) * Pole;
263 H(MyNbVar, 2) = MyLinearForm(0).X() * Vect(kk)
264 + MyLinearForm(0).Y() * Vect(kk+1);
265 }
266
267 if (MyContrOrder2 >= 1) {
268 H(MyNbVar-MyWithAuxValue, 2) = 0; // correct if there are less than 3 nodes
269 if (MyContrOrder2 == 2) {H(MyNbVar-MyWithAuxValue-1, 2) = 0;}
270 }
271 Vk = Vdeb;
272
273 Standard_Real Xaux = (MyLinearForm(0) + 2*Lambda0*MyQuadForm(0)).X(),
274 Yaux = (MyLinearForm(0) + 2*Lambda0*MyQuadForm(0)).Y();
275
276 kk = Indice(Vk, Vdeb-2*MyContrOrder1+2);
277 for (ii=DebH; ii<=FinH; ii++) {
278 H(ii, 2) = MyLinearForm(0).X() * Vect(kk)
279 + MyLinearForm(0).Y() * Vect(kk+1);
280 H(ii, 1) += Xaux * Vect(kk) + Yaux*Vect(kk+1);
281 kk += Vk;
282 Vk++;
283 }
284 }
285
286 // calculate the right lambda line -----------------------
287 if (MyContrOrder2 >= 1) {
288
289 jj = FinH + 1;
290 Vk = Vfin + 2*MyContrOrder2 - 1;
291 kk = Indice(Vk, Vdeb);
292 for (ii=DebH; ii<=FinH; ii++) {
293 H(jj, ii) = MyLinearForm(1).X() * Vect(kk)
294 + MyLinearForm(1).Y() * Vect(kk+Vk);
295 kk++;
296 }
297
298 kk = Indice(Vk, Vk);
299 H(jj, jj) = Cos1 * Vect(kk) + CosSin1 * Vect(kk+Vk) + Sin1 * Vect(kk+Vk+1);
300
301 if (MyContrOrder2 >= 2) {
302 // H(jj,jj) +=
303 gp_XY Laux = (MyLinearForm(1) + 2*Lambda1*MyQuadForm(1));
304 jj = Vfin + 2*MyContrOrder2 - 3;
305 kk=Indice(jj+2, jj); // Xn-1Xn-2
306 ii=Indice(jj+3, jj); //Yn-1Xn-2
307 Standard_Integer ll = Indice(jj, jj);
308
309 H (FinH+1, FinH+1) += 2 * (
310 ( MyQuadForm(1).X() * Vect(jj) + MyQuadForm(1).Y() * Vect(jj+1) )
311 + ( Laux.X() * ( MyLinearForm(1).X()*Vect(kk) + MyLinearForm(1).Y()*Vect(ii))
312 + Laux.Y() * ( MyLinearForm(1).X()*Vect(kk+1) + MyLinearForm(1).Y()*Vect(ii+1)) )
313 + Laux.X() * Laux.Y() * Vect(ll+jj) )
314 + ( Pow(Laux.X(),2) * Vect(ll) + Pow(Laux.Y(),2) * Vect(ll+jj+1));
315
316 H(FinH+2, FinH+1) = Cos1 * Vect(kk) + CosSin1*(Vect(ii)+Vect(kk+1))/2 + Sin1 * Vect(ii+1);
317 gp_XY Aux(Vect(ll), Vect(ll+jj+1));
318 H(FinH+2, FinH+1) += Laux * MyLinearForm(1).Multiplied(Aux)
319 + (Laux.X()*MyLinearForm(1).Y() + Laux.Y()*MyLinearForm(1).X())
320 * Vect(ll+jj);
321 // H(FinH+2, FinH+1) = 0; // No better alternative. Bug in the previous expression
322 }
323 }
324
325 // calculate the right line mu -----------------------
326 if (MyContrOrder2 >= 2) {
327 jj = FinH + 2;
328 Vk = Vfin + 2*MyContrOrder2 - 3;
329 kk = Indice(Vk, Vdeb);
330
331 Standard_Real Xaux = (MyLinearForm(1) + 2*Lambda1*MyQuadForm(1)).X(),
332 Yaux = (MyLinearForm(1) + 2*Lambda1*MyQuadForm(1)).Y();
333
334 for (ii=DebH; ii<=FinH; ii++) {
335 H(jj, ii) = MyLinearForm(1).X() * Vect(kk)
336 + MyLinearForm(1).Y() * Vect(kk+Vk);
337 // update the right line Lambda
338 H(jj-1,ii) += Xaux*Vect(kk) + Yaux*Vect(kk+Vk);
339 kk++;
340 }
341 kk = Indice(Vk, Vk);
342 ii = Indice(Vk+1, Vk);
343 H(jj,jj) = Cos1*Vect(kk) + CosSin1*Vect(ii) + Sin1*Vect(ii+1);
344 }
345
346 // calculate the Auxiliary Variable line -----------------------
347 if (MyWithAuxValue) {
348
349 kk = Indice(Vup, Vdeb);
350 for (ii=DebH; ii<=FinH; ii++) {
351 H(MyNbVar, ii) = Vect(kk);
352 kk++;
353 }
354
355 if (MyContrOrder2 >= 1) {
356 kk = Indice(Vup, Vfin+2*MyContrOrder2-1);
357 H(MyNbVar, FinH+1) =
358 MyLinearForm(1).X() * Vect(kk) + MyLinearForm(1).Y() * Vect(kk+1);
359 }
360 if (MyContrOrder2 >= 2) {
361 kk = Indice(Vup, Vfin+2*MyContrOrder2-3);
362 gp_XY Pole( Vect(kk), Vect(kk+1));
363 H(MyNbVar, FinH+1) += (MyLinearForm(1) + 2*Lambda1*MyQuadForm(1)) * Pole;
364 H(MyNbVar, FinH+2) = MyLinearForm(1) * Pole;
365 }
366 kk = Indice(Vup, Vup);
367 H(H.UpperRow(), H.UpperRow()) = Vect(kk);
368 }
369
370 // recopy the internal block -----------------------------------
371
372 kk = Indice(Vdeb, Vdeb);
373 for (ii = DebH; ii <=FinH; ii++) {
374 for (jj = DebH; jj<=ii; jj++) {
375 H(ii,jj) = Vect(kk);
376 kk++;
377 }
378 kk += Vdeb-1;
379 }
380 // symmetry
381 for (ii = H.LowerRow(); ii <= H.UpperRow(); ii++)
382 for (jj = ii+1; jj <= H.UpperRow(); jj++) H(ii,jj) = H(jj,ii);
383 }
384
385 //=======================================================================
Variable(math_Vector & X) const386 Standard_Boolean FairCurve_Energy::Variable(math_Vector& X) const
387 //=======================================================================
388 {
389 Standard_Integer ii,
390 IndexDeb1 = MyPoles->Lower()+1,
391 IndexDeb2 = X.Lower(),
392 IndexFin1 = MyPoles->Upper()-1,
393 IndexFin2 = X.Upper() - MyWithAuxValue; // decrease by 1 if the sliding is
394 // free as the last value of X is reserved.
395
396
397 // calculation of variables of constraints
398 if (MyContrOrder1 >= 1) {
399 X(IndexDeb2) = MyPoles->Value(MyPoles->Lower())
400 .Distance( MyPoles->Value(MyPoles->Lower()+1) );
401 IndexDeb1 += 1;
402 IndexDeb2 += 1;
403 }
404 if (MyContrOrder1 == 2) {
405 gp_Vec2d b1b2( MyPoles->Value(MyPoles->Lower()+1),
406 MyPoles->Value(MyPoles->Lower()+2) );
407 X(IndexDeb2) = b1b2.XY() * MyLinearForm(0);
408 IndexDeb1 += 1;
409 IndexDeb2 += 1;
410 }
411 if (MyContrOrder2 == 2) {
412 gp_Vec2d bn2bn1( MyPoles->Value(MyPoles->Upper()-2),
413 MyPoles->Value(MyPoles->Upper()-1));
414 X(IndexFin2) = - bn2bn1.XY() * MyLinearForm(1);
415 IndexFin1 -= 1;
416 IndexFin2 -= 1;
417 }
418 if (MyContrOrder2 >= 1) {
419 X(IndexFin2) = MyPoles->Value(MyPoles->Upper())
420 .Distance(MyPoles->Value(MyPoles->Upper()-1) );
421 IndexFin1 -= 1;
422 }
423
424 // Recopy of auxiliary variables is not done in the abstract class
425
426 // copy poles to variables
427 for (ii=IndexDeb1; ii<=IndexFin1; ii++) {
428 X(IndexDeb2) = MyPoles->Value(ii).X();
429 X(IndexDeb2+1) = MyPoles->Value(ii).Y();
430 IndexDeb2 +=2;
431 }
432 return Standard_True;
433 }
434
435 //=======================================================================
ComputePoles(const math_Vector & X)436 void FairCurve_Energy::ComputePoles(const math_Vector& X)
437 //=======================================================================
438 {
439 Standard_Integer ii,
440 IndexDeb1 = MyPoles->Lower()+1,
441 IndexDeb2 = X.Lower(),
442 IndexFin1 = MyPoles->Upper()-1,
443 IndexFin2 = X.Upper() - MyWithAuxValue; // decrease by 1 if the sliding is
444 // is free as the last value of X is reserved.
445 // calculation of pole constraints
446 // for (ii=MyPoles->Lower();ii<=MyPoles->Upper();ii++) {
447 // std::cout << ii << " X = " << MyPoles->Value(ii).X() <<
448 // " Y = " << MyPoles->Value(ii).Y() << std::endl;}
449
450 if (MyContrOrder1 >= 1) {
451 IndexDeb1 += 1;
452 IndexDeb2 += 1;
453 ComputePolesG1( 0, X(1), MyPoles->Value(MyPoles->Lower()),
454 MyPoles->ChangeValue(MyPoles->Lower()+1) );
455 }
456 if (MyContrOrder1 == 2) {
457 IndexDeb1 += 1;
458 IndexDeb2 += 1;
459 ComputePolesG2( 0, X(1), X(2), MyPoles->Value(MyPoles->Lower()),
460 MyPoles->ChangeValue(MyPoles->Lower()+2) );
461 }
462 if (MyContrOrder2 == 2) {
463 IndexFin1 -= 1;
464 IndexFin2 -= 1;
465 ComputePolesG2( 1, X(IndexFin2), X(IndexFin2+1),
466 MyPoles->Value(MyPoles->Upper()),
467 MyPoles->ChangeValue(MyPoles->Upper()-2) );
468 }
469 if (MyContrOrder2 >= 1) {
470 IndexFin1 -= 1;
471 ComputePolesG1( 1, X(IndexFin2), MyPoles->Value(MyPoles->Upper()),
472 MyPoles->ChangeValue(MyPoles->Upper()-1) );
473 }
474
475 // if (MyWithAuxValue) { MyLengthSliding = X(X.Upper()); }
476 // recopy others
477 for (ii=IndexDeb1; ii<=IndexFin1; ii++) {
478 MyPoles -> ChangeValue(ii).SetX( X(IndexDeb2) );
479 MyPoles -> ChangeValue(ii).SetY( X(IndexDeb2+1) );
480 IndexDeb2 += 2;
481 }
482 }
483
484 //=======================================================================
ComputePolesG1(const Standard_Integer Side,const Standard_Real Lambda,const gp_Pnt2d & P1,gp_Pnt2d & P2) const485 void FairCurve_Energy::ComputePolesG1(const Standard_Integer Side,
486 const Standard_Real Lambda,
487 const gp_Pnt2d& P1,
488 gp_Pnt2d& P2) const
489 //=======================================================================
490 { P2.SetXY ( P1.XY() + MyLinearForm(Side) * Lambda ); }
491
492 //=======================================================================
ComputePolesG2(const Standard_Integer Side,const Standard_Real Lambda,const Standard_Real Rho,const gp_Pnt2d & P1,gp_Pnt2d & P2) const493 void FairCurve_Energy::ComputePolesG2(const Standard_Integer Side,
494 const Standard_Real Lambda,
495 const Standard_Real Rho,
496 const gp_Pnt2d& P1,
497 gp_Pnt2d& P2) const
498 //=======================================================================
499 { P2.SetXY ( P1.XY()
500 + MyLinearForm(Side) * (Lambda + Rho )
501 + MyQuadForm(Side) * (Lambda * Lambda) ) ; }
502
503
504
505
506