0024285: Updates of PLib::EvalPolynomial for code acceleration

Functions PLib::EvalPolynomial and PLib::NoDerivativeEvalPolynomial are refactored to allow generation of faster code:
1. Iteration by degree is made in outer loop
2. Avoided pointer arithmetic
3. Recursive templates are used to expand loop by dimension in specific cases (1-15)

diff --git a/src/PLib/PLib.cxx b/src/PLib/PLib.cxx
index 2de1ca6..99330d0 100644 (file)
--- a/src/PLib/PLib.cxx
+++ b/src/PLib/PLib.cxx
@@ -636,17 +636,111 @@ void  PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
   }
 }
 
+//=======================================================================
+// Auxiliary template functions used for optimized evaluation of polynome
+// and its derivatives for smaller dimensions of the polynome
+//=======================================================================
+
+namespace {
+  // recursive template for evaluating value or first derivative
+  template<int dim> 
+  inline void eval_step1 (double* poly, double par, double* coef)
+  {
+    eval_step1<dim - 1> (poly, par, coef);
+    poly[dim] = poly[dim] * par + coef[dim];
+  }
+
+  // recursion end
+  template<>
+  inline void eval_step1<-1> (double*, double, double*)
+  {
+  }
+
+  // recursive template for evaluating second derivative
+  template<int dim> 
+  inline void eval_step2 (double* poly, double par, double* coef)
+  {
+    eval_step2<dim - 1> (poly, par, coef);
+    poly[dim] = poly[dim] * par + coef[dim] * 2.;
+  }
+
+  // recursion end
+  template<>
+  inline void eval_step2<-1> (double*, double, double*)
+  {
+  }
+
+  // evaluation of only value
+  template<int dim>
+  inline void eval_poly0 (double* aRes, double* aCoeffs, int Degree, double Par)
+  {
+    Standard_Real* aRes0 = aRes;
+    memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
+
+    for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
+    {
+      aCoeffs -= dim;
+      // Calculating the value of the polynomial
+      eval_step1<dim-1> (aRes0, Par, aCoeffs);
+    }
+  }
+
+  // evaluation of value and first derivative
+  template<int dim>
+  inline void eval_poly1 (double* aRes, double* aCoeffs, int Degree, double Par)
+  {
+    Standard_Real* aRes0 = aRes;
+    Standard_Real* aRes1 = aRes + dim;
+
+    memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
+    memset(aRes1, 0, sizeof(Standard_Real) * dim);
+
+    for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
+    {
+      aCoeffs -= dim;
+      // Calculating derivatives of the polynomial
+      eval_step1<dim-1> (aRes1, Par, aRes0);
+      // Calculating the value of the polynomial
+      eval_step1<dim-1> (aRes0, Par, aCoeffs);
+    }
+  }
+
+  // evaluation of value and first and second derivatives
+  template<int dim>
+  inline void eval_poly2 (double* aRes, double* aCoeffs, int Degree, double Par)
+  {
+    Standard_Real* aRes0 = aRes;
+    Standard_Real* aRes1 = aRes + dim;
+    Standard_Real* aRes2 = aRes + 2 * dim;
+
+    memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
+    memset(aRes1, 0, sizeof(Standard_Real) * dim);
+    memset(aRes2, 0, sizeof(Standard_Real) * dim);
+
+    for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
+    {
+      aCoeffs -= dim;
+      // Calculating second derivatives of the polynomial
+      eval_step2<dim-1> (aRes2, Par, aRes1);
+      // Calculating derivatives of the polynomial
+      eval_step1<dim-1> (aRes1, Par, aRes0);
+      // Calculating the value of the polynomial
+      eval_step1<dim-1> (aRes0, Par, aCoeffs);
+    }
+  }
+}
+
 //=======================================================================
 //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)
