src/PLib/PLib.cxx

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