0028327: BSplCLib can cause memory corruption in degenerated cases
[occt.git] / src / BSplCLib / BSplCLib.cxx
1 // Created on: 1991-08-09
2 // Created by: JCV
3 // Copyright (c) 1991-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 RLE 9 Sep 1993
18 // pmn : modified 28-01-97  : fixed a mistake in LocateParameter (PRO6973)
19 // pmn : modified 4-11-96   : fixed a mistake in BuildKnots (PRO6124)
20 // pmn : modified 28-Jun-96 : fixed a mistake in AntiBoorScheme
21 // xab : modified 15-Jun-95 : fixed a mistake in IsRational
22 // xab : modified 15-Mar-95 : removed Epsilon comparison in IsRational
23 //                            added RationalDerivatives.
24 // xab : 30-Mar-95 : fixed coupling with lti in RationalDerivatives
25 // xab : 15-Mar-96 : fixed a typo in Eval with extrapolation
26 // jct : 15-Apr-97 : added TangExtendToConstraint
27 // jct : 24-Apr-97 : correction on computation of Tbord and NewFlatKnots
28 //                   in TangExtendToConstraint; Continuity can be equal to 0
29
30 #include <BSplCLib.hxx>
31 #include <ElCLib.hxx>
32 #include <gp_Pnt.hxx>
33 #include <gp_Pnt2d.hxx>
34 #include <gp_Vec.hxx>
35 #include <gp_Vec2d.hxx>
36 #include <math_Matrix.hxx>
37 #include <NCollection_LocalArray.hxx>
38 #include <PLib.hxx>
39 #include <Precision.hxx>
40 #include <Standard_NotImplemented.hxx>
41
42 typedef gp_Pnt Pnt;
43 typedef gp_Vec Vec;
44 typedef TColgp_Array1OfPnt Array1OfPnt;
45 typedef TColStd_Array1OfReal Array1OfReal;
46 typedef TColStd_Array1OfInteger Array1OfInteger;
47
48 //=======================================================================
49 //class : BSplCLib_LocalMatrix
50 //purpose: Auxiliary class optimizing creation of matrix buffer for
51 //         evaluation of bspline (using stack allocation for main matrix)
52 //=======================================================================
53
54 class BSplCLib_LocalMatrix : public math_Matrix 
55 {
56 public:
57   BSplCLib_LocalMatrix (Standard_Integer DerivativeRequest, Standard_Integer Order)
58     : math_Matrix (myBuffer, 1, DerivativeRequest + 1, 1, Order)
59   {
60     Standard_OutOfRange_Raise_if (DerivativeRequest > BSplCLib::MaxDegree() ||
61         Order > BSplCLib::MaxDegree()+1 || BSplCLib::MaxDegree() > 25,
62         "BSplCLib: bspline degree is greater than maximum supported");
63   }
64
65  private:
66   // local buffer, to be sufficient for addressing by index [Degree+1][Degree+1]
67   // (see math_Matrix implementation)
68   Standard_Real myBuffer[27*27];
69 };
70
71 //=======================================================================
72 //function : Hunt
73 //purpose  : 
74 //=======================================================================
75
76 void BSplCLib::Hunt (const Array1OfReal& XX, 
77                      const Standard_Real X,
78                      Standard_Integer&   Ilc)
79 {
80   // replaced by simple dichotomy (RLE)
81   Ilc = XX.Lower();
82   if (XX.Length() <= 1) return;
83   const Standard_Real *px = &XX(Ilc);
84   px -= Ilc;
85
86   if (X < px[Ilc]) {
87     Ilc--;
88     return;
89   }
90   Standard_Integer Ihi = XX.Upper();
91   if (X > px[Ihi]) {
92     Ilc = Ihi + 1;
93     return;
94   }
95   Standard_Integer Im;
96
97   while (Ihi - Ilc != 1) {
98     Im = (Ihi + Ilc) >> 1;
99     if (X > px[Im]) Ilc = Im;
100     else            Ihi = Im;
101   }
102 }
103
104 //=======================================================================
105 //function : FirstUKnotIndex
106 //purpose  : 
107 //=======================================================================
108
109 Standard_Integer BSplCLib::FirstUKnotIndex (const Standard_Integer Degree,
110                                    const TColStd_Array1OfInteger& Mults)
111
112   Standard_Integer Index = Mults.Lower();
113   Standard_Integer SigmaMult = Mults(Index);
114
115   while (SigmaMult <= Degree) {
116     Index++;
117     SigmaMult += Mults (Index);
118   }
119   return Index;
120 }
121
122 //=======================================================================
123 //function : LastUKnotIndex
124 //purpose  : 
125 //=======================================================================
126
127 Standard_Integer BSplCLib::LastUKnotIndex  (const Standard_Integer Degree,
128                                    const Array1OfInteger& Mults) 
129
130    Standard_Integer Index = Mults.Upper();
131    Standard_Integer SigmaMult = Mults(Index);
132
133    while (SigmaMult <= Degree) {
134      Index--;
135      SigmaMult += Mults.Value (Index);
136    }
137    return Index;
138 }
139
140 //=======================================================================
141 //function : FlatIndex
142 //purpose  : 
143 //=======================================================================
144
145 Standard_Integer  BSplCLib::FlatIndex
146   (const Standard_Integer Degree,
147    const Standard_Integer Index,
148    const TColStd_Array1OfInteger& Mults,
149    const Standard_Boolean Periodic)
150 {
151   Standard_Integer i, index = Index;
152   const Standard_Integer MLower = Mults.Lower();
153   const Standard_Integer *pmu = &Mults(MLower);
154   pmu -= MLower;
155
156   for (i = MLower + 1; i <= Index; i++)
157     index += pmu[i] - 1;
158   if ( Periodic)
159     index += Degree;
160   else
161     index += pmu[MLower] - 1;
162   return index;
163 }
164
165 //=======================================================================
166 //function : LocateParameter
167 //purpose  : Processing of nodes with multiplicities
168 //pmn  28-01-97 -> compute eventual of the period.
169 //=======================================================================
170
171 void BSplCLib::LocateParameter
172 (const Standard_Integer          , //Degree,
173  const Array1OfReal&    Knots,
174  const Array1OfInteger& , //Mults,
175  const Standard_Real             U,
176  const Standard_Boolean          IsPeriodic,
177  const Standard_Integer          FromK1,
178  const Standard_Integer          ToK2,
179  Standard_Integer&               KnotIndex,
180  Standard_Real&                  NewU)
181 {
182   Standard_Real uf = 0, ul=1;
183   if (IsPeriodic) {
184     uf = Knots(Knots.Lower());
185     ul = Knots(Knots.Upper());
186   }
187   BSplCLib::LocateParameter(Knots,U,IsPeriodic,FromK1,ToK2,
188                             KnotIndex,NewU, uf, ul);
189 }
190
191 //=======================================================================
192 //function : LocateParameter
193 //purpose  : For plane nodes 
194 //   pmn  28-01-97 -> There is a need of the degre to calculate
195 //   the eventual period
196 //=======================================================================
197
198 void BSplCLib::LocateParameter
199 (const Standard_Integer          Degree,
200  const Array1OfReal&    Knots,
201  const Standard_Real             U,
202  const Standard_Boolean          IsPeriodic,
203  const Standard_Integer          FromK1,
204  const Standard_Integer          ToK2,
205  Standard_Integer&               KnotIndex,
206  Standard_Real&                  NewU)
207
208   if (IsPeriodic)
209     BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
210                               KnotIndex, NewU,
211                               Knots(Knots.Lower() + Degree),
212                               Knots(Knots.Upper() - Degree));
213   else 
214     BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
215                               KnotIndex, NewU,
216                               0.,
217                               1.);
218 }
219
220 //=======================================================================
221 //function : LocateParameter
222 //purpose  : Effective computation
223 // pmn 28-01-97 : Add limits of the period as input argument,  
224 //                as it is imposible to produce them at this level.
225 //=======================================================================
226
227 void BSplCLib::LocateParameter 
228 (const TColStd_Array1OfReal& Knots,
229  const Standard_Real         U,
230  const Standard_Boolean      IsPeriodic,
231  const Standard_Integer      FromK1,
232  const Standard_Integer      ToK2,
233  Standard_Integer&           KnotIndex,
234  Standard_Real&              NewU,
235  const Standard_Real         UFirst,
236  const Standard_Real         ULast)
237 {
238   /*
239   Let Knots are distributed as follows (the array is sorted in ascending order):
240     
241       K1, K1,..., K1, K1, K2, K2,..., K2, K2,..., Kn, Kn,..., Kn
242            M1 times             M2 times             Mn times
243
244   NbKnots = sum(M1+M2+...+Mn)
245   If U <= K1 then KnotIndex should be equal to M1.
246   If U >= Kn then KnotIndex should be equal to NbKnots-Mn-1.
247   If Ki <= U < K(i+1) then KnotIndex should be equal to sum (M1+M2+...+Mi).
248   */
249
250   Standard_Integer First,Last;
251   if (FromK1 < ToK2) {
252     First = FromK1;
253     Last  = ToK2;
254   }
255   else {
256     First = ToK2;
257     Last  = FromK1;
258   }
259   Standard_Integer Last1 = Last - 1;
260   NewU = U;
261   if (IsPeriodic && (NewU < UFirst || NewU > ULast))
262     NewU = ElCLib::InPeriod(NewU, UFirst, ULast);
263   
264   BSplCLib::Hunt (Knots, NewU, KnotIndex);
265   
266   Standard_Real val;
267   const Standard_Integer  KLower = Knots.Lower(),
268                           KUpper = Knots.Upper();
269
270   const Standard_Real Eps = Epsilon(Min(Abs(Knots(KUpper)), Abs(U)));
271
272   const Standard_Real *knots = &Knots(KLower);
273   knots -= KLower;
274   if ( KnotIndex < Knots.Upper()) {
275     val = NewU - knots[KnotIndex + 1];
276     if (val < 0) val = - val;
277     // <= to be coherent with Segment where Eps corresponds to a bit of error.
278     if (val <= Eps) KnotIndex++; 
279   }
280   if (KnotIndex < First) KnotIndex = First;
281   if (KnotIndex > Last1) KnotIndex = Last1;
282   
283   if (KnotIndex != Last1) {
284     Standard_Real K1 = knots[KnotIndex];
285     Standard_Real K2 = knots[KnotIndex + 1];
286     val = K2 - K1;
287     if (val < 0) val = - val;
288
289     while (val <= Eps) {
290       KnotIndex++;
291
292       if(KnotIndex >= Knots.Upper())
293         break;
294
295       K1 = K2;
296       K2 = knots[KnotIndex + 1];
297       val = K2 - K1;
298       if (val < 0) val = - val;
299     }
300   }
301 }
302
303 //=======================================================================
304 //function : LocateParameter
305 //purpose  : the index is recomputed only if out of range
306 //pmn  28-01-97 -> eventual computation of the period.
307 //=======================================================================
308
309 void BSplCLib::LocateParameter 
310 (const Standard_Integer         Degree,
311  const TColStd_Array1OfReal&    Knots,
312  const TColStd_Array1OfInteger* Mults,
313  const Standard_Real            U,
314  const Standard_Boolean         Periodic,
315  Standard_Integer&              KnotIndex,
316  Standard_Real&                 NewU) 
317 {
318   Standard_Integer first,last;
319   if (Mults) {
320     if (Periodic) {
321       first = Knots.Lower();
322       last  = Knots.Upper();
323     }
324     else {
325       first = FirstUKnotIndex(Degree,*Mults);
326       last  = LastUKnotIndex (Degree,*Mults);
327     }
328   }
329   else {
330     first = Knots.Lower() + Degree;
331     last  = Knots.Upper() - Degree;
332   }
333   if ( KnotIndex < first || KnotIndex > last)
334     BSplCLib::LocateParameter(Knots, U, Periodic, first, last,
335                               KnotIndex, NewU, Knots(first), Knots(last));
336   else
337     NewU = U;
338 }
339
340 //=======================================================================
341 //function : MaxKnotMult
342 //purpose  : 
343 //=======================================================================
344
345 Standard_Integer BSplCLib::MaxKnotMult
346 (const Array1OfInteger& Mults,
347  const Standard_Integer          FromK1,
348  const Standard_Integer          ToK2)
349 {
350   Standard_Integer MLower = Mults.Lower();
351   const Standard_Integer *pmu = &Mults(MLower);
352   pmu -= MLower;
353   Standard_Integer MaxMult = pmu[FromK1];
354
355   for (Standard_Integer i = FromK1; i <= ToK2; i++) {
356     if (MaxMult < pmu[i]) MaxMult = pmu[i];
357   }
358   return MaxMult;
359 }
360
361 //=======================================================================
362 //function : MinKnotMult
363 //purpose  : 
364 //=======================================================================
365
366 Standard_Integer BSplCLib::MinKnotMult
367 (const Array1OfInteger& Mults,
368  const Standard_Integer          FromK1,
369  const Standard_Integer          ToK2)
370 {
371   Standard_Integer MLower = Mults.Lower();
372   const Standard_Integer *pmu = &Mults(MLower);
373   pmu -= MLower;
374   Standard_Integer MinMult = pmu[FromK1];
375
376   for (Standard_Integer i = FromK1; i <= ToK2; i++) {
377     if (MinMult > pmu[i]) MinMult = pmu[i];
378   }
379   return MinMult;
380 }
381
382 //=======================================================================
383 //function : NbPoles
384 //purpose  : 
385 //=======================================================================
386
387 Standard_Integer BSplCLib::NbPoles(const Standard_Integer Degree,
388                                    const Standard_Boolean Periodic,
389                                    const TColStd_Array1OfInteger& Mults)
390 {
391   Standard_Integer i,sigma = 0;
392   Standard_Integer f = Mults.Lower();
393   Standard_Integer l = Mults.Upper();
394   const Standard_Integer * pmu = &Mults(f);
395   pmu -= f;
396   Standard_Integer Mf = pmu[f];
397   Standard_Integer Ml = pmu[l];
398   if (Mf <= 0) return 0;
399   if (Ml <= 0) return 0;
400   if (Periodic) {
401     if (Mf > Degree) return 0;
402     if (Ml > Degree) return 0;
403     if (Mf != Ml   ) return 0;
404     sigma = Mf;
405   }
406   else {
407     Standard_Integer Deg1 = Degree + 1;
408     if (Mf > Deg1) return 0;
409     if (Ml > Deg1) return 0;
410     sigma = Mf + Ml - Deg1;
411   }
412     
413   for (i = f + 1; i < l; i++) {
414     if (pmu[i] <= 0    ) return 0;
415     if (pmu[i] > Degree) return 0;
416     sigma += pmu[i];
417   }
418   return sigma;
419 }
420
421 //=======================================================================
422 //function : KnotSequenceLength
423 //purpose  : 
424 //=======================================================================
425
426 Standard_Integer BSplCLib::KnotSequenceLength
427 (const TColStd_Array1OfInteger& Mults,
428  const Standard_Integer         Degree,
429  const Standard_Boolean         Periodic)
430 {
431   Standard_Integer i,l = 0;
432   Standard_Integer MLower = Mults.Lower();
433   Standard_Integer MUpper = Mults.Upper();
434   const Standard_Integer * pmu = &Mults(MLower);
435   pmu -= MLower;
436
437   for (i = MLower; i <= MUpper; i++)
438     l += pmu[i];
439   if (Periodic) l += 2 * (Degree + 1 - pmu[MLower]);
440   return l;
441 }
442
443 //=======================================================================
444 //function : KnotSequence
445 //purpose  : 
446 //=======================================================================
447
448 void BSplCLib::KnotSequence 
449 (const TColStd_Array1OfReal&    Knots,
450  const TColStd_Array1OfInteger& Mults,
451  TColStd_Array1OfReal&          KnotSeq,
452  const Standard_Boolean         Periodic)
453 {
454   BSplCLib::KnotSequence(Knots,Mults,0,Periodic,KnotSeq);
455 }
456
457 //=======================================================================
458 //function : KnotSequence
459 //purpose  : 
460 //=======================================================================
461
462 void BSplCLib::KnotSequence 
463 (const TColStd_Array1OfReal&    Knots,
464  const TColStd_Array1OfInteger& Mults,
465  const Standard_Integer         Degree,
466  const Standard_Boolean         Periodic,
467  TColStd_Array1OfReal&          KnotSeq)
468 {
469   Standard_Real K;
470   Standard_Integer Mult;
471   Standard_Integer MLower = Mults.Lower();
472   const Standard_Integer * pmu = &Mults(MLower);
473   pmu -= MLower;
474   Standard_Integer KLower = Knots.Lower();
475   Standard_Integer KUpper = Knots.Upper();
476   const Standard_Real * pkn = &Knots(KLower);
477   pkn -= KLower;
478   Standard_Integer M1 = Degree + 1 - pmu[MLower];  // for periodic
479   Standard_Integer i,j,index = Periodic ? M1 + 1 : 1;
480
481   for (i = KLower; i <= KUpper; i++) {
482     Mult = pmu[i];
483     K    = pkn[i];
484
485     for (j = 1; j <= Mult; j++) { 
486       KnotSeq (index) = K;   
487       index++;
488     }
489   }
490   if (Periodic) {
491     Standard_Real period = pkn[KUpper] - pkn[KLower];
492     Standard_Integer m;
493     m = 1;
494     j = KUpper - 1;
495
496     for (i = M1; i >= 1; i--) {
497       KnotSeq(i) = pkn[j] - period;
498       m++;
499       if (m > pmu[j]) {
500         j--;
501         m = 1;
502       }
503     }
504     m = 1;
505     j = KLower + 1;
506
507     for (i = index; i <= KnotSeq.Upper(); i++) {
508       KnotSeq(i) = pkn[j] + period;
509       m++;
510       if (m > pmu[j]) {
511         j++;
512         m = 1;
513       }
514     }
515   }
516 }
517
518 //=======================================================================
519 //function : KnotsLength
520 //purpose  : 
521 //=======================================================================
522  Standard_Integer BSplCLib::KnotsLength(const TColStd_Array1OfReal& SeqKnots,
523 //                                      const Standard_Boolean Periodic)
524                                         const Standard_Boolean )
525 {
526   Standard_Integer sizeMult = 1; 
527   Standard_Real val = SeqKnots(1);
528   for (Standard_Integer jj=2;
529        jj<=SeqKnots.Length();jj++)
530     {
531       // test on strict equality on nodes
532       if (SeqKnots(jj)!=val)
533         {
534           val = SeqKnots(jj);
535           sizeMult++;
536         }
537     }
538   return sizeMult;
539 }
540
541 //=======================================================================
542 //function : Knots
543 //purpose  : 
544 //=======================================================================
545 void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots, 
546                      TColStd_Array1OfReal &knots,
547                      TColStd_Array1OfInteger &mult,
548 //                   const Standard_Boolean Periodic)
549                      const Standard_Boolean )
550 {
551   Standard_Real val = SeqKnots(1);
552   Standard_Integer kk=1;
553   knots(kk) = val;
554   mult(kk)  = 1;
555
556   for (Standard_Integer jj=2;jj<=SeqKnots.Length();jj++)
557     {
558       // test on strict equality on nodes
559       if (SeqKnots(jj)!=val)
560         {
561           val = SeqKnots(jj);
562           kk++;
563           knots(kk) = val;
564           mult(kk)  = 1;
565         }
566       else
567         {
568           mult(kk)++;
569         }
570     }
571 }
572
573 //=======================================================================
574 //function : KnotForm
575 //purpose  : 
576 //=======================================================================
577
578 BSplCLib_KnotDistribution BSplCLib::KnotForm
579 (const Array1OfReal& Knots,
580  const Standard_Integer       FromK1,
581  const Standard_Integer       ToK2)
582 {
583    Standard_Real DU0,DU1,Ui,Uj,Eps0,val;
584    BSplCLib_KnotDistribution  KForm = BSplCLib_Uniform;
585
586    if (FromK1 + 1 > Knots.Upper())
587    {
588      return BSplCLib_Uniform;
589    }
590
591    Ui  = Knots(FromK1);
592    if (Ui < 0) Ui = - Ui;
593    Uj  = Knots(FromK1 + 1);
594    if (Uj < 0) Uj = - Uj;
595    DU0 = Uj - Ui;
596    if (DU0 < 0) DU0 = - DU0;
597    Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
598    Standard_Integer i = FromK1 + 1;
599
600    while (KForm != BSplCLib_NonUniform && i < ToK2) {
601      Ui = Knots(i);
602      if (Ui < 0) Ui = - Ui;
603      i++;
604      Uj = Knots(i);
605      if (Uj < 0) Uj = - Uj;
606      DU1 = Uj - Ui;
607      if (DU1 < 0) DU1 = - DU1;
608      val = DU1 - DU0;
609      if (val < 0) val = -val;
610      if (val > Eps0) KForm = BSplCLib_NonUniform;
611      DU0 = DU1;
612      Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
613    }
614    return KForm;
615 }
616
617 //=======================================================================
618 //function : MultForm
619 //purpose  : 
620 //=======================================================================
621
622 BSplCLib_MultDistribution BSplCLib::MultForm
623 (const Array1OfInteger& Mults,
624  const Standard_Integer          FromK1,
625  const Standard_Integer          ToK2)
626 {
627   Standard_Integer First,Last;
628   if (FromK1 < ToK2) {
629     First = FromK1;
630     Last  = ToK2;
631   }
632   else {
633     First = ToK2;
634     Last  = FromK1;
635   }
636   if (First + 1 > Mults.Upper())
637   {
638     return BSplCLib_Constant;
639   }
640
641   Standard_Integer FirstMult = Mults(First);
642   BSplCLib_MultDistribution MForm = BSplCLib_Constant;
643   Standard_Integer i    = First + 1;
644   Standard_Integer Mult = Mults(i);
645   
646 //  while (MForm != BSplCLib_NonUniform && i <= Last) { ???????????JR????????
647   while (MForm != BSplCLib_NonConstant && i <= Last) {
648     if (i == First + 1) {  
649       if (Mult != FirstMult)      MForm = BSplCLib_QuasiConstant;
650     }
651     else if (i == Last)  {
652       if (MForm == BSplCLib_QuasiConstant) {
653         if (FirstMult != Mults(i))  MForm = BSplCLib_NonConstant;
654       }
655       else {
656         if (Mult != Mults(i))       MForm = BSplCLib_NonConstant;
657       }
658     }
659     else {
660       if (Mult != Mults(i))         MForm = BSplCLib_NonConstant;
661       Mult = Mults(i);
662     }
663     i++;
664   }
665   return MForm;
666 }
667
668 //=======================================================================
669 //function : KnotAnalysis
670 //purpose  : 
671 //=======================================================================
672
673 void BSplCLib::KnotAnalysis (const Standard_Integer         Degree,
674                              const Standard_Boolean         Periodic,
675                              const TColStd_Array1OfReal&    CKnots,
676                              const TColStd_Array1OfInteger& CMults,
677                              GeomAbs_BSplKnotDistribution&  KnotForm,
678                              Standard_Integer&              MaxKnotMult)
679 {
680   KnotForm = GeomAbs_NonUniform;
681
682   BSplCLib_KnotDistribution KSet = 
683     BSplCLib::KnotForm (CKnots, 1, CKnots.Length());
684   
685
686   if (KSet == BSplCLib_Uniform) {
687     BSplCLib_MultDistribution MSet =
688       BSplCLib::MultForm (CMults, 1, CMults.Length());
689     switch (MSet) {
690     case BSplCLib_NonConstant   :       
691       break;
692     case BSplCLib_Constant      : 
693       if (CKnots.Length() == 2) {
694         KnotForm = GeomAbs_PiecewiseBezier;
695       }
696       else {
697         if (CMults (1) == 1)  KnotForm = GeomAbs_Uniform;   
698       }
699       break;
700     case BSplCLib_QuasiConstant :   
701       if (CMults (1) == Degree + 1) {
702         Standard_Real M = CMults (2);
703         if (M == Degree )   KnotForm = GeomAbs_PiecewiseBezier;
704         else if  (M == 1)   KnotForm = GeomAbs_QuasiUniform;
705       }
706       break;
707     }
708   }
709
710   Standard_Integer FirstKM = 
711     Periodic ? CKnots.Lower() :  BSplCLib::FirstUKnotIndex (Degree,CMults);
712   Standard_Integer LastKM =
713     Periodic ? CKnots.Upper() :  BSplCLib::LastUKnotIndex (Degree,CMults);
714   MaxKnotMult = 0;
715   if (LastKM - FirstKM != 1) {
716     Standard_Integer Multi;
717     for (Standard_Integer i = FirstKM + 1; i < LastKM; i++) {
718       Multi = CMults (i);
719       MaxKnotMult = Max (MaxKnotMult, Multi);
720     }
721   }
722 }
723
724 //=======================================================================
725 //function : Reparametrize
726 //purpose  : 
727 //=======================================================================
728
729 void BSplCLib::Reparametrize
730 (const Standard_Real      U1,
731  const Standard_Real      U2,
732  Array1OfReal&   Knots)
733 {
734   Standard_Integer Lower  = Knots.Lower();
735   Standard_Integer Upper  = Knots.Upper();
736   Standard_Real UFirst    = Min (U1, U2);
737   Standard_Real ULast     = Max (U1, U2);
738   Standard_Real NewLength = ULast - UFirst;
739   BSplCLib_KnotDistribution KSet = BSplCLib::KnotForm (Knots, Lower, Upper);
740   if (KSet == BSplCLib_Uniform) {
741     Standard_Real DU = NewLength / (Upper - Lower);
742     Knots (Lower) = UFirst;
743
744     for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
745       Knots (i) = Knots (i-1) + DU;
746     }
747   }
748   else {
749     Standard_Real K2;
750     Standard_Real Ratio;
751     Standard_Real K1 = Knots (Lower);
752     Standard_Real Length = Knots (Upper) - Knots (Lower);
753     Knots (Lower) = UFirst;
754
755     for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
756       K2 = Knots (i);
757       Ratio = (K2 - K1) / Length;
758       Knots (i) = Knots (i-1) + (NewLength * Ratio);
759
760       //for CheckCurveData
761       Standard_Real Eps = Epsilon( Abs(Knots(i-1)) );
762       if (Knots(i) - Knots(i-1) <= Eps)
763         Knots(i) = NextAfter (Knots(i-1) + Eps, RealLast());
764
765       K1 = K2;
766     }
767   }
768 }
769
770 //=======================================================================
771 //function : Reverse
772 //purpose  : 
773 //=======================================================================
774
775 void  BSplCLib::Reverse(TColStd_Array1OfReal& Knots)
776 {
777   Standard_Integer first = Knots.Lower();
778   Standard_Integer last  = Knots.Upper();
779   Standard_Real kfirst = Knots(first);
780   Standard_Real klast = Knots(last);
781   Standard_Real tfirst = kfirst;
782   Standard_Real tlast  = klast;
783   first++;
784   last--;
785
786   while (first <= last) {
787     tfirst += klast - Knots(last);
788     tlast  -= Knots(first) - kfirst;
789     kfirst = Knots(first);
790     klast  = Knots(last);
791     Knots(first) = tfirst;
792     Knots(last)  = tlast;
793     first++;
794     last--;
795   }
796 }
797
798 //=======================================================================
799 //function : Reverse
800 //purpose  : 
801 //=======================================================================
802
803 void  BSplCLib::Reverse(TColStd_Array1OfInteger& Mults)
804 {
805   Standard_Integer first = Mults.Lower();
806   Standard_Integer last  = Mults.Upper();
807   Standard_Integer temp;
808
809   while (first < last) {
810     temp = Mults(first);
811     Mults(first) = Mults(last);
812     Mults(last) = temp;
813     first++;
814     last--;
815   }
816 }
817
818 //=======================================================================
819 //function : Reverse
820 //purpose  : 
821 //=======================================================================
822
823 void  BSplCLib::Reverse(TColStd_Array1OfReal& Weights,
824                         const Standard_Integer L)
825 {
826   Standard_Integer i, l = L;
827   l = Weights.Lower()+(l-Weights.Lower())%(Weights.Upper()-Weights.Lower()+1);
828
829   TColStd_Array1OfReal temp(0,Weights.Length()-1);
830
831   for (i = Weights.Lower(); i <= l; i++)
832     temp(l-i) = Weights(i);
833
834   for (i = l+1; i <= Weights.Upper(); i++)
835     temp(l-Weights.Lower()+Weights.Upper()-i+1) = Weights(i);
836
837   for (i = Weights.Lower(); i <= Weights.Upper(); i++)
838     Weights(i) = temp(i-Weights.Lower());
839 }
840
841 //=======================================================================
842 //function : IsRational
843 //purpose  : 
844 //=======================================================================
845
846 Standard_Boolean  BSplCLib::IsRational(const TColStd_Array1OfReal& Weights,
847                                        const Standard_Integer I1,
848                                        const Standard_Integer I2,
849 //                                     const Standard_Real Epsi)
850                                        const Standard_Real )
851 {
852   Standard_Integer i, f = Weights.Lower(), l = Weights.Length();
853   Standard_Integer I3 = I2 - f;
854   const Standard_Real * WG = &Weights(f);
855   WG -= f;
856
857   for (i = I1 - f; i < I3; i++) {
858     if (WG[f + (i % l)] != WG[f + ((i + 1) % l)]) return Standard_True;
859   }
860   return Standard_False ;
861 }
862
863 //=======================================================================
864 //function : Eval
865 //purpose  : evaluate point and derivatives
866 //=======================================================================
867
868 void  BSplCLib::Eval(const Standard_Real U,
869                      const Standard_Integer Degree,
870                      Standard_Real& Knots, 
871                      const Standard_Integer Dimension, 
872                      Standard_Real& Poles)
873 {
874   Standard_Integer step,i,Dms,Dm1,Dpi,Sti;
875   Standard_Real X, Y, *poles, *knots = &Knots;
876   Dm1 = Dms = Degree;
877   Dm1--;
878   Dms++;
879   switch (Dimension) { 
880
881   case 1 : {
882     
883     for (step = - 1; step < Dm1; step++) {
884       Dms--;
885       poles = &Poles;
886       Dpi   = Dm1;
887       Sti   = step;
888       
889       for (i = 0; i < Dms; i++) {
890         Dpi++;
891         Sti++;
892         X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
893         Y = 1 - X;
894         poles[0] *= X; poles[0] += Y * poles[1];
895         poles += 1;
896       }
897     }
898     break;
899   }
900   case 2 : {
901     
902     for (step = - 1; step < Dm1; step++) {
903       Dms--;
904       poles = &Poles;
905       Dpi   = Dm1;
906       Sti   = step;
907       
908       for (i = 0; i < Dms; i++) {
909         Dpi++;
910         Sti++;
911         X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
912         Y = 1 - X;
913         poles[0] *= X; poles[0] += Y * poles[2];
914         poles[1] *= X; poles[1] += Y * poles[3];
915         poles += 2;
916       }
917     }
918     break;
919   }
920   case 3 : {
921     
922     for (step = - 1; step < Dm1; step++) {
923       Dms--;
924       poles = &Poles;
925       Dpi   = Dm1;
926       Sti   = step;
927       
928       for (i = 0; i < Dms; i++) {
929         Dpi++;
930         Sti++;
931         X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
932         Y = 1 - X;
933         poles[0] *= X; poles[0] += Y * poles[3];
934         poles[1] *= X; poles[1] += Y * poles[4];
935         poles[2] *= X; poles[2] += Y * poles[5];
936         poles += 3;
937       }
938     }
939     break;
940   }
941   case 4 : {
942     
943     for (step = - 1; step < Dm1; step++) {
944       Dms--;
945       poles = &Poles;
946       Dpi   = Dm1;
947       Sti   = step;
948       
949       for (i = 0; i < Dms; i++) {
950         Dpi++;
951         Sti++;
952         X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
953         Y = 1 - X;
954         poles[0] *= X; poles[0] += Y * poles[4];
955         poles[1] *= X; poles[1] += Y * poles[5];
956         poles[2] *= X; poles[2] += Y * poles[6];
957         poles[3] *= X; poles[3] += Y * poles[7];
958         poles += 4;
959       }
960     }
961     break;
962   }
963     default : {
964       Standard_Integer k;
965       
966       for (step = - 1; step < Dm1; step++) {
967         Dms--;
968         poles = &Poles;
969         Dpi   = Dm1;
970         Sti   = step;
971         
972         for (i = 0; i < Dms; i++) {
973           Dpi++;
974           Sti++;
975           X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
976           Y = 1 - X;
977           
978           for (k = 0; k < Dimension; k++) {
979             poles[k] *= X;
980             poles[k] += Y * poles[k + Dimension];
981           }
982           poles += Dimension;
983         }
984       }
985     }
986   }
987 }
988
989 //=======================================================================
990 //function : BoorScheme
991 //purpose  : 
992 //=======================================================================
993
994 void  BSplCLib::BoorScheme(const Standard_Real U,
995                            const Standard_Integer Degree,
996                            Standard_Real& Knots, 
997                            const Standard_Integer Dimension, 
998                            Standard_Real& Poles, 
999                            const Standard_Integer Depth, 
1000                            const Standard_Integer Length)
1001 {
1002   //
1003   // Compute the values
1004   //
1005   //  P(i,j) (U).
1006   //
1007   // for i = 0 to Depth, 
1008   // j = 0 to Length - i
1009   //
1010   //  The Boor scheme is :
1011   //
1012   //  P(0,j) = Pole(j)
1013   //  P(i,j) = x * P(i-1,j) + (1-x) * P(i-1,j+1)
1014   //
1015   //    where x = (knot(i+j+Degree) - U) / (knot(i+j+Degree) - knot(i+j))
1016   //
1017   //
1018   //  The values are stored in the array Poles
1019   //  They are alternatively written if the odd and even positions.
1020   //
1021   //  The successives contents of the array are
1022   //   ***** means unitialised, l = Degree + Length
1023   //
1024   //  P(0,0) ****** P(0,1) ...... P(0,l-1) ******** P(0,l)
1025   //  P(0,0) P(1,0) P(0,1) ...... P(0,l-1) P(1,l-1) P(0,l)
1026   //  P(0,0) P(1,0) P(2,0) ...... P(2,l-1) P(1,l-1) P(0,l)
1027   //
1028
1029   Standard_Integer i,k,step;
1030   Standard_Real *knots = &Knots;
1031   Standard_Real *pole, *firstpole = &Poles - 2 * Dimension;
1032   // the steps of recursion
1033
1034   for (step = 0; step < Depth; step++) {
1035     firstpole += Dimension;
1036     pole = firstpole;
1037     // compute the new row of poles
1038
1039     for (i = step; i < Length; i++) {
1040       pole += 2 * Dimension;
1041       // coefficient
1042       Standard_Real X = (knots[i+Degree-step] - U) 
1043         / (knots[i+Degree-step] - knots[i]);
1044       Standard_Real Y = 1. - X;
1045       // Value
1046       // P(i,j) = X * P(i-1,j) + (1-X) * P(i-1,j+1)
1047
1048       for (k = 0; k < Dimension; k++)
1049         pole[k] = X * pole[k - Dimension] + Y * pole[k + Dimension];
1050     }
1051   }
1052 }
1053
1054 //=======================================================================
1055 //function : AntiBoorScheme
1056 //purpose  : 
1057 //=======================================================================
1058
1059 Standard_Boolean  BSplCLib::AntiBoorScheme(const Standard_Real    U,
1060                                            const Standard_Integer Degree,
1061                                            Standard_Real&         Knots, 
1062                                            const Standard_Integer Dimension, 
1063                                            Standard_Real&         Poles, 
1064                                            const Standard_Integer Depth, 
1065                                            const Standard_Integer Length,
1066                                            const Standard_Real    Tolerance)
1067 {
1068   // do the Boor scheme reverted.
1069
1070   Standard_Integer i,k,step, half_length;
1071   Standard_Real *knots = &Knots;
1072   Standard_Real z,X,Y,*pole, *firstpole = &Poles + (Depth-1) * Dimension;
1073
1074   // Test the special case length = 1 
1075   // only verification of the central point
1076
1077   if (Length == 1) {
1078     X = (knots[Degree] - U) / (knots[Degree] - knots[0]);
1079     Y = 1. - X;
1080
1081     for (k = 0; k < Dimension; k++) {
1082       z = X * firstpole[k] + Y * firstpole[k+2*Dimension];
1083       if (Abs(z - firstpole[k+Dimension]) > Tolerance) 
1084         return Standard_False;
1085     }
1086     return Standard_True;
1087   }
1088
1089   // General case
1090   // the steps of recursion
1091
1092   for (step = Depth-1; step >= 0; step--) {
1093     firstpole -= Dimension;
1094     pole = firstpole;
1095
1096     // first step from left to right
1097
1098     for (i = step; i < Length-1; i++) {
1099       pole += 2 * Dimension;
1100
1101       X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1102       Y = 1. - X;
1103
1104       for (k = 0; k < Dimension; k++)
1105         pole[k+Dimension] = (pole[k] - X*pole[k-Dimension]) / Y;
1106
1107     }
1108
1109     // second step from right to left
1110     pole += 4* Dimension;
1111     half_length = (Length - 1 + step) / 2  ;
1112     //
1113     // only do half of the way from right to left 
1114     // otherwise it start degenerating because of 
1115     // overflows
1116     // 
1117
1118     for (i = Length-1; i > half_length ; i--) {
1119       pole -= 2 * Dimension;
1120
1121       // coefficient
1122       X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1123       Y = 1. - X;
1124
1125       for (k = 0; k < Dimension; k++) {
1126         z = (pole[k] - Y * pole[k+Dimension]) / X;
1127         if (Abs(z-pole[k-Dimension]) > Tolerance) 
1128           return Standard_False;
1129         pole[k-Dimension] += z;
1130         pole[k-Dimension] /= 2.;
1131       }
1132     }
1133   }
1134   return Standard_True;
1135 }
1136
1137 //=======================================================================
1138 //function : Derivative
1139 //purpose  : 
1140 //=======================================================================
1141
1142 void  BSplCLib::Derivative(const Standard_Integer Degree, 
1143                            Standard_Real& Knots, 
1144                            const Standard_Integer Dimension, 
1145                            const Standard_Integer Length, 
1146                            const Standard_Integer Order, 
1147                            Standard_Real& Poles)
1148 {
1149   Standard_Integer i,k,step,span = Degree;
1150   Standard_Real *knot = &Knots;
1151
1152   for (step = 1; step <= Order; step++) {
1153     Standard_Real* pole = &Poles;
1154
1155     for (i = step; i < Length; i++) {
1156       Standard_Real coef = - span / (knot[i+span] - knot[i]);
1157
1158       for (k = 0; k < Dimension; k++) {
1159         pole[k] -= pole[k+Dimension];
1160         pole[k] *= coef;
1161       }
1162       pole += Dimension;
1163     }
1164     span--;
1165   }
1166 }
1167
1168 //=======================================================================
1169 //function : Bohm
1170 //purpose  : 
1171 //=======================================================================
1172
1173 void  BSplCLib::Bohm(const Standard_Real U,
1174                      const Standard_Integer Degree,
1175                      const Standard_Integer N,
1176                      Standard_Real& Knots,
1177                      const Standard_Integer Dimension,
1178                      Standard_Real& Poles)
1179 {
1180   // First phase independant of U, compute the poles of the derivatives
1181   Standard_Integer i,j,iDim,min,Dmi,DDmi,jDmi,Degm1;
1182   Standard_Real *knot = &Knots, *pole, coef, *tbis, *psav, *psDD, *psDDmDim;
1183   psav     = &Poles;
1184   if (N < Degree) min = N;
1185   else            min = Degree;
1186   Degm1 = Degree - 1;
1187   DDmi = (Degree << 1) + 1;
1188   switch (Dimension) { 
1189   case 1 : {
1190     psDD     = psav + Degree;
1191     psDDmDim = psDD - 1;
1192     
1193     for (i = 0; i < Degree; i++) {
1194       DDmi--;
1195       pole = psDD;
1196       tbis = psDDmDim;
1197       jDmi = DDmi;
1198       
1199       for (j = Degm1; j >= i; j--) {
1200         jDmi--;
1201         *pole -= *tbis;
1202   *pole = (knot[jDmi] == knot[j]) ? 0.0 :  *pole / (knot[jDmi] - knot[j]);
1203         pole--;
1204         tbis--;
1205       }
1206     }
1207     // Second phase, dependant of U
1208     iDim = - 1;
1209     
1210     for (i = 0; i < Degree; i++) {
1211       iDim += 1;
1212       pole  = psav + iDim;
1213       tbis  = pole + 1;
1214       coef  = U - knot[i];
1215       
1216       for (j = i; j >= 0; j--) {
1217         *pole += coef * (*tbis);
1218         pole--;
1219         tbis--;
1220       }
1221     }
1222     // multiply by the degrees
1223     coef = Degree;
1224     Dmi  = Degree;
1225     pole = psav + 1;
1226     
1227     for (i = 1; i <= min; i++) {
1228       *pole *= coef; pole++;
1229       Dmi--;
1230       coef  *= Dmi;
1231     }
1232     break;
1233   }
1234   case 2 : {
1235     psDD     = psav + (Degree << 1);
1236     psDDmDim = psDD - 2;
1237     
1238     for (i = 0; i < Degree; i++) {
1239       DDmi--;
1240       pole = psDD;
1241       tbis = psDDmDim;
1242       jDmi = DDmi;
1243       
1244       for (j = Degm1; j >= i; j--) {
1245         jDmi--;
1246         coef   = (knot[jDmi] == knot[j]) ? 0.0 : 1. / (knot[jDmi] - knot[j]);
1247         *pole -= *tbis; *pole *= coef; pole++; tbis++;
1248         *pole -= *tbis; *pole *= coef;
1249         pole  -= 3;
1250         tbis  -= 3;
1251       }
1252     }
1253     // Second phase, dependant of U
1254     iDim = - 2;
1255     
1256     for (i = 0; i < Degree; i++) {
1257       iDim += 2;
1258       pole  = psav + iDim;
1259       tbis  = pole + 2;
1260       coef  = U - knot[i];
1261       
1262       for (j = i; j >= 0; j--) {
1263         *pole += coef * (*tbis); pole++; tbis++;
1264         *pole += coef * (*tbis);
1265         pole  -= 3;
1266         tbis  -= 3;
1267       }
1268     }
1269     // multiply by the degrees
1270     coef = Degree;
1271     Dmi  = Degree;
1272     pole = psav + 2;
1273     
1274     for (i = 1; i <= min; i++) {
1275       *pole *= coef; pole++;
1276       *pole *= coef; pole++;
1277       Dmi--;
1278       coef  *= Dmi;
1279     }
1280     break;
1281   }
1282   case 3 : {
1283     psDD     = psav + (Degree << 1) + Degree;
1284     psDDmDim = psDD - 3;
1285     
1286     for (i = 0; i < Degree; i++) {
1287       DDmi--;
1288       pole = psDD;
1289       tbis = psDDmDim;
1290       jDmi = DDmi;
1291       
1292       for (j = Degm1; j >= i; j--) {
1293         jDmi--;
1294         coef   = (knot[jDmi] == knot[j]) ? 0.0 : 1. / (knot[jDmi] - knot[j]);
1295         *pole -= *tbis; *pole *= coef; pole++; tbis++;
1296         *pole -= *tbis; *pole *= coef; pole++; tbis++;
1297         *pole -= *tbis; *pole *= coef;
1298         pole  -= 5;
1299         tbis  -= 5;
1300       }
1301     }
1302     // Second phase, dependant of U
1303     iDim = - 3;
1304     
1305     for (i = 0; i < Degree; i++) {
1306       iDim += 3;
1307       pole  = psav + iDim;
1308       tbis  = pole + 3;
1309       coef  = U - knot[i];
1310       
1311       for (j = i; j >= 0; j--) {
1312         *pole += coef * (*tbis); pole++; tbis++;
1313         *pole += coef * (*tbis); pole++; tbis++;
1314         *pole += coef * (*tbis);
1315         pole  -= 5;
1316         tbis  -= 5;
1317       }
1318     }
1319     // multiply by the degrees
1320     coef = Degree;
1321     Dmi  = Degree;
1322     pole = psav + 3;
1323     
1324     for (i = 1; i <= min; i++) {
1325       *pole *= coef; pole++;
1326       *pole *= coef; pole++;
1327       *pole *= coef; pole++;
1328       Dmi--;
1329       coef  *= Dmi;
1330     }
1331     break;
1332   }
1333   case 4 : {
1334     psDD     = psav + (Degree << 2);
1335     psDDmDim = psDD - 4;
1336     
1337     for (i = 0; i < Degree; i++) {
1338       DDmi--;
1339       pole = psDD;
1340       tbis = psDDmDim;
1341       jDmi = DDmi;
1342       
1343       for (j = Degm1; j >= i; j--) {
1344         jDmi--;
1345         coef   = (knot[jDmi]  == knot[j]) ? 0.0 : 1. /(knot[jDmi] - knot[j]) ;
1346         *pole -= *tbis; *pole *= coef; pole++; tbis++;
1347         *pole -= *tbis; *pole *= coef; pole++; tbis++;
1348         *pole -= *tbis; *pole *= coef; pole++; tbis++;
1349         *pole -= *tbis; *pole *= coef;
1350         pole  -= 7;
1351         tbis  -= 7;
1352       }
1353     }
1354     // Second phase, dependant of U
1355     iDim = - 4;
1356     
1357     for (i = 0; i < Degree; i++) {
1358       iDim += 4;
1359       pole  = psav + iDim;
1360       tbis  = pole + 4;
1361       coef  = U - knot[i];
1362       
1363       for (j = i; j >= 0; j--) {
1364         *pole += coef * (*tbis); pole++; tbis++;
1365         *pole += coef * (*tbis); pole++; tbis++;
1366         *pole += coef * (*tbis); pole++; tbis++;
1367         *pole += coef * (*tbis);
1368         pole  -= 7;
1369         tbis  -= 7;
1370       }
1371     }
1372     // multiply by the degrees
1373     coef = Degree; 
1374     Dmi  = Degree;
1375     pole = psav + 4;
1376    
1377     for (i = 1; i <= min; i++) {
1378       *pole *= coef; pole++;
1379       *pole *= coef; pole++;
1380       *pole *= coef; pole++;
1381       *pole *= coef; pole++;
1382       Dmi--;
1383       coef  *= Dmi;
1384     }
1385     break;
1386   }
1387     default : {
1388       Standard_Integer k;
1389       Standard_Integer Dim2 = Dimension << 1;
1390       psDD     = psav + Degree * Dimension;
1391       psDDmDim = psDD - Dimension;
1392       
1393       for (i = 0; i < Degree; i++) {
1394         DDmi--;
1395         pole = psDD;
1396         tbis = psDDmDim;
1397         jDmi = DDmi;
1398         
1399         for (j = Degm1; j >= i; j--) {
1400           jDmi--;
1401           coef = (knot[jDmi] == knot[j]) ? 0.0 : 1. / (knot[jDmi] - knot[j]);
1402           
1403           for (k = 0; k < Dimension; k++) {
1404             *pole -= *tbis; *pole *= coef; pole++; tbis++;
1405           }
1406           pole -= Dim2;
1407           tbis -= Dim2;
1408         }
1409       }
1410       // Second phase, dependant of U
1411       iDim = - Dimension;
1412       
1413       for (i = 0; i < Degree; i++) {
1414         iDim += Dimension;
1415         pole  = psav + iDim;
1416         tbis  = pole + Dimension;
1417         coef  = U - knot[i];
1418         
1419         for (j = i; j >= 0; j--) {
1420           
1421           for (k = 0; k < Dimension; k++) {
1422             *pole += coef * (*tbis); pole++; tbis++;
1423           }
1424           pole -= Dim2;
1425           tbis -= Dim2;
1426         }
1427       }
1428       // multiply by the degrees
1429       coef = Degree;
1430       Dmi  = Degree;
1431       pole = psav + Dimension;
1432       
1433       for (i = 1; i <= min; i++) {
1434         
1435         for (k = 0; k < Dimension; k++) {
1436           *pole *= coef; pole++;
1437         }
1438         Dmi--;
1439         coef *= Dmi;
1440       }
1441     }
1442   }
1443 }
1444
1445 //=======================================================================
1446 //function : BuildKnots
1447 //purpose  : 
1448 //=======================================================================
1449
1450 void BSplCLib::BuildKnots(const Standard_Integer         Degree,
1451                           const Standard_Integer         Index,
1452                           const Standard_Boolean         Periodic,
1453                           const TColStd_Array1OfReal&    Knots,
1454                           const TColStd_Array1OfInteger* Mults,
1455                           Standard_Real&                 LK)
1456 {
1457   Standard_Integer KLower = Knots.Lower();
1458   const Standard_Real * pkn = &Knots(KLower);
1459   pkn -= KLower;
1460   Standard_Real *knot = &LK;
1461   if (Mults == NULL) {
1462     switch (Degree) {
1463     case 1 : {
1464       Standard_Integer j = Index    ;
1465       knot[0] = pkn[j]; j++;
1466       knot[1] = pkn[j];
1467       break;
1468     }
1469     case 2 : {
1470       Standard_Integer j = Index - 1;
1471       knot[0] = pkn[j]; j++;
1472       knot[1] = pkn[j]; j++;
1473       knot[2] = pkn[j]; j++;
1474       knot[3] = pkn[j];
1475       break;
1476     }
1477     case 3 : {
1478       Standard_Integer j = Index - 2;
1479       knot[0] = pkn[j]; j++;
1480       knot[1] = pkn[j]; j++;
1481       knot[2] = pkn[j]; j++;
1482       knot[3] = pkn[j]; j++;
1483       knot[4] = pkn[j]; j++;
1484       knot[5] = pkn[j];
1485       break;
1486     }
1487     case 4 : {
1488       Standard_Integer j = Index - 3;
1489       knot[0] = pkn[j]; j++;
1490       knot[1] = pkn[j]; j++;
1491       knot[2] = pkn[j]; j++;
1492       knot[3] = pkn[j]; j++;
1493       knot[4] = pkn[j]; j++;
1494       knot[5] = pkn[j]; j++;
1495       knot[6] = pkn[j]; j++;
1496       knot[7] = pkn[j];
1497       break;
1498     }
1499     case 5 : {
1500       Standard_Integer j = Index - 4;
1501       knot[0] = pkn[j]; j++;
1502       knot[1] = pkn[j]; j++;
1503       knot[2] = pkn[j]; j++;
1504       knot[3] = pkn[j]; j++;
1505       knot[4] = pkn[j]; j++;
1506       knot[5] = pkn[j]; j++;
1507       knot[6] = pkn[j]; j++;
1508       knot[7] = pkn[j]; j++;
1509       knot[8] = pkn[j]; j++;
1510       knot[9] = pkn[j];
1511       break;
1512     }
1513     case 6 : {
1514       Standard_Integer j = Index - 5;
1515       knot[ 0] = pkn[j]; j++;
1516       knot[ 1] = pkn[j]; j++;
1517       knot[ 2] = pkn[j]; j++;
1518       knot[ 3] = pkn[j]; j++;
1519       knot[ 4] = pkn[j]; j++;
1520       knot[ 5] = pkn[j]; j++;
1521       knot[ 6] = pkn[j]; j++;
1522       knot[ 7] = pkn[j]; j++;
1523       knot[ 8] = pkn[j]; j++;
1524       knot[ 9] = pkn[j]; j++;
1525       knot[10] = pkn[j]; j++;
1526       knot[11] = pkn[j];
1527       break;
1528     }
1529       default : {
1530         Standard_Integer i,j;
1531         Standard_Integer Deg2 = Degree << 1;
1532         j = Index - Degree;
1533         
1534         for (i = 0; i < Deg2; i++) {
1535           j++;
1536           knot[i] = pkn[j];
1537         }
1538       }
1539     }
1540   }
1541   else {
1542     Standard_Integer i;
1543     Standard_Integer Deg1 = Degree - 1;
1544     Standard_Integer KUpper = Knots.Upper();
1545     Standard_Integer MLower = Mults->Lower();
1546     Standard_Integer MUpper = Mults->Upper();
1547     const Standard_Integer * pmu = &(*Mults)(MLower);
1548     pmu -= MLower;
1549     Standard_Real dknot = 0;
1550     Standard_Integer ilow = Index    , mlow = 0;
1551     Standard_Integer iupp = Index + 1, mupp = 0;
1552     Standard_Real loffset = 0., uoffset = 0.;
1553     Standard_Boolean getlow = Standard_True, getupp = Standard_True;
1554     if (Periodic) {
1555       dknot = pkn[KUpper] - pkn[KLower];
1556       if (iupp > MUpper) {
1557         iupp = MLower + 1;
1558         uoffset = dknot;
1559       }
1560     }
1561     // Find the knots around Index
1562
1563     for (i = 0; i < Degree; i++) {
1564       if (getlow) {
1565         mlow++;
1566         if (mlow > pmu[ilow]) {
1567           mlow = 1;
1568           ilow--;
1569           getlow =  (ilow >= MLower);
1570           if (Periodic && !getlow) {
1571             ilow = MUpper - 1;
1572             loffset = dknot;
1573             getlow = Standard_True;
1574           }
1575         }
1576         if (getlow)
1577           knot[Deg1 - i] = pkn[ilow] - loffset;
1578       }
1579       if (getupp) {
1580         mupp++;
1581         if (mupp > pmu[iupp]) {
1582           mupp = 1;
1583           iupp++;
1584           getupp = (iupp <= MUpper);
1585           if (Periodic && !getupp) {
1586             iupp = MLower + 1;
1587             uoffset = dknot;
1588             getupp = Standard_True;
1589           }
1590         }
1591         if (getupp)
1592           knot[Degree + i] = pkn[iupp] + uoffset;
1593       }
1594     }
1595   } 
1596 }
1597
1598 //=======================================================================
1599 //function : PoleIndex
1600 //purpose  : 
1601 //=======================================================================
1602
1603 Standard_Integer BSplCLib::PoleIndex(const Standard_Integer         Degree,
1604                                      const Standard_Integer         Index,
1605                                      const Standard_Boolean         Periodic,
1606                                      const TColStd_Array1OfInteger& Mults)
1607 {
1608   Standard_Integer i, pindex = 0;
1609
1610   for (i = Mults.Lower(); i <= Index; i++)
1611     pindex += Mults(i);
1612   if (Periodic)
1613     pindex -= Mults(Mults.Lower());
1614   else
1615     pindex -= Degree + 1;
1616
1617   return pindex;
1618 }
1619
1620 //=======================================================================
1621 //function : BuildBoor
1622 //purpose  : builds the local array for boor
1623 //=======================================================================
1624
1625 void  BSplCLib::BuildBoor(const Standard_Integer         Index,
1626                           const Standard_Integer         Length,
1627                           const Standard_Integer         Dimension,
1628                           const TColStd_Array1OfReal&    Poles,
1629                           Standard_Real&                 LP)
1630 {
1631   Standard_Real *poles = &LP;
1632   Standard_Integer i,k, ip = Poles.Lower() + Index * Dimension;
1633   
1634   for (i = 0; i < Length+1; i++) {
1635
1636     for (k = 0; k < Dimension; k++) {
1637       poles[k] = Poles(ip);
1638       ip++;
1639       if (ip > Poles.Upper()) ip = Poles.Lower();
1640     }
1641     poles += 2 * Dimension;
1642   }
1643 }
1644
1645 //=======================================================================
1646 //function : BoorIndex
1647 //purpose  : 
1648 //=======================================================================
1649
1650 Standard_Integer  BSplCLib::BoorIndex(const Standard_Integer Index,
1651                                       const Standard_Integer Length,
1652                                       const Standard_Integer Depth)
1653 {
1654   if (Index <= Depth)  return Index;
1655   if (Index <= Length) return 2 * Index - Depth;
1656   return                      Length + Index - Depth;
1657 }
1658
1659 //=======================================================================
1660 //function : GetPole
1661 //purpose  : 
1662 //=======================================================================
1663
1664 void  BSplCLib::GetPole(const Standard_Integer         Index,
1665                         const Standard_Integer         Length,
1666                         const Standard_Integer         Depth,
1667                         const Standard_Integer         Dimension,
1668                         Standard_Real&                 LP,
1669                         Standard_Integer&              Position,
1670                         TColStd_Array1OfReal&          Pole)
1671 {
1672   Standard_Integer k;
1673   Standard_Real* pole = &LP + BoorIndex(Index,Length,Depth) * Dimension;
1674
1675   for (k = 0; k < Dimension; k++) {
1676     Pole(Position) = pole[k];
1677     Position++;
1678   }
1679   if (Position > Pole.Upper()) Position = Pole.Lower();
1680 }
1681
1682 //=======================================================================
1683 //function : PrepareInsertKnots
1684 //purpose  : 
1685 //=======================================================================
1686
1687 Standard_Boolean  BSplCLib::PrepareInsertKnots
1688 (const Standard_Integer         Degree,
1689  const Standard_Boolean         Periodic, 
1690  const TColStd_Array1OfReal&    Knots,
1691  const TColStd_Array1OfInteger& Mults,
1692  const TColStd_Array1OfReal&    AddKnots,
1693  const TColStd_Array1OfInteger* AddMults,
1694  Standard_Integer&              NbPoles,
1695  Standard_Integer&              NbKnots, 
1696  const Standard_Real            Tolerance,
1697  const Standard_Boolean         Add)
1698 {
1699   Standard_Boolean addflat = AddMults == NULL;
1700   
1701   Standard_Integer first,last;
1702   if (Periodic) {
1703     first = Knots.Lower();
1704     last  = Knots.Upper();
1705   }
1706   else {
1707     first = FirstUKnotIndex(Degree,Mults);
1708     last  = LastUKnotIndex(Degree,Mults);
1709   }
1710   Standard_Real adeltaK1 = Knots(first)-AddKnots(AddKnots.Lower());
1711    Standard_Real adeltaK2 = AddKnots(AddKnots.Upper())-Knots(last);
1712   if (adeltaK1 > Tolerance) return Standard_False;
1713   if (adeltaK2  > Tolerance) return Standard_False;
1714   
1715   Standard_Integer sigma = 0, mult, amult;
1716   NbKnots = 0;
1717   Standard_Integer  k  = Knots.Lower() - 1;
1718   Standard_Integer  ak = AddKnots.Lower();
1719
1720   if(Periodic && AddKnots.Length() > 1)
1721   {
1722     //gka for case when segments was produced on full period only one knot
1723     //was added in the end of curve
1724     if(fabs(adeltaK1) <= gp::Resolution() && 
1725        fabs(adeltaK2) <= gp::Resolution())
1726       ak++;
1727   }
1728   
1729   Standard_Integer aLastKnotMult = Mults (Knots.Upper());
1730   Standard_Real au,oldau = AddKnots(ak),Eps;
1731   
1732   while (ak <= AddKnots.Upper()) {
1733     au = AddKnots(ak);
1734     if (au < oldau) return Standard_False;
1735     oldau = au;
1736
1737     Eps = Max(Tolerance,Epsilon(au));
1738     
1739     while ((k < Knots.Upper()) && (Knots(k+1)  - au <= Eps)) {
1740       k++;
1741       NbKnots++;
1742       sigma += Mults(k);
1743     }
1744
1745     if (addflat) amult = 1;
1746     else         amult = Max(0,(*AddMults)(ak));
1747     
1748     while ((ak < AddKnots.Upper()) &&
1749            (Abs(au - AddKnots(ak+1)) <= Eps)) {
1750       ak++;
1751       if (Add) {
1752         if (addflat) amult++;
1753         else         amult += Max(0,(*AddMults)(ak));
1754       }
1755     }
1756     
1757     
1758     if (Abs(au - Knots(k)) <= Eps) {
1759       // identic to existing knot
1760       mult = Mults(k);
1761       if (Add) {
1762         if (mult + amult > Degree)
1763           amult = Max(0,Degree - mult);
1764         sigma += amult;
1765         
1766       }
1767       else if (amult > mult) {
1768         if (amult > Degree) amult = Degree;
1769         if (k == Knots.Upper () && Periodic)
1770         {
1771           aLastKnotMult = Max (amult, mult);
1772           sigma += 2 * (aLastKnotMult - mult);
1773         }
1774         else
1775         {
1776           sigma += amult - mult;
1777         }
1778       }
1779       /*
1780       // on periodic curves if this is the last knot
1781       // the multiplicity is added twice to account for the first knot
1782       if (k == Knots.Upper() && Periodic) {
1783         if (Add)
1784           sigma += amult;
1785         else
1786           sigma += amult - mult;
1787       }
1788       */
1789     }
1790     else {
1791       // not identic to existing knot
1792       if (amult > 0) {
1793         if (amult > Degree) amult = Degree;
1794         NbKnots++;
1795         sigma += amult;
1796       }
1797     }
1798     
1799     ak++;
1800   }
1801   
1802   // count the last knots
1803   while (k < Knots.Upper()) {
1804     k++;
1805     NbKnots++;
1806     sigma += Mults(k);
1807   }
1808
1809   if (Periodic) {
1810     //for periodic B-Spline the requirement is that multiplicites of the first
1811     //and last knots must be equal (see Geom_BSplineCurve constructor for
1812     //instance);
1813     //respectively AddMults() must meet this requirement if AddKnots() contains
1814     //knot(s) coincident with first or last
1815     NbPoles = sigma - aLastKnotMult;
1816   }
1817   else {
1818     NbPoles = sigma - Degree - 1;
1819   }
1820  
1821   return Standard_True;
1822 }
1823
1824 //=======================================================================
1825 //function : Copy
1826 //purpose  : copy reals from an array to an other
1827 //        
1828 //   NbValues are copied from OldPoles(OldFirst)
1829 //                 to    NewPoles(NewFirst)
1830 //
1831 //   Periodicity is handled.
1832 //   OldFirst and NewFirst are updated 
1833 //   to the position after the last copied pole.
1834 //
1835 //=======================================================================
1836
1837 static void Copy(const Standard_Integer      NbPoles,
1838                  Standard_Integer&           OldFirst,
1839                  const TColStd_Array1OfReal& OldPoles,
1840                  Standard_Integer&           NewFirst,
1841                  TColStd_Array1OfReal&       NewPoles)
1842 {
1843   // reset the index in the range for periodicity
1844
1845   OldFirst = OldPoles.Lower() + 
1846     (OldFirst - OldPoles.Lower()) % (OldPoles.Upper() - OldPoles.Lower() + 1);
1847
1848   NewFirst = NewPoles.Lower() + 
1849     (NewFirst - NewPoles.Lower()) % (NewPoles.Upper() - NewPoles.Lower() + 1);
1850
1851   // copy
1852   Standard_Integer i;
1853
1854   for (i = 1; i <= NbPoles; i++) {
1855     NewPoles(NewFirst) = OldPoles(OldFirst);
1856     OldFirst++;
1857     if (OldFirst > OldPoles.Upper()) OldFirst = OldPoles.Lower();
1858     NewFirst++;
1859     if (NewFirst > NewPoles.Upper()) NewFirst = NewPoles.Lower();
1860   }
1861 }
1862                       
1863 //=======================================================================
1864 //function : InsertKnots
1865 //purpose  : insert an array of knots and multiplicities
1866 //=======================================================================
1867
1868 void BSplCLib::InsertKnots
1869 (const Standard_Integer         Degree, 
1870  const Standard_Boolean         Periodic,
1871  const Standard_Integer         Dimension, 
1872  const TColStd_Array1OfReal&    Poles,  
1873  const TColStd_Array1OfReal&    Knots,    
1874  const TColStd_Array1OfInteger& Mults, 
1875  const TColStd_Array1OfReal&    AddKnots,    
1876  const TColStd_Array1OfInteger* AddMults, 
1877  TColStd_Array1OfReal&          NewPoles,
1878  TColStd_Array1OfReal&          NewKnots,    
1879  TColStd_Array1OfInteger&       NewMults, 
1880  const Standard_Real            Tolerance,
1881  const Standard_Boolean         Add)
1882 {
1883   Standard_Boolean addflat  = AddMults == NULL;
1884   
1885   Standard_Integer i,k,mult,firstmult;
1886   Standard_Integer index,kn,curnk,curk;
1887   Standard_Integer p,np, curp, curnp, length, depth;
1888   Standard_Real u;
1889   Standard_Integer need;
1890   Standard_Real Eps;
1891
1892   // -------------------
1893   // create local arrays
1894   // -------------------
1895
1896   Standard_Real *knots = new Standard_Real[2*Degree];
1897   Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
1898   
1899   //----------------------------
1900   // loop on the knots to insert
1901   //----------------------------
1902   
1903   curk   = Knots.Lower()-1;          // current position in Knots
1904   curnk  = NewKnots.Lower()-1;       // current position in NewKnots
1905   curp   = Poles.Lower();            // current position in Poles
1906   curnp  = NewPoles.Lower();         // current position in NewPoles
1907
1908   // NewKnots, NewMults, NewPoles contains the current state of the curve
1909
1910   // index is the first pole of the current curve for insertion schema
1911
1912   if (Periodic) index = -Mults(Mults.Lower());
1913   else          index = -Degree-1;
1914
1915   // on Periodic curves the first knot and the last knot are inserted later
1916   // (they are the same knot)
1917   firstmult = 0;  // multiplicity of the first-last knot for periodic
1918   
1919
1920   // kn current knot to insert in AddKnots
1921
1922   for (kn = AddKnots.Lower(); kn <= AddKnots.Upper(); kn++) {
1923     
1924     u = AddKnots(kn);
1925     Eps = Max(Tolerance,Epsilon(u));
1926     
1927     //-----------------------------------
1928     // find the position in the old knots
1929     // and copy to the new knots
1930     //-----------------------------------
1931     
1932     while (curk < Knots.Upper() && Knots(curk+1) - u <= Eps) {
1933       curk++; curnk++;
1934       NewKnots(curnk) = Knots(curk);
1935       index += NewMults(curnk) = Mults(curk);
1936     }
1937     
1938     //-----------------------------------
1939     // Slice the knots and the mults
1940     // to the current size of the new curve
1941     //-----------------------------------
1942
1943     i = curnk + Knots.Upper() - curk;
1944     TColStd_Array1OfReal    nknots(NewKnots(NewKnots.Lower()),NewKnots.Lower(),i);
1945     TColStd_Array1OfInteger nmults(NewMults(NewMults.Lower()),NewMults.Lower(),i);
1946
1947     //-----------------------------------
1948     // copy enough knots 
1949     // to compute the insertion schema
1950     //-----------------------------------
1951
1952     k = curk;
1953     i = curnk;
1954     mult = 0;
1955
1956     while (mult < Degree && k < Knots.Upper()) {
1957       k++; i++;
1958       nknots(i) = Knots(k);
1959       mult += nmults(i) = Mults(k);
1960     }
1961
1962     // copy knots at the end for periodic curve
1963     if (Periodic) {
1964       mult = 0;
1965       k = Knots.Upper();
1966       i = nknots.Upper();
1967
1968       while (mult < Degree && i > curnk) {
1969         nknots(i) = Knots(k);
1970         mult += nmults(i) = Mults(k);
1971         k--;
1972         i--;
1973       }
1974       nmults(nmults.Upper()) = nmults(nmults.Lower());
1975     }
1976
1977   
1978
1979     //------------------------------------
1980     // do the boor scheme on the new curve
1981     // to insert the new knot
1982     //------------------------------------
1983     
1984     Standard_Boolean sameknot = (Abs(u-NewKnots(curnk)) <= Eps);
1985     
1986     if (sameknot) length = Max(0,Degree - NewMults(curnk));
1987     else          length = Degree;
1988     
1989     if (addflat) depth = 1;
1990     else         depth = Min(Degree,(*AddMults)(kn));
1991
1992     if (sameknot) {
1993       if (Add) {
1994         if ((NewMults(curnk) + depth) > Degree)
1995           depth = Degree - NewMults(curnk);
1996       }
1997       else {
1998         depth = Max(0,depth-NewMults(curnk));
1999       }
2000
2001       if (Periodic) {
2002         // on periodic curve the first and last knot are delayed to the end
2003         if (curk == Knots.Lower() || (curk == Knots.Upper())) {
2004           if (firstmult == 0) // do that only once
2005             firstmult += depth;
2006           depth = 0;
2007         }
2008       }
2009     }
2010     if (depth <= 0) continue;
2011     
2012     BuildKnots(Degree,curnk,Periodic,nknots,&nmults,*knots);
2013
2014     // copy the poles
2015
2016     need   = NewPoles.Lower()+(index+length+1)*Dimension - curnp;
2017     need = Min(need,Poles.Upper() - curp + 1);
2018
2019     p = curp;
2020     np = curnp;
2021     Copy(need,p,Poles,np,NewPoles);
2022     curp  += need;
2023     curnp += need;
2024
2025     // slice the poles to the current number of poles in case of periodic
2026     TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
2027
2028     BuildBoor(index,length,Dimension,npoles,*poles);
2029     BoorScheme(u,Degree,*knots,Dimension,*poles,depth,length);
2030     
2031     //-------------------
2032     // copy the new poles
2033     //-------------------
2034
2035     curnp += depth * Dimension; // number of poles is increased by depth
2036     TColStd_Array1OfReal ThePoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
2037     np = NewKnots.Lower()+(index+1)*Dimension;
2038
2039     for (i = 1; i <= length + depth; i++)
2040       GetPole(i,length,depth,Dimension,*poles,np,ThePoles);
2041     
2042     //-------------------
2043     // insert the knot
2044     //-------------------
2045
2046     index += depth;
2047     if (sameknot) {
2048       NewMults(curnk) += depth;
2049     }
2050     else {
2051       curnk++;
2052       NewKnots(curnk) = u;
2053       NewMults(curnk) = depth;
2054     }
2055   }
2056   
2057   //------------------------------
2058   // copy the last poles and knots
2059   //------------------------------
2060   
2061   Copy(Poles.Upper() - curp + 1,curp,Poles,curnp,NewPoles);
2062   
2063   while (curk < Knots.Upper()) {
2064     curk++;  curnk++;
2065     NewKnots(curnk) = Knots(curk);
2066     NewMults(curnk) = Mults(curk);
2067   }
2068   
2069   //------------------------------
2070   // process the first-last knot 
2071   // on periodic curves
2072   //------------------------------
2073
2074   if (firstmult > 0) {
2075     curnk = NewKnots.Lower();
2076     if (NewMults(curnk) + firstmult > Degree) {
2077       firstmult = Degree - NewMults(curnk);
2078     }
2079     if (firstmult > 0) {
2080
2081       length = Degree - NewMults(curnk);
2082       depth  = firstmult;
2083
2084       BuildKnots(Degree,curnk,Periodic,NewKnots,&NewMults,*knots);
2085       TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),
2086                                   NewPoles.Lower(),
2087                                   NewPoles.Upper()-depth*Dimension);
2088       BuildBoor(0,length,Dimension,npoles,*poles);
2089       BoorScheme(NewKnots(curnk),Degree,*knots,Dimension,*poles,depth,length);
2090       
2091       //---------------------------
2092       // copy the new poles
2093       // but rotate them with depth
2094       //---------------------------
2095       
2096       np = NewPoles.Lower();
2097
2098       for (i = depth; i < length + depth; i++)
2099         GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2100
2101       np = NewPoles.Upper() - depth*Dimension + 1;
2102
2103       for (i = 0; i < depth; i++)
2104         GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2105       
2106       NewMults(NewMults.Lower()) += depth;
2107       NewMults(NewMults.Upper()) += depth;
2108     }
2109   }
2110   // free local arrays
2111   delete [] knots;
2112   delete [] poles;
2113 }
2114
2115 //=======================================================================
2116 //function : RemoveKnot
2117 //purpose  : 
2118 //=======================================================================
2119
2120 Standard_Boolean BSplCLib::RemoveKnot 
2121 (const Standard_Integer         Index,       
2122  const Standard_Integer         Mult,        
2123  const Standard_Integer         Degree,  
2124  const Standard_Boolean         Periodic,
2125  const Standard_Integer         Dimension,  
2126  const TColStd_Array1OfReal&    Poles,
2127  const TColStd_Array1OfReal&    Knots,  
2128  const TColStd_Array1OfInteger& Mults,
2129  TColStd_Array1OfReal&          NewPoles,
2130  TColStd_Array1OfReal&          NewKnots,  
2131  TColStd_Array1OfInteger&       NewMults,
2132  const Standard_Real            Tolerance)
2133 {
2134   Standard_Integer index,i,j,k,p,np;
2135
2136   Standard_Integer TheIndex = Index;
2137
2138   // protection
2139   Standard_Integer first,last;
2140   if (Periodic) {
2141     first = Knots.Lower();
2142     last  = Knots.Upper();
2143   }
2144   else {
2145     first = BSplCLib::FirstUKnotIndex(Degree,Mults) + 1;
2146     last  = BSplCLib::LastUKnotIndex(Degree,Mults) - 1;
2147   }
2148   if (Index < first) return Standard_False;
2149   if (Index > last)  return Standard_False;
2150
2151   if ( Periodic && (Index == first)) TheIndex = last;
2152
2153   Standard_Integer depth  = Mults(TheIndex) - Mult;
2154   Standard_Integer length = Degree - Mult;
2155
2156   // -------------------
2157   // create local arrays
2158   // -------------------
2159
2160   Standard_Real *knots = new Standard_Real[4*Degree];
2161   Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
2162   
2163
2164   // ------------------------------------
2165   // build the knots for anti Boor Scheme
2166   // ------------------------------------
2167
2168   // the new sequence of knots
2169   // is obtained from the knots at Index-1 and Index
2170   
2171   BSplCLib::BuildKnots(Degree,TheIndex-1,Periodic,Knots,&Mults,*knots);
2172   index = PoleIndex(Degree,TheIndex-1,Periodic,Mults);
2173   BSplCLib::BuildKnots(Degree,TheIndex,Periodic,Knots,&Mults,knots[2*Degree]);
2174
2175   index += Mult;
2176
2177   for (i = 0; i < Degree - Mult; i++)
2178     knots[i] = knots[i+Mult];
2179
2180   for (i = Degree-Mult; i < 2*Degree; i++)
2181     knots[i] = knots[2*Degree+i];
2182
2183
2184   // ------------------------------------
2185   // build the poles for anti Boor Scheme
2186   // ------------------------------------
2187
2188   p = Poles.Lower()+index * Dimension;
2189
2190   for (i = 0; i <= length + depth; i++) {
2191     j = Dimension * BoorIndex(i,length,depth);
2192
2193     for (k = 0; k < Dimension; k++) {
2194       poles[j+k] = Poles(p+k);
2195     }
2196     p += Dimension;
2197     if (p > Poles.Upper()) p = Poles.Lower();
2198   }
2199
2200
2201   // ----------------
2202   // Anti Boor Scheme
2203   // ----------------
2204
2205   Standard_Boolean result = AntiBoorScheme(Knots(TheIndex),Degree,*knots,
2206                                            Dimension,*poles,
2207                                            depth,length,Tolerance);
2208   
2209   // ----------------
2210   // copy the results
2211   // ----------------
2212
2213   if (result) {
2214
2215     // poles
2216
2217     p = Poles.Lower();
2218     np = NewPoles.Lower();
2219     
2220     // unmodified poles before
2221     Copy((index+1)*Dimension,p,Poles,np,NewPoles);
2222     
2223     
2224     // modified
2225
2226     for (i = 1; i <= length; i++)
2227       BSplCLib::GetPole(i,length,0,Dimension,*poles,np,NewPoles);
2228     p += (length + depth) * Dimension ;
2229     
2230     // unmodified poles after
2231     if (p != Poles.Lower()) {
2232       i = Poles.Upper() - p + 1;
2233       Copy(i,p,Poles,np,NewPoles);
2234     }
2235
2236     // knots and mults
2237
2238     if (Mult > 0) {
2239       NewKnots = Knots;
2240       NewMults = Mults;
2241       NewMults(TheIndex) = Mult;
2242       if (Periodic) {
2243         if (TheIndex == first) NewMults(last)  = Mult;
2244         if (TheIndex == last)  NewMults(first) = Mult;
2245       }
2246     }
2247     else {
2248       if (!Periodic || (TheIndex != first && TheIndex != last)) {
2249
2250         for (i = Knots.Lower(); i < TheIndex; i++) {
2251           NewKnots(i) = Knots(i);
2252           NewMults(i) = Mults(i);
2253         }
2254
2255         for (i = TheIndex+1; i <= Knots.Upper(); i++) {
2256           NewKnots(i-1) = Knots(i);
2257           NewMults(i-1) = Mults(i);
2258         }
2259       }
2260       else {
2261         // The interesting case of a Periodic curve 
2262         // where the first and last knot is removed.
2263         
2264         for (i = first; i < last-1; i++) {
2265           NewKnots(i) = Knots(i+1);
2266           NewMults(i) = Mults(i+1);
2267         }
2268         NewKnots(last-1) = NewKnots(first) + Knots(last) - Knots(first);
2269         NewMults(last-1) = NewMults(first);
2270       }
2271     }
2272   }
2273
2274
2275   // free local arrays
2276   delete [] knots;
2277   delete [] poles;
2278   
2279   return result;
2280 }
2281
2282 //=======================================================================
2283 //function : IncreaseDegreeCountKnots
2284 //purpose  : 
2285 //=======================================================================
2286
2287 Standard_Integer  BSplCLib::IncreaseDegreeCountKnots
2288 (const Standard_Integer Degree,
2289  const Standard_Integer NewDegree, 
2290  const Standard_Boolean Periodic, 
2291  const TColStd_Array1OfInteger& Mults)
2292 {
2293   if (Periodic) return Mults.Length();
2294   Standard_Integer f = FirstUKnotIndex(Degree,Mults);
2295   Standard_Integer l = LastUKnotIndex(Degree,Mults);
2296   Standard_Integer m,i,removed = 0, step = NewDegree - Degree;
2297   
2298   i = Mults.Lower();
2299   m = Degree + (f - i + 1) * step + 1;
2300
2301   while (m > NewDegree+1) {
2302     removed++;
2303     m -= Mults(i) + step;
2304     i++;
2305   }
2306   if (m < NewDegree+1) removed--;
2307
2308   i = Mults.Upper();
2309   m = Degree + (i - l + 1) * step + 1;
2310
2311   while (m > NewDegree+1) {
2312     removed++;
2313     m -= Mults(i) + step;
2314     i--;
2315   }
2316   if (m < NewDegree+1) removed--;
2317
2318   return Mults.Length() - removed;
2319 }
2320
2321 //=======================================================================
2322 //function : IncreaseDegree
2323 //purpose  : 
2324 //=======================================================================
2325
2326 void BSplCLib::IncreaseDegree 
2327 (const Standard_Integer         Degree,
2328  const Standard_Integer         NewDegree,
2329  const Standard_Boolean         Periodic,
2330  const Standard_Integer         Dimension,
2331  const TColStd_Array1OfReal&    Poles,
2332  const TColStd_Array1OfReal&    Knots,
2333  const TColStd_Array1OfInteger& Mults,
2334  TColStd_Array1OfReal&          NewPoles,
2335  TColStd_Array1OfReal&          NewKnots,
2336  TColStd_Array1OfInteger&       NewMults)
2337
2338   // Degree elevation of a BSpline Curve
2339
2340   // This algorithms loops on degree incrementation from Degree to NewDegree.
2341   // The variable curDeg is the current degree to increment.
2342
2343   // Before degree incrementations a "working curve" is created.
2344   // It has the same knot, poles and multiplicities.
2345
2346   // If the curve is periodic knots are added on the working curve before
2347   // and after the existing knots to make it a non-periodic curves. 
2348   // The poles are also copied.
2349
2350   // The first and last multiplicity of the working curve are set to Degree+1,
2351   // null poles are  added if necessary.
2352
2353   // Then the degree is incremented on the working curve.
2354   // The knots are unchanged but all multiplicities will be incremented.
2355
2356   // Each degree incrementation is achieved by averaging curDeg+1 curves.
2357
2358   // See : Degree elevation of B-spline curves
2359   //       Hartmut PRAUTZSCH
2360   //       CAGD 1 (1984)
2361
2362
2363   //-------------------------
2364   // create the working curve
2365   //-------------------------
2366
2367   Standard_Integer i,k,f,l,m,pf,pl,firstknot;
2368
2369   pf = 0; // number of null poles added at beginning
2370   pl = 0; // number of null poles added at end
2371
2372   Standard_Integer nbwknots = Knots.Length();
2373   f = FirstUKnotIndex(Degree,Mults);
2374   l = LastUKnotIndex (Degree,Mults);
2375
2376   if (Periodic) {
2377     // Periodic curves are transformed in non-periodic curves
2378
2379     nbwknots += f - Mults.Lower();
2380
2381     pf = -Degree - 1;
2382
2383     for (i = Mults.Lower(); i <= f; i++)
2384       pf += Mults(i);
2385
2386     nbwknots += Mults.Upper() - l;
2387
2388     pl = -Degree - 1;
2389
2390     for (i = l; i <= Mults.Upper(); i++)
2391       pl += Mults(i);
2392   }
2393
2394   // copy the knots and multiplicities 
2395   TColStd_Array1OfReal    wknots(1,nbwknots);
2396   TColStd_Array1OfInteger wmults(1,nbwknots);
2397   if (!Periodic) {
2398     wknots  = Knots;
2399     wmults  = Mults;
2400   }
2401   else {
2402     // copy the knots for a periodic curve
2403     Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2404     i = 0;
2405
2406     for (k = l; k < Knots.Upper(); k++) {
2407       i++; 
2408       wknots(i) = Knots(k) - period;
2409       wmults(i) = Mults(k);
2410     }
2411
2412     for (k = Knots.Lower(); k <= Knots.Upper(); k++) {
2413       i++; 
2414       wknots(i) = Knots(k);
2415       wmults(i) = Mults(k);
2416     }
2417
2418     for (k = Knots.Lower()+1; k <= f; k++) {
2419       i++; 
2420       wknots(i) = Knots(k) + period;
2421       wmults(i) = Mults(k);
2422     }
2423   }
2424
2425   // set the first and last mults to Degree+1
2426   // and add null poles
2427
2428   pf += Degree + 1 - wmults(1);
2429   wmults(1) = Degree + 1;
2430   pl += Degree + 1 - wmults(nbwknots);
2431   wmults(nbwknots) = Degree + 1;
2432
2433   //---------------------------
2434   // poles of the working curve
2435   //---------------------------
2436
2437   Standard_Integer nbwpoles = 0;
2438
2439   for (i = 1; i <= nbwknots; i++) nbwpoles += wmults(i);
2440   nbwpoles -= Degree + 1;
2441
2442   // we provide space for degree elevation
2443   TColStd_Array1OfReal 
2444     wpoles(1,(nbwpoles + (nbwknots-1) * (NewDegree - Degree)) * Dimension);
2445
2446   for (i = 1; i <= pf * Dimension; i++) 
2447     wpoles(i) = 0;
2448
2449   k = Poles.Lower();
2450
2451   for (i = pf * Dimension + 1; i <= (nbwpoles - pl) * Dimension; i++) {
2452     wpoles(i) = Poles(k);
2453     k++;
2454     if (k > Poles.Upper()) k = Poles.Lower();
2455   }
2456
2457   for (i = (nbwpoles-pl)*Dimension+1; i <= nbwpoles*Dimension; i++)
2458     wpoles(i) = 0;
2459   
2460   
2461   //-----------------------------------------------------------
2462   // Declare the temporary arrays used in degree incrementation
2463   //-----------------------------------------------------------
2464
2465   Standard_Integer nbwp = nbwpoles + (nbwknots-1) * (NewDegree - Degree);
2466   // Arrays for storing the temporary curves
2467   TColStd_Array1OfReal tempc1(1,nbwp * Dimension);
2468   TColStd_Array1OfReal tempc2(1,nbwp * Dimension);
2469
2470   // Array for storing the knots to insert
2471   TColStd_Array1OfReal iknots(1,nbwknots);
2472
2473   // Arrays for receiving the knots after insertion
2474   TColStd_Array1OfReal    nknots(1,nbwknots);
2475
2476
2477   
2478   //------------------------------
2479   // Loop on degree incrementation
2480   //------------------------------
2481
2482   Standard_Integer step,curDeg;
2483   Standard_Integer nbp = nbwpoles;
2484   nbwp = nbp;
2485
2486   for (curDeg = Degree; curDeg < NewDegree; curDeg++) {
2487     
2488     nbp  = nbwp;               // current number of poles
2489     nbwp = nbp + nbwknots - 1; // new number of poles
2490
2491     // For the averaging
2492     TColStd_Array1OfReal nwpoles(1,nbwp * Dimension);
2493     nwpoles.Init(0.0e0) ;
2494   
2495     
2496     for (step = 0; step <= curDeg; step++) {
2497     
2498       // Compute the bspline of rank step.
2499
2500       // if not the first time, decrement the multiplicities back
2501       if (step != 0) {
2502         for (i = 1; i <= nbwknots; i++)
2503           wmults(i)--;
2504       }
2505     
2506       // Poles are the current poles 
2507       // but the poles congruent to step are duplicated.
2508       
2509       Standard_Integer offset = 0;
2510
2511       for (i = 0; i < nbp; i++) {
2512         offset++;
2513
2514         for (k = 0; k < Dimension; k++) {
2515           tempc1((offset-1)*Dimension+k+1) = 
2516             wpoles(NewPoles.Lower()+i*Dimension+k);
2517         }
2518         if (i % (curDeg+1) == step) {
2519           offset++;
2520
2521           for (k = 0; k < Dimension; k++) {
2522             tempc1((offset-1)*Dimension+k+1) = 
2523               wpoles(NewPoles.Lower()+i*Dimension+k);
2524           }
2525         }
2526       }
2527         
2528       // Knots multiplicities are increased
2529       // For each knot where the sum of multiplicities is congruent to step
2530       
2531       Standard_Integer stepmult = step+1;
2532       Standard_Integer nbknots = 0;
2533       Standard_Integer smult = 0;
2534
2535       for (k = 1; k <= nbwknots; k++) {
2536         smult += wmults(k);
2537         if (smult  >= stepmult) {
2538           // this knot is increased
2539           stepmult += curDeg+1;
2540           wmults(k)++;
2541         }
2542         else {
2543           // this knot is inserted
2544           nbknots++;
2545           iknots(nbknots) = wknots(k);
2546         }
2547       }
2548       
2549       // the curve is obtained by inserting the knots
2550       // to raise the multiplicities
2551
2552       // we build "slices" of the arrays to set the correct size
2553       if (nbknots > 0) {
2554         TColStd_Array1OfReal aknots(iknots(1),1,nbknots);
2555         TColStd_Array1OfReal curve (tempc1(1),1,offset * Dimension);
2556         TColStd_Array1OfReal ncurve(tempc2(1),1,nbwp   * Dimension);
2557 //      InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2558 //                  aknots,NoMults(),ncurve,nknots,wmults,Epsilon(1.));
2559
2560         InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2561                     aknots,NoMults(),ncurve,nknots,wmults,0.0);
2562         
2563         // add to the average
2564
2565         for (i = 1; i <= nbwp * Dimension; i++)
2566           nwpoles(i) += ncurve(i);
2567       }
2568       else {
2569         // add to the average
2570
2571         for (i = 1; i <= nbwp * Dimension; i++)
2572           nwpoles(i) += tempc1(i);
2573       }
2574     }
2575     
2576     // The result is the average
2577
2578     for (i = 1; i <= nbwp * Dimension; i++) {
2579       wpoles(i) = nwpoles(i) / (curDeg+1);
2580     }
2581   }
2582   
2583   //-----------------
2584   // Copy the results
2585   //-----------------
2586
2587   // index in new knots of the first knot of the curve
2588   if (Periodic)
2589     firstknot = Mults.Upper() - l + 1;
2590   else 
2591     firstknot = f;
2592   
2593   // the new curve starts at index firstknot
2594   // so we must remove knots until the sum of multiplicities
2595   // from the first to the start is NewDegree+1
2596
2597   // m is the current sum of multiplicities
2598   m = 0;
2599
2600   for (k = 1; k <= firstknot; k++)
2601     m += wmults(k);
2602
2603   // compute the new first knot (k), pf will be the index of the first pole
2604   k = 1;
2605   pf = 0;
2606
2607   while (m > NewDegree+1) {
2608     k++;
2609     m  -= wmults(k);
2610     pf += wmults(k);
2611   }
2612   if (m < NewDegree+1) {
2613     k--;
2614     wmults(k) += m - NewDegree - 1;
2615     pf        += m - NewDegree - 1;
2616   }
2617
2618   // on a periodic curve the knots start with firstknot
2619   if (Periodic)
2620     k = firstknot;
2621
2622   // copy knots
2623
2624   for (i = NewKnots.Lower(); i <= NewKnots.Upper(); i++) {
2625     NewKnots(i) = wknots(k);
2626     NewMults(i) = wmults(k);
2627     k++;
2628   }
2629
2630   // copy poles
2631   pf *= Dimension;
2632
2633   for (i = NewPoles.Lower(); i <= NewPoles.Upper(); i++) {
2634     pf++;
2635     NewPoles(i) = wpoles(pf);
2636   }
2637 }
2638
2639 //=======================================================================
2640 //function : PrepareUnperiodize
2641 //purpose  : 
2642 //=======================================================================
2643
2644 void  BSplCLib::PrepareUnperiodize
2645 (const Standard_Integer         Degree, 
2646  const TColStd_Array1OfInteger& Mults, 
2647  Standard_Integer&        NbKnots, 
2648  Standard_Integer&        NbPoles)
2649 {
2650   Standard_Integer i;
2651   // initialize NbKnots and NbPoles
2652   NbKnots = Mults.Length();
2653   NbPoles = - Degree - 1;
2654
2655   for (i = Mults.Lower(); i <= Mults.Upper(); i++) 
2656     NbPoles += Mults(i);
2657
2658   Standard_Integer sigma, k;
2659   // Add knots at the beginning of the curve to raise Multiplicities 
2660   // to Degre + 1;
2661   sigma = Mults(Mults.Lower());
2662   k = Mults.Upper() - 1;
2663
2664   while ( sigma < Degree + 1) {
2665     sigma   += Mults(k);
2666     NbPoles += Mults(k);
2667     k--;
2668     NbKnots++;
2669   }
2670   // We must add exactly until Degree + 1 -> 
2671   //    Supress the excedent.
2672   if ( sigma > Degree + 1)
2673     NbPoles -= sigma - Degree - 1;
2674
2675   // Add knots at the end of the curve to raise Multiplicities 
2676   // to Degre + 1;
2677   sigma = Mults(Mults.Upper());
2678   k = Mults.Lower() + 1;
2679
2680   while ( sigma < Degree + 1) {
2681     sigma   += Mults(k);
2682     NbPoles += Mults(k);
2683     k++;
2684     NbKnots++;
2685   }
2686   // We must add exactly until Degree + 1 -> 
2687   //    Supress the excedent.
2688   if ( sigma > Degree + 1)
2689     NbPoles -= sigma - Degree - 1;
2690 }
2691
2692 //=======================================================================
2693 //function : Unperiodize
2694 //purpose  : 
2695 //=======================================================================
2696
2697 void  BSplCLib::Unperiodize
2698 (const Standard_Integer         Degree,
2699  const Standard_Integer         , // Dimension,
2700  const TColStd_Array1OfInteger& Mults,
2701  const TColStd_Array1OfReal&    Knots,
2702  const TColStd_Array1OfReal&    Poles,
2703  TColStd_Array1OfInteger& NewMults,
2704  TColStd_Array1OfReal&    NewKnots,
2705  TColStd_Array1OfReal&    NewPoles)
2706 {
2707   Standard_Integer sigma, k, index = 0;
2708   // evaluation of index : number of knots to insert before knot(1) to
2709   // raise sum of multiplicities to <Degree + 1>
2710   sigma = Mults(Mults.Lower());
2711   k = Mults.Upper() - 1;
2712
2713   while ( sigma < Degree + 1) {
2714     sigma   += Mults(k);
2715     k--;
2716     index++;
2717   }
2718
2719   Standard_Real    period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2720
2721   // set the 'interior' knots;
2722
2723   for ( k = 1; k <= Knots.Length(); k++) {
2724     NewKnots ( k + index ) = Knots( k);
2725     NewMults ( k + index ) = Mults( k);
2726   }
2727   
2728   // set the 'starting' knots;
2729
2730   for ( k = 1; k <= index; k++) {
2731     NewKnots( k) = NewKnots( k + Knots.Length() - 1) - period;
2732     NewMults( k) = NewMults( k + Knots.Length() - 1);
2733   }
2734   NewMults( 1) -= sigma - Degree -1;
2735   
2736   // set the 'ending' knots;
2737   sigma = NewMults( index + Knots.Length() );
2738
2739   for ( k = Knots.Length() + index + 1; k <= NewKnots.Length(); k++) {
2740     NewKnots( k) = NewKnots( k - Knots.Length() + 1) + period;
2741     NewMults( k) = NewMults( k - Knots.Length() + 1);
2742     sigma += NewMults( k - Knots.Length() + 1);
2743   }
2744   NewMults(NewMults.Length()) -= sigma - Degree - 1;
2745
2746   for ( k = 1 ; k <= NewPoles.Length(); k++) {
2747     NewPoles(k ) = Poles( (k-1) % Poles.Length() + 1);
2748   }
2749 }
2750
2751 //=======================================================================
2752 //function : PrepareTrimming
2753 //purpose  : 
2754 //=======================================================================
2755
2756 void BSplCLib::PrepareTrimming(const Standard_Integer         Degree,
2757                                const Standard_Boolean         Periodic,
2758                                const TColStd_Array1OfReal&    Knots,
2759                                const TColStd_Array1OfInteger& Mults,
2760                                const Standard_Real            U1,
2761                                const Standard_Real            U2,
2762                                      Standard_Integer&        NbKnots,
2763                                      Standard_Integer&        NbPoles)
2764 {
2765   Standard_Integer i;
2766   Standard_Real NewU1, NewU2;
2767   Standard_Integer index1 = 0, index2 = 0;
2768
2769   // Eval index1, index2 : position of U1 and U2 in the Array Knots
2770   // such as Knots(index1-1) <= U1 < Knots(index1)
2771   //         Knots(index2-1) <= U2 < Knots(index2)
2772   LocateParameter( Degree, Knots, Mults, U1, Periodic,
2773                    Knots.Lower(), Knots.Upper(), index1, NewU1);
2774   LocateParameter( Degree, Knots, Mults, U2, Periodic,
2775                    Knots.Lower(), Knots.Upper(), index2, NewU2);
2776   index1++;
2777   if ( Abs(Knots(index2) - U2) <= Epsilon( U1))
2778     index2--;
2779
2780   // eval NbKnots:
2781   NbKnots = index2 - index1 + 3;
2782
2783   // eval NbPoles:
2784   NbPoles = Degree + 1;
2785
2786   for ( i = index1; i <= index2; i++) 
2787     NbPoles += Mults(i);
2788 }
2789
2790 //=======================================================================
2791 //function : Trimming
2792 //purpose  : 
2793 //=======================================================================
2794
2795 void BSplCLib::Trimming(const Standard_Integer         Degree,
2796                         const Standard_Boolean         Periodic,
2797                         const Standard_Integer         Dimension,
2798                         const TColStd_Array1OfReal&    Knots,
2799                         const TColStd_Array1OfInteger& Mults,
2800                         const TColStd_Array1OfReal&    Poles,
2801                         const Standard_Real            U1,
2802                         const Standard_Real            U2,
2803                               TColStd_Array1OfReal&    NewKnots,
2804                               TColStd_Array1OfInteger& NewMults,
2805                               TColStd_Array1OfReal&    NewPoles)
2806 {
2807   Standard_Integer i, nbpoles=0, nbknots=0;
2808   Standard_Real kk[2];
2809   Standard_Integer mm[2];
2810   TColStd_Array1OfReal    K( kk[0], 1, 2 );
2811   TColStd_Array1OfInteger M( mm[0], 1, 2 );
2812
2813   K(1) = U1;  K(2) = U2;
2814   mm[0] = mm[1] = Degree;
2815   if (!PrepareInsertKnots( Degree, Periodic, Knots, Mults, K, &M, 
2816                           nbpoles, nbknots, Epsilon( U1), 0))
2817     throw Standard_OutOfRange();
2818
2819   TColStd_Array1OfReal    TempPoles(1, nbpoles*Dimension);
2820   TColStd_Array1OfReal    TempKnots(1, nbknots);
2821   TColStd_Array1OfInteger TempMults(1, nbknots);
2822
2823 //
2824 // do not allow the multiplicities to Add : they must be less than Degree
2825 //
2826   InsertKnots(Degree, Periodic, Dimension, Poles, Knots, Mults,
2827               K, &M, TempPoles, TempKnots, TempMults, Epsilon(U1),
2828               Standard_False);
2829
2830   // find in TempPoles the index of the pole corresponding to U1
2831   Standard_Integer Kindex = 0, Pindex;
2832   Standard_Real NewU1;
2833   LocateParameter( Degree, TempKnots, TempMults, U1, Periodic,
2834                    TempKnots.Lower(), TempKnots.Upper(), Kindex, NewU1);
2835   Pindex = PoleIndex ( Degree, Kindex, Periodic, TempMults);
2836   Pindex *= Dimension;
2837
2838   for ( i = 1; i <= NewPoles.Length(); i++) NewPoles(i) = TempPoles(Pindex + i);
2839
2840   for ( i = 1; i <= NewKnots.Length(); i++) {
2841     NewKnots(i) = TempKnots( Kindex+i-1);
2842     NewMults(i) = TempMults( Kindex+i-1);
2843   }
2844   NewMults(1) = Min(Degree, NewMults(1)) + 1 ;
2845   NewMults(NewMults.Length())= Min(Degree, NewMults(NewMults.Length())) + 1 ;
2846 }
2847
2848 //=======================================================================
2849 //function : Solves a LU factored Matrix 
2850 //purpose  : 
2851 //=======================================================================
2852
2853 Standard_Integer 
2854 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
2855                             const Standard_Integer UpperBandWidth,
2856                             const Standard_Integer LowerBandWidth,
2857                             const Standard_Integer ArrayDimension,
2858                             Standard_Real&   Array) 
2859 {
2860   Standard_Integer ii,
2861   jj,
2862   kk,
2863   MinIndex,
2864   MaxIndex,
2865   ReturnCode = 0 ;
2866   
2867   Standard_Real   *PolesArray = &Array ;
2868   Standard_Real   Inverse ;
2869   
2870   
2871   if (Matrix.LowerCol() != 1 || 
2872       Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2873     ReturnCode = 1 ;
2874     goto FINISH ;
2875   }
2876   
2877   for (ii = Matrix.LowerRow() + 1; ii <=  Matrix.UpperRow() ; ii++) {
2878     MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ?
2879                 ii - LowerBandWidth : Matrix.LowerRow()) ;
2880     
2881     for ( jj = MinIndex  ; jj < ii  ; jj++) {
2882       
2883       for (kk = 0 ; kk < ArrayDimension ; kk++) {
2884         PolesArray[(ii-1) * ArrayDimension + kk] += 
2885           PolesArray[(jj-1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2886       }
2887     }
2888   }
2889   
2890   for (ii = Matrix.UpperRow() ; ii >=  Matrix.LowerRow() ; ii--) {
2891     MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ? 
2892                 ii + UpperBandWidth : Matrix.UpperRow()) ;
2893     
2894     for (jj = MaxIndex  ; jj > ii ; jj--) {
2895       
2896       for (kk = 0 ; kk < ArrayDimension ; kk++) {
2897         PolesArray[(ii-1)  * ArrayDimension + kk] -=
2898           PolesArray[(jj - 1) * ArrayDimension + kk] * 
2899             Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2900       }
2901     }
2902     
2903     //fixing a bug PRO18577 to avoid divizion by zero
2904     
2905     Standard_Real divizor = Matrix(ii,LowerBandWidth + 1) ;
2906     Standard_Real Toler = 1.0e-16;
2907     if ( Abs(divizor) > Toler )
2908       Inverse = 1.0e0 / divizor ;
2909     else {
2910       Inverse = 1.0e0;
2911 //      cout << "  BSplCLib::SolveBandedSystem() : zero determinant " << endl;
2912       ReturnCode = 1;
2913       goto FINISH;
2914     }
2915         
2916     for (kk = 0 ; kk < ArrayDimension ; kk++) {
2917       PolesArray[(ii-1)  * ArrayDimension + kk] *=  Inverse ; 
2918     }
2919   }
2920   FINISH :
2921     return (ReturnCode) ;
2922 }
2923
2924 //=======================================================================
2925 //function : Solves a LU factored Matrix 
2926 //purpose  : 
2927 //=======================================================================
2928
2929 Standard_Integer 
2930 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
2931                             const Standard_Integer UpperBandWidth,
2932                             const Standard_Integer LowerBandWidth,
2933                             const Standard_Boolean HomogeneousFlag,
2934                             const Standard_Integer ArrayDimension,
2935                             Standard_Real&   Poles,
2936                             Standard_Real&   Weights) 
2937 {
2938   Standard_Integer ii,
2939   kk,
2940   ErrorCode = 0,
2941   ReturnCode = 0 ;
2942   
2943   Standard_Real   Inverse,
2944   *PolesArray   = &Poles,
2945   *WeightsArray = &Weights ;
2946   
2947   if (Matrix.LowerCol() != 1 || 
2948       Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2949     ReturnCode = 1 ;
2950     goto FINISH ;
2951   }
2952   if (HomogeneousFlag == Standard_False) {
2953     
2954     for (ii = 0 ; ii <  Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) {
2955       
2956       for (kk = 0 ; kk < ArrayDimension ; kk++) {
2957         PolesArray[ii * ArrayDimension + kk] *=
2958           WeightsArray[ii] ;
2959       }
2960     }
2961   }
2962   ErrorCode = 
2963     BSplCLib::SolveBandedSystem(Matrix,
2964                                 UpperBandWidth,
2965                                 LowerBandWidth,
2966                                 ArrayDimension,
2967                                 Poles) ;
2968   if (ErrorCode != 0) {
2969     ReturnCode = 2 ;
2970     goto FINISH ;
2971   }
2972   ErrorCode = 
2973     BSplCLib::SolveBandedSystem(Matrix,
2974                                 UpperBandWidth,
2975                                 LowerBandWidth,
2976                                 1,
2977                                 Weights) ;
2978   if (ErrorCode != 0) {
2979     ReturnCode = 3 ;
2980     goto FINISH ;
2981   }
2982   if (HomogeneousFlag == Standard_False) {
2983
2984     for (ii = 0  ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1 ; ii++) {
2985       Inverse = 1.0e0 / WeightsArray[ii] ;
2986       
2987       for (kk = 0  ; kk < ArrayDimension ; kk++) {
2988         PolesArray[ii * ArrayDimension + kk] *= Inverse ;
2989       }
2990     }
2991   }
2992   FINISH : return (ReturnCode) ;
2993 }
2994
2995 //=======================================================================
2996 //function : BuildSchoenbergPoints
2997 //purpose  : 
2998 //=======================================================================
2999
3000 void  BSplCLib::BuildSchoenbergPoints(const Standard_Integer         Degree,
3001                                       const TColStd_Array1OfReal&    FlatKnots,
3002                                       TColStd_Array1OfReal&          Parameters) 
3003 {
3004   Standard_Integer ii,
3005   jj ;
3006   Standard_Real Inverse ;
3007   Inverse = 1.0e0 / (Standard_Real)Degree ;
3008   
3009   for (ii = Parameters.Lower() ;   ii <= Parameters.Upper() ; ii++) {
3010     Parameters(ii) = 0.0e0 ;
3011     
3012     for (jj = 1 ; jj <= Degree ; jj++) {
3013       Parameters(ii) += FlatKnots(jj + ii) ;
3014     } 
3015     Parameters(ii) *= Inverse ; 
3016   }
3017 }
3018
3019 //=======================================================================
3020 //function : Interpolate
3021 //purpose  : 
3022 //=======================================================================
3023
3024 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
3025                             const TColStd_Array1OfReal&    FlatKnots,
3026                             const TColStd_Array1OfReal&    Parameters,
3027                             const TColStd_Array1OfInteger& ContactOrderArray,
3028                             const Standard_Integer         ArrayDimension,
3029                             Standard_Real&                 Poles,
3030                             Standard_Integer&              InversionProblem) 
3031 {
3032   Standard_Integer ErrorCode,
3033   UpperBandWidth,
3034   LowerBandWidth ;
3035 //  Standard_Real *PolesArray = &Poles ;
3036   math_Matrix InterpolationMatrix(1, Parameters.Length(),
3037                                   1, 2 * Degree + 1) ;
3038   ErrorCode =
3039   BSplCLib::BuildBSpMatrix(Parameters,
3040                            ContactOrderArray,
3041                            FlatKnots,
3042                            Degree,
3043                            InterpolationMatrix,
3044                            UpperBandWidth,
3045                            LowerBandWidth) ;
3046   if(ErrorCode)
3047     throw Standard_OutOfRange("BSplCLib::Interpolate");
3048
3049   ErrorCode =
3050   BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3051                            UpperBandWidth,
3052                            LowerBandWidth,
3053                            InversionProblem) ;
3054   if(ErrorCode)
3055     throw Standard_OutOfRange("BSplCLib::Interpolate");
3056
3057   ErrorCode  =
3058   BSplCLib::SolveBandedSystem(InterpolationMatrix,
3059                               UpperBandWidth,
3060                               LowerBandWidth,
3061                               ArrayDimension,
3062                               Poles) ;
3063   if(ErrorCode)
3064     throw Standard_OutOfRange("BSplCLib::Interpolate");
3065 }
3066
3067 //=======================================================================
3068 //function : Interpolate
3069 //purpose  : 
3070 //=======================================================================
3071
3072 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
3073                             const TColStd_Array1OfReal&    FlatKnots,
3074                             const TColStd_Array1OfReal&    Parameters,
3075                             const TColStd_Array1OfInteger& ContactOrderArray,
3076                             const Standard_Integer         ArrayDimension,
3077                             Standard_Real&                 Poles,
3078                             Standard_Real&                 Weights,
3079                             Standard_Integer&              InversionProblem) 
3080 {
3081   Standard_Integer ErrorCode,
3082   UpperBandWidth,
3083   LowerBandWidth ;
3084
3085   math_Matrix InterpolationMatrix(1, Parameters.Length(),
3086                                   1, 2 * Degree + 1) ;
3087   ErrorCode =
3088   BSplCLib::BuildBSpMatrix(Parameters,
3089                            ContactOrderArray,
3090                            FlatKnots,
3091                            Degree,
3092                            InterpolationMatrix,
3093                            UpperBandWidth,
3094                            LowerBandWidth) ;
3095   if(ErrorCode)
3096     throw Standard_OutOfRange("BSplCLib::Interpolate");
3097
3098   ErrorCode =
3099   BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3100                            UpperBandWidth,
3101                            LowerBandWidth,
3102                            InversionProblem) ;
3103   if(ErrorCode)
3104     throw Standard_OutOfRange("BSplCLib::Interpolate");
3105
3106   ErrorCode  =
3107   BSplCLib::SolveBandedSystem(InterpolationMatrix,
3108                               UpperBandWidth,
3109                               LowerBandWidth,
3110                               Standard_False,
3111                               ArrayDimension,
3112                               Poles,
3113                               Weights) ;
3114   if(ErrorCode)
3115     throw Standard_OutOfRange("BSplCLib::Interpolate");
3116 }
3117
3118 //=======================================================================
3119 //function : Evaluates a Bspline function : uses the ExtrapMode 
3120 //purpose  : the function is extrapolated using the Taylor expansion
3121 //           of degree ExtrapMode[0] to the left and the Taylor
3122 //           expansion of degree ExtrapMode[1] to the right 
3123 //  this evaluates the numerator by multiplying by the weights
3124 //  and evaluating it but does not call RationalDerivatives after 
3125 //=======================================================================
3126
3127 void  BSplCLib::Eval
3128 (const Standard_Real                   Parameter,
3129  const Standard_Boolean                PeriodicFlag,
3130  const Standard_Integer                DerivativeRequest,
3131  Standard_Integer&                     ExtrapMode,
3132  const Standard_Integer                Degree,
3133  const  TColStd_Array1OfReal&          FlatKnots, 
3134  const Standard_Integer                ArrayDimension,
3135  Standard_Real&                        Poles,
3136  Standard_Real&                        Weights,
3137  Standard_Real&                        PolesResults,
3138  Standard_Real&                        WeightsResults)
3139 {
3140   Standard_Integer ii,
3141   jj,
3142   kk=0,
3143   Index,
3144   Index1,
3145   Index2,
3146   *ExtrapModeArray,
3147   Modulus,
3148   NewRequest,
3149   ExtrapolatingFlag[2],
3150   ErrorCode,
3151   Order = Degree + 1,
3152   FirstNonZeroBsplineIndex,
3153   LocalRequest = DerivativeRequest ;
3154   Standard_Real  *PResultArray,
3155   *WResultArray,
3156   *PolesArray,
3157   *WeightsArray,
3158   LocalParameter,
3159   Period,
3160   Inverse,
3161   Delta ;
3162   PolesArray = &Poles     ;
3163   WeightsArray = &Weights ;
3164   ExtrapModeArray = &ExtrapMode ;
3165   PResultArray = &PolesResults ;
3166   WResultArray = &WeightsResults ;
3167   LocalParameter = Parameter ;
3168   ExtrapolatingFlag[0] = 
3169     ExtrapolatingFlag[1] = 0 ;
3170   //
3171   // check if we are extrapolating to a degree which is smaller than
3172   // the degree of the Bspline
3173   //
3174   if (PeriodicFlag) {
3175     Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3176
3177     while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3178       LocalParameter -= Period ;
3179     }
3180     
3181     while (LocalParameter < FlatKnots(2)) {
3182       LocalParameter +=  Period ;
3183     }
3184   }
3185   if (Parameter < FlatKnots(2) && 
3186       LocalRequest < ExtrapModeArray[0] &&
3187       ExtrapModeArray[0] < Degree) {
3188     LocalRequest = ExtrapModeArray[0] ;
3189     LocalParameter = FlatKnots(2) ;
3190     ExtrapolatingFlag[0] = 1 ;
3191   }
3192   if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3193       LocalRequest < ExtrapModeArray[1]  &&
3194       ExtrapModeArray[1] < Degree) {
3195     LocalRequest = ExtrapModeArray[1] ;
3196     LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3197     ExtrapolatingFlag[1] = 1 ;
3198   }
3199   Delta = Parameter - LocalParameter ;
3200   if (LocalRequest >= Order) {
3201     LocalRequest = Degree ;
3202   }
3203   if (PeriodicFlag) {
3204     Modulus = FlatKnots.Length() - Degree -1 ;
3205   }
3206   else {
3207     Modulus = FlatKnots.Length() - Degree ;
3208   }
3209
3210   BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3211   ErrorCode =
3212     BSplCLib::EvalBsplineBasis(LocalRequest,
3213                                Order,
3214                                FlatKnots,
3215                                LocalParameter,
3216                                FirstNonZeroBsplineIndex,
3217                                BsplineBasis) ;
3218   if (ErrorCode != 0) {
3219     goto FINISH ;
3220   }
3221   if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3222     Index = 0 ;
3223     Index2 = 0 ;
3224
3225     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3226       Index1 = FirstNonZeroBsplineIndex ;
3227
3228       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3229         PResultArray[Index + kk] = 0.0e0 ;
3230       }
3231       WResultArray[Index] = 0.0e0 ;
3232
3233       for (jj = 1  ; jj <= Order ; jj++) {
3234         
3235         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3236           PResultArray[Index + kk] += 
3237             PolesArray[(Index1-1) * ArrayDimension + kk] 
3238               * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3239         }
3240         WResultArray[Index2]  += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3241         
3242         Index1 = Index1 % Modulus ;
3243         Index1 += 1 ;
3244       }
3245       Index += ArrayDimension ;
3246       Index2 += 1 ;
3247     }
3248   }
3249   else {
3250     // 
3251     //  store Taylor expansion in LocalRealArray
3252     //
3253     NewRequest = DerivativeRequest ;
3254     if (NewRequest > Degree) {
3255       NewRequest = Degree ;
3256     }
3257     NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3258     Index = 0 ;
3259     Inverse = 1.0e0 ;
3260
3261     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3262       Index1 = FirstNonZeroBsplineIndex ;
3263       
3264       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3265         LocalRealArray[Index + kk] = 0.0e0 ;
3266       }
3267
3268       for (jj = 1  ; jj <= Order ; jj++) {
3269
3270         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3271           LocalRealArray[Index + kk] += 
3272             PolesArray[(Index1-1)*ArrayDimension + kk] * 
3273               WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3274         }
3275         Index1 = Index1 % Modulus ;
3276         Index1 += 1 ;
3277       }
3278
3279       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3280         LocalRealArray[Index + kk] *= Inverse ;
3281       }
3282       Index += ArrayDimension ;
3283       Inverse /= (Standard_Real) ii ;
3284     }
3285     PLib::EvalPolynomial(Delta,
3286                          NewRequest,
3287                          Degree,
3288                          ArrayDimension,
3289                          LocalRealArray[0],
3290                          PolesResults) ;
3291     Index = 0 ;
3292     Inverse = 1.0e0 ;
3293
3294     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3295       Index1 = FirstNonZeroBsplineIndex ;
3296       LocalRealArray[Index] = 0.0e0 ;
3297
3298       for (jj = 1  ; jj <= Order ; jj++) {
3299         LocalRealArray[Index] += 
3300           WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3301         Index1 = Index1 % Modulus ;
3302         Index1 += 1 ;
3303       }
3304       LocalRealArray[Index + kk] *= Inverse ;
3305       Index += 1 ;
3306       Inverse /= (Standard_Real) ii ;
3307     }
3308     PLib::EvalPolynomial(Delta,
3309                          NewRequest,
3310                          Degree,
3311                          1,
3312                          LocalRealArray[0],
3313                          WeightsResults) ;
3314   }
3315   FINISH : ;
3316 }
3317
3318 //=======================================================================
3319 //function : Evaluates a Bspline function : uses the ExtrapMode 
3320 //purpose  : the function is extrapolated using the Taylor expansion
3321 //           of degree ExtrapMode[0] to the left and the Taylor
3322 //           expansion of degree ExtrapMode[1] to the right 
3323 // WARNING : the array Results is supposed to have at least 
3324 // (DerivativeRequest + 1) * ArrayDimension slots and the 
3325 // 
3326 //=======================================================================
3327
3328 void  BSplCLib::Eval
3329 (const Standard_Real                   Parameter,
3330  const Standard_Boolean                PeriodicFlag,
3331  const Standard_Integer                DerivativeRequest,
3332  Standard_Integer&                     ExtrapMode,
3333  const Standard_Integer                Degree,
3334  const  TColStd_Array1OfReal&          FlatKnots, 
3335  const Standard_Integer                ArrayDimension,
3336  Standard_Real&                        Poles,
3337  Standard_Real&                        Results) 
3338 {
3339   Standard_Integer ii,
3340   jj,
3341   kk,
3342   Index,
3343   Index1,
3344   *ExtrapModeArray,
3345   Modulus,
3346   NewRequest,
3347   ExtrapolatingFlag[2],
3348   ErrorCode,
3349   Order = Degree + 1,
3350   FirstNonZeroBsplineIndex,
3351   LocalRequest = DerivativeRequest ;
3352
3353   Standard_Real  *ResultArray,
3354   *PolesArray,
3355   LocalParameter,
3356   Period,
3357   Inverse,
3358   Delta ;
3359          
3360   PolesArray = &Poles ;
3361   ExtrapModeArray = &ExtrapMode ;
3362   ResultArray = &Results ;  
3363   LocalParameter = Parameter ;
3364   ExtrapolatingFlag[0] = 
3365     ExtrapolatingFlag[1] = 0 ;
3366   //
3367   // check if we are extrapolating to a degree which is smaller than
3368   // the degree of the Bspline
3369   //
3370   if (PeriodicFlag) {
3371     Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3372
3373     while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3374       LocalParameter -= Period ;
3375     }
3376
3377     while (LocalParameter < FlatKnots(2)) {
3378       LocalParameter +=  Period ;
3379     }
3380   }
3381   if (Parameter < FlatKnots(2) && 
3382       LocalRequest < ExtrapModeArray[0] &&
3383       ExtrapModeArray[0] < Degree) {
3384     LocalRequest = ExtrapModeArray[0] ;
3385     LocalParameter = FlatKnots(2) ;
3386     ExtrapolatingFlag[0] = 1 ;
3387   }
3388   if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3389       LocalRequest < ExtrapModeArray[1]  &&
3390       ExtrapModeArray[1] < Degree) {
3391     LocalRequest = ExtrapModeArray[1] ;
3392     LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3393     ExtrapolatingFlag[1] = 1 ;
3394   }
3395   Delta = Parameter - LocalParameter ;
3396   if (LocalRequest >= Order) {
3397     LocalRequest = Degree ;
3398   }
3399   
3400   if (PeriodicFlag) {
3401     Modulus = FlatKnots.Length() - Degree -1 ;
3402   }
3403   else {
3404     Modulus = FlatKnots.Length() - Degree ;
3405   }
3406   
3407   BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3408   
3409   ErrorCode =
3410     BSplCLib::EvalBsplineBasis(LocalRequest,
3411                                Order,
3412                                FlatKnots,
3413                                LocalParameter,
3414                                FirstNonZeroBsplineIndex,
3415                                BsplineBasis);
3416   if (ErrorCode != 0) {
3417     goto FINISH ;
3418   }
3419   if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3420     Index = 0 ;
3421     
3422     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3423       Index1 = FirstNonZeroBsplineIndex ;
3424       
3425       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3426         ResultArray[Index + kk] = 0.0e0 ;
3427       }
3428
3429       for (jj = 1  ; jj <= Order ; jj++) {
3430         
3431         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3432           ResultArray[Index + kk] += 
3433             PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3434         }
3435         Index1 = Index1 % Modulus ;
3436         Index1 += 1 ;
3437       }
3438       Index += ArrayDimension ;
3439     }
3440   }
3441   else {
3442     // 
3443     //  store Taylor expansion in LocalRealArray
3444     //
3445     NewRequest = DerivativeRequest ;
3446     if (NewRequest > Degree) {
3447       NewRequest = Degree ;
3448     }
3449     NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3450
3451     Index = 0 ;
3452     Inverse = 1.0e0 ;
3453
3454     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3455       Index1 = FirstNonZeroBsplineIndex ;
3456       
3457       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3458         LocalRealArray[Index + kk] = 0.0e0 ;
3459       }
3460
3461       for (jj = 1  ; jj <= Order ; jj++) {
3462
3463         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3464           LocalRealArray[Index + kk] += 
3465             PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3466         }
3467         Index1 = Index1 % Modulus ;
3468         Index1 += 1 ;
3469       }
3470
3471       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3472         LocalRealArray[Index + kk] *= Inverse ;
3473       }
3474       Index += ArrayDimension ;
3475       Inverse /= (Standard_Real) ii ;
3476     }
3477     PLib::EvalPolynomial(Delta,
3478                          NewRequest,
3479                          Degree,
3480                          ArrayDimension,
3481                          LocalRealArray[0],
3482                          Results) ;
3483   }
3484   FINISH : ;
3485 }
3486
3487 //=======================================================================
3488 //function : TangExtendToConstraint 
3489 //purpose  : Extends a Bspline function using the tangency map
3490 // WARNING :  
3491 //  
3492 // 
3493 //=======================================================================
3494
3495 void  BSplCLib::TangExtendToConstraint
3496 (const  TColStd_Array1OfReal&          FlatKnots, 
3497  const Standard_Real                   C1Coefficient,
3498  const Standard_Integer                NumPoles,
3499  Standard_Real&                        Poles,
3500  const Standard_Integer                CDimension,
3501  const Standard_Integer                CDegree,
3502  const  TColStd_Array1OfReal&          ConstraintPoint, 
3503  const Standard_Integer                Continuity,
3504  const Standard_Boolean                After,
3505  Standard_Integer&                     NbPolesResult,
3506  Standard_Integer&                     NbKnotsResult,
3507  Standard_Real&                        KnotsResult, 
3508  Standard_Real&                        PolesResult) 
3509 {
3510 #ifdef OCCT_DEBUG
3511   if (CDegree<Continuity+1) {
3512     cout<<"The BSpline degree must be greater than the order of continuity"<<endl;
3513   }
3514 #endif
3515   Standard_Real * Padr = &Poles ;
3516   Standard_Real * KRadr = &KnotsResult ;
3517   Standard_Real * PRadr = &PolesResult ;
3518
3519 ////////////////////////////////////////////////////////////////////////
3520 //
3521 //    1. calculation of extension nD
3522 //
3523 ////////////////////////////////////////////////////////////////////////
3524
3525 //  Hermite matrix
3526   Standard_Integer Csize = Continuity + 2;
3527   math_Matrix  MatCoefs(1,Csize, 1,Csize);
3528   if (After) {
3529     PLib::HermiteCoefficients(0, 1,           // Limits 
3530                               Continuity, 0,  // Orders of constraints
3531                               MatCoefs);
3532   }
3533   else {
3534     PLib::HermiteCoefficients(0, 1,           // Limits 
3535                               0, Continuity,  // Orders of constraints
3536                               MatCoefs);    
3537   }
3538
3539
3540 //  position at the node of connection
3541   Standard_Real Tbord ;
3542   if (After) {
3543     Tbord = FlatKnots(FlatKnots.Upper()-CDegree);
3544   }
3545   else {
3546     Tbord = FlatKnots(FlatKnots.Lower()+CDegree);
3547   }
3548   Standard_Boolean periodic_flag = Standard_False ;
3549   Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1);
3550   extrap_mode[0] = extrap_mode[1] = CDegree;
3551   TColStd_Array1OfReal  EvalBS(1, CDimension * (derivative_request+1)) ; 
3552   Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
3553   BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
3554                   CDegree,FlatKnots,CDimension,Poles,*Eadr);
3555
3556 //  norm of the tangent at the node of connection
3557   math_Vector Tgte(1,CDimension);
3558
3559   for (ipos=1;ipos<=CDimension;ipos++) {
3560     Tgte(ipos) = EvalBS(ipos+CDimension);
3561   }
3562   Standard_Real L1=Tgte.Norm();
3563
3564
3565 //  matrix of constraints
3566   math_Matrix Contraintes(1,Csize,1,CDimension);
3567   if (After) {
3568
3569     for (ipos=1;ipos<=CDimension;ipos++) {
3570       Contraintes(1,ipos) = EvalBS(ipos);
3571       Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3572       if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3573       if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3574       Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos);
3575     }
3576   }
3577   else {
3578
3579     for (ipos=1;ipos<=CDimension;ipos++) {
3580       Contraintes(1,ipos) = ConstraintPoint(ipos);
3581       Contraintes(2,ipos) = EvalBS(ipos);
3582       if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3583       if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3584       if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3585     }
3586   }
3587
3588 //  calculate the coefficients of extension
3589   Standard_Integer ii, jj, kk;
3590   TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension);
3591   ExtraCoeffs.Init(0.);
3592
3593   for (ii=1; ii<=Csize; ii++) {
3594
3595     for (jj=1; jj<=Csize; jj++) {
3596
3597       for (kk=1; kk<=CDimension; kk++) {
3598         ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk);
3599       }
3600     }
3601   }
3602
3603 //  calculate the poles of extension
3604   TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension);
3605   Standard_Real * EPadr = &ExtrapPoles(1) ;
3606   PLib::CoefficientsPoles(CDimension,
3607                           ExtraCoeffs, PLib::NoWeights(),
3608                           ExtrapPoles, PLib::NoWeights());
3609
3610 //  calculate the nodes of extension with multiplicities
3611   TColStd_Array1OfReal ExtrapNoeuds(1,2);
3612   ExtrapNoeuds(1) = 0.;
3613   ExtrapNoeuds(2) = 1.;
3614   TColStd_Array1OfInteger ExtrapMults(1,2);
3615   ExtrapMults(1) = Csize;
3616   ExtrapMults(2) = Csize;
3617
3618 // flat nodes of extension
3619   TColStd_Array1OfReal FK2(1, Csize*2);
3620   BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2);
3621
3622 //  norm of the tangent at the connection point 
3623   if (After) {
3624     BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0],
3625                   Csize-1,FK2,CDimension,*EPadr,*Eadr);
3626   }
3627   else {
3628     BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0],
3629                   Csize-1,FK2,CDimension,*EPadr,*Eadr);
3630   }
3631
3632   for (ipos=1;ipos<=CDimension;ipos++) {
3633     Tgte(ipos) = EvalBS(ipos+CDimension);
3634   }
3635   Standard_Real L2 = Tgte.Norm();
3636
3637 //  harmonisation of degrees
3638   TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension);
3639   TColStd_Array1OfReal NewK2(1, 2);
3640   TColStd_Array1OfInteger NewM2(1, 2);
3641   if (Csize-1<CDegree) {
3642     BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension,
3643                              ExtrapPoles,ExtrapNoeuds,ExtrapMults,
3644                              NewP2,NewK2,NewM2);
3645   }
3646   else {
3647     NewP2 = ExtrapPoles;
3648     NewK2 = ExtrapNoeuds;
3649     NewM2 = ExtrapMults;
3650   }
3651
3652 //  flat nodes of extension after harmonization of degrees
3653   TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2);
3654   BSplCLib::KnotSequence(NewK2,NewM2,NewFK2);
3655
3656
3657 ////////////////////////////////////////////////////////////////////////
3658 //
3659 //    2.  concatenation C0
3660 //
3661 ////////////////////////////////////////////////////////////////////////
3662
3663 //  ratio of reparametrization
3664   Standard_Real Ratio=1, Delta;
3665   if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) {
3666     Ratio = L2 / L1;
3667   }
3668   if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1;
3669
3670   if (After) {
3671 //    do not touch the first BSpline
3672     Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper());
3673   }
3674   else {