// File:	PLib.cxx
// Created:	Mon Aug 28 16:32:43 1995
// Author:	Laurent BOURESCHE
//		<lbo@phylox>
// Modified: 28/02/1996 by PMN :  HermiteCoefficients added
// Modified: 18/06/1996 by PMN :  NULL reference.
// Modified: 19/02/1997 by JCT :  EvalPoly2Var added

#define No_Standard_RangeError
#define No_Standard_OutOfRange

#include <PLib.ixx>
#include <BSplCLib.hxx>
#include <PLib_LocalArray.hxx>
#include <math_Matrix.hxx> 
#include <math_Gauss.hxx> 
#include <Standard_ConstructionError.hxx>
#include <GeomAbs_Shape.hxx>

// To convert points array into Real ..
// *********************************

#define Dimension_gen 2
#define Array1OfPoints  TColgp_Array1OfPnt2d
#define Point           gp_Pnt2d

#include <PLib_ChangeDim.gxx>

#undef Dimension_gen
#undef Array1OfPoints
#undef Point

#define Dimension_gen 3
#define Array1OfPoints  TColgp_Array1OfPnt
#define Point           gp_Pnt

#include <PLib_ChangeDim.gxx>

#undef Dimension_gen
#undef Array1OfPoints
#undef Point

#include <math_Gauss.hxx>
#include <math.hxx>

//Roman Lygin 19 Sep 2009
//To ensure thread-safety the array of coefficients is initialized in load-time

static Standard_Integer theMaxBinom = -1;
static Standard_Address theBinom;

class BinomAllocator
{
public:
    BinomAllocator()
    {
        PLib::InternalBinomial (BSplCLib::MaxDegree() + 1, theMaxBinom, theBinom);
    }
    virtual ~BinomAllocator()
    {
        if (theBinom)
            delete [] ((Standard_Integer**)theBinom);
        theBinom = NULL;
    }
};
static BinomAllocator theBinomAllocator;


//=======================================================================
//function : Bin
//purpose  : 
//=======================================================================

Standard_Real PLib::Bin(const Standard_Integer N,
                        const Standard_Integer P)
{
  return (Standard_Real)((Standard_Integer**)theBinom)[N][P];
}

//=======================================================================
//function : InternalBinomial
//purpose  : 
//=======================================================================

void PLib::InternalBinomial(const Standard_Integer N,
			    Standard_Integer& maxbinom,
			    Standard_Address& binom)
{
  if (N > maxbinom) {
    Standard_Integer i,im1,ip1,id2,md2,md3,j,k;
    Standard_Integer np1 = N + 1;
    Standard_Integer** nwbin = new Standard_Integer* [np1];
    if (maxbinom >= 0) {

      for (i = 0; i <= maxbinom; i++) nwbin[i] =
	((Standard_Integer**)binom)[i];
      delete [] ((Standard_Integer**)binom);
    }
    else {
      maxbinom = 0;
      nwbin[0] = new Standard_Integer [1];
      nwbin[0][0] = 1;
    }
    binom = nwbin;
    
    for (i = maxbinom + 1; i < np1; i++) {
      im1 = i - 1;
      ip1 = i + 1;
      id2 = i >> 1;
      md2 = im1 >> 1;
      md3 = ip1 >> 1;
      k   = 0;
      ((Standard_Integer**)binom)[i] = new Standard_Integer [ip1];

      for (j = 0; j < id2; j++) {
	((Standard_Integer**)binom)[i][j] =
	    k + ((Standard_Integer**)binom)[im1][j];
	k = ((Standard_Integer**)binom)[im1][j];
      }
      j = id2;
      if (j > md2) j = im1 - j;
      ((Standard_Integer**)binom)[i][id2] =
	  k + ((Standard_Integer**)binom)[im1][j];

      for (j = ip1 - md3; j < ip1; j++) {
	((Standard_Integer**)binom)[i][j] =
	  ((Standard_Integer**)binom)[i][i - j];
      }
    }
    maxbinom = N;
  }
}

//=======================================================================
//function : RationalDerivative
//purpose  : 
//=======================================================================

void  PLib::RationalDerivative(const Standard_Integer Degree,
			       const Standard_Integer DerivativeRequest,
			       const Standard_Integer Dimension,
			       Standard_Real& Ders,
			       Standard_Real& RDers,
			       const Standard_Boolean All)
{
  //
  // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
  // 
  //  We Write  u = fv
  //  Let C(N,P) be the binomial
  //
  //        then we have 
  //
  //         (q)                             (p)   (q-p) 
  //        u    =     SUM          C (q,p) f     v
  //                p = 0 to q
  //
  //
  //        Therefore 
  //        
  //          
  //         (q)         (   (q)                               (p)   (q-p)   )
  //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
  //                     (          p = 0 to q-1                             )
  //
  //
  // make arrays for the binomial since computing it each time could raise a performance
  // issue
  // As oppose to the method below the <Der> array is organized in the following 
  // fashion :
  //
  //  u (1)     u (2) ....   u  (Dimension)     v (1) 
  //
  //   (1)       (1)          (1)                (1) 
  //  u (1)     u (2) ....   u  (Dimension)     v (1)
  //
  //   ............................................
  // 
  //   (Degree)  (Degree)     (Degree)           (Degree) 
  //  u (1)     u (2) ....   u  (Dimension)     v (1) 
  //
  //
  Standard_Real Inverse;
  Standard_Real *PolesArray = &Ders;
  Standard_Real *RationalArray = &RDers;
  Standard_Real Factor ;
  Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
  PLib_LocalArray binomial_array;
  PLib_LocalArray derivative_storage;
  if (Dimension == 3) {
    Standard_Integer DeRequest1 = DerivativeRequest + 1;
    Standard_Integer MinDegRequ = DerivativeRequest;
    if (MinDegRequ > Degree) MinDegRequ = Degree;
    binomial_array.Allocate (DeRequest1);
    
    for (ii = 0 ; ii < DeRequest1 ; ii++) {
      binomial_array[ii] = 1.0e0 ;
    }
    if (!All) {
      Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
      derivative_storage.Allocate (DimDeRequ1);
      RationalArray = derivative_storage ;
    }
    
    Inverse = 1.0e0 / PolesArray[3]  ;
    Index = 0 ;
    Index2 = - 6;
    OtherIndex = 0 ;
    
    for (ii = 0 ; ii <= MinDegRequ ; ii++) {
      Index2 += 3;
      Index1  = Index2;
      RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
      RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
      RationalArray[Index] = PolesArray[OtherIndex];
      Index      -= 2;
      OtherIndex += 2;
      
      for (jj = ii - 1 ; jj >= 0 ; jj--) {
	Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3]; 
	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	RationalArray[Index] -=  Factor * RationalArray[Index1];
	Index  -= 2;
	Index1 -= 5;
      }
      
      for (jj = ii ; jj >=  1 ; jj--) {
	binomial_array[jj] += binomial_array[jj - 1] ;
      }
      RationalArray[Index] *= Inverse ; Index++;
      RationalArray[Index] *= Inverse ; Index++;
      RationalArray[Index] *= Inverse ; Index++;
    } 
    
    for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
      Index2 += 3;
      Index1 = Index2;
      RationalArray[Index] = 0.0e0; Index++;
      RationalArray[Index] = 0.0e0; Index++;
      RationalArray[Index] = 0.0e0;
      Index -= 2;
      
      for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
	Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3]; 
	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	RationalArray[Index] -=  Factor * RationalArray[Index1];
	Index  -= 2;
	Index1 -= 5;
      }
      
      for (jj = ii ; jj >=  1 ; jj--) {
	binomial_array[jj] += binomial_array[jj - 1] ;
      }
      RationalArray[Index] *= Inverse; Index++;
      RationalArray[Index] *= Inverse; Index++;
      RationalArray[Index] *= Inverse; Index++;
    }
    
    if (!All) {
      RationalArray = &RDers ;
      Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
      RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
      RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
      RationalArray[2] = derivative_storage[DimDeRequ];
    }
  }
  else {
    Standard_Integer kk;
    Standard_Integer Dimension1 = Dimension + 1;
    Standard_Integer Dimension2 = Dimension << 1;
    Standard_Integer DeRequest1 = DerivativeRequest + 1;
    Standard_Integer MinDegRequ = DerivativeRequest;
    if (MinDegRequ > Degree) MinDegRequ = Degree;
    binomial_array.Allocate (DeRequest1);
    
    for (ii = 0 ; ii < DeRequest1 ; ii++) {
      binomial_array[ii] = 1.0e0 ;
    }
    if (!All) {
      Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
      derivative_storage.Allocate (DimDeRequ1);
      RationalArray = derivative_storage ;
    }
    
    Inverse = 1.0e0 / PolesArray[Dimension]  ;
    Index = 0 ;
    Index2 = - Dimension2;
    OtherIndex = 0 ;
    
    for (ii = 0 ; ii <= MinDegRequ ; ii++) {
      Index2 += Dimension;
      Index1  = Index2;
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
      }
      Index -= Dimension;
      OtherIndex ++;;
      
      for (jj = ii - 1 ; jj >= 0 ; jj--) {
	Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension]; 
	
	for (kk = 0 ; kk < Dimension ; kk++) {
	  RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	}
	Index  -= Dimension ;
	Index1 -= Dimension2 ;
      }
      
      for (jj = ii ; jj >=  1 ; jj--) {
	binomial_array[jj] += binomial_array[jj - 1] ;
      }
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[Index] *= Inverse ; Index++;
      }
    } 
    
    for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
      Index2 += Dimension;
      Index1  = Index2;
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[Index] = 0.0e0 ; Index++;
      }
      Index -= Dimension;
      
      for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
	Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension]; 
	
	for (kk = 0 ; kk < Dimension ; kk++) {
	  RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	}
	Index  -= Dimension ;
	Index1 -= Dimension2 ;
      }
      
      for (jj = ii ; jj >=  1 ; jj--) {
	binomial_array[jj] += binomial_array[jj - 1] ;
      }
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[Index] *= Inverse; Index++;
      }
    }
    
    if (!All) {
      RationalArray = &RDers ;
      Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
      }
    }
  }
} 

