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