9bcd5889ac00a77732ea116899af46921e3eb2a0
[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-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21 // Modified: 28/02/1996 by PMN :  HermiteCoefficients added
22 // Modified: 18/06/1996 by PMN :  NULL reference.
23 // Modified: 19/02/1997 by JCT :  EvalPoly2Var added
24
25 #include <PLib.ixx>
26 #include <NCollection_LocalArray.hxx>
27 #include <math_Matrix.hxx> 
28 #include <math_Gauss.hxx> 
29 #include <Standard_ConstructionError.hxx>
30 #include <GeomAbs_Shape.hxx>
31
32 #include <math_Gauss.hxx>
33 #include <math.hxx>
34
35 // To convert points array into Real ..
36 // *********************************
37
38 //=======================================================================
39 //function : SetPoles
40 //purpose  : 
41 //=======================================================================
42
43 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
44                     TColStd_Array1OfReal& FP)
45 {
46   Standard_Integer j      = FP   .Lower();
47   Standard_Integer PLower = Poles.Lower();
48   Standard_Integer PUpper = Poles.Upper();
49     
50   for (Standard_Integer i = PLower; i <= PUpper; i++) {
51     const gp_Pnt2d& P = Poles(i);
52     FP(j) = P.Coord(1); j++;
53     FP(j) = P.Coord(2); j++;
54   }
55 }
56
57 //=======================================================================
58 //function : SetPoles
59 //purpose  : 
60 //=======================================================================
61
62 void PLib::SetPoles(const TColgp_Array1OfPnt2d&       Poles,
63                     const TColStd_Array1OfReal& Weights,
64                     TColStd_Array1OfReal&       FP)
65 {
66   Standard_Integer j      = FP   .Lower();
67   Standard_Integer PLower = Poles.Lower();
68   Standard_Integer PUpper = Poles.Upper();
69     
70   for (Standard_Integer i = PLower; i <= PUpper; i++) {
71     Standard_Real w = Weights(i);
72     const gp_Pnt2d& P = Poles(i);
73     FP(j) = P.Coord(1) * w; j++;
74     FP(j) = P.Coord(2) * w; j++;
75     FP(j) =              w; j++;
76   }
77 }
78
79 //=======================================================================
80 //function : GetPoles
81 //purpose  : 
82 //=======================================================================
83
84 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
85                     TColgp_Array1OfPnt2d& Poles)
86 {
87   Standard_Integer j      = FP   .Lower();
88   Standard_Integer PLower = Poles.Lower();
89   Standard_Integer PUpper = Poles.Upper();
90     
91   for (Standard_Integer i = PLower; i <= PUpper; i++) {
92     gp_Pnt2d& P = Poles(i);
93     P.SetCoord(1,FP(j)); j++;
94     P.SetCoord(2,FP(j)); j++;
95   }
96 }
97
98 //=======================================================================
99 //function : GetPoles
100 //purpose  : 
101 //=======================================================================
102
103 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
104                     TColgp_Array1OfPnt2d&       Poles,
105                     TColStd_Array1OfReal&       Weights)
106 {
107   Standard_Integer j      = FP   .Lower();
108   Standard_Integer PLower = Poles.Lower();
109   Standard_Integer PUpper = Poles.Upper();
110     
111   for (Standard_Integer i = PLower; i <= PUpper; i++) {
112     Standard_Real w = FP(j + 2);
113     Weights(i) = w;
114     gp_Pnt2d& P = Poles(i);
115     P.SetCoord(1,FP(j) / w); j++;
116     P.SetCoord(2,FP(j) / w); j++;
117     j++;
118   }
119 }
120
121 //=======================================================================
122 //function : SetPoles
123 //purpose  : 
124 //=======================================================================
125
126 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
127                     TColStd_Array1OfReal& FP)
128 {
129   Standard_Integer j      = FP   .Lower();
130   Standard_Integer PLower = Poles.Lower();
131   Standard_Integer PUpper = Poles.Upper();
132
133   for (Standard_Integer i = PLower; i <= PUpper; i++) {
134     const gp_Pnt& P = Poles(i);
135     FP(j) = P.Coord(1); j++;
136     FP(j) = P.Coord(2); j++;
137     FP(j) = P.Coord(3); j++;
138   }
139 }
140
141 //=======================================================================
142 //function : SetPoles
143 //purpose  : 
144 //=======================================================================
145
146 void PLib::SetPoles(const TColgp_Array1OfPnt&   Poles,
147                     const TColStd_Array1OfReal& Weights,
148                     TColStd_Array1OfReal&       FP)
149 {
150   Standard_Integer j      = FP   .Lower();
151   Standard_Integer PLower = Poles.Lower();
152   Standard_Integer PUpper = Poles.Upper();
153     
154   for (Standard_Integer i = PLower; i <= PUpper; i++) {
155     Standard_Real w = Weights(i);
156     const gp_Pnt& P = Poles(i);
157     FP(j) = P.Coord(1) * w; j++;
158     FP(j) = P.Coord(2) * w; j++;
159     FP(j) = P.Coord(3) * w; j++;
160     FP(j) =              w; j++;
161   }
162 }
163
164 //=======================================================================
165 //function : GetPoles
166 //purpose  : 
167 //=======================================================================
168
169 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
170                     TColgp_Array1OfPnt&         Poles)
171 {
172   Standard_Integer j      = FP   .Lower();
173   Standard_Integer PLower = Poles.Lower();
174   Standard_Integer PUpper = Poles.Upper();
175     
176   for (Standard_Integer i = PLower; i <= PUpper; i++) {
177     gp_Pnt& P = Poles(i);
178     P.SetCoord(1,FP(j)); j++;
179     P.SetCoord(2,FP(j)); j++;
180     P.SetCoord(3,FP(j)); j++;
181   }
182 }
183
184 //=======================================================================
185 //function : GetPoles
186 //purpose  : 
187 //=======================================================================
188
189 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
190                     TColgp_Array1OfPnt&         Poles,
191                     TColStd_Array1OfReal&       Weights)
192 {
193   Standard_Integer j      = FP   .Lower();
194   Standard_Integer PLower = Poles.Lower();
195   Standard_Integer PUpper = Poles.Upper();
196     
197   for (Standard_Integer i = PLower; i <= PUpper; i++) {
198     Standard_Real w = FP(j + 3);
199     Weights(i) = w;
200     gp_Pnt& P = Poles(i);
201     P.SetCoord(1,FP(j) / w); j++;
202     P.SetCoord(2,FP(j) / w); j++;
203     P.SetCoord(3,FP(j) / w); j++;
204     j++;
205   }
206 }
207
208 // specialized allocator
209 namespace
210 {
211
212 class BinomAllocator
213 {
214 public:
215
216   //! Main constructor
217   BinomAllocator (const Standard_Integer theMaxBinom)
218   : myBinom (NULL),
219     myMaxBinom (theMaxBinom)
220   {
221     Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
222     Standard_Integer np1 = myMaxBinom + 1;
223     myBinom = new Standard_Integer*[np1];
224     myBinom[0] = new Standard_Integer[1];
225     myBinom[0][0] = 1;
226     for (i = 1; i < np1; ++i)
227     {
228       im1 = i - 1;
229       ip1 = i + 1;
230       id2 = i >> 1;
231       md2 = im1 >> 1;
232       md3 = ip1 >> 1;
233       k   = 0;
234       myBinom[i] = new Standard_Integer[ip1];
235
236       for (j = 0; j < id2; ++j)
237       {
238         myBinom[i][j] = k + myBinom[im1][j];
239         k = myBinom[im1][j];
240       }
241       j = id2;
242       if (j > md2) j = im1 - j;
243       myBinom[i][id2] = k + myBinom[im1][j];
244
245       for (j = ip1 - md3; j < ip1; j++)
246       {
247         myBinom[i][j] = myBinom[i][i - j];
248       }
249     }
250   }
251
252   //! Destructor
253   ~BinomAllocator()
254   {
255     // free memory
256     for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
257     {
258       delete[] myBinom[i];
259     }
260     delete[] myBinom;
261   }
262
263   Standard_Real Value (const Standard_Integer N,
264                        const Standard_Integer P) const
265   {
266     Standard_OutOfRange_Raise_if (N > myMaxBinom,
267       "PLib, BinomAllocator: requested degree is greater than maximum supported");
268     return Standard_Real (myBinom[N][P]);
269   }
270
271 private:
272   Standard_Integer**  myBinom;
273   Standard_Integer myMaxBinom;
274
275 };
276
277   // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
278   //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
279   static BinomAllocator THE_BINOM (25 + 1);
280 };
281
282 //=======================================================================
283 //function : Bin
284 //purpose  : 
285 //=======================================================================
286
287 Standard_Real PLib::Bin(const Standard_Integer N,
288                         const Standard_Integer P)
289 {
290   return THE_BINOM.Value (N, P);
291 }
292
293 //=======================================================================
294 //function : RationalDerivative
295 //purpose  : 
296 //=======================================================================
297
298 void  PLib::RationalDerivative(const Standard_Integer Degree,
299                                const Standard_Integer DerivativeRequest,
300                                const Standard_Integer Dimension,
301                                Standard_Real& Ders,
302                                Standard_Real& RDers,
303                                const Standard_Boolean All)
304 {
305   //
306   // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
307   // 
308   //  We Write  u = fv
309   //  Let C(N,P) be the binomial
310   //
311   //        then we have 
312   //
313   //         (q)                             (p)   (q-p) 
314   //        u    =     SUM          C (q,p) f     v
315   //                p = 0 to q
316   //
317   //
318   //        Therefore 
319   //        
320   //          
321   //         (q)         (   (q)                               (p)   (q-p)   )
322   //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
323   //                     (          p = 0 to q-1                             )
324   //
325   //
326   // make arrays for the binomial since computing it each time could raise a performance
327   // issue
328   // As oppose to the method below the <Der> array is organized in the following 
329   // fashion :
330   //
331   //  u (1)     u (2) ....   u  (Dimension)     v (1) 
332   //
333   //   (1)       (1)          (1)                (1) 
334   //  u (1)     u (2) ....   u  (Dimension)     v (1)
335   //
336   //   ............................................
337   // 
338   //   (Degree)  (Degree)     (Degree)           (Degree) 
339   //  u (1)     u (2) ....   u  (Dimension)     v (1) 
340   //
341   //
342   Standard_Real Inverse;
343   Standard_Real *PolesArray = &Ders;
344   Standard_Real *RationalArray = &RDers;
345   Standard_Real Factor ;
346   Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
347   NCollection_LocalArray<Standard_Real> binomial_array;
348   NCollection_LocalArray<Standard_Real> derivative_storage;
349   if (Dimension == 3) {
350     Standard_Integer DeRequest1 = DerivativeRequest + 1;
351     Standard_Integer MinDegRequ = DerivativeRequest;
352     if (MinDegRequ > Degree) MinDegRequ = Degree;
353     binomial_array.Allocate (DeRequest1);
354     
355     for (ii = 0 ; ii < DeRequest1 ; ii++) {
356       binomial_array[ii] = 1.0e0 ;
357     }
358     if (!All) {
359       Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
360       derivative_storage.Allocate (DimDeRequ1);
361       RationalArray = derivative_storage ;
362     }
363     
364     Inverse = 1.0e0 / PolesArray[3]  ;
365     Index = 0 ;
366     Index2 = - 6;
367     OtherIndex = 0 ;
368     
369     for (ii = 0 ; ii <= MinDegRequ ; ii++) {
370       Index2 += 3;
371       Index1  = Index2;
372       RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
373       RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
374       RationalArray[Index] = PolesArray[OtherIndex];
375       Index      -= 2;
376       OtherIndex += 2;
377       
378       for (jj = ii - 1 ; jj >= 0 ; jj--) {
379         Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3]; 
380         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
381         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
382         RationalArray[Index] -=  Factor * RationalArray[Index1];
383         Index  -= 2;
384         Index1 -= 5;
385       }
386       
387       for (jj = ii ; jj >=  1 ; jj--) {
388         binomial_array[jj] += binomial_array[jj - 1] ;
389       }
390       RationalArray[Index] *= Inverse ; Index++;
391       RationalArray[Index] *= Inverse ; Index++;
392       RationalArray[Index] *= Inverse ; Index++;
393     } 
394     
395     for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
396       Index2 += 3;
397       Index1 = Index2;
398       RationalArray[Index] = 0.0e0; Index++;
399       RationalArray[Index] = 0.0e0; Index++;
400       RationalArray[Index] = 0.0e0;
401       Index -= 2;
402       
403       for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
404         Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3]; 
405         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
406         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
407         RationalArray[Index] -=  Factor * RationalArray[Index1];
408         Index  -= 2;
409         Index1 -= 5;
410       }
411       
412       for (jj = ii ; jj >=  1 ; jj--) {
413         binomial_array[jj] += binomial_array[jj - 1] ;
414       }
415       RationalArray[Index] *= Inverse; Index++;
416       RationalArray[Index] *= Inverse; Index++;
417       RationalArray[Index] *= Inverse; Index++;
418     }
419     
420     if (!All) {
421       RationalArray = &RDers ;
422       Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
423       RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
424       RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
425       RationalArray[2] = derivative_storage[DimDeRequ];
426     }
427   }
428   else {
429     Standard_Integer kk;
430     Standard_Integer Dimension1 = Dimension + 1;
431     Standard_Integer Dimension2 = Dimension << 1;
432     Standard_Integer DeRequest1 = DerivativeRequest + 1;
433     Standard_Integer MinDegRequ = DerivativeRequest;
434     if (MinDegRequ > Degree) MinDegRequ = Degree;
435     binomial_array.Allocate (DeRequest1);
436     
437     for (ii = 0 ; ii < DeRequest1 ; ii++) {
438       binomial_array[ii] = 1.0e0 ;
439     }
440     if (!All) {
441       Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
442       derivative_storage.Allocate (DimDeRequ1);
443       RationalArray = derivative_storage ;
444     }
445     
446     Inverse = 1.0e0 / PolesArray[Dimension]  ;
447     Index = 0 ;
448     Index2 = - Dimension2;
449     OtherIndex = 0 ;
450     
451     for (ii = 0 ; ii <= MinDegRequ ; ii++) {
452       Index2 += Dimension;
453       Index1  = Index2;
454       
455       for (kk = 0 ; kk < Dimension ; kk++) {
456         RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
457       }
458       Index -= Dimension;
459       OtherIndex ++;;
460       
461       for (jj = ii - 1 ; jj >= 0 ; jj--) {
462         Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension]; 
463         
464         for (kk = 0 ; kk < Dimension ; kk++) {
465           RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
466         }
467         Index  -= Dimension ;
468         Index1 -= Dimension2 ;
469       }
470       
471       for (jj = ii ; jj >=  1 ; jj--) {
472         binomial_array[jj] += binomial_array[jj - 1] ;
473       }
474       
475       for (kk = 0 ; kk < Dimension ; kk++) {
476         RationalArray[Index] *= Inverse ; Index++;
477       }
478     } 
479     
480     for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
481       Index2 += Dimension;
482       Index1  = Index2;
483       
484       for (kk = 0 ; kk < Dimension ; kk++) {
485         RationalArray[Index] = 0.0e0 ; Index++;
486       }
487       Index -= Dimension;
488       
489       for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
490         Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension]; 
491         
492         for (kk = 0 ; kk < Dimension ; kk++) {
493           RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
494         }
495         Index  -= Dimension ;
496         Index1 -= Dimension2 ;
497       }
498       
499       for (jj = ii ; jj >=  1 ; jj--) {
500         binomial_array[jj] += binomial_array[jj - 1] ;
501       }
502       
503       for (kk = 0 ; kk < Dimension ; kk++) {
504         RationalArray[Index] *= Inverse; Index++;
505       }
506     }
507     
508     if (!All) {
509       RationalArray = &RDers ;
510       Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
511       
512       for (kk = 0 ; kk < Dimension ; kk++) {
513         RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
514       }
515     }
516   }
517
518
519 //=======================================================================
520 //function : RationalDerivatives
521 //purpose  : Uses Homogeneous Poles Derivatives and Deivatives of Weights
522 //=======================================================================
523
524 void  PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
525                                 const Standard_Integer  Dimension,
526                                 Standard_Real& PolesDerivates, 
527                                 // must be an array with 
528                                 // (DerivativeRequest + 1) * Dimension slots
529                                 Standard_Real& WeightsDerivates,
530                                 // must be an array with 
531                                 // (DerivativeRequest + 1) slots
532                                 Standard_Real& RationalDerivates) 
533 {
534   //
535   // Our purpose is to compute f = (u/v) derivated N times
536   // 
537   //  We Write  u = fv
538   //  Let C(N,P) be the binomial
539   //
540   //        then we have 
541   //
542   //         (q)                             (p)   (q-p) 
543   //        u    =     SUM          C (q,p) f     v
544   //                p = 0 to q
545   //
546   //
547   //        Therefore 
548   //        
549   //          
550   //         (q)         (   (q)                               (p)   (q-p)   )
551   //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
552   //                     (          p = 0 to q-1                             )
553   //
554   //
555   // make arrays for the binomial since computing it each time could 
556   // raize a performance issue
557   // 
558   Standard_Real Inverse;
559   Standard_Real *PolesArray = &PolesDerivates;
560   Standard_Real *WeightsArray = &WeightsDerivates;
561   Standard_Real *RationalArray = &RationalDerivates;
562   Standard_Real Factor ;
563   
564   Standard_Integer ii, Index, Index1, Index2, jj;
565   Standard_Integer DeRequest1 = DerivativeRequest + 1;
566   
567   NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
568   NCollection_LocalArray<Standard_Real> derivative_storage;
569
570   for (ii = 0 ; ii < DeRequest1 ; ii++) {
571     binomial_array[ii] = 1.0e0 ;
572   }
573   Inverse = 1.0e0 / WeightsArray[0]  ;
574   if (Dimension == 3) {
575     Index  = 0 ;
576     Index2 = - 6 ;
577
578     for (ii = 0 ; ii < DeRequest1 ; ii++) {
579       Index2 += 3;
580       Index1 = Index2;
581       RationalArray[Index] = PolesArray[Index] ; Index++;
582       RationalArray[Index] = PolesArray[Index] ; Index++;
583       RationalArray[Index] = PolesArray[Index] ;
584       Index -= 2;
585       
586       for (jj = ii - 1 ; jj >= 0 ; jj--) {
587         Factor = binomial_array[jj] * WeightsArray[ii - jj] ; 
588         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
589         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
590         RationalArray[Index] -=  Factor * RationalArray[Index1];
591         Index  -= 2;
592         Index1 -= 5;
593       }
594       
595       for (jj = ii ; jj >=  1 ; jj--) {
596         binomial_array[jj] += binomial_array[jj - 1] ;
597       }
598       RationalArray[Index] *= Inverse ; Index++;
599       RationalArray[Index] *= Inverse ; Index++;
600       RationalArray[Index] *= Inverse ; Index++;
601     }
602   }
603   else {
604     Standard_Integer kk;
605     Standard_Integer Dimension2 = Dimension << 1;
606     Index  = 0 ;
607     Index2 = - Dimension2;
608
609     for (ii = 0 ; ii < DeRequest1 ; ii++) {
610       Index2 += Dimension;
611       Index1  = Index2;
612       
613       for (kk = 0 ; kk < Dimension ; kk++) {
614         RationalArray[Index] = PolesArray[Index]; Index++;
615       }
616       Index  -= Dimension;
617       
618       for (jj = ii - 1 ; jj >= 0 ; jj--) {
619         Factor = binomial_array[jj] * WeightsArray[ii - jj] ; 
620         
621         for (kk = 0 ; kk < Dimension ; kk++) {
622           RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
623         }
624         Index  -= Dimension;
625         Index1 -= Dimension2;
626       }
627       
628       for (jj = ii ; jj >=  1 ; jj--) {
629         binomial_array[jj] += binomial_array[jj - 1] ;
630       }
631       
632       for (kk = 0 ; kk < Dimension ; kk++) {
633         RationalArray[Index] *= Inverse ; Index++;
634       }
635     }
636   }
637 }
638
639 //=======================================================================
640 //function : This evaluates a polynomial and its derivatives 
641 //purpose  : up to the requested order 
642 //=======================================================================
643
644 void  PLib::EvalPolynomial(const Standard_Real    Par,
645                            const Standard_Integer DerivativeRequest,
646                            const Standard_Integer Degree,
647                            const Standard_Integer Dimension, 
648                            Standard_Real&         PolynomialCoeff,
649                            Standard_Real&         Results)
650      //
651      // the polynomial coefficients are assumed to be stored as follows :
652      //                                                      0    
653      // [0]                  [Dimension -1]                 X     coefficient
654      //                                                      1 
655      // [Dimension]          [Dimension + Dimension -1]     X     coefficient
656      //                                                      2 
657      // [2 * Dimension]      [2 * Dimension + Dimension-1]  X     coefficient
658      //
659      //   ...................................................
660      //
661      //  
662      //                                                      d 
663      // [d * Dimension]      [d * Dimension + Dimension-1]  X     coefficient
664      //
665      //  where d is the Degree
666      //
667 {
668   Standard_Integer DegreeDimension = Degree * Dimension;
669
670   Standard_Integer jj;
671   Standard_Real *RA = &Results ;  
672   Standard_Real *PA = &PolynomialCoeff ;
673   Standard_Real *tmpRA = RA;
674   Standard_Real *tmpPA = PA + DegreeDimension;
675   
676   switch (Dimension) {
677     
678   case 1 : {
679     *tmpRA = *tmpPA; 
680     if (DerivativeRequest > 0 ) {
681       tmpRA++ ;
682       Standard_Real *valRA;
683       Standard_Integer ii, LocalRequest;
684       Standard_Integer Index1, Index2;
685       Standard_Integer MaxIndex1, MaxIndex2;
686       if      (DerivativeRequest < Degree) {
687         LocalRequest = DerivativeRequest;
688         MaxIndex2 = MaxIndex1 = LocalRequest;
689       }
690       else {
691         LocalRequest = Degree;
692         MaxIndex2 = MaxIndex1 = Degree;
693       }
694       MaxIndex2 --;
695       
696       for (ii = 1; ii <= LocalRequest; ii++) {
697         *tmpRA = 0.0e0;  tmpRA ++ ;          
698       }
699       
700       for (jj = Degree ; jj >  0 ; jj--) {
701         tmpPA --;
702         Index1 = MaxIndex1;
703         Index2 = MaxIndex2;
704         
705         for (ii = LocalRequest ; ii > 0 ; ii--) {
706           valRA  = &RA[Index1];
707           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
708           Index1 --;
709           Index2 --;
710         }
711         valRA  = &RA[Index1];
712         *valRA = Par * (*valRA) + (*tmpPA);
713       }
714     }
715     else {
716       
717       for (jj = Degree  ; jj >  0 ; jj--) {
718         tmpPA --;
719         *tmpRA = Par * (*tmpRA) + (*tmpPA);
720       }
721     }
722     break;
723   }
724   case 2 : {
725     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
726     *tmpRA = *tmpPA; tmpRA++;
727     tmpPA --;
728     if (DerivativeRequest > 0 ) {
729       Standard_Real *valRA;
730       Standard_Integer ii, LocalRequest;
731       Standard_Integer Index1, Index2;
732       Standard_Integer MaxIndex1, MaxIndex2;
733       if      (DerivativeRequest < Degree) {
734         LocalRequest = DerivativeRequest;
735         MaxIndex2 = MaxIndex1 = LocalRequest << 1;
736       }
737       else {
738         LocalRequest = Degree;
739         MaxIndex2 = MaxIndex1 = DegreeDimension;
740       }
741       MaxIndex2 -= 2;
742       
743       for (ii = 1; ii <= LocalRequest; ii++) {
744         *tmpRA = 0.0e0; tmpRA++;           
745         *tmpRA = 0.0e0; tmpRA++;           
746       }
747       
748       for (jj = Degree ; jj >  0 ; jj--) {
749         tmpPA -= 2;
750         
751         Index1 = MaxIndex1;
752         Index2 = MaxIndex2;
753         
754         for (ii = LocalRequest ; ii > 0 ; ii--) {
755           valRA  = &RA[Index1];
756           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
757           Index1 -= 2;
758           Index2 -= 2;
759         }
760         valRA  = &RA[Index1];
761         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
762         
763         Index1 = MaxIndex1 + 1;
764         Index2 = MaxIndex2 + 1;
765         
766         for (ii = LocalRequest ; ii > 0 ; ii--) {
767           valRA  = &RA[Index1];
768           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
769           Index1 -= 2;
770           Index2 -= 2;
771         }
772         valRA  = &RA[Index1];
773         *valRA = Par * (*valRA) + (*tmpPA);
774         
775         tmpPA --;
776       }
777     }
778     else {
779       
780       for (jj = Degree  ; jj >  0 ; jj--) {
781         tmpPA -= 2;
782         tmpRA  = RA;
783         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
784         *tmpRA = Par * (*tmpRA) + (*tmpPA);
785         tmpPA --;
786       }
787     }
788     break;
789   }
790   case 3 : {
791     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
792     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
793     *tmpRA = *tmpPA; tmpRA++;
794     tmpPA -= 2;
795     if (DerivativeRequest > 0 ) {
796       Standard_Real *valRA;
797       Standard_Integer ii, LocalRequest;
798       Standard_Integer Index1, Index2;
799       Standard_Integer MaxIndex1, MaxIndex2;
800       if      (DerivativeRequest < Degree) {
801         LocalRequest = DerivativeRequest;
802         MaxIndex2 = MaxIndex1 = (LocalRequest << 1) + LocalRequest;
803       }
804       else {
805         LocalRequest = Degree;
806         MaxIndex2 = MaxIndex1 = DegreeDimension;
807       }
808       MaxIndex2 -= 3;
809       
810       for (ii = 1; ii <= LocalRequest; ii++) {
811         *tmpRA = 0.0e0; tmpRA++;           
812         *tmpRA = 0.0e0; tmpRA++;           
813         *tmpRA = 0.0e0; tmpRA++;           
814       }
815       
816       for (jj = Degree ; jj >  0 ; jj--) {
817         tmpPA -= 3;
818         
819         Index1 = MaxIndex1;
820         Index2 = MaxIndex2;
821         
822         for (ii = LocalRequest ; ii > 0 ; ii--) {
823           valRA  = &RA[Index1];
824           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
825           Index1 -= 3;
826           Index2 -= 3;
827         }
828         valRA  = &RA[Index1];
829         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
830         
831         Index1 = MaxIndex1 + 1;
832         Index2 = MaxIndex2 + 1;
833         
834         for (ii = LocalRequest ; ii > 0 ; ii--) {
835           valRA  = &RA[Index1];
836           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
837           Index1 -= 3;
838           Index2 -= 3;
839         }
840         valRA  = &RA[Index1];
841         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
842         
843         Index1 = MaxIndex1 + 2;
844         Index2 = MaxIndex2 + 2;
845         
846         for (ii = LocalRequest ; ii > 0 ; ii--) {
847           valRA  = &RA[Index1];
848           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
849           Index1 -= 3;
850           Index2 -= 3;
851         }
852         valRA  = &RA[Index1];
853         *valRA = Par * (*valRA) + (*tmpPA);
854         
855         tmpPA -= 2;
856       }
857     }
858     else {
859       
860       for (jj = Degree  ; jj >  0 ; jj--) {
861         tmpPA -= 3;
862         tmpRA  = RA;
863         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
864         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
865         *tmpRA = Par * (*tmpRA) + (*tmpPA);
866         tmpPA -= 2;
867       }
868     }
869     break;
870   }
871   case 6 : {
872     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
873     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
874     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
875
876     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
877     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
878     *tmpRA = *tmpPA; tmpRA++;
879     tmpPA -= 5;
880     if (DerivativeRequest > 0 ) {
881       Standard_Real *valRA;
882       Standard_Integer ii, LocalRequest;
883       Standard_Integer Index1, Index2;
884       Standard_Integer MaxIndex1, MaxIndex2;
885       if      (DerivativeRequest < Degree) {
886         LocalRequest = DerivativeRequest;
887         MaxIndex2 = MaxIndex1 = (LocalRequest << 2) + (LocalRequest << 1);
888       }
889       else {
890         LocalRequest = Degree;
891         MaxIndex2 = MaxIndex1 = DegreeDimension;
892       }
893       MaxIndex2 -= 6;
894       
895       for (ii = 1; ii <= LocalRequest; ii++) {
896         *tmpRA = 0.0e0; tmpRA++;           
897         *tmpRA = 0.0e0; tmpRA++;           
898         *tmpRA = 0.0e0; tmpRA++;           
899
900         *tmpRA = 0.0e0; tmpRA++;           
901         *tmpRA = 0.0e0; tmpRA++;           
902         *tmpRA = 0.0e0; tmpRA++;           
903       }
904       
905       for (jj = Degree ; jj >  0 ; jj--) {
906         tmpPA -= 6;
907         
908         Index1 = MaxIndex1;
909         Index2 = MaxIndex2;
910         
911         for (ii = LocalRequest ; ii > 0 ; ii--) {
912           valRA  = &RA[Index1];
913           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
914           Index1 -= 6;
915           Index2 -= 6;
916         }
917         valRA  = &RA[Index1];
918         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
919         
920         Index1 = MaxIndex1 + 1;
921         Index2 = MaxIndex2 + 1;
922         
923         for (ii = LocalRequest ; ii > 0 ; ii--) {
924           valRA  = &RA[Index1];
925           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
926           Index1 -= 6;
927           Index2 -= 6;
928         }
929         valRA  = &RA[Index1];
930         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
931         
932         Index1 = MaxIndex1 + 2;
933         Index2 = MaxIndex2 + 2;
934         
935         for (ii = LocalRequest ; ii > 0 ; ii--) {
936           valRA  = &RA[Index1];
937           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
938           Index1 -= 6;
939           Index2 -= 6;
940         }
941         valRA  = &RA[Index1];
942         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
943         
944         Index1 = MaxIndex1 + 3;
945         Index2 = MaxIndex2 + 3;
946         
947         for (ii = LocalRequest ; ii > 0 ; ii--) {
948           valRA  = &RA[Index1];
949           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
950           Index1 -= 6;
951           Index2 -= 6;
952         }
953         valRA  = &RA[Index1];
954         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
955         
956         Index1 = MaxIndex1 + 4;
957         Index2 = MaxIndex2 + 4;
958         
959         for (ii = LocalRequest ; ii > 0 ; ii--) {
960           valRA  = &RA[Index1];
961           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
962           Index1 -= 6;
963           Index2 -= 6;
964         }
965         valRA  = &RA[Index1];
966         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
967         
968         Index1 = MaxIndex1 + 5;
969         Index2 = MaxIndex2 + 5;
970         
971         for (ii = LocalRequest ; ii > 0 ; ii--) {
972           valRA  = &RA[Index1];
973           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
974           Index1 -= 6;
975           Index2 -= 6;
976         }
977         valRA  = &RA[Index1];
978         *valRA = Par * (*valRA) + (*tmpPA);
979         
980         tmpPA -= 5;
981       }
982     }
983     else {
984       
985       for (jj = Degree  ; jj >  0 ; jj--) {
986         tmpPA -= 6;
987         tmpRA  = RA;
988
989         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
990         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
991         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
992
993         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
994         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
995         *tmpRA = Par * (*tmpRA) + (*tmpPA);
996         tmpPA -= 5;
997       }
998     }
999     break;
1000   }
1001   case 9 : {
1002     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1003     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1004     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1005
1006     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1007     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1008     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1009
1010     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1011     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1012     *tmpRA = *tmpPA; tmpRA++;
1013     tmpPA -= 8;
1014     if (DerivativeRequest > 0 ) {
1015       Standard_Real *valRA;
1016       Standard_Integer ii, LocalRequest;
1017       Standard_Integer Index1, Index2;
1018       Standard_Integer MaxIndex1, MaxIndex2;
1019       if      (DerivativeRequest < Degree) {
1020         LocalRequest = DerivativeRequest;
1021         MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + LocalRequest;
1022       }
1023       else {
1024         LocalRequest = Degree;
1025         MaxIndex2 = MaxIndex1 = DegreeDimension;
1026       }
1027       MaxIndex2 -= 9;
1028       
1029       for (ii = 1; ii <= LocalRequest; ii++) {
1030         *tmpRA = 0.0e0; tmpRA++;           
1031         *tmpRA = 0.0e0; tmpRA++;           
1032         *tmpRA = 0.0e0; tmpRA++;           
1033
1034         *tmpRA = 0.0e0; tmpRA++;           
1035         *tmpRA = 0.0e0; tmpRA++;           
1036         *tmpRA = 0.0e0; tmpRA++;           
1037
1038         *tmpRA = 0.0e0; tmpRA++;           
1039         *tmpRA = 0.0e0; tmpRA++;           
1040         *tmpRA = 0.0e0; tmpRA++;           
1041       }
1042       
1043       for (jj = Degree ; jj >  0 ; jj--) {
1044         tmpPA -= 9;
1045         
1046         Index1 = MaxIndex1;
1047         Index2 = MaxIndex2;
1048         
1049         for (ii = LocalRequest ; ii > 0 ; ii--) {
1050           valRA  = &RA[Index1];
1051           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1052           Index1 -= 9;
1053           Index2 -= 9;
1054         }
1055         valRA  = &RA[Index1];
1056         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1057         
1058         Index1 = MaxIndex1 + 1;
1059         Index2 = MaxIndex2 + 1;
1060         
1061         for (ii = LocalRequest ; ii > 0 ; ii--) {
1062           valRA  = &RA[Index1];
1063           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1064           Index1 -= 9;
1065           Index2 -= 9;
1066         }
1067         valRA  = &RA[Index1];
1068         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1069         
1070         Index1 = MaxIndex1 + 2;
1071         Index2 = MaxIndex2 + 2;
1072         
1073         for (ii = LocalRequest ; ii > 0 ; ii--) {
1074           valRA  = &RA[Index1];
1075           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1076           Index1 -= 9;
1077           Index2 -= 9;
1078         }
1079         valRA  = &RA[Index1];
1080         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1081         
1082         Index1 = MaxIndex1 + 3;
1083         Index2 = MaxIndex2 + 3;
1084         
1085         for (ii = LocalRequest ; ii > 0 ; ii--) {
1086           valRA  = &RA[Index1];
1087           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1088           Index1 -= 9;
1089           Index2 -= 9;
1090         }
1091         valRA  = &RA[Index1];
1092         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1093         
1094         Index1 = MaxIndex1 + 4;
1095         Index2 = MaxIndex2 + 4;
1096         
1097         for (ii = LocalRequest ; ii > 0 ; ii--) {
1098           valRA  = &RA[Index1];
1099           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1100           Index1 -= 9;
1101           Index2 -= 9;
1102         }
1103         valRA  = &RA[Index1];
1104         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1105         
1106         Index1 = MaxIndex1 + 5;
1107         Index2 = MaxIndex2 + 5;
1108         
1109         for (ii = LocalRequest ; ii > 0 ; ii--) {
1110           valRA  = &RA[Index1];
1111           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1112           Index1 -= 9;
1113           Index2 -= 9;
1114         }
1115         valRA  = &RA[Index1];
1116         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1117         
1118         Index1 = MaxIndex1 + 6;
1119         Index2 = MaxIndex2 + 6;
1120         
1121         for (ii = LocalRequest ; ii > 0 ; ii--) {
1122           valRA  = &RA[Index1];
1123           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1124           Index1 -= 9;
1125           Index2 -= 9;
1126         }
1127         valRA  = &RA[Index1];
1128         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1129         
1130         Index1 = MaxIndex1 + 7;
1131         Index2 = MaxIndex2 + 7;
1132         
1133         for (ii = LocalRequest ; ii > 0 ; ii--) {
1134           valRA  = &RA[Index1];
1135           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1136           Index1 -= 9;
1137           Index2 -= 9;
1138         }
1139         valRA  = &RA[Index1];
1140         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1141         
1142         Index1 = MaxIndex1 + 8;
1143         Index2 = MaxIndex2 + 8;
1144         
1145         for (ii = LocalRequest ; ii > 0 ; ii--) {
1146           valRA  = &RA[Index1];
1147           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1148           Index1 -= 9;
1149           Index2 -= 9;
1150         }
1151         valRA  = &RA[Index1];
1152         *valRA = Par * (*valRA) + (*tmpPA);
1153         
1154         tmpPA -= 8;
1155       }
1156     }
1157     else {
1158       
1159       for (jj = Degree  ; jj >  0 ; jj--) {
1160         tmpPA -= 9;
1161         tmpRA  = RA;
1162
1163         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1164         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1165         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1166
1167         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1168         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1169         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1170
1171         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1172         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1173         *tmpRA = Par * (*tmpRA) + (*tmpPA);
1174         tmpPA -= 8;
1175       }
1176     }
1177     break;
1178   }
1179   case 12 : {
1180     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1181     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1182     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1183     
1184     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1185     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1186     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1187     
1188     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1189     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1190     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1191     
1192     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1193     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1194     *tmpRA = *tmpPA; tmpRA++;
1195     tmpPA -= 11;
1196     if (DerivativeRequest > 0 ) {
1197       Standard_Real *valRA;
1198       Standard_Integer ii, LocalRequest;
1199       Standard_Integer Index1, Index2;
1200       Standard_Integer MaxIndex1, MaxIndex2;
1201       if      (DerivativeRequest < Degree) {
1202         LocalRequest = DerivativeRequest;
1203         MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + (LocalRequest << 2);
1204       }
1205       else {
1206         LocalRequest = Degree;
1207         MaxIndex2 = MaxIndex1 = DegreeDimension;
1208       }
1209       MaxIndex2 -= 12;
1210       
1211       for (ii = 1; ii <= LocalRequest; ii++) {
1212         *tmpRA = 0.0e0; tmpRA++;           
1213         *tmpRA = 0.0e0; tmpRA++;           
1214         *tmpRA = 0.0e0; tmpRA++;           
1215         
1216         *tmpRA = 0.0e0; tmpRA++;           
1217         *tmpRA = 0.0e0; tmpRA++;           
1218         *tmpRA = 0.0e0; tmpRA++;           
1219         
1220         *tmpRA = 0.0e0; tmpRA++;           
1221         *tmpRA = 0.0e0; tmpRA++;           
1222         *tmpRA = 0.0e0; tmpRA++;           
1223         
1224         *tmpRA = 0.0e0; tmpRA++;           
1225         *tmpRA = 0.0e0; tmpRA++;           
1226         *tmpRA = 0.0e0; tmpRA++;           
1227       }
1228       
1229       for (jj = Degree ; jj >  0 ; jj--) {
1230         tmpPA -= 12;
1231         
1232         Index1 = MaxIndex1;
1233         Index2 = MaxIndex2;
1234         
1235         for (ii = LocalRequest ; ii > 0 ; ii--) {
1236           valRA  = &RA[Index1];
1237           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1238           Index1 -= 12;
1239           Index2 -= 12;
1240         }
1241         valRA  = &RA[Index1];
1242         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1243         
1244         Index1 = MaxIndex1 + 1;
1245         Index2 = MaxIndex2 + 1;
1246         
1247         for (ii = LocalRequest ; ii > 0 ; ii--) {
1248           valRA  = &RA[Index1];
1249           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1250           Index1 -= 12;
1251           Index2 -= 12;
1252         }
1253         valRA  = &RA[Index1];
1254         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1255         
1256         Index1 = MaxIndex1 + 2;
1257         Index2 = MaxIndex2 + 2;
1258         
1259         for (ii = LocalRequest ; ii > 0 ; ii--) {
1260           valRA  = &RA[Index1];
1261           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1262           Index1 -= 12;
1263           Index2 -= 12;
1264         }
1265         valRA  = &RA[Index1];
1266         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1267         
1268         Index1 = MaxIndex1 + 3;
1269         Index2 = MaxIndex2 + 3;
1270         
1271         for (ii = LocalRequest ; ii > 0 ; ii--) {
1272           valRA  = &RA[Index1];
1273           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1274           Index1 -= 12;
1275           Index2 -= 12;
1276         }
1277         valRA  = &RA[Index1];
1278         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1279         
1280         Index1 = MaxIndex1 + 4;
1281         Index2 = MaxIndex2 + 4;
1282         
1283         for (ii = LocalRequest ; ii > 0 ; ii--) {
1284           valRA  = &RA[Index1];
1285           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1286           Index1 -= 12;
1287           Index2 -= 12;
1288         }
1289         valRA  = &RA[Index1];
1290         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1291         
1292         Index1 = MaxIndex1 + 5;
1293         Index2 = MaxIndex2 + 5;
1294         
1295         for (ii = LocalRequest ; ii > 0 ; ii--) {
1296           valRA  = &RA[Index1];
1297           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1298           Index1 -= 12;
1299           Index2 -= 12;
1300         }
1301         valRA  = &RA[Index1];
1302         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1303         
1304         Index1 = MaxIndex1 + 6;
1305         Index2 = MaxIndex2 + 6;
1306         
1307         for (ii = LocalRequest ; ii > 0 ; ii--) {
1308           valRA  = &RA[Index1];
1309           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1310           Index1 -= 12;
1311           Index2 -= 12;
1312         }
1313         valRA  = &RA[Index1];
1314         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1315         
1316         Index1 = MaxIndex1 + 7;
1317         Index2 = MaxIndex2 + 7;
1318         
1319         for (ii = LocalRequest ; ii > 0 ; ii--) {
1320           valRA  = &RA[Index1];
1321           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1322           Index1 -= 12;
1323           Index2 -= 12;
1324         }
1325         valRA  = &RA[Index1];
1326         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1327         
1328         Index1 = MaxIndex1 + 8;
1329         Index2 = MaxIndex2 + 8;
1330         
1331         for (ii = LocalRequest ; ii > 0 ; ii--) {
1332           valRA  = &RA[Index1];
1333           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1334           Index1 -= 12;
1335           Index2 -= 12;
1336         }
1337         valRA  = &RA[Index1];
1338         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1339         
1340         Index1 = MaxIndex1 + 9;
1341         Index2 = MaxIndex2 + 9;
1342         
1343         for (ii = LocalRequest ; ii > 0 ; ii--) {
1344           valRA  = &RA[Index1];
1345           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1346           Index1 -= 12;
1347           Index2 -= 12;
1348         }
1349         valRA  = &RA[Index1];
1350         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1351         
1352         Index1 = MaxIndex1 + 10;
1353         Index2 = MaxIndex2 + 10;
1354         
1355         for (ii = LocalRequest ; ii > 0 ; ii--) {
1356           valRA  = &RA[Index1];
1357           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1358           Index1 -= 12;
1359           Index2 -= 12;
1360         }
1361         valRA  = &RA[Index1];
1362         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1363         
1364         Index1 = MaxIndex1 + 11;
1365         Index2 = MaxIndex2 + 11;
1366         
1367         for (ii = LocalRequest ; ii > 0 ; ii--) {
1368           valRA  = &RA[Index1];
1369           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1370           Index1 -= 12;
1371           Index2 -= 12;
1372         }
1373         valRA  = &RA[Index1];
1374         *valRA = Par * (*valRA) + (*tmpPA);
1375         
1376         tmpPA -= 11;
1377       }
1378     }
1379     else {
1380       
1381       for (jj = Degree  ; jj >  0 ; jj--) {
1382         tmpPA -= 12;
1383         tmpRA  = RA;
1384
1385         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1386         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1387         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1388
1389         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1390         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1391         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1392
1393         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1394         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1395         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1396
1397         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1398         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1399         *tmpRA = Par * (*tmpRA) + (*tmpPA);
1400         tmpPA -= 11;
1401       }
1402     }
1403     break;
1404     break;
1405   }
1406   case 15 : {
1407     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1408     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1409     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1410
1411     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1412     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1413     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1414
1415     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1416     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1417     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1418
1419     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1420     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1421     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1422
1423     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1424     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1425     *tmpRA = *tmpPA; tmpRA++;
1426     tmpPA -= 14;
1427     if (DerivativeRequest > 0 ) {
1428       Standard_Real *valRA;
1429       Standard_Integer ii, LocalRequest;
1430       Standard_Integer Index1, Index2;
1431       Standard_Integer MaxIndex1, MaxIndex2;
1432       if      (DerivativeRequest < Degree) {
1433         LocalRequest = DerivativeRequest;
1434         MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest;
1435       }
1436       else {
1437         LocalRequest = Degree;
1438         MaxIndex2 = MaxIndex1 = DegreeDimension;
1439       }
1440       MaxIndex2 -= 15;
1441       
1442       for (ii = 1; ii <= LocalRequest; ii++) {
1443         *tmpRA = 0.0e0; tmpRA++;           
1444         *tmpRA = 0.0e0; tmpRA++;           
1445         *tmpRA = 0.0e0; tmpRA++;           
1446
1447         *tmpRA = 0.0e0; tmpRA++;           
1448         *tmpRA = 0.0e0; tmpRA++;           
1449         *tmpRA = 0.0e0; tmpRA++;           
1450
1451         *tmpRA = 0.0e0; tmpRA++;           
1452         *tmpRA = 0.0e0; tmpRA++;           
1453         *tmpRA = 0.0e0; tmpRA++;           
1454
1455         *tmpRA = 0.0e0; tmpRA++;           
1456         *tmpRA = 0.0e0; tmpRA++;           
1457         *tmpRA = 0.0e0; tmpRA++;           
1458
1459         *tmpRA = 0.0e0; tmpRA++;           
1460         *tmpRA = 0.0e0; tmpRA++;           
1461         *tmpRA = 0.0e0; tmpRA++;           
1462       }
1463       
1464       for (jj = Degree ; jj >  0 ; jj--) {
1465         tmpPA -= 15;
1466         
1467         Index1 = MaxIndex1;
1468         Index2 = MaxIndex2;
1469         
1470         for (ii = LocalRequest ; ii > 0 ; ii--) {
1471           valRA  = &RA[Index1];
1472           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1473           Index1 -= 15;
1474           Index2 -= 15;
1475         }
1476         valRA  = &RA[Index1];
1477         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1478         
1479         Index1 = MaxIndex1 + 1;
1480         Index2 = MaxIndex2 + 1;
1481         
1482         for (ii = LocalRequest ; ii > 0 ; ii--) {
1483           valRA  = &RA[Index1];
1484           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1485           Index1 -= 15;
1486           Index2 -= 15;
1487         }
1488         valRA  = &RA[Index1];
1489         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1490         
1491         Index1 = MaxIndex1 + 2;
1492         Index2 = MaxIndex2 + 2;
1493         
1494         for (ii = LocalRequest ; ii > 0 ; ii--) {
1495           valRA  = &RA[Index1];
1496           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1497           Index1 -= 15;
1498           Index2 -= 15;
1499         }
1500         valRA  = &RA[Index1];
1501         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1502         
1503         Index1 = MaxIndex1 + 3;
1504         Index2 = MaxIndex2 + 3;
1505         
1506         for (ii = LocalRequest ; ii > 0 ; ii--) {
1507           valRA  = &RA[Index1];
1508           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1509           Index1 -= 15;
1510           Index2 -= 15;
1511         }
1512         valRA  = &RA[Index1];
1513         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1514         
1515         Index1 = MaxIndex1 + 4;
1516         Index2 = MaxIndex2 + 4;
1517         
1518         for (ii = LocalRequest ; ii > 0 ; ii--) {
1519           valRA  = &RA[Index1];
1520           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1521           Index1 -= 15;
1522           Index2 -= 15;
1523         }
1524         valRA  = &RA[Index1];
1525         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1526         
1527         Index1 = MaxIndex1 + 5;
1528         Index2 = MaxIndex2 + 5;
1529         
1530         for (ii = LocalRequest ; ii > 0 ; ii--) {
1531           valRA  = &RA[Index1];
1532           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1533           Index1 -= 15;
1534           Index2 -= 15;
1535         }
1536         valRA  = &RA[Index1];
1537         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1538         
1539         Index1 = MaxIndex1 + 6;
1540         Index2 = MaxIndex2 + 6;
1541         
1542         for (ii = LocalRequest ; ii > 0 ; ii--) {
1543           valRA  = &RA[Index1];
1544           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1545           Index1 -= 15;
1546           Index2 -= 15;
1547         }
1548         valRA  = &RA[Index1];
1549         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1550         
1551         Index1 = MaxIndex1 + 7;
1552         Index2 = MaxIndex2 + 7;
1553         
1554         for (ii = LocalRequest ; ii > 0 ; ii--) {
1555           valRA  = &RA[Index1];
1556           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1557           Index1 -= 15;
1558           Index2 -= 15;
1559         }
1560         valRA  = &RA[Index1];
1561         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1562         
1563         Index1 = MaxIndex1 + 8;
1564         Index2 = MaxIndex2 + 8;
1565         
1566         for (ii = LocalRequest ; ii > 0 ; ii--) {
1567           valRA  = &RA[Index1];
1568           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1569           Index1 -= 15;
1570           Index2 -= 15;
1571         }
1572         valRA  = &RA[Index1];
1573         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1574         
1575         Index1 = MaxIndex1 + 9;
1576         Index2 = MaxIndex2 + 9;
1577         
1578         for (ii = LocalRequest ; ii > 0 ; ii--) {
1579           valRA  = &RA[Index1];
1580           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1581           Index1 -= 15;
1582           Index2 -= 15;
1583         }
1584         valRA  = &RA[Index1];
1585         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1586         
1587         Index1 = MaxIndex1 + 10;
1588         Index2 = MaxIndex2 + 10;
1589         
1590         for (ii = LocalRequest ; ii > 0 ; ii--) {
1591           valRA  = &RA[Index1];
1592           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1593           Index1 -= 15;
1594           Index2 -= 15;
1595         }
1596         valRA  = &RA[Index1];
1597         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1598         
1599         Index1 = MaxIndex1 + 11;
1600         Index2 = MaxIndex2 + 11;
1601         
1602         for (ii = LocalRequest ; ii > 0 ; ii--) {
1603           valRA  = &RA[Index1];
1604           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1605           Index1 -= 15;
1606           Index2 -= 15;
1607         }
1608         valRA  = &RA[Index1];
1609         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1610         
1611         Index1 = MaxIndex1 + 12;
1612         Index2 = MaxIndex2 + 12;
1613         
1614         for (ii = LocalRequest ; ii > 0 ; ii--) {
1615           valRA  = &RA[Index1];
1616           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1617           Index1 -= 15;
1618           Index2 -= 15;
1619         }
1620         valRA  = &RA[Index1];
1621         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1622         
1623         Index1 = MaxIndex1 + 13;
1624         Index2 = MaxIndex2 + 13;
1625         
1626         for (ii = LocalRequest ; ii > 0 ; ii--) {
1627           valRA  = &RA[Index1];
1628           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1629           Index1 -= 15;
1630           Index2 -= 15;
1631         }
1632         valRA  = &RA[Index1];
1633         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1634         
1635         Index1 = MaxIndex1 + 14;
1636         Index2 = MaxIndex2 + 14;
1637         
1638         for (ii = LocalRequest ; ii > 0 ; ii--) {
1639           valRA  = &RA[Index1];
1640           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1641           Index1 -= 15;
1642           Index2 -= 15;
1643         }
1644         valRA  = &RA[Index1];
1645         *valRA = Par * (*valRA) + (*tmpPA);
1646         
1647         tmpPA -= 14;
1648       }
1649     }
1650     else {
1651       
1652       for (jj = Degree  ; jj >  0 ; jj--) {
1653         tmpPA -= 15;
1654         tmpRA  = RA;
1655
1656         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1657         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1658         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1659
1660         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1661         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1662         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1663
1664         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1665         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1666         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1667
1668         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1669         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1670         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1671
1672         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1673         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1674         *tmpRA = Par * (*tmpRA) + (*tmpPA);
1675         tmpPA -= 14;
1676       }
1677     }
1678     break;
1679   }
1680     default : {
1681       Standard_Integer kk ;
1682       for ( kk = 0 ; kk < Dimension ; kk++) {
1683         *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1684       }
1685       tmpPA -= Dimension;
1686       if (DerivativeRequest > 0 ) {
1687         Standard_Real *valRA;
1688         Standard_Integer ii, LocalRequest;
1689         Standard_Integer Index1, Index2;
1690         Standard_Integer MaxIndex1, MaxIndex2;
1691         if      (DerivativeRequest < Degree) {
1692           LocalRequest = DerivativeRequest;
1693           MaxIndex2 = MaxIndex1 = LocalRequest * Dimension;
1694         }
1695         else {
1696           LocalRequest = Degree;
1697           MaxIndex2 = MaxIndex1 = DegreeDimension;
1698         }
1699         MaxIndex2 -= Dimension;
1700         
1701         for (ii = 1; ii <= MaxIndex1; ii++) {
1702           *tmpRA = 0.0e0; tmpRA++;           
1703         }
1704         
1705         for (jj = Degree ; jj >  0 ; jj--) {
1706           tmpPA -= Dimension;
1707           
1708           for (kk = 0 ; kk < Dimension ; kk++) {
1709             Index1 = MaxIndex1 + kk;
1710             Index2 = MaxIndex2 + kk;
1711             
1712             for (ii = LocalRequest ; ii > 0 ; ii--) {
1713               valRA  = &RA[Index1];
1714               *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1715               Index1 -= Dimension;
1716               Index2 -= Dimension;
1717             }
1718             valRA  = &RA[Index1];
1719             *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1720           }
1721           tmpPA -= Dimension;
1722         }
1723       }
1724       else {
1725         
1726         for (jj = Degree  ; jj >  0 ; jj--) {
1727           tmpPA -= Dimension;
1728           tmpRA  = RA;
1729           
1730           for (kk = 0 ; kk < Dimension ; kk++) {
1731             *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1732           }
1733           tmpPA -= Dimension;
1734         }
1735       }
1736     }
1737   }
1738 }
1739
1740 //=======================================================================
1741 //function : This evaluates a polynomial without derivative 
1742 //purpose  :
1743 //=======================================================================
1744
1745 void  PLib::NoDerivativeEvalPolynomial(const Standard_Real    Par,
1746                                        const Standard_Integer Degree,
1747                                        const Standard_Integer Dimension, 
1748                                        const Standard_Integer DegreeDimension, 
1749                                        Standard_Real&         PolynomialCoeff,
1750                                        Standard_Real&         Results)
1751 {
1752   Standard_Integer jj;
1753   Standard_Real *RA = &Results ;  
1754   Standard_Real *PA = &PolynomialCoeff ;
1755   Standard_Real *tmpRA = RA;
1756   Standard_Real *tmpPA = PA + DegreeDimension;
1757
1758   switch (Dimension) {
1759   
1760   case 1 : {
1761     *tmpRA = *tmpPA;
1762     
1763     for (jj = Degree  ; jj >  0 ; jj--) {
1764       tmpPA--;
1765
1766       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1767     }
1768     break;
1769   }
1770   case 2 : {
1771     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1772     *tmpRA = *tmpPA;
1773     tmpPA--;
1774     
1775     for (jj = Degree  ; jj >  0 ; jj--) {
1776       tmpPA -= 2;
1777       tmpRA  = RA;
1778
1779       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1780       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1781       tmpPA--;
1782     }
1783     break;
1784   }
1785   case 3 : {
1786     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1787     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1788     *tmpRA = *tmpPA;
1789     tmpPA -= 2;
1790     
1791     for (jj = Degree  ; jj >  0 ; jj--) {
1792       tmpPA -= 3;
1793       tmpRA  = RA;
1794
1795       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1796       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1797       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1798       tmpPA -= 2;
1799     }
1800     break;
1801   }
1802   case 6 : {
1803     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1804     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1805     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1806
1807     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1808     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1809     *tmpRA = *tmpPA;
1810     tmpPA -= 5;
1811     
1812     for (jj = Degree  ; jj >  0 ; jj--) {
1813       tmpPA -= 6;
1814       tmpRA  = RA;
1815
1816       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1817       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1818       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1819
1820       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1821       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1822       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1823       tmpPA -= 5;
1824     }
1825     break;
1826   }
1827   case 9 : {
1828     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1829     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1830     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1831
1832     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1833     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1834     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1835
1836     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1837     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1838     *tmpRA = *tmpPA;
1839     tmpPA -= 8;
1840     
1841     for (jj = Degree  ; jj >  0 ; jj--) {
1842       tmpPA -= 9;
1843       tmpRA  = RA;
1844
1845       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1846       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1847       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1848
1849       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1850       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1851       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1852
1853       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1854       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1855       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1856       tmpPA -= 8;
1857     }
1858     break;
1859   }
1860   case 12 : {
1861     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1862     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1863     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1864
1865     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1866     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1867     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1868
1869     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1870     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1871     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1872
1873     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1874     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1875     *tmpRA = *tmpPA;
1876     tmpPA -= 11;
1877     
1878     for (jj = Degree  ; jj >  0 ; jj--) {
1879       tmpPA -= 12;
1880       tmpRA  = RA;
1881
1882       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1883       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1884       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1885
1886       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1887       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1888       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1889
1890       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1891       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1892       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1893
1894       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1895       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1896       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1897       tmpPA -= 11;
1898     }
1899     break;
1900   }
1901   case 15 : {
1902     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1903     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1904     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1905
1906     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1907     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1908     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1909
1910     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1911     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1912     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1913
1914     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1915     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1916     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1917
1918     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1919     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1920     *tmpRA = *tmpPA;
1921     tmpPA -= 14;
1922     
1923     for (jj = Degree  ; jj >  0 ; jj--) {
1924       tmpPA -= 15;
1925       tmpRA  = RA;
1926
1927       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1928       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1929       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1930
1931       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1932       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1933       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1934
1935       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1936       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1937       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1938
1939       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1940       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1941       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1942
1943       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1944       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1945       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1946       tmpPA -= 14;
1947     }
1948     break;
1949   }
1950     default : {
1951       Standard_Integer kk ;
1952       for ( kk = 0 ; kk < Dimension ; kk++) {
1953         *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1954       }
1955       tmpPA -= Dimension;
1956       
1957       for (jj = Degree  ; jj >  0 ; jj--) {
1958         tmpPA -= Dimension;
1959         tmpRA  = RA;
1960         
1961         for (kk = 0 ; kk < Dimension ; kk++) {
1962           *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1963         }
1964         tmpPA -= Dimension;
1965       }
1966     }
1967   }
1968 }
1969
1970 //=======================================================================
1971 //function : This evaluates a polynomial of 2 variables 
1972 //purpose  : or its derivative at the requested orders 
1973 //=======================================================================
1974
1975 void  PLib::EvalPoly2Var(const Standard_Real    UParameter,
1976                          const Standard_Real    VParameter,
1977                          const Standard_Integer UDerivativeRequest,
1978                          const Standard_Integer VDerivativeRequest,
1979                          const Standard_Integer UDegree,
1980                          const Standard_Integer VDegree,
1981                          const Standard_Integer Dimension, 
1982                          Standard_Real&         PolynomialCoeff,
1983                          Standard_Real&         Results)
1984      //
1985      // the polynomial coefficients are assumed to be stored as follows :
1986      //                                                      0 0   
1987      // [0]                  [Dimension -1]                 U V    coefficient
1988      //                                                      1 0
1989      // [Dimension]          [Dimension + Dimension -1]     U V    coefficient
1990      //                                                      2 0
1991      // [2 * Dimension]      [2 * Dimension + Dimension-1]  U V    coefficient
1992      //
1993      //   ...................................................
1994      //
1995      //  
1996      //                                                      m 0
1997      // [m * Dimension]      [m * Dimension + Dimension-1]  U V    coefficient
1998      //
1999      //  where m = UDegree
2000      //
2001      //                                                           0 1
2002      // [(m+1) * Dimension]  [(m+1) * Dimension + Dimension-1]   U V   coefficient
2003      //
2004      //   ...................................................
2005      //
2006      //                                                          m 1
2007      // [2*m * Dimension]      [2*m * Dimension + Dimension-1]  U V    coefficient
2008      //  
2009      //   ...................................................
2010      //
2011      //                                                          m n
2012      // [m*n * Dimension]      [m*n * Dimension + Dimension-1]  U V    coefficient
2013      //
2014      //  where n = VDegree
2015 {
2016   Standard_Integer Udim = (VDegree+1)*Dimension,
2017   index = Udim*UDerivativeRequest;
2018   TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
2019   TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1)); 
2020   Standard_Real * Result =  (Standard_Real *) &Curve.ChangeValue(1);
2021   Standard_Real * Digit  =  (Standard_Real *) &Point.ChangeValue(1);
2022   Standard_Real * ResultArray ;
2023   ResultArray = &Results ;
2024   
2025   PLib::EvalPolynomial(UParameter,
2026                        UDerivativeRequest,
2027                        UDegree,
2028                        Udim,
2029                        PolynomialCoeff,
2030                        Result[0]);
2031   
2032   PLib::EvalPolynomial(VParameter,
2033                        VDerivativeRequest,
2034                        VDegree,  
2035                        Dimension,
2036                        Result[index],
2037                        Digit[0]);
2038
2039   index = Dimension*VDerivativeRequest;
2040
2041   for (Standard_Integer i=0;i<Dimension;i++) {
2042     ResultArray[i] = Digit[index+i];
2043   }
2044 }
2045
2046
2047
2048 //=======================================================================
2049 //function : This evaluates the lagrange polynomial and its derivatives 
2050 //purpose  : up to the requested order that interpolates a series of
2051 //points of dimension <Dimension> with given assigned parameters
2052 //=======================================================================
2053
2054 Standard_Integer  
2055 PLib::EvalLagrange(const Standard_Real                   Parameter,
2056                    const Standard_Integer                DerivativeRequest,
2057                    const Standard_Integer                Degree,
2058                    const Standard_Integer                Dimension, 
2059                    Standard_Real&                        Values,
2060                    Standard_Real&                        Parameters,
2061                    Standard_Real&                        Results)
2062 {
2063   //
2064   // the points  are assumed to be stored as follows in the Values array :
2065   //                                                      
2066   // [0]              [Dimension -1]                first  point   coefficients
2067   //                                                       
2068   // [Dimension]      [Dimension + Dimension -1]    second point   coefficients
2069   //                                                     
2070   // [2 * Dimension]  [2 * Dimension + Dimension-1] third  point   coefficients
2071   //
2072   //   ...................................................
2073   //
2074   //  
2075   //                                                      
2076   // [d * Dimension]  [d * Dimension + Dimension-1] d + 1 point   coefficients
2077   //
2078   //  where d is the Degree
2079   // 
2080   //  The ParameterArray stores the parameter value assign to each point in
2081   //  order described above, that is 
2082   //  [0]   is assign to first  point
2083   //  [1]   is assign to second point
2084   //
2085   Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
2086   Standard_Integer local_request = DerivativeRequest;
2087   Standard_Real  *ParameterArray;
2088   Standard_Real  difference;
2089   Standard_Real  *PointsArray;
2090   Standard_Real  *ResultArray ;
2091   
2092   PointsArray    = &Values ;
2093   ParameterArray = &Parameters ;
2094   ResultArray = &Results ;
2095   if (local_request >= Degree) {
2096     local_request = Degree ;
2097   }
2098   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
2099   //
2100   //  Build the divided differences array
2101   //
2102   
2103   for (ii = 0 ; ii <  (Degree + 1) * Dimension  ; ii++) {
2104     divided_differences_array[ii] = PointsArray[ii] ;  
2105   }
2106
2107   for (ii = Degree ; ii >= 0   ; ii--) {
2108
2109     for (jj = Degree  ; jj > Degree - ii  ; jj--) {
2110       Index = jj * Dimension ;
2111       Index1 = Index - Dimension ;
2112
2113       for (kk = 0 ; kk < Dimension ; kk++) {
2114         divided_differences_array[Index + kk] -=
2115           divided_differences_array[Index1 + kk] ;
2116       }
2117       difference = 
2118         ParameterArray[jj] - ParameterArray[jj - Degree -1 +  ii] ; 
2119       if (Abs(difference) < RealSmall()) {
2120         ReturnCode = 1 ;
2121         goto FINISH ;
2122       }
2123       difference = 1.0e0 / difference ;
2124
2125       for (kk = 0 ; kk < Dimension ; kk++) {
2126         divided_differences_array[Index + kk] *= difference ;
2127       }
2128     }
2129   }
2130   //
2131   //
2132   // Evaluate the divided difference array polynomial which expresses as 
2133   //
2134   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2135   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2136   //
2137   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2138   //
2139   // 
2140   Index = Degree * Dimension ;
2141
2142   for (kk = 0 ; kk < Dimension ; kk++) {
2143     ResultArray[kk] = divided_differences_array[Index + kk] ;
2144   }  
2145
2146   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
2147     ResultArray[ii] = 0.0e0 ;
2148   }
2149
2150   for (ii = Degree ; ii >= 1 ; ii--) {
2151     difference =  Parameter - ParameterArray[ii - 1] ;
2152
2153     for (jj = local_request ; jj > 0 ; jj--) {
2154       Index = jj * Dimension ;
2155       Index1 = Index - Dimension ;
2156       
2157       for (kk = 0 ; kk < Dimension ; kk++) {
2158         ResultArray[Index + kk] *= difference ;
2159         ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
2160       }
2161     }
2162     Index = (ii -1) * Dimension ;
2163
2164     for (kk = 0 ; kk < Dimension ; kk++) {
2165       ResultArray[kk] *= difference ;
2166       ResultArray[kk] += divided_differences_array[Index+kk] ;
2167     }
2168   }
2169   FINISH : 
2170     return (ReturnCode) ;
2171 }
2172
2173 //=======================================================================
2174 //function : This evaluates the hermite polynomial and its derivatives 
2175 //purpose  : up to the requested order that interpolates a series of
2176 //points of dimension <Dimension> with given assigned parameters
2177 //=======================================================================
2178
2179 Standard_Integer PLib::EvalCubicHermite
2180 (const Standard_Real          Parameter,
2181  const Standard_Integer       DerivativeRequest,
2182  const Standard_Integer       Dimension, 
2183  Standard_Real&               Values,
2184  Standard_Real&               Derivatives,
2185  Standard_Real&               theParameters,
2186  Standard_Real&               Results)
2187 {
2188   //
2189   // the points  are assumed to be stored as follows in the Values array :
2190   //                                                      
2191   // [0]            [Dimension -1]             first  point   coefficients
2192   //                                                       
2193   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
2194   // 
2195   //
2196   // the derivatives  are assumed to be stored as follows in 
2197   // the Derivatives array :
2198   //                                                      
2199   // [0]            [Dimension -1]             first  point   coefficients
2200   //                                                       
2201   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
2202   // 
2203   //  The ParameterArray stores the parameter value assign to each point in
2204   //  order described above, that is 
2205   //  [0]   is assign to first  point
2206   //  [1]   is assign to last   point
2207   //
2208   Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
2209   Standard_Integer local_request = DerivativeRequest ;
2210   
2211   ReturnCode = 0 ;
2212   Degree = 3 ;
2213   Standard_Real  ParametersArray[4];
2214   Standard_Real  difference;
2215   Standard_Real  inverse;
2216   Standard_Real  *FirstLast;
2217   Standard_Real  *PointsArray;
2218   Standard_Real  *DerivativesArray;
2219   Standard_Real  *ResultArray ;
2220
2221   DerivativesArray = &Derivatives ;
2222   PointsArray    = &Values ;
2223   FirstLast = &theParameters ;
2224   ResultArray = &Results ;
2225   if (local_request >= Degree) {
2226     local_request = Degree ;
2227   }   
2228   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
2229
2230   for (ii = 0, jj = 0  ; ii < 2 ; ii++, jj+= 2) {
2231     ParametersArray[jj] =
2232       ParametersArray[jj+1] = FirstLast[ii] ;
2233   }
2234   //
2235   //  Build the divided differences array
2236   //
2237   //
2238   // initialise it at the stage 2 of the building algorithm
2239   // for devided differences
2240   //
2241   inverse = FirstLast[1] - FirstLast[0] ;
2242   inverse = 1.0e0 / inverse ;
2243
2244   for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension  ; 
2245        ii <  Dimension  ; 
2246        ii++, jj++, kk++, pp++) {
2247     divided_differences_array[ii] = PointsArray[ii] ;  
2248     divided_differences_array[kk] = inverse * 
2249       (PointsArray[jj] - PointsArray[ii]) ; 
2250     divided_differences_array[jj] = DerivativesArray[ii] ;
2251     divided_differences_array[pp] = DerivativesArray[jj] ;
2252   }
2253
2254   for (ii = 1 ; ii <= Degree   ; ii++) {
2255
2256     for (jj = Degree  ; jj >=  ii+1  ; jj--) {
2257       Index = jj * Dimension ;
2258       Index1 = Index - Dimension ;
2259
2260       for (kk = 0 ; kk < Dimension ; kk++) {
2261         divided_differences_array[Index + kk] -=
2262           divided_differences_array[Index1 + kk] ;
2263       }
2264
2265       for (kk = 0 ; kk < Dimension ; kk++) {
2266         divided_differences_array[Index + kk] *= inverse ;
2267       }
2268     }
2269   }
2270   //
2271   //
2272   // Evaluate the divided difference array polynomial which expresses as 
2273   //
2274   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2275   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2276   //
2277   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2278   //
2279   // 
2280   Index = Degree * Dimension ;
2281
2282   for (kk = 0 ; kk < Dimension ; kk++) {
2283     ResultArray[kk] = divided_differences_array[Index + kk] ;
2284   }  
2285
2286   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
2287     ResultArray[ii] = 0.0e0 ;
2288   }
2289
2290   for (ii = Degree ; ii >= 1 ; ii--) {
2291     difference =  Parameter - ParametersArray[ii - 1] ;
2292
2293     for (jj = local_request ; jj > 0 ; jj--) {
2294       Index = jj * Dimension ;
2295       Index1 = Index - Dimension ;
2296
2297       for (kk = 0 ; kk < Dimension ; kk++) {
2298         ResultArray[Index + kk] *= difference ;
2299         ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
2300       }
2301     }
2302     Index = (ii -1) * Dimension ;
2303
2304     for (kk = 0 ; kk < Dimension ; kk++) {
2305       ResultArray[kk] *= difference ;
2306       ResultArray[kk] += divided_differences_array[Index+kk] ;
2307     }
2308   }
2309 //  FINISH : 
2310     return (ReturnCode) ;
2311 }
2312
2313 //=======================================================================
2314 //function : HermiteCoefficients
2315 //purpose  : calcul des polynomes d'Hermite
2316 //=======================================================================
2317
2318 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter, 
2319                                            const Standard_Real LastParameter,
2320                                            const Standard_Integer FirstOrder,
2321                                            const Standard_Integer LastOrder, 
2322                                            math_Matrix& MatrixCoefs)
2323 {
2324   Standard_Integer NbCoeff =  FirstOrder +  LastOrder + 2, Ordre[2];
2325   Standard_Integer ii, jj, pp, cote, iof=0;
2326   Standard_Real Prod, TBorne = FirstParameter;
2327   math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
2328   math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
2329
2330   // Test de validites
2331
2332   if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
2333   Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
2334   if (D1 > 100 || D2 > 100) return Standard_False;
2335   D2 += D1;
2336   if (D2 < 0.01) return Standard_False;
2337   if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False; 
2338
2339   // Calcul de la matrice a inverser (MAT)
2340
2341   Ordre[0] = FirstOrder+1;
2342   Ordre[1] = LastOrder+1;
2343
2344   for (cote=0; cote<=1; cote++) {
2345     Coeff.Init(1);
2346
2347     for (pp=1; pp<=Ordre[cote]; pp++) {
2348       ii = pp + iof;
2349       Prod = 1;
2350
2351       for (jj=pp; jj<=NbCoeff; jj++) {
2352         //        tout se passe dans les 3 lignes suivantes
2353         MAT(ii, jj) = Coeff(jj) * Prod;
2354         Coeff(jj) *= jj - pp;
2355         Prod      *= TBorne;
2356       }
2357     }
2358     TBorne = LastParameter;
2359     iof = Ordre[0];
2360   }
2361
2362   // resolution du systemes
2363   math_Gauss ResolCoeff(MAT, 1.0e-10);
2364   if (!ResolCoeff.IsDone()) return Standard_False;
2365   
2366   for (ii=1; ii<=NbCoeff; ii++) {
2367     B(ii) = 1;
2368     ResolCoeff.Solve(B, Coeff);
2369     MatrixCoefs.SetRow( ii, Coeff);
2370     B(ii) = 0;
2371   }
2372   return Standard_True;      
2373 }
2374
2375 //=======================================================================
2376 //function : CoefficientsPoles
2377 //purpose  : 
2378 //=======================================================================
2379
2380 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt&   Coefs,
2381                               const TColStd_Array1OfReal& WCoefs,
2382                               TColgp_Array1OfPnt&         Poles, 
2383                               TColStd_Array1OfReal&       Weights)
2384 {
2385   TColStd_Array1OfReal tempC(1,3*Coefs.Length());
2386   PLib::SetPoles(Coefs,tempC);
2387   TColStd_Array1OfReal tempP(1,3*Poles.Length());
2388   PLib::SetPoles(Coefs,tempP);
2389   PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
2390   PLib::GetPoles(tempP,Poles);
2391 }
2392
2393 //=======================================================================
2394 //function : CoefficientsPoles
2395 //purpose  : 
2396 //=======================================================================
2397
2398 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
2399                               const TColStd_Array1OfReal& WCoefs,
2400                               TColgp_Array1OfPnt2d&       Poles, 
2401                               TColStd_Array1OfReal&       Weights)
2402 {
2403   TColStd_Array1OfReal tempC(1,2*Coefs.Length());
2404   PLib::SetPoles(Coefs,tempC);
2405   TColStd_Array1OfReal tempP(1,2*Poles.Length());
2406   PLib::SetPoles(Coefs,tempP);
2407   PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
2408   PLib::GetPoles(tempP,Poles);
2409 }
2410
2411 //=======================================================================
2412 //function : CoefficientsPoles
2413 //purpose  : 
2414 //=======================================================================
2415
2416 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
2417                               const TColStd_Array1OfReal& WCoefs,
2418                               TColStd_Array1OfReal&       Poles, 
2419                               TColStd_Array1OfReal&       Weights)
2420 {
2421   PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
2422 }
2423
2424 //=======================================================================
2425 //function : CoefficientsPoles
2426 //purpose  : 
2427 //=======================================================================
2428
2429 void PLib::CoefficientsPoles (const Standard_Integer      dim,
2430                               const TColStd_Array1OfReal& Coefs,
2431                               const TColStd_Array1OfReal& WCoefs,
2432                               TColStd_Array1OfReal&       Poles, 
2433                               TColStd_Array1OfReal&       Weights)
2434 {
2435   Standard_Boolean rat = &WCoefs != NULL;
2436   Standard_Integer loc = Coefs.Lower();
2437   Standard_Integer lop = Poles.Lower();
2438   Standard_Integer lowc=0;
2439   Standard_Integer lowp=0;
2440   Standard_Integer upc = Coefs.Upper();
2441   Standard_Integer upp = Poles.Upper();
2442   Standard_Integer upwc=0;
2443   Standard_Integer upwp=0;
2444   Standard_Integer reflen = Coefs.Length()/dim;
2445   Standard_Integer i,j,k; 
2446   //Les Extremites.
2447   if (rat) { 
2448     lowc = WCoefs.Lower(); lowp = Weights.Lower();
2449     upwc = WCoefs.Upper(); upwp = Weights.Upper();
2450   }
2451
2452   for (i = 0; i < dim; i++){
2453     Poles (lop + i) = Coefs (loc + i);
2454     Poles (upp - i) = Coefs (upc - i);
2455   }
2456   if (rat) {
2457     Weights (lowp) = WCoefs (lowc);
2458     Weights (upwp) = WCoefs (upwc);
2459   }
2460   
2461   Standard_Real Cnp;
2462   for (i = 2; i < reflen; i++ ) {
2463     Cnp = PLib::Bin(reflen - 1, i - 1);
2464     if (rat) Weights (lowp + i - 1) = WCoefs (lowc + i - 1) / Cnp;
2465
2466     for(j = 0; j < dim; j++){
2467       Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
2468     }
2469   }
2470   
2471   for (i = 1; i <= reflen - 1; i++) {
2472
2473     for (j = reflen - 1; j >= i; j--) {
2474       if (rat) Weights (lowp + j) += Weights (lowp + j -1);
2475
2476       for(k = 0; k < dim; k++){
2477         Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
2478       }
2479     }
2480   }
2481   if (rat) {
2482
2483     for (i = 1; i <= reflen; i++) {
2484
2485       for(j = 0; j < dim; j++){
2486         Poles(lop + dim * (i-1) + j) /= Weights(lowp + i -1);
2487       }
2488     }
2489   }
2490 }
2491
2492 //=======================================================================
2493 //function : Trimming
2494 //purpose  : 
2495 //=======================================================================
2496
2497 void PLib::Trimming(const Standard_Real   U1, 
2498                     const Standard_Real   U2, 
2499                     TColgp_Array1OfPnt&   Coefs,
2500                     TColStd_Array1OfReal& WCoefs)
2501 {
2502   TColStd_Array1OfReal temp(1,3*Coefs.Length());
2503   PLib::SetPoles(Coefs,temp);
2504   PLib::Trimming(U1,U2,3,temp,WCoefs);
2505   PLib::GetPoles(temp,Coefs);
2506 }
2507
2508 //=======================================================================
2509 //function : Trimming
2510 //purpose  : 
2511 //=======================================================================
2512
2513 void PLib::Trimming(const Standard_Real   U1, 
2514                     const Standard_Real   U2, 
2515                     TColgp_Array1OfPnt2d& Coefs,
2516                     TColStd_Array1OfReal& WCoefs)
2517 {
2518   TColStd_Array1OfReal temp(1,2*Coefs.Length());
2519   PLib::SetPoles(Coefs,temp);
2520   PLib::Trimming(U1,U2,2,temp,WCoefs);
2521   PLib::GetPoles(temp,Coefs);
2522 }
2523
2524 //=======================================================================
2525 //function : Trimming
2526 //purpose  : 
2527 //=======================================================================
2528
2529 void PLib::Trimming(const Standard_Real   U1, 
2530                     const Standard_Real   U2, 
2531                     TColStd_Array1OfReal& Coefs,
2532                     TColStd_Array1OfReal& WCoefs)
2533 {
2534   PLib::Trimming(U1,U2,1,Coefs,WCoefs);
2535 }
2536
2537 //=======================================================================
2538 //function : Trimming
2539 //purpose  : 
2540 //=======================================================================
2541
2542 void PLib::Trimming(const Standard_Real U1, 
2543                     const Standard_Real U2, 
2544                     const Standard_Integer dim, 
2545                     TColStd_Array1OfReal& Coefs,
2546                     TColStd_Array1OfReal& WCoefs)
2547 {
2548
2549   // principe :
2550   // on fait le changement de variable v = (u-U1) / (U2-U1)
2551   // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
2552   // decomposee sous la forme du schema iteratif de horner.
2553
2554   Standard_Real lsp = U2 - U1;
2555   Standard_Integer indc, indw=0;
2556   Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
2557   Standard_Integer len = Coefs.Length()/dim;
2558   Standard_Boolean rat = &WCoefs != NULL;
2559
2560   if (rat) {
2561     if(len != WCoefs.Length())
2562       Standard_Failure::Raise("PLib::Trimming : nbcoefs/dim != nbweights !!!");
2563     upw = WCoefs.Upper();
2564   }
2565   len --;
2566
2567   for (Standard_Integer i = 1; i <= len; i++) {
2568     Standard_Integer j ;
2569     indc = upc - dim*(i-1);
2570     if (rat) indw = upw - i + 1;
2571     //calcul du coefficient de degre le plus faible a l'iteration i
2572
2573     for( j = 0; j < dim; j++){
2574       Coefs(indc - dim + j) += U1 * Coefs(indc + j);
2575     }
2576     if (rat) WCoefs(indw - 1) += U1 * WCoefs(indw);
2577     
2578     //calcul des coefficients intermediaires :
2579
2580     while (indc < upc){
2581       indc += dim;
2582
2583       for(Standard_Integer k = 0; k < dim; k++){
2584         Coefs(indc - dim + k) = 
2585           U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
2586       }
2587       if (rat) {
2588         indw ++;
2589         WCoefs(indw - 1) = U1 * WCoefs(indw) + lsp * WCoefs(indw - 1);
2590       }
2591     }
2592
2593     //calcul du coefficient de degre le plus eleve :
2594
2595     for(j = 0; j < dim; j++){
2596       Coefs(upc + j) *= lsp;
2597     }
2598     if (rat) WCoefs(upw) *= lsp;
2599   }
2600 }
2601
2602 //=======================================================================
2603 //function : CoefficientsPoles
2604 //purpose  : 
2605 // Modified: 21/10/1996 by PMN :  PolesCoefficient (PRO5852).
2606 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
2607 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
2608 //=======================================================================
2609
2610 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt&   Coefs,
2611                               const TColStd_Array2OfReal& WCoefs,
2612                               TColgp_Array2OfPnt&         Poles,
2613                               TColStd_Array2OfReal&       Weights) 
2614 {
2615   Standard_Boolean rat = (&WCoefs != NULL);
2616   Standard_Integer LowerRow  = Poles.LowerRow();
2617   Standard_Integer UpperRow  = Poles.UpperRow();
2618   Standard_Integer LowerCol  = Poles.LowerCol();
2619   Standard_Integer UpperCol  = Poles.UpperCol();
2620   Standard_Integer ColLength = Poles.ColLength();
2621   Standard_Integer RowLength = Poles.RowLength();
2622
2623   // Bidouille pour retablir u et v pour les coefs calcules 
2624   // par buildcache
2625 //  Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
2626
2627   Standard_Integer Row, Col;
2628   Standard_Real W, Cnp;
2629
2630   Standard_Integer I1, I2;
2631   Standard_Integer NPoleu , NPolev;
2632   gp_XYZ Temp;
2633
2634   for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
2635     Poles (NPoleu, LowerCol) =  Coefs (NPoleu, LowerCol);
2636     if (rat) {
2637       Weights (NPoleu, LowerCol) =  WCoefs (NPoleu, LowerCol);
2638     }
2639
2640     for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
2641       Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
2642       Temp = Coefs (NPoleu, Col).XYZ();
2643       Temp.Divide (Cnp);
2644       Poles (NPoleu, Col).SetXYZ (Temp);
2645       if (rat) {
2646         Weights (NPoleu, Col) = WCoefs (NPoleu, Col) /  Cnp;
2647       }
2648     }
2649     Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
2650     if (rat) {
2651       Weights (NPoleu, UpperCol) = WCoefs (NPoleu, UpperCol);
2652     }
2653
2654     for (I1 = 1; I1 <= RowLength - 1; I1++) {
2655
2656       for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
2657         Temp.SetLinearForm 
2658           (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
2659         Poles (NPoleu, I2).SetXYZ (Temp);
2660         if (rat) Weights(NPoleu, I2) += Weights(NPoleu, I2-1);
2661       }
2662     } 
2663   }
2664   
2665   for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
2666
2667     for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
2668       Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
2669       Temp = Poles (Row, NPolev).XYZ();
2670       Temp.Divide (Cnp);
2671       Poles (Row, NPolev).SetXYZ (Temp);
2672       if (rat) Weights(Row, NPolev) /= Cnp;
2673     }
2674
2675     for (I1 = 1; I1 <= ColLength - 1; I1++) {
2676
2677       for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
2678         Temp.SetLinearForm 
2679           (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
2680         Poles (I2, NPolev).SetXYZ (Temp);
2681         if (rat) Weights(I2, NPolev) += Weights(I2-1, NPolev);
2682       }
2683     }
2684   }
2685   if (rat) {
2686
2687     for (Row = LowerRow; Row <= UpperRow; Row++) {
2688
2689       for (Col = LowerCol; Col <= UpperCol; Col++) {
2690         W = Weights (Row, Col);
2691         Temp = Poles(Row, Col).XYZ();
2692         Temp.Divide (W);
2693         Poles(Row, Col).SetXYZ (Temp);
2694       }
2695     }
2696   }
2697 }  
2698
2699 //=======================================================================
2700 //function : UTrimming
2701 //purpose  : 
2702 //=======================================================================
2703
2704 void PLib::UTrimming(const Standard_Real U1, 
2705                      const Standard_Real U2, 
2706                      TColgp_Array2OfPnt& Coeffs, 
2707                      TColStd_Array2OfReal& WCoeffs)
2708 {
2709   Standard_Boolean rat = &WCoeffs != NULL;
2710   Standard_Integer lr = Coeffs.LowerRow();
2711   Standard_Integer ur = Coeffs.UpperRow();
2712   Standard_Integer lc = Coeffs.LowerCol();
2713   Standard_Integer uc = Coeffs.UpperCol();
2714   TColgp_Array1OfPnt  Temp (lr,ur);
2715   TColStd_Array1OfReal Temw (lr,ur);
2716
2717   for (Standard_Integer icol = lc; icol <= uc; icol++) {
2718     Standard_Integer irow ;
2719     for ( irow = lr; irow <= ur; irow++) {
2720       Temp (irow) = Coeffs  (irow, icol);
2721       if (rat) Temw (irow) = WCoeffs (irow, icol);
2722     }
2723     if (rat) PLib::Trimming (U1, U2, Temp, Temw);
2724     else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
2725
2726     for (irow = lr; irow <= ur; irow++) {
2727       Coeffs  (irow, icol) = Temp (irow);
2728       if (rat) WCoeffs (irow, icol) = Temw (irow);
2729     }
2730   }
2731 }
2732
2733 //=======================================================================
2734 //function : VTrimming
2735 //purpose  : 
2736 //=======================================================================
2737
2738 void PLib::VTrimming(const Standard_Real V1, 
2739                      const Standard_Real V2, 
2740                      TColgp_Array2OfPnt& Coeffs, 
2741                      TColStd_Array2OfReal& WCoeffs)
2742 {
2743   Standard_Boolean rat = &WCoeffs != NULL;
2744   Standard_Integer lr = Coeffs.LowerRow();
2745   Standard_Integer ur = Coeffs.UpperRow();
2746   Standard_Integer lc = Coeffs.LowerCol();
2747   Standard_Integer uc = Coeffs.UpperCol();
2748   TColgp_Array1OfPnt  Temp (lc,uc);
2749   TColStd_Array1OfReal Temw (lc,uc);
2750
2751   for (Standard_Integer irow = lr; irow <= ur; irow++) {
2752     Standard_Integer icol ;
2753     for ( icol = lc; icol <= uc; icol++) {
2754       Temp (icol) = Coeffs  (irow, icol);
2755       if (rat) Temw (icol) = WCoeffs (irow, icol);
2756     }
2757     if (rat) PLib::Trimming (V1, V2, Temp, Temw);
2758     else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
2759
2760     for (icol = lc; icol <= uc; icol++) {
2761       Coeffs  (irow, icol) = Temp (icol);
2762       if (rat) WCoeffs (irow, icol) = Temw (icol);
2763     }
2764   }
2765 }
2766
2767 //=======================================================================
2768 //function : HermiteInterpolate
2769 //purpose  : 
2770 //=======================================================================
2771
2772 Standard_Boolean PLib::HermiteInterpolate
2773 (const Standard_Integer Dimension, 
2774  const Standard_Real FirstParameter, 
2775  const Standard_Real LastParameter, 
2776  const Standard_Integer FirstOrder, 
2777  const Standard_Integer LastOrder, 
2778  const TColStd_Array2OfReal& FirstConstr,
2779  const TColStd_Array2OfReal& LastConstr,
2780  TColStd_Array1OfReal& Coefficients)
2781 {
2782   Standard_Real Pattern[3][6];
2783
2784   // portage HP : il faut les initialiser 1 par 1
2785
2786   Pattern[0][0] = 1;
2787   Pattern[0][1] = 1;
2788   Pattern[0][2] = 1;
2789   Pattern[0][3] = 1;
2790   Pattern[0][4] = 1;
2791   Pattern[0][5] = 1;
2792   Pattern[1][0] = 0;
2793   Pattern[1][1] = 1;
2794   Pattern[1][2] = 2;
2795   Pattern[1][3] = 3;
2796   Pattern[1][4] = 4;
2797   Pattern[1][5] = 5;
2798   Pattern[2][0] = 0;
2799   Pattern[2][1] = 0;
2800   Pattern[2][2] = 2;
2801   Pattern[2][3] = 6;
2802   Pattern[2][4] = 12;
2803   Pattern[2][5] = 20;
2804
2805   math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
2806   //  The initialisation of the matrix A
2807   Standard_Integer irow ;
2808   for ( irow=0; irow<=FirstOrder; irow++) {
2809     Standard_Real FirstVal = 1.;
2810
2811     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2812       A(irow,icol) = Pattern[irow][icol]*FirstVal;
2813       if (irow <= icol) FirstVal *= FirstParameter;
2814     }
2815   }
2816
2817   for (irow=0; irow<=LastOrder; irow++) {
2818     Standard_Real LastVal  = 1.;
2819
2820     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2821       A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
2822       if (irow <= icol) LastVal *= LastParameter;
2823     }
2824   }
2825   //  
2826   //  The filled matrix A for FirstOrder=LastOrder=2 is:  
2827   // 
2828   //      1   FP  FP**2   FP**3    FP**4     FP**5  
2829   //      0   1   2*FP    3*FP**2  4*FP**3   5*FP**4        FP - FirstParameter
2830   //      0   0   2       6*FP     12*FP**2  20*FP**3
2831   //      1   LP  LP**2   LP**3    LP**4     LP**5  
2832   //      0   1   2*LP    3*LP**2  4*LP**3   5*LP**4        LP - LastParameter 
2833   //      0   0   2       6*LP     12*LP**2  20*LP**3
2834   //  
2835   //  If FirstOrder or LastOrder <=2 then some rows and columns are missing. 
2836   //  For example:
2837   //  If FirstOrder=1 then 3th row and 6th column are missing
2838   //  If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
2839
2840   math_Gauss Equations(A);
2841   //  cout << "A=" << A << endl;
2842
2843   for (Standard_Integer idim=1; idim<=Dimension; idim++) {
2844     //  cout << "idim=" << idim << endl;
2845
2846     math_Vector B(0,FirstOrder+LastOrder+1);
2847     Standard_Integer icol ;
2848     for ( icol=0; icol<=FirstOrder; icol++) 
2849       B(icol) = FirstConstr(idim,icol);
2850
2851     for (icol=0; icol<=LastOrder; icol++) 
2852       B(FirstOrder+1+icol) = LastConstr(idim,icol);
2853     //  cout << "B=" << B << endl;
2854
2855     //  The solving of equations system A * X = B. Then B = X
2856     Equations.Solve(B);         
2857     //  cout << "After Solving" << endl << "B=" << B << endl;
2858
2859     if (Equations.IsDone()==Standard_False) return Standard_False;
2860     
2861     //  the filling of the Coefficients
2862
2863     for (icol=0; icol<=FirstOrder+LastOrder+1; icol++) 
2864       Coefficients(Dimension*icol+idim-1) = B(icol);
2865   } 
2866   return Standard_True;
2867 }
2868
2869 //=======================================================================
2870 //function : JacobiParameters
2871 //purpose  : 
2872 //=======================================================================
2873
2874 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder, 
2875                             const Standard_Integer MaxDegree, 
2876                             const Standard_Integer Code, 
2877                             Standard_Integer& NbGaussPoints, 
2878                             Standard_Integer& WorkDegree)
2879 {
2880   // ConstraintOrder: Ordre de contrainte aux extremites :
2881   //            C0 = contraintes de passage aux bornes;
2882   //            C1 = C0 + contraintes de derivees 1eres;
2883   //            C2 = C1 + contraintes de derivees 2ndes.
2884   // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
2885   //            d' approximation (doit etre superieur ou egal a
2886   //            2*NivConstr+2 et inferieur ou egal a 50).
2887   // Code:      Code d' init. des parametres de discretisation.
2888   //            (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
2889   //            = -5 Calcul tres rapide mais peu precis (8pts)
2890   //            = -4    '     '    '      '   '    '    (10pts)
2891   //            = -3    '     '    '      '   '    '    (15pts)
2892   //            = -2    '     '    '      '   '    '    (20pts)
2893   //            = -1    '     '    '      '   '    '    (25pts)
2894   //            = 1 calcul rapide avec precision moyenne (30pts).
2895   //            = 2 calcul rapide avec meilleure precision (40pts).
2896   //            = 3 calcul un peu plus lent avec bonne precision (50 pts).
2897   //            = 4 calcul lent avec la meilleure precision possible
2898   //             (61pts).
2899
2900   // The possible values of NbGaussPoints
2901
2902   const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25, 
2903   NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
2904
2905   Standard_Integer NivConstr=0;
2906   switch (ConstraintOrder) {
2907   case GeomAbs_C0: NivConstr = 0; break;
2908   case GeomAbs_C1: NivConstr = 1; break;
2909   case GeomAbs_C2: NivConstr = 2; break;
2910   default: 
2911     Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2912   }
2913   if (MaxDegree < 2*NivConstr+1)
2914     Standard_ConstructionError::Raise("Invalid MaxDegree");
2915   
2916   if (Code >= 1)
2917     WorkDegree = MaxDegree + 9;
2918   else 
2919     WorkDegree = MaxDegree + 6;
2920   
2921   //---> Nbre mini de points necessaires.
2922   Standard_Integer IPMIN=0;
2923   if (WorkDegree < NDEG8) 
2924     IPMIN=NDEG8;
2925   else if (WorkDegree < NDEG10)
2926     IPMIN=NDEG10;
2927   else if (WorkDegree < NDEG15) 
2928     IPMIN=NDEG15;
2929   else if (WorkDegree < NDEG20) 
2930     IPMIN=NDEG20;
2931   else if (WorkDegree < NDEG25)
2932     IPMIN=NDEG25;
2933   else if (WorkDegree < NDEG30) 
2934     IPMIN=NDEG30;
2935   else if (WorkDegree < NDEG40) 
2936     IPMIN=NDEG40;
2937   else if (WorkDegree < NDEG50) 
2938     IPMIN=NDEG50;
2939   else if (WorkDegree < NDEG61) 
2940     IPMIN=NDEG61;
2941   else
2942     Standard_ConstructionError::Raise("Invalid MaxDegree");
2943   // ---> Nbre de points voulus.
2944   Standard_Integer IWANT=0;
2945   switch (Code) {
2946   case -5: IWANT=NDEG8;  break;
2947   case -4: IWANT=NDEG10; break;
2948   case -3: IWANT=NDEG15; break;
2949   case -2: IWANT=NDEG20; break;
2950   case -1: IWANT=NDEG25; break;
2951   case  1: IWANT=NDEG30; break;
2952   case  2: IWANT=NDEG40; break;
2953   case  3: IWANT=NDEG50; break;
2954   case  4: IWANT=NDEG61; break;
2955   default: 
2956     Standard_ConstructionError::Raise("Invalid Code");
2957   }      
2958   //-->  NbGaussPoints est le nombre de points de discretisation de la fonction,
2959   //     il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
2960   //     NbGaussPoints doit etre superieur strictement a WorkDegree.
2961   NbGaussPoints = Max(IPMIN,IWANT);
2962   //  NbGaussPoints +=2;
2963 }
2964
2965 //=======================================================================
2966 //function : NivConstr
2967 //purpose  : translates from GeomAbs_Shape to Integer
2968 //=======================================================================
2969
2970  Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder) 
2971 {
2972   Standard_Integer NivConstr=0;
2973   switch (ConstraintOrder) {
2974     case GeomAbs_C0: NivConstr = 0; break;
2975     case GeomAbs_C1: NivConstr = 1; break;
2976     case GeomAbs_C2: NivConstr = 2; break;
2977     default: 
2978       Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2979   }
2980   return NivConstr;
2981 }
2982
2983 //=======================================================================
2984 //function : ConstraintOrder
2985 //purpose  : translates from Integer to GeomAbs_Shape
2986 //=======================================================================
2987
2988  GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr) 
2989 {
2990   GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
2991   switch (NivConstr) {
2992     case 0: ConstraintOrder = GeomAbs_C0; break;
2993     case 1: ConstraintOrder = GeomAbs_C1; break;
2994     case 2: ConstraintOrder = GeomAbs_C2; break;
2995     default: 
2996       Standard_ConstructionError::Raise("Invalid NivConstr");
2997   }
2998   return ConstraintOrder;
2999 }
3000
3001 //=======================================================================
3002 //function : EvalLength
3003 //purpose  : 
3004 //=======================================================================
3005
3006  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension, 
3007                       Standard_Real& PolynomialCoeff,
3008                       const Standard_Real U1, const Standard_Real U2,
3009                       Standard_Real& Length)
3010 {
3011   Standard_Integer i,j,idim, degdim;
3012   Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
3013
3014   Standard_Real *PolynomialArray = &PolynomialCoeff ;
3015
3016   Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
3017
3018   math_Vector GaussPoints(1,NbGaussPoints);
3019   math::GaussPoints(NbGaussPoints,GaussPoints);
3020
3021   math_Vector GaussWeights(1,NbGaussPoints);
3022   math::GaussWeights(NbGaussPoints,GaussWeights);
3023
3024   C1 = (U2 + U1) / 2.;
3025   C2 = (U2 - U1) / 2.;
3026
3027 //-----------------------------------------------------------
3028 //****** Integration - Boucle sur les intervalles de GAUSS **
3029 //-----------------------------------------------------------
3030
3031   Sum = 0;
3032
3033   for (j=1; j<=NbGaussPoints/2; j++) {
3034     // Integration en tenant compte de la symetrie 
3035     Tran  = C2 * GaussPoints(j);
3036     X1 = C1 + Tran;
3037     X2 = C1 - Tran;
3038
3039     //****** Derivation sur la dimension de l'espace **
3040
3041     degdim =  Degree*Dimension;
3042     Der1 = Der2 = 0.;
3043     for (idim=0; idim<Dimension; idim++) {
3044       D1 = D2 = Degree * PolynomialArray [idim + degdim];
3045       for (i=Degree-1; i>=1; i--) {
3046         DD = i * PolynomialArray [idim + i*Dimension];
3047         D1 = D1 * X1 + DD;
3048         D2 = D2 * X2 + DD;
3049       }
3050       Der1 += D1 * D1;
3051       Der2 += D2 * D2;
3052     }
3053
3054     //****** Integration **
3055
3056     Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
3057
3058 //****** Fin de boucle dur les intervalles de GAUSS **
3059   }
3060   Length = Sum;
3061 }
3062
3063
3064 //=======================================================================
3065 //function : EvalLength
3066 //purpose  : 
3067 //=======================================================================
3068
3069  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension, 
3070                       Standard_Real& PolynomialCoeff,
3071                       const Standard_Real U1, const Standard_Real U2,
3072                       const Standard_Real Tol, 
3073                       Standard_Real& Length, Standard_Real& Error)
3074 {
3075   Standard_Integer i;
3076   Standard_Integer NbSubInt = 1,    // Current number of subintervals
3077                    MaxNbIter = 13,  // Max number of iterations
3078                    NbIter    = 1;   // Current number of iterations
3079   Standard_Real    dU,OldLen,LenI;
3080
3081   PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
3082
3083   do {
3084     OldLen = Length;
3085     Length = 0.;
3086     NbSubInt *= 2;
3087     dU = (U2-U1)/NbSubInt;    
3088     for (i=1; i<=NbSubInt; i++) {
3089       PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
3090       Length += LenI;
3091     }
3092     NbIter++;
3093     Error = Abs(OldLen-Length);
3094   }    
3095   while (Error > Tol && NbIter <= MaxNbIter);
3096 }