//=======================================================================
//function : RationalDerivatives
//purpose  : Uses Homogeneous Poles Derivatives and Deivatives of Weights
//=======================================================================

void  PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
				const Standard_Integer  Dimension,
				Standard_Real& PolesDerivates, 
				// must be an array with 
				// (DerivativeRequest + 1) * Dimension slots
				Standard_Real& WeightsDerivates,
				// must be an array with 
				// (DerivativeRequest + 1) slots
				Standard_Real& RationalDerivates) 
{
  //
  // Our purpose is to compute f = (u/v) derivated N times
  // 
  //  We Write  u = fv
  //  Let C(N,P) be the binomial
  //
  //        then we have 
  //
  //         (q)                             (p)   (q-p) 
  //        u    =     SUM          C (q,p) f     v
  //                p = 0 to q
  //
  //
  //        Therefore 
  //        
  //          
  //         (q)         (   (q)                               (p)   (q-p)   )
  //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
  //                     (          p = 0 to q-1                             )
  //
  //
  // make arrays for the binomial since computing it each time could 
  // raize a performance issue
  // 
  Standard_Real Inverse;
  Standard_Real *PolesArray = &PolesDerivates;
  Standard_Real *WeightsArray = &WeightsDerivates;
  Standard_Real *RationalArray = &RationalDerivates;
  Standard_Real Factor ;
  
  Standard_Integer ii, Index, Index1, Index2, jj;
  Standard_Integer DeRequest1 = DerivativeRequest + 1;
  
  PLib_LocalArray binomial_array (DeRequest1);
  PLib_LocalArray derivative_storage;

  for (ii = 0 ; ii < DeRequest1 ; ii++) {
    binomial_array[ii] = 1.0e0 ;
  }
  Inverse = 1.0e0 / WeightsArray[0]  ;
  if (Dimension == 3) {
    Index  = 0 ;
    Index2 = - 6 ;

    for (ii = 0 ; ii < DeRequest1 ; ii++) {
      Index2 += 3;
      Index1 = Index2;
      RationalArray[Index] = PolesArray[Index] ; Index++;
      RationalArray[Index] = PolesArray[Index] ; Index++;
      RationalArray[Index] = PolesArray[Index] ;
      Index -= 2;
      
      for (jj = ii - 1 ; jj >= 0 ; jj--) {
	Factor = binomial_array[jj] * WeightsArray[ii - jj] ; 
	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	RationalArray[Index] -=  Factor * RationalArray[Index1];
	Index  -= 2;
	Index1 -= 5;
      }
      
      for (jj = ii ; jj >=  1 ; jj--) {
	binomial_array[jj] += binomial_array[jj - 1] ;
      }
      RationalArray[Index] *= Inverse ; Index++;
      RationalArray[Index] *= Inverse ; Index++;
      RationalArray[Index] *= Inverse ; Index++;
    }
  }
  else {
    Standard_Integer kk;
    Standard_Integer Dimension2 = Dimension << 1;
    Index  = 0 ;
    Index2 = - Dimension2;

    for (ii = 0 ; ii < DeRequest1 ; ii++) {
      Index2 += Dimension;
      Index1  = Index2;
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[Index] = PolesArray[Index]; Index++;
      }
      Index  -= Dimension;
      
      for (jj = ii - 1 ; jj >= 0 ; jj--) {
	Factor = binomial_array[jj] * WeightsArray[ii - jj] ; 
	
	for (kk = 0 ; kk < Dimension ; kk++) {
	  RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
	}
	Index  -= Dimension;
	Index1 -= Dimension2;
      }
      
      for (jj = ii ; jj >=  1 ; jj--) {
	binomial_array[jj] += binomial_array[jj - 1] ;
      }
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	RationalArray[Index] *= Inverse ; Index++;
      }
    }
  }
}

//=======================================================================
//function : This evaluates a polynomial and its derivatives 
//purpose  : up to the requested order 
//=======================================================================

