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