0031007: Coding - eliminate warnings issued while compiling with -pedantic flag
[occt.git] / src / PLib / PLib.cxx
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 //=======================================================================
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
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
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
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
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
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
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
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
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
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
255   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
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
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
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> 
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<>
651   inline void eval_step1<-1> (double*, double, double*)
652   {
653   }
654
655   // recursive template for evaluating second derivative
656   template<int dim> 
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<>
665   inline void eval_step2<-1> (double*, double, double*)
666   {
667   }
668
669   // evaluation of only value
670   template<int dim>
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>
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>
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
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
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
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  
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
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 devided 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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 }