void  PLib::EvalPolynomial(const Standard_Real    Par,
			   const Standard_Integer DerivativeRequest,
			   const Standard_Integer Degree,
			   const Standard_Integer Dimension, 
			   Standard_Real&         PolynomialCoeff,
			   Standard_Real&         Results)
     //
     // the polynomial coefficients are assumed to be stored as follows :
     //                                                      0    
     // [0]                  [Dimension -1]                 X     coefficient
     //                                                      1 
     // [Dimension]          [Dimension + Dimension -1]     X     coefficient
     //                                                      2 
     // [2 * Dimension]      [2 * Dimension + Dimension-1]  X     coefficient
     //
     //   ...................................................
     //
     //  
     //                                                      d 
     // [d * Dimension]      [d * Dimension + Dimension-1]  X     coefficient
     //
     //  where d is the Degree
     //
{
  Standard_Integer DegreeDimension = Degree * Dimension;

  Standard_Integer jj;
  Standard_Real *RA = &Results ;  
  Standard_Real *PA = &PolynomialCoeff ;
  Standard_Real *tmpRA = RA;
  Standard_Real *tmpPA = PA + DegreeDimension;
  
  switch (Dimension) {
    
  case 1 : {
    *tmpRA = *tmpPA; 
    if (DerivativeRequest > 0 ) {
      tmpRA++ ;
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = LocalRequest;
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = Degree;
      }
      MaxIndex2 --;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0;  tmpRA ++ ;          
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA --;
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 --;
	  Index2 --;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA --;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
      }
    }
    break;
  }
  case 2 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++;
    tmpPA --;
    if (DerivativeRequest > 0 ) {
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = LocalRequest << 1;
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = DegreeDimension;
      }
      MaxIndex2 -= 2;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA -= 2;
	
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 2;
	  Index2 -= 2;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 1;
	Index2 = MaxIndex2 + 1;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 2;
	  Index2 -= 2;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
	
	tmpPA --;
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= 2;
	tmpRA  = RA;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
	tmpPA --;
      }
    }
    break;
  }
  case 3 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++;
    tmpPA -= 2;
    if (DerivativeRequest > 0 ) {
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = (LocalRequest << 1) + LocalRequest;
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = DegreeDimension;
      }
      MaxIndex2 -= 3;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA -= 3;
	
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 3;
	  Index2 -= 3;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 1;
	Index2 = MaxIndex2 + 1;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 3;
	  Index2 -= 3;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 2;
	Index2 = MaxIndex2 + 2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 3;
	  Index2 -= 3;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
	
	tmpPA -= 2;
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= 3;
	tmpRA  = RA;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
	tmpPA -= 2;
      }
    }
    break;
  }
  case 6 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++;
    tmpPA -= 5;
    if (DerivativeRequest > 0 ) {
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = (LocalRequest << 2) + (LocalRequest << 1);
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = DegreeDimension;
      }
      MaxIndex2 -= 6;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA -= 6;
	
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 6;
	  Index2 -= 6;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 1;
	Index2 = MaxIndex2 + 1;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 6;
	  Index2 -= 6;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 2;
	Index2 = MaxIndex2 + 2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 6;
	  Index2 -= 6;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 3;
	Index2 = MaxIndex2 + 3;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 6;
	  Index2 -= 6;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 4;
	Index2 = MaxIndex2 + 4;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 6;
	  Index2 -= 6;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 5;
	Index2 = MaxIndex2 + 5;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 6;
	  Index2 -= 6;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
	
	tmpPA -= 5;
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= 6;
	tmpRA  = RA;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
	tmpPA -= 5;
      }
    }
    break;
  }
  case 9 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++;
    tmpPA -= 8;
    if (DerivativeRequest > 0 ) {
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + LocalRequest;
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = DegreeDimension;
      }
      MaxIndex2 -= 9;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA -= 9;
	
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 1;
	Index2 = MaxIndex2 + 1;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 2;
	Index2 = MaxIndex2 + 2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 3;
	Index2 = MaxIndex2 + 3;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 4;
	Index2 = MaxIndex2 + 4;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 5;
	Index2 = MaxIndex2 + 5;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 6;
	Index2 = MaxIndex2 + 6;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 7;
	Index2 = MaxIndex2 + 7;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 8;
	Index2 = MaxIndex2 + 8;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 9;
	  Index2 -= 9;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
	
	tmpPA -= 8;
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= 9;
	tmpRA  = RA;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
	tmpPA -= 8;
      }
    }
    break;
  }
  case 12 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++;
    tmpPA -= 11;
    if (DerivativeRequest > 0 ) {
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + (LocalRequest << 2);
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = DegreeDimension;
      }
      MaxIndex2 -= 12;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA -= 12;
	
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 1;
	Index2 = MaxIndex2 + 1;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 2;
	Index2 = MaxIndex2 + 2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 3;
	Index2 = MaxIndex2 + 3;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 4;
	Index2 = MaxIndex2 + 4;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 5;
	Index2 = MaxIndex2 + 5;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 6;
	Index2 = MaxIndex2 + 6;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 7;
	Index2 = MaxIndex2 + 7;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 8;
	Index2 = MaxIndex2 + 8;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 9;
	Index2 = MaxIndex2 + 9;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 10;
	Index2 = MaxIndex2 + 10;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 11;
	Index2 = MaxIndex2 + 11;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 12;
	  Index2 -= 12;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
	
	tmpPA -= 11;
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= 12;
	tmpRA  = RA;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
	tmpPA -= 11;
      }
    }
    break;
    break;
  }
  case 15 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++;
    tmpPA -= 14;
    if (DerivativeRequest > 0 ) {
      Standard_Real *valRA;
      Standard_Integer ii, LocalRequest;
      Standard_Integer Index1, Index2;
      Standard_Integer MaxIndex1, MaxIndex2;
      if      (DerivativeRequest < Degree) {
	LocalRequest = DerivativeRequest;
	MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest;
      }
      else {
	LocalRequest = Degree;
	MaxIndex2 = MaxIndex1 = DegreeDimension;
      }
      MaxIndex2 -= 15;
      
      for (ii = 1; ii <= LocalRequest; ii++) {
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           

	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
	*tmpRA = 0.0e0; tmpRA++;           
      }
      
      for (jj = Degree ; jj >  0 ; jj--) {
	tmpPA -= 15;
	
	Index1 = MaxIndex1;
	Index2 = MaxIndex2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 1;
	Index2 = MaxIndex2 + 1;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 2;
	Index2 = MaxIndex2 + 2;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 3;
	Index2 = MaxIndex2 + 3;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 4;
	Index2 = MaxIndex2 + 4;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 5;
	Index2 = MaxIndex2 + 5;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 6;
	Index2 = MaxIndex2 + 6;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 7;
	Index2 = MaxIndex2 + 7;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 8;
	Index2 = MaxIndex2 + 8;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 9;
	Index2 = MaxIndex2 + 9;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 10;
	Index2 = MaxIndex2 + 10;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 11;
	Index2 = MaxIndex2 + 11;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 12;
	Index2 = MaxIndex2 + 12;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 13;
	Index2 = MaxIndex2 + 13;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	
	Index1 = MaxIndex1 + 14;
	Index2 = MaxIndex2 + 14;
	
	for (ii = LocalRequest ; ii > 0 ; ii--) {
	  valRA  = &RA[Index1];
	  *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	  Index1 -= 15;
	  Index2 -= 15;
	}
	valRA  = &RA[Index1];
	*valRA = Par * (*valRA) + (*tmpPA);
	
	tmpPA -= 14;
      }
    }
    else {
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= 15;
	tmpRA  = RA;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	*tmpRA = Par * (*tmpRA) + (*tmpPA);
	tmpPA -= 14;
      }
    }
    break;
  }
    default : {
      Standard_Integer kk ;
      for ( kk = 0 ; kk < Dimension ; kk++) {
	*tmpRA = *tmpPA; tmpRA++; tmpPA++;
      }
      tmpPA -= Dimension;
      if (DerivativeRequest > 0 ) {
	Standard_Real *valRA;
	Standard_Integer ii, LocalRequest;
	Standard_Integer Index1, Index2;
	Standard_Integer MaxIndex1, MaxIndex2;
	if      (DerivativeRequest < Degree) {
	  LocalRequest = DerivativeRequest;
	  MaxIndex2 = MaxIndex1 = LocalRequest * Dimension;
	}
	else {
	  LocalRequest = Degree;
	  MaxIndex2 = MaxIndex1 = DegreeDimension;
	}
	MaxIndex2 -= Dimension;
	
	for (ii = 1; ii <= MaxIndex1; ii++) {
	  *tmpRA = 0.0e0; tmpRA++;           
	}
	
	for (jj = Degree ; jj >  0 ; jj--) {
	  tmpPA -= Dimension;
	  
	  for (kk = 0 ; kk < Dimension ; kk++) {
	    Index1 = MaxIndex1 + kk;
	    Index2 = MaxIndex2 + kk;
	    
	    for (ii = LocalRequest ; ii > 0 ; ii--) {
	      valRA  = &RA[Index1];
	      *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
	      Index1 -= Dimension;
	      Index2 -= Dimension;
	    }
	    valRA  = &RA[Index1];
	    *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
	  }
	  tmpPA -= Dimension;
	}
      }
      else {
	
	for (jj = Degree  ; jj >  0 ; jj--) {
	  tmpPA -= Dimension;
	  tmpRA  = RA;
	  
	  for (kk = 0 ; kk < Dimension ; kk++) {
	    *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	  }
	  tmpPA -= Dimension;
	}
      }
    }
  }
}