+                           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    
@@ -665,1072 +759,121 @@ void  PLib::EvalPolynomial(const Standard_Real    Par,
      //  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;
+  Standard_Real* aCoeffs = &PolynomialCoeff + Degree * Dimension;
+  Standard_Real* aRes = &Results;
+  Standard_Real* anOriginal;
+  Standard_Integer ind = 0;
+  switch (DerivativeRequest)
+  {
+  case 1:
+  {
+    switch (Dimension)
+    {
+    case 1:  eval_poly1<1>  (aRes, aCoeffs, Degree, Par); break;
+    case 2:  eval_poly1<2>  (aRes, aCoeffs, Degree, Par); break;
+    case 3:  eval_poly1<3>  (aRes, aCoeffs, Degree, Par); break;
+    case 4:  eval_poly1<4>  (aRes, aCoeffs, Degree, Par); break;
+    case 5:  eval_poly1<5>  (aRes, aCoeffs, Degree, Par); break;
+    case 6:  eval_poly1<6>  (aRes, aCoeffs, Degree, Par); break;
+    case 7:  eval_poly1<7>  (aRes, aCoeffs, Degree, Par); break;
+    case 8:  eval_poly1<8>  (aRes, aCoeffs, Degree, Par); break;
+    case 9:  eval_poly1<9>  (aRes, aCoeffs, Degree, Par); break;
+    case 10: eval_poly1<10> (aRes, aCoeffs, Degree, Par); break;
+    case 11: eval_poly1<11> (aRes, aCoeffs, Degree, Par); break;
+    case 12: eval_poly1<12> (aRes, aCoeffs, Degree, Par); break;
+    case 13: eval_poly1<13> (aRes, aCoeffs, Degree, Par); break;
+    case 14: eval_poly1<14> (aRes, aCoeffs, Degree, Par); break;
+    case 15: eval_poly1<15> (aRes, aCoeffs, Degree, Par); break;
+    default:
+      {
+        Standard_Real* aRes0 = aRes;
+        Standard_Real* aRes1 = aRes + Dimension;
+
+        memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * Dimension);
+        memset(aRes1, 0, sizeof(Standard_Real) * Dimension);
+
+        for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
+        {
+          aCoeffs -= Dimension;
+          // Calculating derivatives of the polynomial
+          for (ind = 0; ind < Dimension; ind++)
+            aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
+          // Calculating the value of the polynomial
+          for (ind = 0; ind < Dimension; ind++)
+            aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
+        }
       }
     }
     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;
+  case 2:
+  {
+    switch (Dimension)
+    {
+    case 1:  eval_poly2<1>  (aRes, aCoeffs, Degree, Par); break;
+    case 2:  eval_poly2<2>  (aRes, aCoeffs, Degree, Par); break;
+    case 3:  eval_poly2<3>  (aRes, aCoeffs, Degree, Par); break;
+    case 4:  eval_poly2<4>  (aRes, aCoeffs, Degree, Par); break;
+    case 5:  eval_poly2<5>  (aRes, aCoeffs, Degree, Par); break;
+    case 6:  eval_poly2<6>  (aRes, aCoeffs, Degree, Par); break;
+    case 7:  eval_poly2<7>  (aRes, aCoeffs, Degree, Par); break;
+    case 8:  eval_poly2<8>  (aRes, aCoeffs, Degree, Par); break;
+    case 9:  eval_poly2<9>  (aRes, aCoeffs, Degree, Par); break;
+    case 10: eval_poly2<10> (aRes, aCoeffs, Degree, Par); break;
+    case 11: eval_poly2<11> (aRes, aCoeffs, Degree, Par); break;
+    case 12: eval_poly2<12> (aRes, aCoeffs, Degree, Par); break;
+    case 13: eval_poly2<13> (aRes, aCoeffs, Degree, Par); break;
+    case 14: eval_poly2<14> (aRes, aCoeffs, Degree, Par); break;
+    case 15: eval_poly2<15> (aRes, aCoeffs, Degree, Par); break;
+    default: 
+      {
+        Standard_Real* aRes0 = aRes;
+        Standard_Real* aRes1 = aRes + Dimension;
+        Standard_Real* aRes2 = aRes1 + Dimension;
+
+        // Nullify the results
+        Standard_Integer aSize = 2 * Dimension;
+        memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
+        memset(aRes1, 0, sizeof(Standard_Real) * aSize);
+
+        for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
+        {
+          aCoeffs -= Dimension;
+          // Calculating derivatives of the polynomial
+          for (ind = 0; ind < Dimension; ind++)
+            aRes2[ind] = aRes2[ind] * Par + aRes1[ind] * 2.0;
+          for (ind = 0; ind < Dimension; ind++)
+            aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
+          // Calculating the value of the polynomial
+          for (ind = 0; ind < Dimension; ind++)
+            aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
+        }
+        break;
       }
     }
     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;
-       }
+  default:
+    {
+      // Nullify the results
+      Standard_Integer aResSize = (1 + DerivativeRequest) * Dimension;
+      memset(aRes, 0, sizeof(Standard_Real) * aResSize);
+
+      for (Standard_Integer aDeg = 0; aDeg <= Degree; aDeg++)
+      {
+        aRes = &Results + aResSize - Dimension;
+        // Calculating derivatives of the polynomial
+        for (Standard_Integer aDeriv = DerivativeRequest; aDeriv > 0; aDeriv--)
+        {
+          anOriginal = aRes - Dimension; // pointer to the derivative minus 1
+          for (ind = 0; ind < Dimension; ind++)
+            aRes[ind] = aRes[ind] * Par + anOriginal[ind] * aDeriv;
+          aRes = anOriginal;
+        }
+        // Calculating the value of the polynomial
+        for (ind = 0; ind < Dimension; ind++)
+          aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
+        aCoeffs -= Dimension;
       }
     }
   }
@@ -1742,225 +885,40 @@ void  PLib::EvalPolynomial(const Standard_Real    Par,
 //=======================================================================
 
 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)
+                                       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--;
+  Standard_Real* aCoeffs = &PolynomialCoeff + DegreeDimension;
+  Standard_Real* aRes = &Results;
 
-      *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;
+  switch (Dimension)
+  {
+  case 1:  eval_poly0<1>  (aRes, aCoeffs, Degree, Par); break;
+  case 2:  eval_poly0<2>  (aRes, aCoeffs, Degree, Par); break;
+  case 3:  eval_poly0<3>  (aRes, aCoeffs, Degree, Par); break;
+  case 4:  eval_poly0<4>  (aRes, aCoeffs, Degree, Par); break;
+  case 5:  eval_poly0<5>  (aRes, aCoeffs, Degree, Par); break;
+  case 6:  eval_poly0<6>  (aRes, aCoeffs, Degree, Par); break;
+  case 7:  eval_poly0<7>  (aRes, aCoeffs, Degree, Par); break;
+  case 8:  eval_poly0<8>  (aRes, aCoeffs, Degree, Par); break;
+  case 9:  eval_poly0<9>  (aRes, aCoeffs, Degree, Par); break;
+  case 10: eval_poly0<10> (aRes, aCoeffs, Degree, Par); break;
+  case 11: eval_poly0<11> (aRes, aCoeffs, Degree, Par); break;
+  case 12: eval_poly0<12> (aRes, aCoeffs, Degree, Par); break;
+  case 13: eval_poly0<13> (aRes, aCoeffs, Degree, Par); break;
+  case 14: eval_poly0<14> (aRes, aCoeffs, Degree, Par); break;
+  case 15: eval_poly0<15> (aRes, aCoeffs, Degree, Par); break;
+  default:
+    {
+      memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
+      for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
+      {
+        aCoeffs -= Dimension;
+        for (Standard_Integer ind = 0; ind < Dimension; ind++)
+          aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
       }
     }
   }

diff --git a/tests/perf/bspline/intersect b/tests/perf/bspline/intersect
new file mode 100644 (file)
index 0000000..513fb16
--- /dev/null
+++ b/tests/perf/bspline/intersect
@@ -0,0 +1,213 @@
+# Test performance of intersection of several NURBS surfaces
+
+if { [regexp {Debug mode} [dversion]] } {
+   cpulimit 1000
+} else {
+   cpulimit 500
+}
+
+bsplinesurf surf1 \
+3 4 0 4 1 1 2 1 3 4 \
+3 4 0 4 1 1 2 1 3 4 \
+0  0  0 1   2  0  0 1   3  0 15 1   5  0 15 1   7  0  0 1   10  0  0 1 \
+0  2  0 1   1  3  0 1   4  2 15 1   6  3 15 1   8  2  0 1   10  3  0 1 \
+0  4  0 1   3  4  0 1   4  3 15 1   5  3 15 1   7  4  0 1   10  5  0 1 \
+0  6  0 1   3  6  0 1   4  6 15 1   5  6 15 1   8  5  0 1   10  7  0 1 \
+0  8  0 1   2  8  0 1   4  8 15 1   6  8 15 1   7  7  0 1   10  8  0 1 \
+0 10  0 1   2 10  0 1   4 10 15 1   6 10 15 1   7 10  0 1   10 10  0 1
+
+bsplinesurf surf2 \
+3 4 0 4 1 1 2 1 3 4 \
+3 4 0 4 1 1 2 1 3 4 \
+0  0 10  1   2  0  5  1   3  0  4  1   5  0  6  1   7  0 10  1   10  0  5  1 \
+0  2  5  1   1  3  7  1   4  2  7  1   6  3  4  1   8  2  4  4   10  3  7  1 \
+0  4  8  1   3  4 10  1   4  3  6  1   5  3  8  1   7  4  7  4   10  5  5  1 \
+0  6  8  1   3  6 10  1   4  6  6  1   5  6  4  1   8  5  7  4   10  7 10  1 \
+0  8  6  1   2  8  5  1   4  8  8  1   6  8  8  1   7  7  3  4   10  8  5  1 \
+0 10  8  1   2 10 10  1   4 10  6  1   6 10  5  1   7 10  3  1   10 10 10  1
+
+bsplinesurf surf3 \
+2 7 0 3 1 1 2 1 3 1 4 1 5 1 6 3 \
+2 3 0 3 1 1 2 3 \
+-20  20  10  1   20  20  10  1    20 -20  10  1   -20 -20  10  1   -10  10  10  1   10  10  10  1   10 -10  10  1   -10 -10  10  1 \
+-10  10   5  1   15  15   5  1    10 -10   5  1   -15 -15   5  1    -5   5   5  1    5   5   5  1    5  -5   5  1    -5  -5   5  1 \
+-15  15  -5  1   10  10  -5  1    15 -15  -5  1   -10 -10  -5  1    -5   5  -5  1    5   5  -5  1    5  -5  -5  1    -5  -5  -5  1 \
+-20  20 -10  1   20  20 -10  1    20 -20 -10  1   -20 -20 -10  1   -10  10 -10  1   10  10 -10  1   10 -10 -10  1   -10 -10 -10  1
+
+bsplinesurf surf4 \
+4 5 0 5 2 1 4 1 5 1 6 5 \
+4 5 0 5 1 1 4 1 6 1 7 5 \
+ 10  20 -20 1     9  19 -19 1     8  18 -18 1     7  17 -17 1     7  17  17 1     8  18  18 1     9  19  19 1    10  20  20 1 \
+ 10 -20  20 1     9 -19  19 1     8 -18  18 1     7 -17  17 1     7 -17 -17 1     8 -18 -18 1     9 -19 -19 1    10 -20 -20 1 \
+ 10  10 -10 1     9   9  -9 1     8   8  -8 1     7   7  -7 1     7   7   7 1     8   8   8 1     9   9   9 1    10  10  10 1 \
+ 10 -10  10 1     9  -9   9 1     8  -8   8 1     7  -7   7 1     7  -7  -7 1     8  -8  -8 1     9  -9  -9 1    10 -10 -10 1 \
+-10  20 -20 1    -9  19 -19 1    -8  18 -18 1    -7  17 -17 1    -7  17  17 1    -8  18  18 1    -9  19  19 1   -10  20  20 1 \
+-10 -20  20 1    -9 -19  19 1    -8 -18  18 1    -7 -17  17 1    -7 -17 -17 1    -8 -18 -18 1    -9 -19 -19 1   -10 -20 -20 1 \
+-10  10 -10 1    -9   9  -9 1    -8   8  -8 1    -7   7  -7 1    -7   7   7 1    -8   8   8 1    -9   9   9 1   -10  10  10 1 \
+-10 -10  10 1    -9  -9   9 1    -8  -8   8 1    -7  -7   7 1    -7  -7  -7 1    -8  -8  -8 1    -9  -9  -9 1   -10 -10 -10 1
+
+bsplinesurf surf5 \
+1 2 0 2 1 2 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0 0 0 1 10 0 5 2 \
+0 10  2  1   10 10  3  1 \
+0 20 10  1   10 20 20  1 \
+0 30  0  1   10 30  0  1 \
+0 40 -1  1   10 40  5  1 \
+0 50  5  1   10 50  5  1 \
+0 60  4  1   10 60  4  1 \
+0 70 -5  1   10 70 -3  1 \
+0 80  7  1   10 80  0  1
+
+
+bsplinesurf surf6 \
+2 3 0 3 1 1 2 3 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0 0 0 1 2 0 0 1 5 0 -1 1 10 0 5 2 \
+0 10  2  1   3 10  0  1   8 10  5  1   10 10  3  1 \
+0 20 10  1   4 20  4  1   7 20  4  1   10 20 20  1 \
+0 30  0  1   2 30  0  1   8 30  0  1   10 30  0  1 \
+0 40 -1  1   4 40  5  1   9 40  1  1   10 40  5  1 \
+0 50  5  1   4 50 10  1   6 50 10  1   10 50  5  1 \
+0 60  4  1   3 60 -3  1   7 60 -4  1   10 60  4  1 \
+0 70 -5  1   3 70  0  1   5 70  0  1   10 70 -3  1 \
+0 80  7  1   3 80  1  1   7 80  3  1   10 80  0  1
+
+bsplinesurf surf7 \
+5 5 0 6 1 1 4 1 5 1 8 6 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0  0  0  1   2  0  0  1   5  0 -1  1   10  0  5  1   12  0  1  1   15  0 -3  1   16  0 -3  1   19  0 -4 1   24  0  0  1 \
+0 10  2  1   3 10  0  1   8 10  5  1   10 10  3  1   12 10  2  1   15 10  0  1   20 10  5  1   21 10  3 1   24 10  0  1 \
+0 20 10  1   4 20  4  1   7 20  4  1   10 20 20  1   12 20 10  1   16 20  4  1   19 20  4  1   20 20 10 1   24 20  0  1 \
+0 30  0  1   2 30  0  1   8 30  0  1   10 30  0  1   12 30  0  1   14 30  0  1   20 30  0  1   22 30  0 1   24 30  0  1 \
+0 40 -1  1   4 40  5  1   9 40  1  1   10 40  5  1   12 40 -1  1   16 40  5  1   21 40  1  1   22 40  5 1   24 40  0  1 \
+0 50  5  1   4 50 10  1   6 50 10  1   10 50  5  1   12 50  5  1   16 50 10  1   18 50 10  1   20 50  5 1   24 50  0  1 \
+0 60  4  1   3 60 -3  1   7 60 -4  1   10 60  4  1   12 60  4  1   15 60 -3  1   19 60 -4  1   20 60  4 1   24 60  0  1 \
+0 70 -5  1   3 70  0  1   5 70  0  1   10 70 -3  1   12 70 -5  1   15 70  0  1   17 70  0  1   20 70 -3 1   24 70  0  1 \
+0 80  7  1   3 80  1  1   7 80  3  1   10 80  0  1   12 80  7  1   15 80  1  1   19 80  3  1   21 80  0 1   24 80  0  1
+
+bsplinesurf surf8 \
+5 5 0 6 1 1 4 1 5 1 8 6 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0  0 -8  1   2  0 -7  1   5  0 -9  1   10  0 -8  1   12  0 -9  1   15  0 -7  1   16  0 -8  1   19  0 -8 1   24  0 -8  1 \
+0 10 -8  1   3 10 -7  1   8 10 -9  1   10 10 -8  1   12 10 -9  1   15 10 -7  1   20 10 -8  1   21 10 -8 1   24 10 -8  1 \
+0 20  8  1   4 20  7  1   7 20  9  1   10 20  8  1   12 20  9  1   16 20  7  1   19 20  8  1   20 20  8 1   24 20  8  1 \
+0 30  8  1   2 30  7  1   8 30  9  1   10 30  8  1   12 30  9  1   14 30  7  1   20 30  8  1   22 30  8 1   24 30  8  1 \
+0 40 -8  1   4 40 -7  1   9 40 -9  1   10 40 -8  1   12 40 -9  1   16 40 -7  1   21 40 -8  1   22 40 -8 1   24 40 -8  1 \
+0 50 -8  1   4 50 -7  1   6 50 -9  1   10 50 -8  1   12 50 -9  1   16 50 -7  1   18 50 -8  1   20 50 -8 1   24 50 -8  1 \
+0 60  8  1   3 60  7  1   7 60  9  1   10 60  8  1   12 60  9  1   15 60  7  1   19 60  8  1   20 60  8 1   24 60  8  1 \
+0 70  8  1   3 70  7  1   5 70  9  1   10 70  8  1   12 70  9  1   15 70  7  1   17 70  8  1   20 70  8 1   24 70  8  1 \
+0 80  8  1   3 80  7  1   7 80  9  1   10 80  8  1   12 80  9  1   15 80  7  1   19 80  8  1   21 80  8 1   24 80  8  1
+
+bsplinesurf rsurf1 \
+3 4 0 4 1 1 2 1 3 4 \
+3 4 0 4 1 1 2 1 3 4 \
+0  0  0 1   2  0  0 1   3  0 15 1   5  0 15 1   7  0  0 1   10  0  0 1 \
+0  2  0 1   1  3  0 2   4  2 15 2   6  3 15 2   8  2  0 2   10  3  0 1 \
+0  4  0 1   3  4  0 2   4  3 15 3   5  3 15 3   7  4  0 2   10  5  0 1 \
+0  6  0 1   3  6  0 2   4  6 15 3   5  6 15 3   8  5  0 2   10  7  0 1 \
+0  8  0 1   2  8  0 2   4  8 15 2   6  8 15 2   7  7  0 2   10  8  0 1 \
+0 10  0 1   2 10  0 1   4 10 15 1   6 10 15 1   7 10  0 1   10 10  0 1
+
+bsplinesurf rsurf2 \
+3 4 0 4 1 1 2 1 3 4 \
+3 4 0 4 1 1 2 1 3 4 \
+0  0 10  1   2  0  5  1   3  0  4  1   5  0  6  1   7  0 10  1   10  0  5  1 \
+0  2  5  1   1  3  7  4   4  2  7  4   6  3  4  4   8  2  4  4   10  3  7  1 \
+0  4  8  1   3  4 10  4   4  3  6  2   5  3  8  2   7  4  7  4   10  5  5  1 \
+0  6  8  1   3  6 10  4   4  6  6  2   5  6  4  2   8  5  7  4   10  7 10  1 \
+0  8  6  1   2  8  5  4   4  8  8  4   6  8  8  4   7  7  3  4   10  8  5  1 \
+0 10  8  1   2 10 10  1   4 10  6  1   6 10  5  1   7 10  3  1   10 10 10  1
+
+bsplinesurf rsurf3 \
+2 7 0 3 1 1 2 1 3 1 4 1 5 1 6 3 \
+2 3 0 3 1 1 2 3 \
+-20  20  10  1   20  20  10  1    20 -20  10  2   -20 -20  10 2   -10  10  10  3   10  10  10  4   10 -10  10  4   -10 -10  10  2 \
+-10  10   5  1   15  15   5  2    10 -10   5  3   -15 -15   5 1    -5   5   5  1    5   5   5  3    5  -5   5  4    -5  -5   5  3 \
+-15  15  -5  3   10  10  -5  1    15 -15  -5  2   -10 -10  -5 2    -5   5  -5  3    5   5  -5  3    5  -5  -5  1    -5  -5  -5  4 \
+-20  20 -10  2   20  20 -10  4    20 -20 -10  3   -20 -20 -10 3   -10  10 -10  2   10  10 -10  2   10 -10 -10  1   -10 -10 -10  1
+
+bsplinesurf rsurf4 \
+4 5 0 5 2 1 4 1 5 1 6 5 \
+4 5 0 5 1 1 4 1 6 1 7 5 \
+ 10  20 -20 1     9  19 -19 1     8  18 -18 2     7  17 -17 2     7  17  17 1     8  18  18 1     9  19  19 4    10  20  20 1 \
+ 10 -20  20 2     9 -19  19 2     8 -18  18 1     7 -17  17 2     7 -17 -17 4     8 -18 -18 3     9 -19 -19 3    10 -20 -20 2 \
+ 10  10 -10 3     9   9  -9 4     8   8  -8 1     7   7  -7 2     7   7   7 4     8   8   8 1     9   9   9 2    10  10  10 3 \
+ 10 -10  10 4     9  -9   9 2     8  -8   8 3     7  -7   7 2     7  -7  -7 1     8  -8  -8 1     9  -9  -9 1    10 -10 -10 4 \
+-10  20 -20 4    -9  19 -19 3    -8  18 -18 3    -7  17 -17 2    -7  17  17 1    -8  18  18 2    -9  19  19 4   -10  20  20 4 \
+-10 -20  20 3    -9 -19  19 2    -8 -18  18 3    -7 -17  17 1    -7 -17 -17 4    -8 -18 -18 2    -9 -19 -19 3   -10 -20 -20 3 \
+-10  10 -10 2    -9   9  -9 1    -8   8  -8 3    -7   7  -7 1    -7   7   7 1    -8   8   8 1    -9   9   9 2   -10  10  10 2 \
+-10 -10  10 1    -9  -9   9 1    -8  -8   8 2    -7  -7   7 1    -7  -7  -7 4    -8  -8  -8 1    -9  -9  -9 1   -10 -10 -10 1
+
+bsplinesurf rsurf5 \
+1 2 0 2 1 2 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0 0 0 1 10 0 5 2 \
+0 10  2  3   10 10  3  1 \
+0 20 10  1   10 20 20  4 \
+0 30  0  2   10 30  0  2 \
+0 40 -1  1   10 40  5  2 \
+0 50  5  1   10 50  5  1 \
+0 60  4  3   10 60  4  1 \
+0 70 -5  5   10 70 -3  2 \
+0 80  7  1   10 80  0  1
+
+bsplinesurf rsurf6 \
+2 3 0 3 1 1 2 3 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0 0 0 1 2 0 0 1 5 0 -1 1 10 0 5 2 \
+0 10  2  3   3 10  0  1   8 10  5  3   10 10  3  1 \
+0 20 10  1   4 20  4  3   7 20  4  2   10 20 20  4 \
+0 30  0  2   2 30  0  1   8 30  0  1   10 30  0  2 \
+0 40 -1  1   4 40  5  2   9 40  1  3   10 40  5  2 \
+0 50  5  1   4 50 10  5   6 50 10  5   10 50  5  1 \
+0 60  4  3   3 60 -3  1   7 60 -4  1   10 60  4  1 \
+0 70 -5  5   3 70  0  1   5 70  0  1   10 70 -3  2 \
+0 80  7  1   3 80  1  2   7 80  3  1   10 80  0  1
+
+bsplinesurf rsurf7 \
+5 5 0 6 1 1 4 1 5 1 8 6 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0  0  0  1   2  0  0  1   5  0 -1  1   10  0  5  2   12  0  1  1   15  0 -3  2   16  0 -3  1   19  0 -4 1   24  0  0  1 \
+0 10  2  3   3 10  0  1   8 10  5  3   10 10  3  1   12 10  2  3   15 10  0  1   20 10  5  3   21 10  3 1   24 10  0  2 \
+0 20 10  1   4 20  4  3   7 20  4  2   10 20 20  4   12 20 10  1   16 20  4  3   19 20  4  2   20 20 10 4   24 20  0  1 \
+0 30  0  2   2 30  0  1   8 30  0  1   10 30  0  2   12 30  0  2   14 30  0  1   20 30  0  1   22 30  0 2   24 30  0  3 \
+0 40 -1  1   4 40  5  2   9 40  1  3   10 40  5  2   12 40 -1  1   16 40  5  2   21 40  1  3   22 40  5 2   24 40  0  1 \
+0 50  5  1   4 50 10  5   6 50 10  5   10 50  5  1   12 50  5  1   16 50 10  5   18 50 10  5   20 50  5 1   24 50  0  2 \
+0 60  4  3   3 60 -3  1   7 60 -4  1   10 60  4  1   12 60  4  3   15 60 -3  1   19 60 -4  1   20 60  4 1   24 60  0  1 \
+0 70 -5  5   3 70  0  1   5 70  0  1   10 70 -3  2   12 70 -5  5   15 70  0  1   17 70  0  1   20 70 -3 2   24 70  0  3 \
+0 80  7  1   3 80  1  2   7 80  3  1   10 80  0  1   12 80  7  1   15 80  1  2   19 80  3  1   21 80  0 1   24 80  0  1
+
+bsplinesurf rsurf8 \
+5 5 0 6 1 1 4 1 5 1 8 6 \
+5 5 0 6 2 1 3 1 6 1 7 6 \
+0  0 -8  1   2  0 -7  1   5  0 -9  1   10  0 -8  2   12  0 -9  1   15  0 -7  2   16  0 -8  1   19  0 -8 1   24  0 -8  1 \
+0 10 -8  3   3 10 -7  1   8 10 -9  3   10 10 -8  1   12 10 -9  3   15 10 -7  1   20 10 -8  3   21 10 -8 1   24 10 -8  2 \
+0 20  8  1   4 20  7  3   7 20  9  2   10 20  8  4   12 20  9  1   16 20  7  3   19 20  8  2   20 20  8 4   24 20  8  1 \
+0 30  8  2   2 30  7  1   8 30  9  1   10 30  8  2   12 30  9  2   14 30  7  1   20 30  8  1   22 30  8 2   24 30  8  3 \
+0 40 -8  1   4 40 -7  2   9 40 -9  3   10 40 -8  2   12 40 -9  1   16 40 -7  2   21 40 -8  3   22 40 -8 2   24 40 -8  1 \
+0 50 -8  1   4 50 -7  5   6 50 -9  5   10 50 -8  1   12 50 -9  1   16 50 -7  5   18 50 -8  5   20 50 -8 1   24 50 -8  2 \
+0 60  8  3   3 60  7  1   7 60  9  1   10 60  8  1   12 60  9  3   15 60  7  1   19 60  8  1   20 60  8 1   24 60  8  1 \
+0 70  8  5   3 70  7  1   5 70  9  1   10 70  8  2   12 70  9  5   15 70  7  1   17 70  8  1   20 70  8 2   24 70  8  3 \
+0 80  8  1   3 80  7  2   7 80  9  1   10 80  8  1   12 80  9  1   15 80  7  2   19 80  8  1   21 80  8 1   24 80  8  1
+
+# intersect all surfaces
+set surfaces [list surf1 surf2 rsurf1 surf3 rsurf2 surf4 rsurf3 surf5 rsurf4 surf6 rsurf5 surf7 rsurf6 surf8 rsurf7 rsurf8]
+#set i 1
+#set shapes {}
+foreach s1 $surfaces {
+#    mkface f_$s1 $s1
+#    lappend shapes f_$s1
+    foreach s2 $surfaces {
+        if { $s1 != $s2 } {
+            intersect r $s1 $s2
+#            if { [regexp {a 3d curve} [whatis r]] } {
+#                mkedge e_$i r
+#                lappend shapes e_$i
+#                incr i
+#            }
+        }
+    }
+}
+
+#eval vdisplay $shapes; vfit

diff --git a/tests/perf/grids.list b/tests/perf/grids.list
index 714fd38..c7c4c77 100644 (file)
--- a/tests/perf/grids.list
+++ b/tests/perf/grids.list
@@ -1,2 +1,3 @@
 001 bop
 002 ncollection
+003 bspline