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