//=======================================================================
//function : This evaluates a polynomial without derivative 
//purpose  :
//=======================================================================

void  PLib::NoDerivativeEvalPolynomial(const Standard_Real    Par,
				       const Standard_Integer Degree,
				       const Standard_Integer Dimension, 
				       const Standard_Integer DegreeDimension, 
				       Standard_Real&         PolynomialCoeff,
				       Standard_Real&         Results)
{
  Standard_Integer jj;
  Standard_Real *RA = &Results ;  
  Standard_Real *PA = &PolynomialCoeff ;
  Standard_Real *tmpRA = RA;
  Standard_Real *tmpPA = PA + DegreeDimension;

  switch (Dimension) {
  
  case 1 : {
    *tmpRA = *tmpPA;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA--;

      *tmpRA = Par * (*tmpRA) + (*tmpPA);
    }
    break;
  }
  case 2 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA;
    tmpPA--;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA -= 2;
      tmpRA  = RA;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA);
      tmpPA--;
    }
    break;
  }
  case 3 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA;
    tmpPA -= 2;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA -= 3;
      tmpRA  = RA;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA);
      tmpPA -= 2;
    }
    break;
  }
  case 6 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA;
    tmpPA -= 5;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA -= 6;
      tmpRA  = RA;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA);
      tmpPA -= 5;
    }
    break;
  }
  case 9 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA;
    tmpPA -= 8;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA -= 9;
      tmpRA  = RA;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA);
      tmpPA -= 8;
    }
    break;
  }
  case 12 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA;
    tmpPA -= 11;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA -= 12;
      tmpRA  = RA;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA);
      tmpPA -= 11;
    }
    break;
  }
  case 15 : {
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;

    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA; tmpRA++; tmpPA++;
    *tmpRA = *tmpPA;
    tmpPA -= 14;
    
    for (jj = Degree  ; jj >  0 ; jj--) {
      tmpPA -= 15;
      tmpRA  = RA;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;

      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
      *tmpRA = Par * (*tmpRA) + (*tmpPA);
      tmpPA -= 14;
    }
    break;
  }
    default : {
      Standard_Integer kk ;
      for ( kk = 0 ; kk < Dimension ; kk++) {
	*tmpRA = *tmpPA; tmpRA++; tmpPA++;
      }
      tmpPA -= Dimension;
      
      for (jj = Degree  ; jj >  0 ; jj--) {
	tmpPA -= Dimension;
	tmpRA  = RA;
	
	for (kk = 0 ; kk < Dimension ; kk++) {
	  *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
	}
	tmpPA -= Dimension;
      }
    }
  }
}

//=======================================================================
//function : This evaluates a polynomial of 2 variables 
//purpose  : or its derivative at the requested orders 
//=======================================================================

void  PLib::EvalPoly2Var(const Standard_Real    UParameter,
			 const Standard_Real    VParameter,
			 const Standard_Integer UDerivativeRequest,
			 const Standard_Integer VDerivativeRequest,
			 const Standard_Integer UDegree,
			 const Standard_Integer VDegree,
			 const Standard_Integer Dimension, 
			 Standard_Real&         PolynomialCoeff,
			 Standard_Real&         Results)
     //
     // the polynomial coefficients are assumed to be stored as follows :
     //                                                      0 0   
     // [0]                  [Dimension -1]                 U V    coefficient
     //                                                      1 0
     // [Dimension]          [Dimension + Dimension -1]     U V    coefficient
     //                                                      2 0
     // [2 * Dimension]      [2 * Dimension + Dimension-1]  U V    coefficient
     //
     //   ...................................................
     //
     //  
     //                                                      m 0
     // [m * Dimension]      [m * Dimension + Dimension-1]  U V    coefficient
     //
     //  where m = UDegree
     //
     //                                                           0 1
     // [(m+1) * Dimension]  [(m+1) * Dimension + Dimension-1]   U V   coefficient
     //
     //   ...................................................
     //
     //                                                          m 1
     // [2*m * Dimension]      [2*m * Dimension + Dimension-1]  U V    coefficient
     //  
     //   ...................................................
     //
     //                                                          m n
     // [m*n * Dimension]      [m*n * Dimension + Dimension-1]  U V    coefficient
     //
     //  where n = VDegree
{
  Standard_Integer Udim = (VDegree+1)*Dimension,
  index = Udim*UDerivativeRequest;
  TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
  TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1)); 
  Standard_Real * Result =  (Standard_Real *) &Curve.ChangeValue(1);
  Standard_Real * Digit  =  (Standard_Real *) &Point.ChangeValue(1);
  Standard_Real * ResultArray ;
  ResultArray = &Results ;
  
  PLib::EvalPolynomial(UParameter,
		       UDerivativeRequest,
		       UDegree,
		       Udim,
		       PolynomialCoeff,
		       Result[0]);
  
  PLib::EvalPolynomial(VParameter,
		       VDerivativeRequest,
		       VDegree,  
		       Dimension,
		       Result[index],
		       Digit[0]);

  index = Dimension*VDerivativeRequest;

  for (Standard_Integer i=0;i<Dimension;i++) {
    ResultArray[i] = Digit[index+i];
  }
}


static Standard_Integer  storage_divided = 0 ;
static Standard_Real     *divided_differences_array = NULL;

//=======================================================================
//function : This evaluates the lagrange polynomial and its derivatives 
//purpose  : up to the requested order that interpolates a series of
//points of dimension <Dimension> with given assigned parameters
//=======================================================================

Standard_Integer  
PLib::EvalLagrange(const Standard_Real                   Parameter,
		   const Standard_Integer                DerivativeRequest,
		   const Standard_Integer                Degree,
		   const Standard_Integer                Dimension, 
		   Standard_Real&                        Values,
		   Standard_Real&                        Parameters,
		   Standard_Real&                        Results)
{
  //
  // the points  are assumed to be stored as follows in the Values array :
  //                                                      
  // [0]              [Dimension -1]                first  point   coefficients
  //                                                       
  // [Dimension]      [Dimension + Dimension -1]    second point   coefficients
  //                                                     
  // [2 * Dimension]  [2 * Dimension + Dimension-1] third  point   coefficients
  //
  //   ...................................................
  //
  //  
  //                                                      
  // [d * Dimension]  [d * Dimension + Dimension-1] d + 1 point   coefficients
  //
  //  where d is the Degree
  // 
  //  The ParameterArray stores the parameter value assign to each point in
  //  order described above, that is 
  //  [0]   is assign to first  point
  //  [1]   is assign to second point
  //
  Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
  Standard_Integer local_request = DerivativeRequest;
  Standard_Real  *ParameterArray;
  Standard_Real  difference;
  Standard_Real  *PointsArray;
  Standard_Real  *ResultArray ;
  
  PointsArray    = &Values ;
  ParameterArray = &Parameters ;
  ResultArray = &Results ;
  if (local_request >= Degree) {
    local_request = Degree ;
  }
  PLib_LocalArray divided_differences_array ((Degree + 1) * Dimension);
  //
  //  Build the divided differences array
  //
  
  for (ii = 0 ; ii <  (Degree + 1) * Dimension  ; ii++) {
    divided_differences_array[ii] = PointsArray[ii] ;  
  }

  for (ii = Degree ; ii >= 0   ; ii--) {

    for (jj = Degree  ; jj > Degree - ii  ; jj--) {
      Index = jj * Dimension ;
      Index1 = Index - Dimension ;

      for (kk = 0 ; kk < Dimension ; kk++) {
	divided_differences_array[Index + kk] -=
	  divided_differences_array[Index1 + kk] ;
      }
      difference = 
	ParameterArray[jj] - ParameterArray[jj - Degree -1 +  ii] ; 
      if (Abs(difference) < RealSmall()) {
	ReturnCode = 1 ;
	goto FINISH ;
      }
      difference = 1.0e0 / difference ;

      for (kk = 0 ; kk < Dimension ; kk++) {
	divided_differences_array[Index + kk] *= difference ;
      }
    }
  }
  //
  //
  // Evaluate the divided difference array polynomial which expresses as 
  //
  //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
  //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
  //
  // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
  //
  // 
  Index = Degree * Dimension ;

  for (kk = 0 ; kk < Dimension ; kk++) {
    ResultArray[kk] = divided_differences_array[Index + kk] ;
  }  

  for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
    ResultArray[ii] = 0.0e0 ;
  }

  for (ii = Degree ; ii >= 1 ; ii--) {
    difference =  Parameter - ParameterArray[ii - 1] ;

    for (jj = local_request ; jj > 0 ; jj--) {
      Index = jj * Dimension ;
      Index1 = Index - Dimension ;
      
      for (kk = 0 ; kk < Dimension ; kk++) {
	ResultArray[Index + kk] *= difference ;
	ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
      }
    }
    Index = (ii -1) * Dimension ;

    for (kk = 0 ; kk < Dimension ; kk++) {
      ResultArray[kk] *= difference ;
      ResultArray[kk] += divided_differences_array[Index+kk] ;
    }
  }
  FINISH : 
    return (ReturnCode) ;
}

