1 // Created on: 1995-08-28
2 // Created by: Laurent BOURESCHE
3 // Copyright (c) 1995-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 // Modified: 28/02/1996 by PMN :  HermiteCoefficients added
18 // Modified: 18/06/1996 by PMN :  NULL reference.
19 // Modified: 19/02/1997 by JCT :  EvalPoly2Var added
20 
21 #include <GeomAbs_Shape.hxx>
22 #include <math.hxx>
23 #include <math_Gauss.hxx>
24 #include <math_Matrix.hxx>
25 #include <NCollection_LocalArray.hxx>
26 #include <PLib.hxx>
27 #include <Standard_ConstructionError.hxx>
28 
29 // To convert points array into Real ..
30 // *********************************
31 //=======================================================================
32 //function : SetPoles
33 //purpose  :
34 //=======================================================================
SetPoles(const TColgp_Array1OfPnt2d & Poles,TColStd_Array1OfReal & FP)35 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
36 		    TColStd_Array1OfReal& FP)
37 {
38   Standard_Integer j      = FP   .Lower();
39   Standard_Integer PLower = Poles.Lower();
40   Standard_Integer PUpper = Poles.Upper();
41 
42   for (Standard_Integer i = PLower; i <= PUpper; i++) {
43     const gp_Pnt2d& P = Poles(i);
44     FP(j) = P.Coord(1); j++;
45     FP(j) = P.Coord(2); j++;
46   }
47 }
48 
49 //=======================================================================
50 //function : SetPoles
51 //purpose  :
52 //=======================================================================
53 
SetPoles(const TColgp_Array1OfPnt2d & Poles,const TColStd_Array1OfReal & Weights,TColStd_Array1OfReal & FP)54 void PLib::SetPoles(const TColgp_Array1OfPnt2d&       Poles,
55 		    const TColStd_Array1OfReal& Weights,
56 		    TColStd_Array1OfReal&       FP)
57 {
58   Standard_Integer j      = FP   .Lower();
59   Standard_Integer PLower = Poles.Lower();
60   Standard_Integer PUpper = Poles.Upper();
61 
62   for (Standard_Integer i = PLower; i <= PUpper; i++) {
63     Standard_Real w = Weights(i);
64     const gp_Pnt2d& P = Poles(i);
65     FP(j) = P.Coord(1) * w; j++;
66     FP(j) = P.Coord(2) * w; j++;
67     FP(j) =              w; j++;
68   }
69 }
70 
71 //=======================================================================
72 //function : GetPoles
73 //purpose  :
74 //=======================================================================
75 
GetPoles(const TColStd_Array1OfReal & FP,TColgp_Array1OfPnt2d & Poles)76 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
77 		    TColgp_Array1OfPnt2d& Poles)
78 {
79   Standard_Integer j      = FP   .Lower();
80   Standard_Integer PLower = Poles.Lower();
81   Standard_Integer PUpper = Poles.Upper();
82 
83   for (Standard_Integer i = PLower; i <= PUpper; i++) {
84     gp_Pnt2d& P = Poles(i);
85     P.SetCoord(1,FP(j)); j++;
86     P.SetCoord(2,FP(j)); j++;
87   }
88 }
89 
90 //=======================================================================
91 //function : GetPoles
92 //purpose  :
93 //=======================================================================
94 
GetPoles(const TColStd_Array1OfReal & FP,TColgp_Array1OfPnt2d & Poles,TColStd_Array1OfReal & Weights)95 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
96 		    TColgp_Array1OfPnt2d&       Poles,
97 		    TColStd_Array1OfReal&       Weights)
98 {
99   Standard_Integer j      = FP   .Lower();
100   Standard_Integer PLower = Poles.Lower();
101   Standard_Integer PUpper = Poles.Upper();
102 
103   for (Standard_Integer i = PLower; i <= PUpper; i++) {
104     Standard_Real w = FP(j + 2);
105     Weights(i) = w;
106     gp_Pnt2d& P = Poles(i);
107     P.SetCoord(1,FP(j) / w); j++;
108     P.SetCoord(2,FP(j) / w); j++;
109     j++;
110   }
111 }
112 
113 //=======================================================================
114 //function : SetPoles
115 //purpose  :
116 //=======================================================================
117 
SetPoles(const TColgp_Array1OfPnt & Poles,TColStd_Array1OfReal & FP)118 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
119 		    TColStd_Array1OfReal& FP)
120 {
121   Standard_Integer j      = FP   .Lower();
122   Standard_Integer PLower = Poles.Lower();
123   Standard_Integer PUpper = Poles.Upper();
124 
125   for (Standard_Integer i = PLower; i <= PUpper; i++) {
126     const gp_Pnt& P = Poles(i);
127     FP(j) = P.Coord(1); j++;
128     FP(j) = P.Coord(2); j++;
129     FP(j) = P.Coord(3); j++;
130   }
131 }
132 
133 //=======================================================================
134 //function : SetPoles
135 //purpose  :
136 //=======================================================================
137 
SetPoles(const TColgp_Array1OfPnt & Poles,const TColStd_Array1OfReal & Weights,TColStd_Array1OfReal & FP)138 void PLib::SetPoles(const TColgp_Array1OfPnt&   Poles,
139 		    const TColStd_Array1OfReal& Weights,
140 		    TColStd_Array1OfReal&       FP)
141 {
142   Standard_Integer j      = FP   .Lower();
143   Standard_Integer PLower = Poles.Lower();
144   Standard_Integer PUpper = Poles.Upper();
145 
146   for (Standard_Integer i = PLower; i <= PUpper; i++) {
147     Standard_Real w = Weights(i);
148     const gp_Pnt& P = Poles(i);
149     FP(j) = P.Coord(1) * w; j++;
150     FP(j) = P.Coord(2) * w; j++;
151     FP(j) = P.Coord(3) * w; j++;
152     FP(j) =              w; j++;
153   }
154 }
155 
156 //=======================================================================
157 //function : GetPoles
158 //purpose  :
159 //=======================================================================
160 
GetPoles(const TColStd_Array1OfReal & FP,TColgp_Array1OfPnt & Poles)161 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
162 		    TColgp_Array1OfPnt&         Poles)
163 {
164   Standard_Integer j      = FP   .Lower();
165   Standard_Integer PLower = Poles.Lower();
166   Standard_Integer PUpper = Poles.Upper();
167 
168   for (Standard_Integer i = PLower; i <= PUpper; i++) {
169     gp_Pnt& P = Poles(i);
170     P.SetCoord(1,FP(j)); j++;
171     P.SetCoord(2,FP(j)); j++;
172     P.SetCoord(3,FP(j)); j++;
173   }
174 }
175 
176 //=======================================================================
177 //function : GetPoles
178 //purpose  :
179 //=======================================================================
180 
GetPoles(const TColStd_Array1OfReal & FP,TColgp_Array1OfPnt & Poles,TColStd_Array1OfReal & Weights)181 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
182 		    TColgp_Array1OfPnt&         Poles,
183 		    TColStd_Array1OfReal&       Weights)
184 {
185   Standard_Integer j      = FP   .Lower();
186   Standard_Integer PLower = Poles.Lower();
187   Standard_Integer PUpper = Poles.Upper();
188 
189   for (Standard_Integer i = PLower; i <= PUpper; i++) {
190     Standard_Real w = FP(j + 3);
191     Weights(i) = w;
192     gp_Pnt& P = Poles(i);
193     P.SetCoord(1,FP(j) / w); j++;
194     P.SetCoord(2,FP(j) / w); j++;
195     P.SetCoord(3,FP(j) / w); j++;
196     j++;
197   }
198 }
199 
200 // specialized allocator
201 namespace
202 {
203 
204 class BinomAllocator
205 {
206 public:
207 
208   //! Main constructor
BinomAllocator(const Standard_Integer theMaxBinom)209   BinomAllocator (const Standard_Integer theMaxBinom)
210   : myBinom (NULL),
211     myMaxBinom (theMaxBinom)
212   {
213     Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
214     Standard_Integer np1 = myMaxBinom + 1;
215     myBinom = new Standard_Integer*[np1];
216     myBinom[0] = new Standard_Integer[1];
217     myBinom[0][0] = 1;
218     for (i = 1; i < np1; ++i)
219     {
220       im1 = i - 1;
221       ip1 = i + 1;
222       id2 = i >> 1;
223       md2 = im1 >> 1;
224       md3 = ip1 >> 1;
225       k   = 0;
226       myBinom[i] = new Standard_Integer[ip1];
227 
228       for (j = 0; j < id2; ++j)
229       {
230         myBinom[i][j] = k + myBinom[im1][j];
231         k = myBinom[im1][j];
232       }
233       j = id2;
234       if (j > md2) j = im1 - j;
235       myBinom[i][id2] = k + myBinom[im1][j];
236 
237       for (j = ip1 - md3; j < ip1; j++)
238       {
239         myBinom[i][j] = myBinom[i][i - j];
240       }
241     }
242   }
243 
244   //! Destructor
~BinomAllocator()245   ~BinomAllocator()
246   {
247     // free memory
248     for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
249     {
250       delete[] myBinom[i];
251     }
252     delete[] myBinom;
253   }
254 
Value(const Standard_Integer N,const Standard_Integer P) const255   Standard_Real Value (const Standard_Integer N,
256                        const Standard_Integer P) const
257   {
258     Standard_OutOfRange_Raise_if (N > myMaxBinom,
259       "PLib, BinomAllocator: requested degree is greater than maximum supported");
260     return Standard_Real (myBinom[N][P]);
261   }
262 
263 private:
264   BinomAllocator (const BinomAllocator&);
265   BinomAllocator& operator= (const BinomAllocator&);
266 
267 private:
268   Standard_Integer**  myBinom;
269   Standard_Integer myMaxBinom;
270 
271 };
272 
273   // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
274   //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
275   static BinomAllocator THE_BINOM (25 + 1);
276 }
277 
278 //=======================================================================
279 //function : Bin
280 //purpose  :
281 //=======================================================================
282 
Bin(const Standard_Integer N,const Standard_Integer P)283 Standard_Real PLib::Bin(const Standard_Integer N,
284                         const Standard_Integer P)
285 {
286   return THE_BINOM.Value (N, P);
287 }
288 
289 //=======================================================================
290 //function : RationalDerivative
291 //purpose  :
292 //=======================================================================
293 
RationalDerivative(const Standard_Integer Degree,const Standard_Integer DerivativeRequest,const Standard_Integer Dimension,Standard_Real & Ders,Standard_Real & RDers,const Standard_Boolean All)294 void  PLib::RationalDerivative(const Standard_Integer Degree,
295 			       const Standard_Integer DerivativeRequest,
296 			       const Standard_Integer Dimension,
297 			       Standard_Real& Ders,
298 			       Standard_Real& RDers,
299 			       const Standard_Boolean All)
300 {
301   //
302   // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
303   //
304   //  We Write  u = fv
305   //  Let C(N,P) be the binomial
306   //
307   //        then we have
308   //
309   //         (q)                             (p)   (q-p)
310   //        u    =     SUM          C (q,p) f     v
311   //                p = 0 to q
312   //
313   //
314   //        Therefore
315   //
316   //
317   //         (q)         (   (q)                               (p)   (q-p)   )
318   //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
319   //                     (          p = 0 to q-1                             )
320   //
321   //
322   // make arrays for the binomial since computing it each time could raise a performance
323   // issue
324   // As oppose to the method below the <Der> array is organized in the following
325   // fashion :
326   //
327   //  u (1)     u (2) ....   u  (Dimension)     v (1)
328   //
329   //   (1)       (1)          (1)                (1)
330   //  u (1)     u (2) ....   u  (Dimension)     v (1)
331   //
332   //   ............................................
333   //
334   //   (Degree)  (Degree)     (Degree)           (Degree)
335   //  u (1)     u (2) ....   u  (Dimension)     v (1)
336   //
337   //
338   Standard_Real Inverse;
339   Standard_Real *PolesArray = &Ders;
340   Standard_Real *RationalArray = &RDers;
341   Standard_Real Factor ;
342   Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
343   NCollection_LocalArray<Standard_Real> binomial_array;
344   NCollection_LocalArray<Standard_Real> derivative_storage;
345   if (Dimension == 3) {
346     Standard_Integer DeRequest1 = DerivativeRequest + 1;
347     Standard_Integer MinDegRequ = DerivativeRequest;
348     if (MinDegRequ > Degree) MinDegRequ = Degree;
349     binomial_array.Allocate (DeRequest1);
350 
351     for (ii = 0 ; ii < DeRequest1 ; ii++) {
352       binomial_array[ii] = 1.0e0 ;
353     }
354     if (!All) {
355       Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
356       derivative_storage.Allocate (DimDeRequ1);
357       RationalArray = derivative_storage ;
358     }
359 
360     Inverse = 1.0e0 / PolesArray[3]  ;
361     Index = 0 ;
362     Index2 = - 6;
363     OtherIndex = 0 ;
364 
365     for (ii = 0 ; ii <= MinDegRequ ; ii++) {
366       Index2 += 3;
367       Index1  = Index2;
368       RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
369       RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
370       RationalArray[Index] = PolesArray[OtherIndex];
371       Index      -= 2;
372       OtherIndex += 2;
373 
374       for (jj = ii - 1 ; jj >= 0 ; jj--) {
375 	Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
376 	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
377 	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
378 	RationalArray[Index] -=  Factor * RationalArray[Index1];
379 	Index  -= 2;
380 	Index1 -= 5;
381       }
382 
383       for (jj = ii ; jj >=  1 ; jj--) {
384 	binomial_array[jj] += binomial_array[jj - 1] ;
385       }
386       RationalArray[Index] *= Inverse ; Index++;
387       RationalArray[Index] *= Inverse ; Index++;
388       RationalArray[Index] *= Inverse ; Index++;
389     }
390 
391     for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
392       Index2 += 3;
393       Index1 = Index2;
394       RationalArray[Index] = 0.0e0; Index++;
395       RationalArray[Index] = 0.0e0; Index++;
396       RationalArray[Index] = 0.0e0;
397       Index -= 2;
398 
399       for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
400 	Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
401 	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
402 	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
403 	RationalArray[Index] -=  Factor * RationalArray[Index1];
404 	Index  -= 2;
405 	Index1 -= 5;
406       }
407 
408       for (jj = ii ; jj >=  1 ; jj--) {
409 	binomial_array[jj] += binomial_array[jj - 1] ;
410       }
411       RationalArray[Index] *= Inverse; Index++;
412       RationalArray[Index] *= Inverse; Index++;
413       RationalArray[Index] *= Inverse; Index++;
414     }
415 
416     if (!All) {
417       RationalArray = &RDers ;
418       Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
419       RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
420       RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
421       RationalArray[2] = derivative_storage[DimDeRequ];
422     }
423   }
424   else {
425     Standard_Integer kk;
426     Standard_Integer Dimension1 = Dimension + 1;
427     Standard_Integer Dimension2 = Dimension << 1;
428     Standard_Integer DeRequest1 = DerivativeRequest + 1;
429     Standard_Integer MinDegRequ = DerivativeRequest;
430     if (MinDegRequ > Degree) MinDegRequ = Degree;
431     binomial_array.Allocate (DeRequest1);
432 
433     for (ii = 0 ; ii < DeRequest1 ; ii++) {
434       binomial_array[ii] = 1.0e0 ;
435     }
436     if (!All) {
437       Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
438       derivative_storage.Allocate (DimDeRequ1);
439       RationalArray = derivative_storage ;
440     }
441 
442     Inverse = 1.0e0 / PolesArray[Dimension]  ;
443     Index = 0 ;
444     Index2 = - Dimension2;
445     OtherIndex = 0 ;
446 
447     for (ii = 0 ; ii <= MinDegRequ ; ii++) {
448       Index2 += Dimension;
449       Index1  = Index2;
450 
451       for (kk = 0 ; kk < Dimension ; kk++) {
452 	RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
453       }
454       Index -= Dimension;
455       ++OtherIndex;
456 
457       for (jj = ii - 1 ; jj >= 0 ; jj--) {
458 	Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
459 
460 	for (kk = 0 ; kk < Dimension ; kk++) {
461 	  RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
462 	}
463 	Index  -= Dimension ;
464 	Index1 -= Dimension2 ;
465       }
466 
467       for (jj = ii ; jj >=  1 ; jj--) {
468 	binomial_array[jj] += binomial_array[jj - 1] ;
469       }
470 
471       for (kk = 0 ; kk < Dimension ; kk++) {
472 	RationalArray[Index] *= Inverse ; Index++;
473       }
474     }
475 
476     for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
477       Index2 += Dimension;
478       Index1  = Index2;
479 
480       for (kk = 0 ; kk < Dimension ; kk++) {
481 	RationalArray[Index] = 0.0e0 ; Index++;
482       }
483       Index -= Dimension;
484 
485       for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
486 	Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
487 
488 	for (kk = 0 ; kk < Dimension ; kk++) {
489 	  RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
490 	}
491 	Index  -= Dimension ;
492 	Index1 -= Dimension2 ;
493       }
494 
495       for (jj = ii ; jj >=  1 ; jj--) {
496 	binomial_array[jj] += binomial_array[jj - 1] ;
497       }
498 
499       for (kk = 0 ; kk < Dimension ; kk++) {
500 	RationalArray[Index] *= Inverse; Index++;
501       }
502     }
503 
504     if (!All) {
505       RationalArray = &RDers ;
506       Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
507 
508       for (kk = 0 ; kk < Dimension ; kk++) {
509 	RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
510       }
511     }
512   }
513 }
514 
515 //=======================================================================
516 //function : RationalDerivatives
517 //purpose  : Uses Homogeneous Poles Derivatives and Deivatives of Weights
518 //=======================================================================
519 
RationalDerivatives(const Standard_Integer DerivativeRequest,const Standard_Integer Dimension,Standard_Real & PolesDerivates,Standard_Real & WeightsDerivates,Standard_Real & RationalDerivates)520 void  PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
521 				const Standard_Integer  Dimension,
522 				Standard_Real& PolesDerivates,
523 				// must be an array with
524 				// (DerivativeRequest + 1) * Dimension slots
525 				Standard_Real& WeightsDerivates,
526 				// must be an array with
527 				// (DerivativeRequest + 1) slots
528 				Standard_Real& RationalDerivates)
529 {
530   //
531   // Our purpose is to compute f = (u/v) derivated N times
532   //
533   //  We Write  u = fv
534   //  Let C(N,P) be the binomial
535   //
536   //        then we have
537   //
538   //         (q)                             (p)   (q-p)
539   //        u    =     SUM          C (q,p) f     v
540   //                p = 0 to q
541   //
542   //
543   //        Therefore
544   //
545   //
546   //         (q)         (   (q)                               (p)   (q-p)   )
547   //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
548   //                     (          p = 0 to q-1                             )
549   //
550   //
551   // make arrays for the binomial since computing it each time could
552   // raize a performance issue
553   //
554   Standard_Real Inverse;
555   Standard_Real *PolesArray = &PolesDerivates;
556   Standard_Real *WeightsArray = &WeightsDerivates;
557   Standard_Real *RationalArray = &RationalDerivates;
558   Standard_Real Factor ;
559 
560   Standard_Integer ii, Index, Index1, Index2, jj;
561   Standard_Integer DeRequest1 = DerivativeRequest + 1;
562 
563   NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
564   NCollection_LocalArray<Standard_Real> derivative_storage;
565 
566   for (ii = 0 ; ii < DeRequest1 ; ii++) {
567     binomial_array[ii] = 1.0e0 ;
568   }
569   Inverse = 1.0e0 / WeightsArray[0]  ;
570   if (Dimension == 3) {
571     Index  = 0 ;
572     Index2 = - 6 ;
573 
574     for (ii = 0 ; ii < DeRequest1 ; ii++) {
575       Index2 += 3;
576       Index1 = Index2;
577       RationalArray[Index] = PolesArray[Index] ; Index++;
578       RationalArray[Index] = PolesArray[Index] ; Index++;
579       RationalArray[Index] = PolesArray[Index] ;
580       Index -= 2;
581 
582       for (jj = ii - 1 ; jj >= 0 ; jj--) {
583 	Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
584 	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
585 	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
586 	RationalArray[Index] -=  Factor * RationalArray[Index1];
587 	Index  -= 2;
588 	Index1 -= 5;
589       }
590 
591       for (jj = ii ; jj >=  1 ; jj--) {
592 	binomial_array[jj] += binomial_array[jj - 1] ;
593       }
594       RationalArray[Index] *= Inverse ; Index++;
595       RationalArray[Index] *= Inverse ; Index++;
596       RationalArray[Index] *= Inverse ; Index++;
597     }
598   }
599   else {
600     Standard_Integer kk;
601     Standard_Integer Dimension2 = Dimension << 1;
602     Index  = 0 ;
603     Index2 = - Dimension2;
604 
605     for (ii = 0 ; ii < DeRequest1 ; ii++) {
606       Index2 += Dimension;
607       Index1  = Index2;
608 
609       for (kk = 0 ; kk < Dimension ; kk++) {
610 	RationalArray[Index] = PolesArray[Index]; Index++;
611       }
612       Index  -= Dimension;
613 
614       for (jj = ii - 1 ; jj >= 0 ; jj--) {
615 	Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
616 
617 	for (kk = 0 ; kk < Dimension ; kk++) {
618 	  RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
619 	}
620 	Index  -= Dimension;
621 	Index1 -= Dimension2;
622       }
623 
624       for (jj = ii ; jj >=  1 ; jj--) {
625 	binomial_array[jj] += binomial_array[jj - 1] ;
626       }
627 
628       for (kk = 0 ; kk < Dimension ; kk++) {
629 	RationalArray[Index] *= Inverse ; Index++;
630       }
631     }
632   }
633 }
634 
635 //=======================================================================
636 // Auxiliary template functions used for optimized evaluation of polynome
637 // and its derivatives for smaller dimensions of the polynome
638 //=======================================================================
639 
640 namespace {
641   // recursive template for evaluating value or first derivative
642   template<int dim>
eval_step1(double * poly,double par,double * coef)643   inline void eval_step1 (double* poly, double par, double* coef)
644   {
645     eval_step1<dim - 1> (poly, par, coef);
646     poly[dim] = poly[dim] * par + coef[dim];
647   }
648 
649   // recursion end
650   template<>
eval_step1(double *,double,double *)651   inline void eval_step1<-1> (double*, double, double*)
652   {
653   }
654 
655   // recursive template for evaluating second derivative
656   template<int dim>
eval_step2(double * poly,double par,double * coef)657   inline void eval_step2 (double* poly, double par, double* coef)
658   {
659     eval_step2<dim - 1> (poly, par, coef);
660     poly[dim] = poly[dim] * par + coef[dim] * 2.;
661   }
662 
663   // recursion end
664   template<>
eval_step2(double *,double,double *)665   inline void eval_step2<-1> (double*, double, double*)
666   {
667   }
668 
669   // evaluation of only value
670   template<int dim>
eval_poly0(double * aRes,double * aCoeffs,int Degree,double Par)671   inline void eval_poly0 (double* aRes, double* aCoeffs, int Degree, double Par)
672   {
673     Standard_Real* aRes0 = aRes;
674     memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
675 
676     for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
677     {
678       aCoeffs -= dim;
679       // Calculating the value of the polynomial
680       eval_step1<dim-1> (aRes0, Par, aCoeffs);
681     }
682   }
683 
684   // evaluation of value and first derivative
685   template<int dim>
eval_poly1(double * aRes,double * aCoeffs,int Degree,double Par)686   inline void eval_poly1 (double* aRes, double* aCoeffs, int Degree, double Par)
687   {
688     Standard_Real* aRes0 = aRes;
689     Standard_Real* aRes1 = aRes + dim;
690 
691     memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
692     memset(aRes1, 0, sizeof(Standard_Real) * dim);
693 
694     for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
695     {
696       aCoeffs -= dim;
697       // Calculating derivatives of the polynomial
698       eval_step1<dim-1> (aRes1, Par, aRes0);
699       // Calculating the value of the polynomial
700       eval_step1<dim-1> (aRes0, Par, aCoeffs);
701     }
702   }
703 
704   // evaluation of value and first and second derivatives
705   template<int dim>
eval_poly2(double * aRes,double * aCoeffs,int Degree,double Par)706   inline void eval_poly2 (double* aRes, double* aCoeffs, int Degree, double Par)
707   {
708     Standard_Real* aRes0 = aRes;
709     Standard_Real* aRes1 = aRes + dim;
710     Standard_Real* aRes2 = aRes + 2 * dim;
711 
712     memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
713     memset(aRes1, 0, sizeof(Standard_Real) * dim);
714     memset(aRes2, 0, sizeof(Standard_Real) * dim);
715 
716     for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
717     {
718       aCoeffs -= dim;
719       // Calculating second derivatives of the polynomial
720       eval_step2<dim-1> (aRes2, Par, aRes1);
721       // Calculating derivatives of the polynomial
722       eval_step1<dim-1> (aRes1, Par, aRes0);
723       // Calculating the value of the polynomial
724       eval_step1<dim-1> (aRes0, Par, aCoeffs);
725     }
726   }
727 }
728 
729 //=======================================================================
730 //function : This evaluates a polynomial and its derivatives
731 //purpose  : up to the requested order
732 //=======================================================================
733 
EvalPolynomial(const Standard_Real Par,const Standard_Integer DerivativeRequest,const Standard_Integer Degree,const Standard_Integer Dimension,Standard_Real & PolynomialCoeff,Standard_Real & Results)734 void  PLib::EvalPolynomial(const Standard_Real    Par,
735                            const Standard_Integer DerivativeRequest,
736                            const Standard_Integer Degree,
737                            const Standard_Integer Dimension,
738                                  Standard_Real&   PolynomialCoeff,
739                                  Standard_Real&   Results)
740      //
741      // the polynomial coefficients are assumed to be stored as follows :
742      //                                                      0
743      // [0]                  [Dimension -1]                 X     coefficient
744      //                                                      1
745      // [Dimension]          [Dimension + Dimension -1]     X     coefficient
746      //                                                      2
747      // [2 * Dimension]      [2 * Dimension + Dimension-1]  X     coefficient
748      //
749      //   ...................................................
750      //
751      //
752      //                                                      d
753      // [d * Dimension]      [d * Dimension + Dimension-1]  X     coefficient
754      //
755      //  where d is the Degree
756      //
757 {
758   Standard_Real* aCoeffs = &PolynomialCoeff + Degree * Dimension;
759   Standard_Real* aRes = &Results;
760   Standard_Real* anOriginal;
761   Standard_Integer ind = 0;
762   switch (DerivativeRequest)
763   {
764   case 1:
765   {
766     switch (Dimension)
767     {
768     case 1:  eval_poly1<1>  (aRes, aCoeffs, Degree, Par); break;
769     case 2:  eval_poly1<2>  (aRes, aCoeffs, Degree, Par); break;
770     case 3:  eval_poly1<3>  (aRes, aCoeffs, Degree, Par); break;
771     case 4:  eval_poly1<4>  (aRes, aCoeffs, Degree, Par); break;
772     case 5:  eval_poly1<5>  (aRes, aCoeffs, Degree, Par); break;
773     case 6:  eval_poly1<6>  (aRes, aCoeffs, Degree, Par); break;
774     case 7:  eval_poly1<7>  (aRes, aCoeffs, Degree, Par); break;
775     case 8:  eval_poly1<8>  (aRes, aCoeffs, Degree, Par); break;
776     case 9:  eval_poly1<9>  (aRes, aCoeffs, Degree, Par); break;
777     case 10: eval_poly1<10> (aRes, aCoeffs, Degree, Par); break;
778     case 11: eval_poly1<11> (aRes, aCoeffs, Degree, Par); break;
779     case 12: eval_poly1<12> (aRes, aCoeffs, Degree, Par); break;
780     case 13: eval_poly1<13> (aRes, aCoeffs, Degree, Par); break;
781     case 14: eval_poly1<14> (aRes, aCoeffs, Degree, Par); break;
782     case 15: eval_poly1<15> (aRes, aCoeffs, Degree, Par); break;
783     default:
784       {
785         Standard_Real* aRes0 = aRes;
786         Standard_Real* aRes1 = aRes + Dimension;
787 
788         memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * Dimension);
789         memset(aRes1, 0, sizeof(Standard_Real) * Dimension);
790 
791         for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
792         {
793           aCoeffs -= Dimension;
794           // Calculating derivatives of the polynomial
795           for (ind = 0; ind < Dimension; ind++)
796             aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
797           // Calculating the value of the polynomial
798           for (ind = 0; ind < Dimension; ind++)
799             aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
800         }
801       }
802     }
803     break;
804   }
805   case 2:
806   {
807     switch (Dimension)
808     {
809     case 1:  eval_poly2<1>  (aRes, aCoeffs, Degree, Par); break;
810     case 2:  eval_poly2<2>  (aRes, aCoeffs, Degree, Par); break;
811     case 3:  eval_poly2<3>  (aRes, aCoeffs, Degree, Par); break;
812     case 4:  eval_poly2<4>  (aRes, aCoeffs, Degree, Par); break;
813     case 5:  eval_poly2<5>  (aRes, aCoeffs, Degree, Par); break;
814     case 6:  eval_poly2<6>  (aRes, aCoeffs, Degree, Par); break;
815     case 7:  eval_poly2<7>  (aRes, aCoeffs, Degree, Par); break;
816     case 8:  eval_poly2<8>  (aRes, aCoeffs, Degree, Par); break;
817     case 9:  eval_poly2<9>  (aRes, aCoeffs, Degree, Par); break;
818     case 10: eval_poly2<10> (aRes, aCoeffs, Degree, Par); break;
819     case 11: eval_poly2<11> (aRes, aCoeffs, Degree, Par); break;
820     case 12: eval_poly2<12> (aRes, aCoeffs, Degree, Par); break;
821     case 13: eval_poly2<13> (aRes, aCoeffs, Degree, Par); break;
822     case 14: eval_poly2<14> (aRes, aCoeffs, Degree, Par); break;
823     case 15: eval_poly2<15> (aRes, aCoeffs, Degree, Par); break;
824     default:
825       {
826         Standard_Real* aRes0 = aRes;
827         Standard_Real* aRes1 = aRes + Dimension;
828         Standard_Real* aRes2 = aRes1 + Dimension;
829 
830         // Nullify the results
831         Standard_Integer aSize = 2 * Dimension;
832         memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
833         memset(aRes1, 0, sizeof(Standard_Real) * aSize);
834 
835         for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
836         {
837           aCoeffs -= Dimension;
838           // Calculating derivatives of the polynomial
839           for (ind = 0; ind < Dimension; ind++)
840             aRes2[ind] = aRes2[ind] * Par + aRes1[ind] * 2.0;
841           for (ind = 0; ind < Dimension; ind++)
842             aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
843           // Calculating the value of the polynomial
844           for (ind = 0; ind < Dimension; ind++)
845             aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
846         }
847         break;
848       }
849     }
850     break;
851   }
852   default:
853     {
854       // Nullify the results
855       Standard_Integer aResSize = (1 + DerivativeRequest) * Dimension;
856       memset(aRes, 0, sizeof(Standard_Real) * aResSize);
857 
858       for (Standard_Integer aDeg = 0; aDeg <= Degree; aDeg++)
859       {
860         aRes = &Results + aResSize - Dimension;
861         // Calculating derivatives of the polynomial
862         for (Standard_Integer aDeriv = DerivativeRequest; aDeriv > 0; aDeriv--)
863         {
864           anOriginal = aRes - Dimension; // pointer to the derivative minus 1
865           for (ind = 0; ind < Dimension; ind++)
866             aRes[ind] = aRes[ind] * Par + anOriginal[ind] * aDeriv;
867           aRes = anOriginal;
868         }
869         // Calculating the value of the polynomial
870         for (ind = 0; ind < Dimension; ind++)
871           aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
872         aCoeffs -= Dimension;
873       }
874     }
875   }
876 }
877 
878 //=======================================================================
879 //function : This evaluates a polynomial without derivative
880 //purpose  :
881 //=======================================================================
882 
NoDerivativeEvalPolynomial(const Standard_Real Par,const Standard_Integer Degree,const Standard_Integer Dimension,const Standard_Integer DegreeDimension,Standard_Real & PolynomialCoeff,Standard_Real & Results)883 void  PLib::NoDerivativeEvalPolynomial(const Standard_Real    Par,
884                                        const Standard_Integer Degree,
885                                        const Standard_Integer Dimension,
886                                        const Standard_Integer DegreeDimension,
887                                        Standard_Real&         PolynomialCoeff,
888                                        Standard_Real&         Results)
889 {
890   Standard_Real* aCoeffs = &PolynomialCoeff + DegreeDimension;
891   Standard_Real* aRes = &Results;
892 
893   switch (Dimension)
894   {
895   case 1:  eval_poly0<1>  (aRes, aCoeffs, Degree, Par); break;
896   case 2:  eval_poly0<2>  (aRes, aCoeffs, Degree, Par); break;
897   case 3:  eval_poly0<3>  (aRes, aCoeffs, Degree, Par); break;
898   case 4:  eval_poly0<4>  (aRes, aCoeffs, Degree, Par); break;
899   case 5:  eval_poly0<5>  (aRes, aCoeffs, Degree, Par); break;
900   case 6:  eval_poly0<6>  (aRes, aCoeffs, Degree, Par); break;
901   case 7:  eval_poly0<7>  (aRes, aCoeffs, Degree, Par); break;
902   case 8:  eval_poly0<8>  (aRes, aCoeffs, Degree, Par); break;
903   case 9:  eval_poly0<9>  (aRes, aCoeffs, Degree, Par); break;
904   case 10: eval_poly0<10> (aRes, aCoeffs, Degree, Par); break;
905   case 11: eval_poly0<11> (aRes, aCoeffs, Degree, Par); break;
906   case 12: eval_poly0<12> (aRes, aCoeffs, Degree, Par); break;
907   case 13: eval_poly0<13> (aRes, aCoeffs, Degree, Par); break;
908   case 14: eval_poly0<14> (aRes, aCoeffs, Degree, Par); break;
909   case 15: eval_poly0<15> (aRes, aCoeffs, Degree, Par); break;
910   default:
911     {
912       memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
913       for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
914       {
915         aCoeffs -= Dimension;
916         for (Standard_Integer ind = 0; ind < Dimension; ind++)
917           aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
918       }
919     }
920   }
921 }
922 
923 //=======================================================================
924 //function : This evaluates a polynomial of 2 variables
925 //purpose  : or its derivative at the requested orders
926 //=======================================================================
927 
EvalPoly2Var(const Standard_Real UParameter,const Standard_Real VParameter,const Standard_Integer UDerivativeRequest,const Standard_Integer VDerivativeRequest,const Standard_Integer UDegree,const Standard_Integer VDegree,const Standard_Integer Dimension,Standard_Real & PolynomialCoeff,Standard_Real & Results)928 void  PLib::EvalPoly2Var(const Standard_Real    UParameter,
929 			 const Standard_Real    VParameter,
930 			 const Standard_Integer UDerivativeRequest,
931 			 const Standard_Integer VDerivativeRequest,
932 			 const Standard_Integer UDegree,
933 			 const Standard_Integer VDegree,
934 			 const Standard_Integer Dimension,
935 			 Standard_Real&         PolynomialCoeff,
936 			 Standard_Real&         Results)
937      //
938      // the polynomial coefficients are assumed to be stored as follows :
939      //                                                      0 0
940      // [0]                  [Dimension -1]                 U V    coefficient
941      //                                                      1 0
942      // [Dimension]          [Dimension + Dimension -1]     U V    coefficient
943      //                                                      2 0
944      // [2 * Dimension]      [2 * Dimension + Dimension-1]  U V    coefficient
945      //
946      //   ...................................................
947      //
948      //
949      //                                                      m 0
950      // [m * Dimension]      [m * Dimension + Dimension-1]  U V    coefficient
951      //
952      //  where m = UDegree
953      //
954      //                                                           0 1
955      // [(m+1) * Dimension]  [(m+1) * Dimension + Dimension-1]   U V   coefficient
956      //
957      //   ...................................................
958      //
959      //                                                          m 1
960      // [2*m * Dimension]      [2*m * Dimension + Dimension-1]  U V    coefficient
961      //
962      //   ...................................................
963      //
964      //                                                          m n
965      // [m*n * Dimension]      [m*n * Dimension + Dimension-1]  U V    coefficient
966      //
967      //  where n = VDegree
968 {
969   Standard_Integer Udim = (VDegree+1)*Dimension,
970   index = Udim*UDerivativeRequest;
971   TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
972   TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
973   Standard_Real * Result =  (Standard_Real *) &Curve.ChangeValue(1);
974   Standard_Real * Digit  =  (Standard_Real *) &Point.ChangeValue(1);
975   Standard_Real * ResultArray ;
976   ResultArray = &Results ;
977 
978   PLib::EvalPolynomial(UParameter,
979 		       UDerivativeRequest,
980 		       UDegree,
981 		       Udim,
982 		       PolynomialCoeff,
983 		       Result[0]);
984 
985   PLib::EvalPolynomial(VParameter,
986 		       VDerivativeRequest,
987 		       VDegree,
988 		       Dimension,
989 		       Result[index],
990 		       Digit[0]);
991 
992   index = Dimension*VDerivativeRequest;
993 
994   for (Standard_Integer i=0;i<Dimension;i++) {
995     ResultArray[i] = Digit[index+i];
996   }
997 }
998 
999 
1000 
1001 //=======================================================================
1002 //function : This evaluates the lagrange polynomial and its derivatives
1003 //purpose  : up to the requested order that interpolates a series of
1004 //points of dimension <Dimension> with given assigned parameters
1005 //=======================================================================
1006 
1007 Standard_Integer
EvalLagrange(const Standard_Real Parameter,const Standard_Integer DerivativeRequest,const Standard_Integer Degree,const Standard_Integer Dimension,Standard_Real & Values,Standard_Real & Parameters,Standard_Real & Results)1008 PLib::EvalLagrange(const Standard_Real                   Parameter,
1009 		   const Standard_Integer                DerivativeRequest,
1010 		   const Standard_Integer                Degree,
1011 		   const Standard_Integer                Dimension,
1012 		   Standard_Real&                        Values,
1013 		   Standard_Real&                        Parameters,
1014 		   Standard_Real&                        Results)
1015 {
1016   //
1017   // the points  are assumed to be stored as follows in the Values array :
1018   //
1019   // [0]              [Dimension -1]                first  point   coefficients
1020   //
1021   // [Dimension]      [Dimension + Dimension -1]    second point   coefficients
1022   //
1023   // [2 * Dimension]  [2 * Dimension + Dimension-1] third  point   coefficients
1024   //
1025   //   ...................................................
1026   //
1027   //
1028   //
1029   // [d * Dimension]  [d * Dimension + Dimension-1] d + 1 point   coefficients
1030   //
1031   //  where d is the Degree
1032   //
1033   //  The ParameterArray stores the parameter value assign to each point in
1034   //  order described above, that is
1035   //  [0]   is assign to first  point
1036   //  [1]   is assign to second point
1037   //
1038   Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
1039   Standard_Integer local_request = DerivativeRequest;
1040   Standard_Real  *ParameterArray;
1041   Standard_Real  difference;
1042   Standard_Real  *PointsArray;
1043   Standard_Real  *ResultArray ;
1044 
1045   PointsArray    = &Values ;
1046   ParameterArray = &Parameters ;
1047   ResultArray = &Results ;
1048   if (local_request >= Degree) {
1049     local_request = Degree ;
1050   }
1051   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1052   //
1053   //  Build the divided differences array
1054   //
1055 
1056   for (ii = 0 ; ii <  (Degree + 1) * Dimension  ; ii++) {
1057     divided_differences_array[ii] = PointsArray[ii] ;
1058   }
1059 
1060   for (ii = Degree ; ii >= 0   ; ii--) {
1061 
1062     for (jj = Degree  ; jj > Degree - ii  ; jj--) {
1063       Index = jj * Dimension ;
1064       Index1 = Index - Dimension ;
1065 
1066       for (kk = 0 ; kk < Dimension ; kk++) {
1067 	divided_differences_array[Index + kk] -=
1068 	  divided_differences_array[Index1 + kk] ;
1069       }
1070       difference =
1071 	ParameterArray[jj] - ParameterArray[jj - Degree -1 +  ii] ;
1072       if (Abs(difference) < RealSmall()) {
1073 	ReturnCode = 1 ;
1074 	goto FINISH ;
1075       }
1076       difference = 1.0e0 / difference ;
1077 
1078       for (kk = 0 ; kk < Dimension ; kk++) {
1079 	divided_differences_array[Index + kk] *= difference ;
1080       }
1081     }
1082   }
1083   //
1084   //
1085   // Evaluate the divided difference array polynomial which expresses as
1086   //
1087   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1088   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1089   //
1090   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1091   //
1092   //
1093   Index = Degree * Dimension ;
1094 
1095   for (kk = 0 ; kk < Dimension ; kk++) {
1096     ResultArray[kk] = divided_differences_array[Index + kk] ;
1097   }
1098 
1099   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
1100     ResultArray[ii] = 0.0e0 ;
1101   }
1102 
1103   for (ii = Degree ; ii >= 1 ; ii--) {
1104     difference =  Parameter - ParameterArray[ii - 1] ;
1105 
1106     for (jj = local_request ; jj > 0 ; jj--) {
1107       Index = jj * Dimension ;
1108       Index1 = Index - Dimension ;
1109 
1110       for (kk = 0 ; kk < Dimension ; kk++) {
1111 	ResultArray[Index + kk] *= difference ;
1112 	ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
1113       }
1114     }
1115     Index = (ii -1) * Dimension ;
1116 
1117     for (kk = 0 ; kk < Dimension ; kk++) {
1118       ResultArray[kk] *= difference ;
1119       ResultArray[kk] += divided_differences_array[Index+kk] ;
1120     }
1121   }
1122   FINISH :
1123     return (ReturnCode) ;
1124 }
1125 
1126 //=======================================================================
1127 //function : This evaluates the hermite polynomial and its derivatives
1128 //purpose  : up to the requested order that interpolates a series of
1129 //points of dimension <Dimension> with given assigned parameters
1130 //=======================================================================
1131 
EvalCubicHermite(const Standard_Real Parameter,const Standard_Integer DerivativeRequest,const Standard_Integer Dimension,Standard_Real & Values,Standard_Real & Derivatives,Standard_Real & theParameters,Standard_Real & Results)1132 Standard_Integer PLib::EvalCubicHermite
1133 (const Standard_Real          Parameter,
1134  const Standard_Integer       DerivativeRequest,
1135  const Standard_Integer       Dimension,
1136  Standard_Real&               Values,
1137  Standard_Real&               Derivatives,
1138  Standard_Real&               theParameters,
1139  Standard_Real&               Results)
1140 {
1141   //
1142   // the points  are assumed to be stored as follows in the Values array :
1143   //
1144   // [0]            [Dimension -1]             first  point   coefficients
1145   //
1146   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
1147   //
1148   //
1149   // the derivatives  are assumed to be stored as follows in
1150   // the Derivatives array :
1151   //
1152   // [0]            [Dimension -1]             first  point   coefficients
1153   //
1154   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
1155   //
1156   //  The ParameterArray stores the parameter value assign to each point in
1157   //  order described above, that is
1158   //  [0]   is assign to first  point
1159   //  [1]   is assign to last   point
1160   //
1161   Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
1162   Standard_Integer local_request = DerivativeRequest ;
1163 
1164   ReturnCode = 0 ;
1165   Degree = 3 ;
1166   Standard_Real  ParametersArray[4];
1167   Standard_Real  difference;
1168   Standard_Real  inverse;
1169   Standard_Real  *FirstLast;
1170   Standard_Real  *PointsArray;
1171   Standard_Real  *DerivativesArray;
1172   Standard_Real  *ResultArray ;
1173 
1174   DerivativesArray = &Derivatives ;
1175   PointsArray    = &Values ;
1176   FirstLast = &theParameters ;
1177   ResultArray = &Results ;
1178   if (local_request >= Degree) {
1179     local_request = Degree ;
1180   }
1181   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1182 
1183   for (ii = 0, jj = 0  ; ii < 2 ; ii++, jj+= 2) {
1184     ParametersArray[jj] =
1185       ParametersArray[jj+1] = FirstLast[ii] ;
1186   }
1187   //
1188   //  Build the divided differences array
1189   //
1190   //
1191   // initialise it at the stage 2 of the building algorithm
1192   // for divided differences
1193   //
1194   inverse = FirstLast[1] - FirstLast[0] ;
1195   inverse = 1.0e0 / inverse ;
1196 
1197   for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension  ;
1198        ii <  Dimension  ;
1199        ii++, jj++, kk++, pp++) {
1200     divided_differences_array[ii] = PointsArray[ii] ;
1201     divided_differences_array[kk] = inverse *
1202       (PointsArray[jj] - PointsArray[ii]) ;
1203     divided_differences_array[jj] = DerivativesArray[ii] ;
1204     divided_differences_array[pp] = DerivativesArray[jj] ;
1205   }
1206 
1207   for (ii = 1 ; ii <= Degree   ; ii++) {
1208 
1209     for (jj = Degree  ; jj >=  ii+1  ; jj--) {
1210       Index = jj * Dimension ;
1211       Index1 = Index - Dimension ;
1212 
1213       for (kk = 0 ; kk < Dimension ; kk++) {
1214 	divided_differences_array[Index + kk] -=
1215 	  divided_differences_array[Index1 + kk] ;
1216       }
1217 
1218       for (kk = 0 ; kk < Dimension ; kk++) {
1219 	divided_differences_array[Index + kk] *= inverse ;
1220       }
1221     }
1222   }
1223   //
1224   //
1225   // Evaluate the divided difference array polynomial which expresses as
1226   //
1227   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1228   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1229   //
1230   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1231   //
1232   //
1233   Index = Degree * Dimension ;
1234 
1235   for (kk = 0 ; kk < Dimension ; kk++) {
1236     ResultArray[kk] = divided_differences_array[Index + kk] ;
1237   }
1238 
1239   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
1240     ResultArray[ii] = 0.0e0 ;
1241   }
1242 
1243   for (ii = Degree ; ii >= 1 ; ii--) {
1244     difference =  Parameter - ParametersArray[ii - 1] ;
1245 
1246     for (jj = local_request ; jj > 0 ; jj--) {
1247       Index = jj * Dimension ;
1248       Index1 = Index - Dimension ;
1249 
1250       for (kk = 0 ; kk < Dimension ; kk++) {
1251 	ResultArray[Index + kk] *= difference ;
1252 	ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
1253       }
1254     }
1255     Index = (ii -1) * Dimension ;
1256 
1257     for (kk = 0 ; kk < Dimension ; kk++) {
1258       ResultArray[kk] *= difference ;
1259       ResultArray[kk] += divided_differences_array[Index+kk] ;
1260     }
1261   }
1262 //  FINISH :
1263     return (ReturnCode) ;
1264 }
1265 
1266 //=======================================================================
1267 //function : HermiteCoefficients
1268 //purpose  : calcul des polynomes d'Hermite
1269 //=======================================================================
1270 
HermiteCoefficients(const Standard_Real FirstParameter,const Standard_Real LastParameter,const Standard_Integer FirstOrder,const Standard_Integer LastOrder,math_Matrix & MatrixCoefs)1271 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
1272 					   const Standard_Real LastParameter,
1273 					   const Standard_Integer FirstOrder,
1274 					   const Standard_Integer LastOrder,
1275 					   math_Matrix& MatrixCoefs)
1276 {
1277   Standard_Integer NbCoeff =  FirstOrder +  LastOrder + 2, Ordre[2];
1278   Standard_Integer ii, jj, pp, cote, iof=0;
1279   Standard_Real Prod, TBorne = FirstParameter;
1280   math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
1281   math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
1282 
1283   // Test de validites
1284 
1285   if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
1286   Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
1287   if (D1 > 100 || D2 > 100) return Standard_False;
1288   D2 += D1;
1289   if (D2 < 0.01) return Standard_False;
1290   if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
1291 
1292   // Calcul de la matrice a inverser (MAT)
1293 
1294   Ordre[0] = FirstOrder+1;
1295   Ordre[1] = LastOrder+1;
1296 
1297   for (cote=0; cote<=1; cote++) {
1298     Coeff.Init(1);
1299 
1300     for (pp=1; pp<=Ordre[cote]; pp++) {
1301       ii = pp + iof;
1302       Prod = 1;
1303 
1304       for (jj=pp; jj<=NbCoeff; jj++) {
1305 	//        tout se passe dans les 3 lignes suivantes
1306 	MAT(ii, jj) = Coeff(jj) * Prod;
1307 	Coeff(jj) *= jj - pp;
1308 	Prod      *= TBorne;
1309       }
1310     }
1311     TBorne = LastParameter;
1312     iof = Ordre[0];
1313   }
1314 
1315   // resolution du systemes
1316   math_Gauss ResolCoeff(MAT, 1.0e-10);
1317   if (!ResolCoeff.IsDone()) return Standard_False;
1318 
1319   for (ii=1; ii<=NbCoeff; ii++) {
1320     B(ii) = 1;
1321     ResolCoeff.Solve(B, Coeff);
1322     MatrixCoefs.SetRow( ii, Coeff);
1323     B(ii) = 0;
1324   }
1325   return Standard_True;
1326 }
1327 
1328 //=======================================================================
1329 //function : CoefficientsPoles
1330 //purpose  :
1331 //=======================================================================
1332 
CoefficientsPoles(const TColgp_Array1OfPnt & Coefs,const TColStd_Array1OfReal * WCoefs,TColgp_Array1OfPnt & Poles,TColStd_Array1OfReal * Weights)1333 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt&   Coefs,
1334 			      const TColStd_Array1OfReal* WCoefs,
1335 			      TColgp_Array1OfPnt&         Poles,
1336 			      TColStd_Array1OfReal*       Weights)
1337 {
1338   TColStd_Array1OfReal tempC(1,3*Coefs.Length());
1339   PLib::SetPoles(Coefs,tempC);
1340   TColStd_Array1OfReal tempP(1,3*Poles.Length());
1341   PLib::SetPoles(Coefs,tempP);
1342   PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
1343   PLib::GetPoles(tempP,Poles);
1344 }
1345 
1346 //=======================================================================
1347 //function : CoefficientsPoles
1348 //purpose  :
1349 //=======================================================================
1350 
CoefficientsPoles(const TColgp_Array1OfPnt2d & Coefs,const TColStd_Array1OfReal * WCoefs,TColgp_Array1OfPnt2d & Poles,TColStd_Array1OfReal * Weights)1351 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
1352 			      const TColStd_Array1OfReal* WCoefs,
1353 			      TColgp_Array1OfPnt2d&       Poles,
1354 			      TColStd_Array1OfReal*       Weights)
1355 {
1356   TColStd_Array1OfReal tempC(1,2*Coefs.Length());
1357   PLib::SetPoles(Coefs,tempC);
1358   TColStd_Array1OfReal tempP(1,2*Poles.Length());
1359   PLib::SetPoles(Coefs,tempP);
1360   PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
1361   PLib::GetPoles(tempP,Poles);
1362 }
1363 
1364 //=======================================================================
1365 //function : CoefficientsPoles
1366 //purpose  :
1367 //=======================================================================
1368 
CoefficientsPoles(const TColStd_Array1OfReal & Coefs,const TColStd_Array1OfReal * WCoefs,TColStd_Array1OfReal & Poles,TColStd_Array1OfReal * Weights)1369 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
1370 			      const TColStd_Array1OfReal* WCoefs,
1371 			      TColStd_Array1OfReal&       Poles,
1372 			      TColStd_Array1OfReal*       Weights)
1373 {
1374   PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
1375 }
1376 
1377 //=======================================================================
1378 //function : CoefficientsPoles
1379 //purpose  :
1380 //=======================================================================
1381 
CoefficientsPoles(const Standard_Integer dim,const TColStd_Array1OfReal & Coefs,const TColStd_Array1OfReal * WCoefs,TColStd_Array1OfReal & Poles,TColStd_Array1OfReal * Weights)1382 void PLib::CoefficientsPoles (const Standard_Integer      dim,
1383 			      const TColStd_Array1OfReal& Coefs,
1384 			      const TColStd_Array1OfReal* WCoefs,
1385 			      TColStd_Array1OfReal&       Poles,
1386 			      TColStd_Array1OfReal*       Weights)
1387 {
1388   Standard_Boolean rat = WCoefs != NULL;
1389   Standard_Integer loc = Coefs.Lower();
1390   Standard_Integer lop = Poles.Lower();
1391   Standard_Integer lowc=0;
1392   Standard_Integer lowp=0;
1393   Standard_Integer upc = Coefs.Upper();
1394   Standard_Integer upp = Poles.Upper();
1395   Standard_Integer upwc=0;
1396   Standard_Integer upwp=0;
1397   Standard_Integer reflen = Coefs.Length()/dim;
1398   Standard_Integer i,j,k;
1399   //Les Extremites.
1400   if (rat) {
1401     lowc = WCoefs->Lower(); lowp = Weights->Lower();
1402     upwc = WCoefs->Upper(); upwp = Weights->Upper();
1403   }
1404 
1405   for (i = 0; i < dim; i++){
1406     Poles (lop + i) = Coefs (loc + i);
1407     Poles (upp - i) = Coefs (upc - i);
1408   }
1409   if (rat) {
1410     (*Weights) (lowp) = (*WCoefs) (lowc);
1411     (*Weights) (upwp) = (*WCoefs) (upwc);
1412   }
1413 
1414   Standard_Real Cnp;
1415   for (i = 2; i < reflen; i++ ) {
1416     Cnp = PLib::Bin(reflen - 1, i - 1);
1417     if (rat) (*Weights)(lowp + i - 1) = (*WCoefs)(lowc + i - 1) / Cnp;
1418 
1419     for(j = 0; j < dim; j++){
1420       Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
1421     }
1422   }
1423 
1424   for (i = 1; i <= reflen - 1; i++) {
1425 
1426     for (j = reflen - 1; j >= i; j--) {
1427       if (rat) (*Weights)(lowp + j) += (*Weights)(lowp + j - 1);
1428 
1429       for(k = 0; k < dim; k++){
1430 	Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
1431       }
1432     }
1433   }
1434   if (rat) {
1435 
1436     for (i = 1; i <= reflen; i++) {
1437 
1438       for(j = 0; j < dim; j++){
1439 	Poles(lop + dim * (i-1) + j) /= (*Weights)(lowp + i -1);
1440       }
1441     }
1442   }
1443 }
1444 
1445 //=======================================================================
1446 //function : Trimming
1447 //purpose  :
1448 //=======================================================================
1449 
Trimming(const Standard_Real U1,const Standard_Real U2,TColgp_Array1OfPnt & Coefs,TColStd_Array1OfReal * WCoefs)1450 void PLib::Trimming(const Standard_Real   U1,
1451 		    const Standard_Real   U2,
1452 		    TColgp_Array1OfPnt&   Coefs,
1453 		    TColStd_Array1OfReal* WCoefs)
1454 {
1455   TColStd_Array1OfReal temp(1,3*Coefs.Length());
1456   PLib::SetPoles(Coefs,temp);
1457   PLib::Trimming(U1,U2,3,temp,WCoefs);
1458   PLib::GetPoles(temp,Coefs);
1459 }
1460 
1461 //=======================================================================
1462 //function : Trimming
1463 //purpose  :
1464 //=======================================================================
1465 
Trimming(const Standard_Real U1,const Standard_Real U2,TColgp_Array1OfPnt2d & Coefs,TColStd_Array1OfReal * WCoefs)1466 void PLib::Trimming(const Standard_Real   U1,
1467 		    const Standard_Real   U2,
1468 		    TColgp_Array1OfPnt2d& Coefs,
1469 		    TColStd_Array1OfReal* WCoefs)
1470 {
1471   TColStd_Array1OfReal temp(1,2*Coefs.Length());
1472   PLib::SetPoles(Coefs,temp);
1473   PLib::Trimming(U1,U2,2,temp,WCoefs);
1474   PLib::GetPoles(temp,Coefs);
1475 }
1476 
1477 //=======================================================================
1478 //function : Trimming
1479 //purpose  :
1480 //=======================================================================
1481 
Trimming(const Standard_Real U1,const Standard_Real U2,TColStd_Array1OfReal & Coefs,TColStd_Array1OfReal * WCoefs)1482 void PLib::Trimming(const Standard_Real   U1,
1483 		    const Standard_Real   U2,
1484 		    TColStd_Array1OfReal& Coefs,
1485 		    TColStd_Array1OfReal* WCoefs)
1486 {
1487   PLib::Trimming(U1,U2,1,Coefs,WCoefs);
1488 }
1489 
1490 //=======================================================================
1491 //function : Trimming
1492 //purpose  :
1493 //=======================================================================
1494 
Trimming(const Standard_Real U1,const Standard_Real U2,const Standard_Integer dim,TColStd_Array1OfReal & Coefs,TColStd_Array1OfReal * WCoefs)1495 void PLib::Trimming(const Standard_Real U1,
1496 		    const Standard_Real U2,
1497 		    const Standard_Integer dim,
1498 		    TColStd_Array1OfReal& Coefs,
1499 		    TColStd_Array1OfReal* WCoefs)
1500 {
1501 
1502   // principe :
1503   // on fait le changement de variable v = (u-U1) / (U2-U1)
1504   // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
1505   // decomposee sous la forme du schema iteratif de horner.
1506 
1507   Standard_Real lsp = U2 - U1;
1508   Standard_Integer indc, indw=0;
1509   Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
1510   Standard_Integer len = Coefs.Length()/dim;
1511   Standard_Boolean rat = WCoefs != NULL;
1512 
1513   if (rat) {
1514     if(len != WCoefs->Length())
1515       throw Standard_Failure("PLib::Trimming : nbcoefs/dim != nbweights !!!");
1516     upw = WCoefs->Upper();
1517   }
1518   len --;
1519 
1520   for (Standard_Integer i = 1; i <= len; i++) {
1521     Standard_Integer j ;
1522     indc = upc - dim*(i-1);
1523     if (rat) indw = upw - i + 1;
1524     //calcul du coefficient de degre le plus faible a l'iteration i
1525 
1526     for( j = 0; j < dim; j++){
1527       Coefs(indc - dim + j) += U1 * Coefs(indc + j);
1528     }
1529     if (rat) (*WCoefs)(indw - 1) += U1 * (*WCoefs)(indw);
1530 
1531     //calcul des coefficients intermediaires :
1532 
1533     while (indc < upc){
1534       indc += dim;
1535 
1536       for(Standard_Integer k = 0; k < dim; k++){
1537 	Coefs(indc - dim + k) =
1538 	  U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
1539       }
1540       if (rat) {
1541 	indw ++;
1542         (*WCoefs)(indw - 1) = U1 * (*WCoefs)(indw) + lsp * (*WCoefs)(indw - 1);
1543       }
1544     }
1545 
1546     //calcul du coefficient de degre le plus eleve :
1547 
1548     for(j = 0; j < dim; j++){
1549       Coefs(upc + j) *= lsp;
1550     }
1551     if (rat) (*WCoefs)(upw) *= lsp;
1552   }
1553 }
1554 
1555 //=======================================================================
1556 //function : CoefficientsPoles
1557 //purpose  :
1558 // Modified: 21/10/1996 by PMN :  PolesCoefficient (PRO5852).
1559 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
1560 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
1561 //=======================================================================
1562 
CoefficientsPoles(const TColgp_Array2OfPnt & Coefs,const TColStd_Array2OfReal * WCoefs,TColgp_Array2OfPnt & Poles,TColStd_Array2OfReal * Weights)1563 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt&   Coefs,
1564 			      const TColStd_Array2OfReal* WCoefs,
1565 			      TColgp_Array2OfPnt&         Poles,
1566 			      TColStd_Array2OfReal*       Weights)
1567 {
1568   Standard_Boolean rat = (WCoefs != NULL);
1569   Standard_Integer LowerRow  = Poles.LowerRow();
1570   Standard_Integer UpperRow  = Poles.UpperRow();
1571   Standard_Integer LowerCol  = Poles.LowerCol();
1572   Standard_Integer UpperCol  = Poles.UpperCol();
1573   Standard_Integer ColLength = Poles.ColLength();
1574   Standard_Integer RowLength = Poles.RowLength();
1575 
1576   // Bidouille pour retablir u et v pour les coefs calcules
1577   // par buildcache
1578 //  Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
1579 
1580   Standard_Integer Row, Col;
1581   Standard_Real W, Cnp;
1582 
1583   Standard_Integer I1, I2;
1584   Standard_Integer NPoleu , NPolev;
1585   gp_XYZ Temp;
1586 
1587   for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
1588     Poles (NPoleu, LowerCol) =  Coefs (NPoleu, LowerCol);
1589     if (rat) {
1590       (*Weights) (NPoleu, LowerCol) =  (*WCoefs) (NPoleu, LowerCol);
1591     }
1592 
1593     for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
1594       Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
1595       Temp = Coefs (NPoleu, Col).XYZ();
1596       Temp.Divide (Cnp);
1597       Poles (NPoleu, Col).SetXYZ (Temp);
1598       if (rat) {
1599 	(*Weights) (NPoleu, Col) = (*WCoefs) (NPoleu, Col) /  Cnp;
1600       }
1601     }
1602     Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
1603     if (rat) {
1604       (*Weights) (NPoleu, UpperCol) = (*WCoefs) (NPoleu, UpperCol);
1605     }
1606 
1607     for (I1 = 1; I1 <= RowLength - 1; I1++) {
1608 
1609       for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
1610         Temp.SetLinearForm
1611 	  (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
1612 	Poles (NPoleu, I2).SetXYZ (Temp);
1613 	if (rat) (*Weights)(NPoleu, I2) += (*Weights)(NPoleu, I2-1);
1614       }
1615     }
1616   }
1617 
1618   for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
1619 
1620     for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
1621       Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
1622       Temp = Poles (Row, NPolev).XYZ();
1623       Temp.Divide (Cnp);
1624       Poles (Row, NPolev).SetXYZ (Temp);
1625       if (rat) (*Weights)(Row, NPolev) /= Cnp;
1626     }
1627 
1628     for (I1 = 1; I1 <= ColLength - 1; I1++) {
1629 
1630       for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
1631         Temp.SetLinearForm
1632 	  (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
1633 	Poles (I2, NPolev).SetXYZ (Temp);
1634 	if (rat) (*Weights)(I2, NPolev) += (*Weights)(I2-1, NPolev);
1635       }
1636     }
1637   }
1638   if (rat) {
1639 
1640     for (Row = LowerRow; Row <= UpperRow; Row++) {
1641 
1642       for (Col = LowerCol; Col <= UpperCol; Col++) {
1643 	W = (*Weights) (Row, Col);
1644 	Temp = Poles(Row, Col).XYZ();
1645 	Temp.Divide (W);
1646 	Poles(Row, Col).SetXYZ (Temp);
1647       }
1648     }
1649   }
1650 }
1651 
1652 //=======================================================================
1653 //function : UTrimming
1654 //purpose  :
1655 //=======================================================================
1656 
UTrimming(const Standard_Real U1,const Standard_Real U2,TColgp_Array2OfPnt & Coeffs,TColStd_Array2OfReal * WCoeffs)1657 void PLib::UTrimming(const Standard_Real U1,
1658 		     const Standard_Real U2,
1659 		     TColgp_Array2OfPnt& Coeffs,
1660 		     TColStd_Array2OfReal* WCoeffs)
1661 {
1662   Standard_Boolean rat = WCoeffs != NULL;
1663   Standard_Integer lr = Coeffs.LowerRow();
1664   Standard_Integer ur = Coeffs.UpperRow();
1665   Standard_Integer lc = Coeffs.LowerCol();
1666   Standard_Integer uc = Coeffs.UpperCol();
1667   TColgp_Array1OfPnt  Temp (lr,ur);
1668   TColStd_Array1OfReal Temw (lr,ur);
1669 
1670   for (Standard_Integer icol = lc; icol <= uc; icol++) {
1671     Standard_Integer irow ;
1672     for ( irow = lr; irow <= ur; irow++) {
1673       Temp (irow) = Coeffs  (irow, icol);
1674       if (rat) Temw (irow) = (*WCoeffs) (irow, icol);
1675     }
1676     if (rat) PLib::Trimming (U1, U2, Temp, &Temw);
1677     else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
1678 
1679     for (irow = lr; irow <= ur; irow++) {
1680       Coeffs  (irow, icol) = Temp (irow);
1681       if (rat) (*WCoeffs) (irow, icol) = Temw (irow);
1682     }
1683   }
1684 }
1685 
1686 //=======================================================================
1687 //function : VTrimming
1688 //purpose  :
1689 //=======================================================================
1690 
VTrimming(const Standard_Real V1,const Standard_Real V2,TColgp_Array2OfPnt & Coeffs,TColStd_Array2OfReal * WCoeffs)1691 void PLib::VTrimming(const Standard_Real V1,
1692 		     const Standard_Real V2,
1693 		     TColgp_Array2OfPnt& Coeffs,
1694 		     TColStd_Array2OfReal* WCoeffs)
1695 {
1696   Standard_Boolean rat = WCoeffs != NULL;
1697   Standard_Integer lr = Coeffs.LowerRow();
1698   Standard_Integer ur = Coeffs.UpperRow();
1699   Standard_Integer lc = Coeffs.LowerCol();
1700   Standard_Integer uc = Coeffs.UpperCol();
1701   TColgp_Array1OfPnt  Temp (lc,uc);
1702   TColStd_Array1OfReal Temw (lc,uc);
1703 
1704   for (Standard_Integer irow = lr; irow <= ur; irow++) {
1705     Standard_Integer icol ;
1706     for ( icol = lc; icol <= uc; icol++) {
1707       Temp (icol) = Coeffs  (irow, icol);
1708       if (rat) Temw (icol) = (*WCoeffs) (irow, icol);
1709     }
1710     if (rat) PLib::Trimming (V1, V2, Temp, &Temw);
1711     else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
1712 
1713     for (icol = lc; icol <= uc; icol++) {
1714       Coeffs  (irow, icol) = Temp (icol);
1715       if (rat) (*WCoeffs) (irow, icol) = Temw (icol);
1716     }
1717   }
1718 }
1719 
1720 //=======================================================================
1721 //function : HermiteInterpolate
1722 //purpose  :
1723 //=======================================================================
1724 
HermiteInterpolate(const Standard_Integer Dimension,const Standard_Real FirstParameter,const Standard_Real LastParameter,const Standard_Integer FirstOrder,const Standard_Integer LastOrder,const TColStd_Array2OfReal & FirstConstr,const TColStd_Array2OfReal & LastConstr,TColStd_Array1OfReal & Coefficients)1725 Standard_Boolean PLib::HermiteInterpolate
1726 (const Standard_Integer Dimension,
1727  const Standard_Real FirstParameter,
1728  const Standard_Real LastParameter,
1729  const Standard_Integer FirstOrder,
1730  const Standard_Integer LastOrder,
1731  const TColStd_Array2OfReal& FirstConstr,
1732  const TColStd_Array2OfReal& LastConstr,
1733  TColStd_Array1OfReal& Coefficients)
1734 {
1735   Standard_Real Pattern[3][6];
1736 
1737   // portage HP : il faut les initialiser 1 par 1
1738 
1739   Pattern[0][0] = 1;
1740   Pattern[0][1] = 1;
1741   Pattern[0][2] = 1;
1742   Pattern[0][3] = 1;
1743   Pattern[0][4] = 1;
1744   Pattern[0][5] = 1;
1745   Pattern[1][0] = 0;
1746   Pattern[1][1] = 1;
1747   Pattern[1][2] = 2;
1748   Pattern[1][3] = 3;
1749   Pattern[1][4] = 4;
1750   Pattern[1][5] = 5;
1751   Pattern[2][0] = 0;
1752   Pattern[2][1] = 0;
1753   Pattern[2][2] = 2;
1754   Pattern[2][3] = 6;
1755   Pattern[2][4] = 12;
1756   Pattern[2][5] = 20;
1757 
1758   math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
1759   //  The initialisation of the matrix A
1760   Standard_Integer irow ;
1761   for ( irow=0; irow<=FirstOrder; irow++) {
1762     Standard_Real FirstVal = 1.;
1763 
1764     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1765       A(irow,icol) = Pattern[irow][icol]*FirstVal;
1766       if (irow <= icol) FirstVal *= FirstParameter;
1767     }
1768   }
1769 
1770   for (irow=0; irow<=LastOrder; irow++) {
1771     Standard_Real LastVal  = 1.;
1772 
1773     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1774       A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
1775       if (irow <= icol) LastVal *= LastParameter;
1776     }
1777   }
1778   //
1779   //  The filled matrix A for FirstOrder=LastOrder=2 is:
1780   //
1781   //      1   FP  FP**2   FP**3    FP**4     FP**5
1782   //      0   1   2*FP    3*FP**2  4*FP**3   5*FP**4        FP - FirstParameter
1783   //      0   0   2       6*FP     12*FP**2  20*FP**3
1784   //      1   LP  LP**2   LP**3    LP**4     LP**5
1785   //      0   1   2*LP    3*LP**2  4*LP**3   5*LP**4        LP - LastParameter
1786   //      0   0   2       6*LP     12*LP**2  20*LP**3
1787   //
1788   //  If FirstOrder or LastOrder <=2 then some rows and columns are missing.
1789   //  For example:
1790   //  If FirstOrder=1 then 3th row and 6th column are missing
1791   //  If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
1792 
1793   math_Gauss Equations(A);
1794   //  std::cout << "A=" << A << std::endl;
1795 
1796   for (Standard_Integer idim=1; idim<=Dimension; idim++) {
1797     //  std::cout << "idim=" << idim << std::endl;
1798 
1799     math_Vector B(0,FirstOrder+LastOrder+1);
1800     Standard_Integer icol ;
1801     for ( icol=0; icol<=FirstOrder; icol++)
1802       B(icol) = FirstConstr(idim,icol);
1803 
1804     for (icol=0; icol<=LastOrder; icol++)
1805       B(FirstOrder+1+icol) = LastConstr(idim,icol);
1806     //  std::cout << "B=" << B << std::endl;
1807 
1808     //  The solving of equations system A * X = B. Then B = X
1809     Equations.Solve(B);
1810     //  std::cout << "After Solving" << std::endl << "B=" << B << std::endl;
1811 
1812     if (Equations.IsDone()==Standard_False) return Standard_False;
1813 
1814     //  the filling of the Coefficients
1815 
1816     for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
1817       Coefficients(Dimension*icol+idim-1) = B(icol);
1818   }
1819   return Standard_True;
1820 }
1821 
1822 //=======================================================================
1823 //function : JacobiParameters
1824 //purpose  :
1825 //=======================================================================
1826 
JacobiParameters(const GeomAbs_Shape ConstraintOrder,const Standard_Integer MaxDegree,const Standard_Integer Code,Standard_Integer & NbGaussPoints,Standard_Integer & WorkDegree)1827 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
1828 			    const Standard_Integer MaxDegree,
1829 			    const Standard_Integer Code,
1830 			    Standard_Integer& NbGaussPoints,
1831 			    Standard_Integer& WorkDegree)
1832 {
1833   // ConstraintOrder: Ordre de contrainte aux extremites :
1834   //            C0 = contraintes de passage aux bornes;
1835   //            C1 = C0 + contraintes de derivees 1eres;
1836   //            C2 = C1 + contraintes de derivees 2ndes.
1837   // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
1838   //            d' approximation (doit etre superieur ou egal a
1839   //            2*NivConstr+2 et inferieur ou egal a 50).
1840   // Code:      Code d' init. des parametres de discretisation.
1841   //            (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
1842   //            = -5 Calcul tres rapide mais peu precis (8pts)
1843   //            = -4    '     '    '      '   '    '    (10pts)
1844   //            = -3    '     '    '      '   '    '    (15pts)
1845   //            = -2    '     '    '      '   '    '    (20pts)
1846   //            = -1    '     '    '      '   '    '    (25pts)
1847   //            = 1 calcul rapide avec precision moyenne (30pts).
1848   //            = 2 calcul rapide avec meilleure precision (40pts).
1849   //            = 3 calcul un peu plus lent avec bonne precision (50 pts).
1850   //            = 4 calcul lent avec la meilleure precision possible
1851   //             (61pts).
1852 
1853   // The possible values of NbGaussPoints
1854 
1855   const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
1856   NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
1857 
1858   Standard_Integer NivConstr=0;
1859   switch (ConstraintOrder) {
1860   case GeomAbs_C0: NivConstr = 0; break;
1861   case GeomAbs_C1: NivConstr = 1; break;
1862   case GeomAbs_C2: NivConstr = 2; break;
1863   default:
1864     throw Standard_ConstructionError("Invalid ConstraintOrder");
1865   }
1866   if (MaxDegree < 2*NivConstr+1)
1867     throw Standard_ConstructionError("Invalid MaxDegree");
1868 
1869   if (Code >= 1)
1870     WorkDegree = MaxDegree + 9;
1871   else
1872     WorkDegree = MaxDegree + 6;
1873 
1874   //---> Nbre mini de points necessaires.
1875   Standard_Integer IPMIN=0;
1876   if (WorkDegree < NDEG8)
1877     IPMIN=NDEG8;
1878   else if (WorkDegree < NDEG10)
1879     IPMIN=NDEG10;
1880   else if (WorkDegree < NDEG15)
1881     IPMIN=NDEG15;
1882   else if (WorkDegree < NDEG20)
1883     IPMIN=NDEG20;
1884   else if (WorkDegree < NDEG25)
1885     IPMIN=NDEG25;
1886   else if (WorkDegree < NDEG30)
1887     IPMIN=NDEG30;
1888   else if (WorkDegree < NDEG40)
1889     IPMIN=NDEG40;
1890   else if (WorkDegree < NDEG50)
1891     IPMIN=NDEG50;
1892   else if (WorkDegree < NDEG61)
1893     IPMIN=NDEG61;
1894   else
1895     throw Standard_ConstructionError("Invalid MaxDegree");
1896   // ---> Nbre de points voulus.
1897   Standard_Integer IWANT=0;
1898   switch (Code) {
1899   case -5: IWANT=NDEG8;  break;
1900   case -4: IWANT=NDEG10; break;
1901   case -3: IWANT=NDEG15; break;
1902   case -2: IWANT=NDEG20; break;
1903   case -1: IWANT=NDEG25; break;
1904   case  1: IWANT=NDEG30; break;
1905   case  2: IWANT=NDEG40; break;
1906   case  3: IWANT=NDEG50; break;
1907   case  4: IWANT=NDEG61; break;
1908   default:
1909     throw Standard_ConstructionError("Invalid Code");
1910   }
1911   //-->  NbGaussPoints est le nombre de points de discretisation de la fonction,
1912   //     il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
1913   //     NbGaussPoints doit etre superieur strictement a WorkDegree.
1914   NbGaussPoints = Max(IPMIN,IWANT);
1915   //  NbGaussPoints +=2;
1916 }
1917 
1918 //=======================================================================
1919 //function : NivConstr
1920 //purpose  : translates from GeomAbs_Shape to Integer
1921 //=======================================================================
1922 
NivConstr(const GeomAbs_Shape ConstraintOrder)1923  Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
1924 {
1925   Standard_Integer NivConstr=0;
1926   switch (ConstraintOrder) {
1927     case GeomAbs_C0: NivConstr = 0; break;
1928     case GeomAbs_C1: NivConstr = 1; break;
1929     case GeomAbs_C2: NivConstr = 2; break;
1930     default:
1931       throw Standard_ConstructionError("Invalid ConstraintOrder");
1932   }
1933   return NivConstr;
1934 }
1935 
1936 //=======================================================================
1937 //function : ConstraintOrder
1938 //purpose  : translates from Integer to GeomAbs_Shape
1939 //=======================================================================
1940 
ConstraintOrder(const Standard_Integer NivConstr)1941  GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
1942 {
1943   GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
1944   switch (NivConstr) {
1945     case 0: ConstraintOrder = GeomAbs_C0; break;
1946     case 1: ConstraintOrder = GeomAbs_C1; break;
1947     case 2: ConstraintOrder = GeomAbs_C2; break;
1948     default:
1949       throw Standard_ConstructionError("Invalid NivConstr");
1950   }
1951   return ConstraintOrder;
1952 }
1953 
1954 //=======================================================================
1955 //function : EvalLength
1956 //purpose  :
1957 //=======================================================================
1958 
EvalLength(const Standard_Integer Degree,const Standard_Integer Dimension,Standard_Real & PolynomialCoeff,const Standard_Real U1,const Standard_Real U2,Standard_Real & Length)1959  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
1960                       Standard_Real& PolynomialCoeff,
1961 		      const Standard_Real U1, const Standard_Real U2,
1962 		      Standard_Real& Length)
1963 {
1964   Standard_Integer i,j,idim, degdim;
1965   Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
1966 
1967   Standard_Real *PolynomialArray = &PolynomialCoeff ;
1968 
1969   Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
1970 
1971   math_Vector GaussPoints(1,NbGaussPoints);
1972   math::GaussPoints(NbGaussPoints,GaussPoints);
1973 
1974   math_Vector GaussWeights(1,NbGaussPoints);
1975   math::GaussWeights(NbGaussPoints,GaussWeights);
1976 
1977   C1 = (U2 + U1) / 2.;
1978   C2 = (U2 - U1) / 2.;
1979 
1980 //-----------------------------------------------------------
1981 //****** Integration - Boucle sur les intervalles de GAUSS **
1982 //-----------------------------------------------------------
1983 
1984   Sum = 0;
1985 
1986   for (j=1; j<=NbGaussPoints/2; j++) {
1987     // Integration en tenant compte de la symetrie
1988     Tran  = C2 * GaussPoints(j);
1989     X1 = C1 + Tran;
1990     X2 = C1 - Tran;
1991 
1992     //****** Derivation sur la dimension de l'espace **
1993 
1994     degdim =  Degree*Dimension;
1995     Der1 = Der2 = 0.;
1996     for (idim=0; idim<Dimension; idim++) {
1997       D1 = D2 = Degree * PolynomialArray [idim + degdim];
1998       for (i=Degree-1; i>=1; i--) {
1999 	DD = i * PolynomialArray [idim + i*Dimension];
2000 	D1 = D1 * X1 + DD;
2001 	D2 = D2 * X2 + DD;
2002       }
2003       Der1 += D1 * D1;
2004       Der2 += D2 * D2;
2005     }
2006 
2007     //****** Integration **
2008 
2009     Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
2010 
2011 //****** Fin de boucle dur les intervalles de GAUSS **
2012   }
2013   Length = Sum;
2014 }
2015 
2016 
2017 //=======================================================================
2018 //function : EvalLength
2019 //purpose  :
2020 //=======================================================================
2021 
EvalLength(const Standard_Integer Degree,const Standard_Integer Dimension,Standard_Real & PolynomialCoeff,const Standard_Real U1,const Standard_Real U2,const Standard_Real Tol,Standard_Real & Length,Standard_Real & Error)2022  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2023                       Standard_Real& PolynomialCoeff,
2024 		      const Standard_Real U1, const Standard_Real U2,
2025 		      const Standard_Real Tol,
2026 		      Standard_Real& Length, Standard_Real& Error)
2027 {
2028   Standard_Integer i;
2029   Standard_Integer NbSubInt = 1,    // Current number of subintervals
2030                    MaxNbIter = 13,  // Max number of iterations
2031                    NbIter    = 1;   // Current number of iterations
2032   Standard_Real    dU,OldLen,LenI;
2033 
2034   PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
2035 
2036   do {
2037     OldLen = Length;
2038     Length = 0.;
2039     NbSubInt *= 2;
2040     dU = (U2-U1)/NbSubInt;
2041     for (i=1; i<=NbSubInt; i++) {
2042       PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
2043       Length += LenI;
2044     }
2045     NbIter++;
2046     Error = Abs(OldLen-Length);
2047   }
2048   while (Error > Tol && NbIter <= MaxNbIter);
2049 }
2050