//=======================================================================
//function : This evaluates the hermite polynomial and its derivatives 
//purpose  : up to the requested order that interpolates a series of
//points of dimension <Dimension> with given assigned parameters
//=======================================================================

Standard_Integer PLib::EvalCubicHermite
(const Standard_Real          Parameter,
 const Standard_Integer       DerivativeRequest,
 const Standard_Integer       Dimension, 
 Standard_Real&               Values,
 Standard_Real&               Derivatives,
 Standard_Real&               theParameters,
 Standard_Real&               Results)
{
  //
  // the points  are assumed to be stored as follows in the Values array :
  //                                                      
  // [0]            [Dimension -1]             first  point   coefficients
  //                                                       
  // [Dimension]    [Dimension + Dimension -1] last point   coefficients
  // 
  //
  // the derivatives  are assumed to be stored as follows in 
  // the Derivatives array :
  //                                                      
  // [0]            [Dimension -1]             first  point   coefficients
  //                                                       
  // [Dimension]    [Dimension + Dimension -1] last point   coefficients
  // 
  //  The ParameterArray stores the parameter value assign to each point in
  //  order described above, that is 
  //  [0]   is assign to first  point
  //  [1]   is assign to last   point
  //
  Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
  Standard_Integer local_request = DerivativeRequest ;
  
  ReturnCode = 0 ;
  Degree = 3 ;
  Standard_Real  ParametersArray[4];
  Standard_Real  difference;
  Standard_Real  inverse;
  Standard_Real  *FirstLast;
  Standard_Real  *PointsArray;
  Standard_Real  *DerivativesArray;
  Standard_Real  *ResultArray ;

  DerivativesArray = &Derivatives ;
  PointsArray    = &Values ;
  FirstLast = &theParameters ;
  ResultArray = &Results ;
  if (local_request >= Degree) {
    local_request = Degree ;
  }   
  PLib_LocalArray divided_differences_array ((Degree + 1) * Dimension);

  for (ii = 0, jj = 0  ; ii < 2 ; ii++, jj+= 2) {
    ParametersArray[jj] =
      ParametersArray[jj+1] = FirstLast[ii] ;
  }
  //
  //  Build the divided differences array
  //
  //
  // initialise it at the stage 2 of the building algorithm
  // for devided differences
  //
  inverse = FirstLast[1] - FirstLast[0] ;
  inverse = 1.0e0 / inverse ;

  for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension  ; 
       ii <  Dimension  ; 
       ii++, jj++, kk++, pp++) {
    divided_differences_array[ii] = PointsArray[ii] ;  
    divided_differences_array[kk] = inverse * 
      (PointsArray[jj] - PointsArray[ii]) ; 
    divided_differences_array[jj] = DerivativesArray[ii] ;
    divided_differences_array[pp] = DerivativesArray[jj] ;
  }

  for (ii = 1 ; ii <= Degree   ; ii++) {

    for (jj = Degree  ; jj >=  ii+1  ; jj--) {
      Index = jj * Dimension ;
      Index1 = Index - Dimension ;

      for (kk = 0 ; kk < Dimension ; kk++) {
	divided_differences_array[Index + kk] -=
	  divided_differences_array[Index1 + kk] ;
      }

      for (kk = 0 ; kk < Dimension ; kk++) {
	divided_differences_array[Index + kk] *= inverse ;
      }
    }
  }
  //
  //
  // Evaluate the divided difference array polynomial which expresses as 
  //
  //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
  //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
  //
  // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
  //
  // 
  Index = Degree * Dimension ;

  for (kk = 0 ; kk < Dimension ; kk++) {
    ResultArray[kk] = divided_differences_array[Index + kk] ;
  }  

  for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
    ResultArray[ii] = 0.0e0 ;
  }

  for (ii = Degree ; ii >= 1 ; ii--) {
    difference =  Parameter - ParametersArray[ii - 1] ;

    for (jj = local_request ; jj > 0 ; jj--) {
      Index = jj * Dimension ;
      Index1 = Index - Dimension ;

      for (kk = 0 ; kk < Dimension ; kk++) {
	ResultArray[Index + kk] *= difference ;
	ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
      }
    }
    Index = (ii -1) * Dimension ;

    for (kk = 0 ; kk < Dimension ; kk++) {
      ResultArray[kk] *= difference ;
      ResultArray[kk] += divided_differences_array[Index+kk] ;
    }
  }
//  FINISH : 
    return (ReturnCode) ;
}

//=======================================================================
//function : HermiteCoefficients
//purpose  : calcul des polynomes d'Hermite
//=======================================================================

Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter, 
					   const Standard_Real LastParameter,
					   const Standard_Integer FirstOrder,
					   const Standard_Integer LastOrder, 
					   math_Matrix& MatrixCoefs)
{
  Standard_Integer NbCoeff =  FirstOrder +  LastOrder + 2, Ordre[2];
  Standard_Integer ii, jj, pp, cote, iof=0;
  Standard_Real Prod, TBorne = FirstParameter;
  math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
  math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);

  // Test de validites

  if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
  Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
  if (D1 > 100 || D2 > 100) return Standard_False;
  D2 += D1;
  if (D2 < 0.01) return Standard_False;
  if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False; 

  // Calcul de la matrice a inverser (MAT)

  Ordre[0] = FirstOrder+1;
  Ordre[1] = LastOrder+1;

  for (cote=0; cote<=1; cote++) {
    Coeff.Init(1);

    for (pp=1; pp<=Ordre[cote]; pp++) {
      ii = pp + iof;
      Prod = 1;

      for (jj=pp; jj<=NbCoeff; jj++) {
	//        tout se passe dans les 3 lignes suivantes
	MAT(ii, jj) = Coeff(jj) * Prod;
	Coeff(jj) *= jj - pp;
	Prod      *= TBorne;
      }
    }
    TBorne = LastParameter;
    iof = Ordre[0];
  }

  // resolution du systemes
  math_Gauss ResolCoeff(MAT, 1.0e-10);
  if (!ResolCoeff.IsDone()) return Standard_False;
  
  for (ii=1; ii<=NbCoeff; ii++) {
    B(ii) = 1;
    ResolCoeff.Solve(B, Coeff);
    MatrixCoefs.SetRow( ii, Coeff);
    B(ii) = 0;
  }
  return Standard_True;      
}

//=======================================================================
//function : CoefficientsPoles
//purpose  : 
//=======================================================================

void PLib::CoefficientsPoles (const TColgp_Array1OfPnt&   Coefs,
			      const TColStd_Array1OfReal& WCoefs,
			      TColgp_Array1OfPnt&         Poles, 
			      TColStd_Array1OfReal&       Weights)
{
  TColStd_Array1OfReal tempC(1,3*Coefs.Length());
  PLib::SetPoles(Coefs,tempC);
  TColStd_Array1OfReal tempP(1,3*Poles.Length());
  PLib::SetPoles(Coefs,tempP);
  PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
  PLib::GetPoles(tempP,Poles);
}

//=======================================================================
//function : CoefficientsPoles
//purpose  : 
//=======================================================================

void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
			      const TColStd_Array1OfReal& WCoefs,
			      TColgp_Array1OfPnt2d&       Poles, 
			      TColStd_Array1OfReal&       Weights)
{
  TColStd_Array1OfReal tempC(1,2*Coefs.Length());
  PLib::SetPoles(Coefs,tempC);
  TColStd_Array1OfReal tempP(1,2*Poles.Length());
  PLib::SetPoles(Coefs,tempP);
  PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
  PLib::GetPoles(tempP,Poles);
}

//=======================================================================
//function : CoefficientsPoles
//purpose  : 
//=======================================================================

void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
			      const TColStd_Array1OfReal& WCoefs,
			      TColStd_Array1OfReal&       Poles, 
			      TColStd_Array1OfReal&       Weights)
{
  PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
}

//=======================================================================
//function : CoefficientsPoles
//purpose  : 
//=======================================================================

void PLib::CoefficientsPoles (const Standard_Integer      dim,
			      const TColStd_Array1OfReal& Coefs,
			      const TColStd_Array1OfReal& WCoefs,
			      TColStd_Array1OfReal&       Poles, 
			      TColStd_Array1OfReal&       Weights)
{
  Standard_Boolean rat = &WCoefs != NULL;
  Standard_Integer loc = Coefs.Lower();
  Standard_Integer lop = Poles.Lower();
  Standard_Integer lowc=0;
  Standard_Integer lowp=0;
  Standard_Integer upc = Coefs.Upper();
  Standard_Integer upp = Poles.Upper();
  Standard_Integer upwc=0;
  Standard_Integer upwp=0;
  Standard_Integer reflen = Coefs.Length()/dim;
  Standard_Integer i,j,k; 
  //Les Extremites.
  if (rat) { 
    lowc = WCoefs.Lower(); lowp = Weights.Lower();
    upwc = WCoefs.Upper(); upwp = Weights.Upper();
  }

  for (i = 0; i < dim; i++){
    Poles (lop + i) = Coefs (loc + i);
    Poles (upp - i) = Coefs (upc - i);
  }
  if (rat) {
    Weights (lowp) = WCoefs (lowc);
    Weights (upwp) = WCoefs (upwc);
  }
  
  Standard_Real Cnp;
  PLib::Binomial(reflen - 1);

  for (i = 2; i < reflen; i++ ) {
    Cnp = PLib::Bin(reflen - 1, i - 1);
    if (rat) Weights (lowp + i - 1) = WCoefs (lowc + i - 1) / Cnp;

    for(j = 0; j < dim; j++){
      Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
    }
  }
  
  for (i = 1; i <= reflen - 1; i++) {

    for (j = reflen - 1; j >= i; j--) {
      if (rat) Weights (lowp + j) += Weights (lowp + j -1);

      for(k = 0; k < dim; k++){
	Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
      }
    }
  }
  if (rat) {

    for (i = 1; i <= reflen; i++) {

      for(j = 0; j < dim; j++){
	Poles(lop + dim * (i-1) + j) /= Weights(lowp + i -1);
      }
    }
  }
}

//=======================================================================
//function : Trimming
//purpose  : 
//=======================================================================

void PLib::Trimming(const Standard_Real   U1, 
		    const Standard_Real   U2, 
		    TColgp_Array1OfPnt&   Coefs,
		    TColStd_Array1OfReal& WCoefs)
{
  TColStd_Array1OfReal temp(1,3*Coefs.Length());
  PLib::SetPoles(Coefs,temp);
  PLib::Trimming(U1,U2,3,temp,WCoefs);
  PLib::GetPoles(temp,Coefs);
}

//=======================================================================
//function : Trimming
//purpose  : 
//=======================================================================

void PLib::Trimming(const Standard_Real   U1, 
		    const Standard_Real   U2, 
		    TColgp_Array1OfPnt2d& Coefs,
		    TColStd_Array1OfReal& WCoefs)
{
  TColStd_Array1OfReal temp(1,2*Coefs.Length());
  PLib::SetPoles(Coefs,temp);
  PLib::Trimming(U1,U2,2,temp,WCoefs);
  PLib::GetPoles(temp,Coefs);
}

//=======================================================================
//function : Trimming
//purpose  : 
//=======================================================================

void PLib::Trimming(const Standard_Real   U1, 
		    const Standard_Real   U2, 
		    TColStd_Array1OfReal& Coefs,
		    TColStd_Array1OfReal& WCoefs)
{
  PLib::Trimming(U1,U2,1,Coefs,WCoefs);
}

//=======================================================================
//function : Trimming
//purpose  : 
//=======================================================================

void PLib::Trimming(const Standard_Real U1, 
		    const Standard_Real U2, 
		    const Standard_Integer dim, 
		    TColStd_Array1OfReal& Coefs,
		    TColStd_Array1OfReal& WCoefs)
{

  // principe :
  // on fait le changement de variable v = (u-U1) / (U2-U1)
  // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
  // decomposee sous la forme du schema iteratif de horner.

  Standard_Real lsp = U2 - U1;
  Standard_Integer indc, indw=0;
  Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
  Standard_Integer len = Coefs.Length()/dim;
  Standard_Boolean rat = &WCoefs != NULL;

  if (rat) {
    if(len != WCoefs.Length())
      Standard_Failure::Raise("PLib::Trimming : nbcoefs/dim != nbweights !!!");
    upw = WCoefs.Upper();
  }
  len --;

  for (Standard_Integer i = 1; i <= len; i++) {
    Standard_Integer j ;
    indc = upc - dim*(i-1);
    if (rat) indw = upw - i + 1;
    //calcul du coefficient de degre le plus faible a l'iteration i

    for( j = 0; j < dim; j++){
      Coefs(indc - dim + j) += U1 * Coefs(indc + j);
    }
    if (rat) WCoefs(indw - 1) += U1 * WCoefs(indw);
    
    //calcul des coefficients intermediaires :

    while (indc < upc){
      indc += dim;

      for(Standard_Integer k = 0; k < dim; k++){
	Coefs(indc - dim + k) = 
	  U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
      }
      if (rat) {
	indw ++;
        WCoefs(indw - 1) = U1 * WCoefs(indw) + lsp * WCoefs(indw - 1);
      }
    }

    //calcul du coefficient de degre le plus eleve :

    for(j = 0; j < dim; j++){
      Coefs(upc + j) *= lsp;
    }
    if (rat) WCoefs(upw) *= lsp;
  }
}

//=======================================================================
//function : CoefficientsPoles
//purpose  : 
// Modified: 21/10/1996 by PMN :  PolesCoefficient (PRO5852).
// on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
// fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
//=======================================================================

void PLib::CoefficientsPoles (const TColgp_Array2OfPnt&   Coefs,
			      const TColStd_Array2OfReal& WCoefs,
			      TColgp_Array2OfPnt&         Poles,
			      TColStd_Array2OfReal&       Weights) 
{
  Standard_Boolean rat = (&WCoefs != NULL);
  Standard_Integer LowerRow  = Poles.LowerRow();
  Standard_Integer UpperRow  = Poles.UpperRow();
  Standard_Integer LowerCol  = Poles.LowerCol();
  Standard_Integer UpperCol  = Poles.UpperCol();
  Standard_Integer ColLength = Poles.ColLength();
  Standard_Integer RowLength = Poles.RowLength();

  // Bidouille pour retablir u et v pour les coefs calcules 
  // par buildcache
//  Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();

  Standard_Integer Row, Col;
  Standard_Real W, Cnp;

  Standard_Integer I1, I2;
  Standard_Integer NPoleu , NPolev;
  gp_XYZ Temp;
  PLib::Binomial(RowLength - 1);

  for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
    Poles (NPoleu, LowerCol) =  Coefs (NPoleu, LowerCol);
    if (rat) {
      Weights (NPoleu, LowerCol) =  WCoefs (NPoleu, LowerCol);
    }

    for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
      Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
      Temp = Coefs (NPoleu, Col).XYZ();
      Temp.Divide (Cnp);
      Poles (NPoleu, Col).SetXYZ (Temp);
      if (rat) {
	Weights (NPoleu, Col) = WCoefs (NPoleu, Col) /  Cnp;
      }
    }
    Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
    if (rat) {
      Weights (NPoleu, UpperCol) = WCoefs (NPoleu, UpperCol);
    }

    for (I1 = 1; I1 <= RowLength - 1; I1++) {

      for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
        Temp.SetLinearForm 
	  (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
	Poles (NPoleu, I2).SetXYZ (Temp);
	if (rat) Weights(NPoleu, I2) += Weights(NPoleu, I2-1);
      }
    } 
  }
  PLib::Binomial(ColLength - 1);
  
  for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){

    for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
      Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
      Temp = Poles (Row, NPolev).XYZ();
      Temp.Divide (Cnp);
      Poles (Row, NPolev).SetXYZ (Temp);
      if (rat) Weights(Row, NPolev) /= Cnp;
    }

    for (I1 = 1; I1 <= ColLength - 1; I1++) {

      for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
        Temp.SetLinearForm 
	  (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
	Poles (I2, NPolev).SetXYZ (Temp);
	if (rat) Weights(I2, NPolev) += Weights(I2-1, NPolev);
      }
    }
  }
  if (rat) {

    for (Row = LowerRow; Row <= UpperRow; Row++) {

      for (Col = LowerCol; Col <= UpperCol; Col++) {
	W = Weights (Row, Col);
	Temp = Poles(Row, Col).XYZ();
	Temp.Divide (W);
	Poles(Row, Col).SetXYZ (Temp);
      }
    }
  }
}  

//=======================================================================
//function : UTrimming
//purpose  : 
//=======================================================================

void PLib::UTrimming(const Standard_Real U1, 
		     const Standard_Real U2, 
		     TColgp_Array2OfPnt& Coeffs, 
		     TColStd_Array2OfReal& WCoeffs)
{
  Standard_Boolean rat = &WCoeffs != NULL;
  Standard_Integer lr = Coeffs.LowerRow();
  Standard_Integer ur = Coeffs.UpperRow();
  Standard_Integer lc = Coeffs.LowerCol();
  Standard_Integer uc = Coeffs.UpperCol();
  TColgp_Array1OfPnt  Temp (lr,ur);
  TColStd_Array1OfReal Temw (lr,ur);

  for (Standard_Integer icol = lc; icol <= uc; icol++) {
    Standard_Integer irow ;
    for ( irow = lr; irow <= ur; irow++) {
      Temp (irow) = Coeffs  (irow, icol);
      if (rat) Temw (irow) = WCoeffs (irow, icol);
    }
    if (rat) PLib::Trimming (U1, U2, Temp, Temw);
    else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());

    for (irow = lr; irow <= ur; irow++) {
      Coeffs  (irow, icol) = Temp (irow);
      if (rat) WCoeffs (irow, icol) = Temw (irow);
    }
  }
}

//=======================================================================
//function : VTrimming
//purpose  : 
//=======================================================================

void PLib::VTrimming(const Standard_Real V1, 
		     const Standard_Real V2, 
		     TColgp_Array2OfPnt& Coeffs, 
		     TColStd_Array2OfReal& WCoeffs)
{
  Standard_Boolean rat = &WCoeffs != NULL;
  Standard_Integer lr = Coeffs.LowerRow();
  Standard_Integer ur = Coeffs.UpperRow();
  Standard_Integer lc = Coeffs.LowerCol();
  Standard_Integer uc = Coeffs.UpperCol();
  TColgp_Array1OfPnt  Temp (lc,uc);
  TColStd_Array1OfReal Temw (lc,uc);

  for (Standard_Integer irow = lr; irow <= ur; irow++) {
    Standard_Integer icol ;
    for ( icol = lc; icol <= uc; icol++) {
      Temp (icol) = Coeffs  (irow, icol);
      if (rat) Temw (icol) = WCoeffs (irow, icol);
    }
    if (rat) PLib::Trimming (V1, V2, Temp, Temw);
    else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());

    for (icol = lc; icol <= uc; icol++) {
      Coeffs  (irow, icol) = Temp (icol);
      if (rat) WCoeffs (irow, icol) = Temw (icol);
    }
  }
}

//=======================================================================
//function : HermiteInterpolate
//purpose  : 
//=======================================================================

Standard_Boolean PLib::HermiteInterpolate
(const Standard_Integer Dimension, 
 const Standard_Real FirstParameter, 
 const Standard_Real LastParameter, 
 const Standard_Integer FirstOrder, 
 const Standard_Integer LastOrder, 
 const TColStd_Array2OfReal& FirstConstr,
 const TColStd_Array2OfReal& LastConstr,
 TColStd_Array1OfReal& Coefficients)
{
  Standard_Real Pattern[3][6];

  // portage HP : il faut les initialiser 1 par 1

  Pattern[0][0] = 1;
  Pattern[0][1] = 1;
  Pattern[0][2] = 1;
  Pattern[0][3] = 1;
  Pattern[0][4] = 1;
  Pattern[0][5] = 1;
  Pattern[1][0] = 0;
  Pattern[1][1] = 1;
  Pattern[1][2] = 2;
  Pattern[1][3] = 3;
  Pattern[1][4] = 4;
  Pattern[1][5] = 5;
  Pattern[2][0] = 0;
  Pattern[2][1] = 0;
  Pattern[2][2] = 2;
  Pattern[2][3] = 6;
  Pattern[2][4] = 12;
  Pattern[2][5] = 20;

  math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
  //  The initialisation of the matrix A
  Standard_Integer irow ;
  for ( irow=0; irow<=FirstOrder; irow++) {
    Standard_Real FirstVal = 1.;

    for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
      A(irow,icol) = Pattern[irow][icol]*FirstVal;
      if (irow <= icol) FirstVal *= FirstParameter;
    }
  }

  for (irow=0; irow<=LastOrder; irow++) {
    Standard_Real LastVal  = 1.;

    for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
      A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
      if (irow <= icol) LastVal *= LastParameter;
    }
  }
  //  
  //  The filled matrix A for FirstOrder=LastOrder=2 is:  
  // 
  //      1   FP  FP**2   FP**3    FP**4     FP**5  
  //      0   1   2*FP    3*FP**2  4*FP**3   5*FP**4        FP - FirstParameter
  //      0   0   2       6*FP     12*FP**2  20*FP**3
  //      1   LP  LP**2   LP**3    LP**4     LP**5  
  //      0   1   2*LP    3*LP**2  4*LP**3   5*LP**4        LP - LastParameter 
  //      0   0   2       6*LP     12*LP**2  20*LP**3
  //  
  //  If FirstOrder or LastOrder <=2 then some rows and columns are missing. 
  //  For example:
  //  If FirstOrder=1 then 3th row and 6th column are missing
  //  If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing

  math_Gauss Equations(A);
  //  cout << "A=" << A << endl;

  for (Standard_Integer idim=1; idim<=Dimension; idim++) {
    //  cout << "idim=" << idim << endl;

    math_Vector B(0,FirstOrder+LastOrder+1);
    Standard_Integer icol ;
    for ( icol=0; icol<=FirstOrder; icol++) 
      B(icol) = FirstConstr(idim,icol);

    for (icol=0; icol<=LastOrder; icol++) 
      B(FirstOrder+1+icol) = LastConstr(idim,icol);
    //  cout << "B=" << B << endl;

    //  The solving of equations system A * X = B. Then B = X
    Equations.Solve(B);         
    //  cout << "After Solving" << endl << "B=" << B << endl;

    if (Equations.IsDone()==Standard_False) return Standard_False;
    
    //  the filling of the Coefficients

    for (icol=0; icol<=FirstOrder+LastOrder+1; icol++) 
      Coefficients(Dimension*icol+idim-1) = B(icol);
  } 
  return Standard_True;
}

//=======================================================================
//function : JacobiParameters
//purpose  : 
//=======================================================================

void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder, 
			    const Standard_Integer MaxDegree, 
			    const Standard_Integer Code, 
			    Standard_Integer& NbGaussPoints, 
			    Standard_Integer& WorkDegree)
{
  // ConstraintOrder: Ordre de contrainte aux extremites :
  //            C0 = contraintes de passage aux bornes;
  //            C1 = C0 + contraintes de derivees 1eres;
  //            C2 = C1 + contraintes de derivees 2ndes.
  // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
  //            d' approximation (doit etre superieur ou egal a
  //            2*NivConstr+2 et inferieur ou egal a 50).
  // Code:      Code d' init. des parametres de discretisation.
  //            (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
  //            = -5 Calcul tres rapide mais peu precis (8pts)
  //            = -4    '     '    '      '   '    '    (10pts)
  //            = -3    '     '    '      '   '    '    (15pts)
  //            = -2    '     '    '      '   '    '    (20pts)
  //            = -1    '     '    '      '   '    '    (25pts)
  //            = 1 calcul rapide avec precision moyenne (30pts).
  //            = 2 calcul rapide avec meilleure precision (40pts).
  //            = 3 calcul un peu plus lent avec bonne precision (50 pts).
  //            = 4 calcul lent avec la meilleure precision possible
  //             (61pts).

  // The possible values of NbGaussPoints

  const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25, 
  NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;

  Standard_Integer NivConstr=0;
  switch (ConstraintOrder) {
  case GeomAbs_C0: NivConstr = 0; break;
  case GeomAbs_C1: NivConstr = 1; break;
  case GeomAbs_C2: NivConstr = 2; break;
  default: 
    Standard_ConstructionError::Raise("Invalid ConstraintOrder");
  }
  if (MaxDegree < 2*NivConstr+1)
    Standard_ConstructionError::Raise("Invalid MaxDegree");
  
  if (Code >= 1)
    WorkDegree = MaxDegree + 9;
  else 
    WorkDegree = MaxDegree + 6;
  
  //---> Nbre mini de points necessaires.
  Standard_Integer IPMIN=0;
  if (WorkDegree < NDEG8) 
    IPMIN=NDEG8;
  else if (WorkDegree < NDEG10)
    IPMIN=NDEG10;
  else if (WorkDegree < NDEG15) 
    IPMIN=NDEG15;
  else if (WorkDegree < NDEG20) 
    IPMIN=NDEG20;
  else if (WorkDegree < NDEG25)
    IPMIN=NDEG25;
  else if (WorkDegree < NDEG30) 
    IPMIN=NDEG30;
  else if (WorkDegree < NDEG40) 
    IPMIN=NDEG40;
  else if (WorkDegree < NDEG50) 
    IPMIN=NDEG50;
  else if (WorkDegree < NDEG61) 
    IPMIN=NDEG61;
  else
    Standard_ConstructionError::Raise("Invalid MaxDegree");
  // ---> Nbre de points voulus.
  Standard_Integer IWANT=0;
  switch (Code) {
  case -5: IWANT=NDEG8;  break;
  case -4: IWANT=NDEG10; break;
  case -3: IWANT=NDEG15; break;
  case -2: IWANT=NDEG20; break;
  case -1: IWANT=NDEG25; break;
  case  1: IWANT=NDEG30; break;
  case  2: IWANT=NDEG40; break;
  case  3: IWANT=NDEG50; break;
  case  4: IWANT=NDEG61; break;
  default: 
    Standard_ConstructionError::Raise("Invalid Code");
  }      
  //-->  NbGaussPoints est le nombre de points de discretisation de la fonction,
  //     il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
  //     NbGaussPoints doit etre superieur strictement a WorkDegree.
  NbGaussPoints = Max(IPMIN,IWANT);
  //  NbGaussPoints +=2;
}

//=======================================================================
//function : NivConstr
//purpose  : translates from GeomAbs_Shape to Integer
//=======================================================================

 Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder) 
{
  Standard_Integer NivConstr=0;
  switch (ConstraintOrder) {
    case GeomAbs_C0: NivConstr = 0; break;
    case GeomAbs_C1: NivConstr = 1; break;
    case GeomAbs_C2: NivConstr = 2; break;
    default: 
      Standard_ConstructionError::Raise("Invalid ConstraintOrder");
  }
  return NivConstr;
}

//=======================================================================
//function : ConstraintOrder
//purpose  : translates from Integer to GeomAbs_Shape
//=======================================================================

 GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr) 
{
  GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
  switch (NivConstr) {
    case 0: ConstraintOrder = GeomAbs_C0; break;
    case 1: ConstraintOrder = GeomAbs_C1; break;
    case 2: ConstraintOrder = GeomAbs_C2; break;
    default: 
      Standard_ConstructionError::Raise("Invalid NivConstr");
  }
  return ConstraintOrder;
}

//=======================================================================
//function : EvalLength
//purpose  : 
//=======================================================================

 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension, 
                      Standard_Real& PolynomialCoeff,
		      const Standard_Real U1, const Standard_Real U2,
		      Standard_Real& Length)
{
  Standard_Integer i,j,idim, degdim;
  Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;

  Standard_Real *PolynomialArray = &PolynomialCoeff ;

  Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);

  math_Vector GaussPoints(1,NbGaussPoints);
  math::GaussPoints(NbGaussPoints,GaussPoints);

  math_Vector GaussWeights(1,NbGaussPoints);
  math::GaussWeights(NbGaussPoints,GaussWeights);

  C1 = (U2 + U1) / 2.;
  C2 = (U2 - U1) / 2.;

//-----------------------------------------------------------
//****** Integration - Boucle sur les intervalles de GAUSS **
//-----------------------------------------------------------

  Sum = 0;

  for (j=1; j<=NbGaussPoints/2; j++) {
    // Integration en tenant compte de la symetrie 
    Tran  = C2 * GaussPoints(j);
    X1 = C1 + Tran;
    X2 = C1 - Tran;

    //****** Derivation sur la dimension de l'espace **

    degdim =  Degree*Dimension;
    Der1 = Der2 = 0.;
    for (idim=0; idim<Dimension; idim++) {
      D1 = D2 = Degree * PolynomialArray [idim + degdim];
      for (i=Degree-1; i>=1; i--) {
	DD = i * PolynomialArray [idim + i*Dimension];
	D1 = D1 * X1 + DD;
	D2 = D2 * X2 + DD;
      }
      Der1 += D1 * D1;
      Der2 += D2 * D2;
    }

    //****** Integration **

    Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));

//****** Fin de boucle dur les intervalles de GAUSS **
  }
  Length = Sum;
}


//=======================================================================
//function : EvalLength
//purpose  : 
//=======================================================================

 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension, 
                      Standard_Real& PolynomialCoeff,
		      const Standard_Real U1, const Standard_Real U2,
		      const Standard_Real Tol, 
		      Standard_Real& Length, Standard_Real& Error)
{
  Standard_Integer i;
  Standard_Integer NbSubInt = 1,    // Current number of subintervals
                   MaxNbIter = 13,  // Max number of iterations
                   NbIter    = 1;   // Current number of iterations
  Standard_Real    dU,OldLen,LenI;

  PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);

  do {
    OldLen = Length;
    Length = 0.;
    NbSubInt *= 2;
    dU = (U2-U1)/NbSubInt;    
    for (i=1; i<=NbSubInt; i++) {
      PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
      Length += LenI;
    }
    NbIter++;
    Error = Abs(OldLen-Length);
  }    
  while (Error > Tol && NbIter <= MaxNbIter);
}  
