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-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 // Modified: 28/02/1996 by PMN :  HermiteCoefficients added
  18 // Modified: 18/06/1996 by PMN :  NULL reference.
  19 // Modified: 19/02/1997 by JCT :  EvalPoly2Var added
  20
  21 #include <PLib.ixx>
  22 #include <NCollection_LocalArray.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <math_Gauss.hxx>
  25 #include <Standard_ConstructionError.hxx>
  26 #include <GeomAbs_Shape.hxx>
  27
  28 #include <math_Gauss.hxx>
  29 #include <math.hxx>
  30
  31 // To convert points array into Real ..
  32 // *********************************
  33
  34 //=======================================================================
  35 //function : SetPoles
  36 //purpose  :
  37 //=======================================================================
  38
  39 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
  40                     TColStd_Array1OfReal& FP)
  41 {
  42   Standard_Integer j      = FP   .Lower();
  43   Standard_Integer PLower = Poles.Lower();
  44   Standard_Integer PUpper = Poles.Upper();
  45
  46   for (Standard_Integer i = PLower; i <= PUpper; i++) {
  47     const gp_Pnt2d& P = Poles(i);
  48     FP(j) = P.Coord(1); j++;
  49     FP(j) = P.Coord(2); j++;
  50   }
  51 }
  52
  53 //=======================================================================
  54 //function : SetPoles
  55 //purpose  :
  56 //=======================================================================
  57
  58 void PLib::SetPoles(const TColgp_Array1OfPnt2d&       Poles,
  59                     const TColStd_Array1OfReal& Weights,
  60                     TColStd_Array1OfReal&       FP)
  61 {
  62   Standard_Integer j      = FP   .Lower();
  63   Standard_Integer PLower = Poles.Lower();
  64   Standard_Integer PUpper = Poles.Upper();
  65
  66   for (Standard_Integer i = PLower; i <= PUpper; i++) {
  67     Standard_Real w = Weights(i);
  68     const gp_Pnt2d& P = Poles(i);
  69     FP(j) = P.Coord(1) * w; j++;
  70     FP(j) = P.Coord(2) * w; j++;
  71     FP(j) =              w; j++;
  72   }
  73 }
  74
  75 //=======================================================================
  76 //function : GetPoles
  77 //purpose  :
  78 //=======================================================================
  79
  80 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
  81                     TColgp_Array1OfPnt2d& Poles)
  82 {
  83   Standard_Integer j      = FP   .Lower();
  84   Standard_Integer PLower = Poles.Lower();
  85   Standard_Integer PUpper = Poles.Upper();
  86
  87   for (Standard_Integer i = PLower; i <= PUpper; i++) {
  88     gp_Pnt2d& P = Poles(i);
  89     P.SetCoord(1,FP(j)); j++;
  90     P.SetCoord(2,FP(j)); j++;
  91   }
  92 }
  93
  94 //=======================================================================
  95 //function : GetPoles
  96 //purpose  :
  97 //=======================================================================
  98
  99 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
 100                     TColgp_Array1OfPnt2d&       Poles,
 101                     TColStd_Array1OfReal&       Weights)
 102 {
 103   Standard_Integer j      = FP   .Lower();
 104   Standard_Integer PLower = Poles.Lower();
 105   Standard_Integer PUpper = Poles.Upper();
 106
 107   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 108     Standard_Real w = FP(j + 2);
 109     Weights(i) = w;
 110     gp_Pnt2d& P = Poles(i);
 111     P.SetCoord(1,FP(j) / w); j++;
 112     P.SetCoord(2,FP(j) / w); j++;
 113     j++;
 114   }
 115 }
 116
 117 //=======================================================================
 118 //function : SetPoles
 119 //purpose  :
 120 //=======================================================================
 121
 122 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
 123                     TColStd_Array1OfReal& FP)
 124 {
 125   Standard_Integer j      = FP   .Lower();
 126   Standard_Integer PLower = Poles.Lower();
 127   Standard_Integer PUpper = Poles.Upper();
 128
 129   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 130     const gp_Pnt& P = Poles(i);
 131     FP(j) = P.Coord(1); j++;
 132     FP(j) = P.Coord(2); j++;
 133     FP(j) = P.Coord(3); j++;
 134   }
 135 }
 136
 137 //=======================================================================
 138 //function : SetPoles
 139 //purpose  :
 140 //=======================================================================
 141
 142 void PLib::SetPoles(const TColgp_Array1OfPnt&   Poles,
 143                     const TColStd_Array1OfReal& Weights,
 144                     TColStd_Array1OfReal&       FP)
 145 {
 146   Standard_Integer j      = FP   .Lower();
 147   Standard_Integer PLower = Poles.Lower();
 148   Standard_Integer PUpper = Poles.Upper();
 149
 150   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 151     Standard_Real w = Weights(i);
 152     const gp_Pnt& P = Poles(i);
 153     FP(j) = P.Coord(1) * w; j++;
 154     FP(j) = P.Coord(2) * w; j++;
 155     FP(j) = P.Coord(3) * w; j++;
 156     FP(j) =              w; j++;
 157   }
 158 }
 159
 160 //=======================================================================
 161 //function : GetPoles
 162 //purpose  :
 163 //=======================================================================
 164
 165 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
 166                     TColgp_Array1OfPnt&         Poles)
 167 {
 168   Standard_Integer j      = FP   .Lower();
 169   Standard_Integer PLower = Poles.Lower();
 170   Standard_Integer PUpper = Poles.Upper();
 171
 172   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 173     gp_Pnt& P = Poles(i);
 174     P.SetCoord(1,FP(j)); j++;
 175     P.SetCoord(2,FP(j)); j++;
 176     P.SetCoord(3,FP(j)); j++;
 177   }
 178 }
 179
 180 //=======================================================================
 181 //function : GetPoles
 182 //purpose  :
 183 //=======================================================================
 184
 185 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
 186                     TColgp_Array1OfPnt&         Poles,
 187                     TColStd_Array1OfReal&       Weights)
 188 {
 189   Standard_Integer j      = FP   .Lower();
 190   Standard_Integer PLower = Poles.Lower();
 191   Standard_Integer PUpper = Poles.Upper();
 192
 193   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 194     Standard_Real w = FP(j + 3);
 195     Weights(i) = w;
 196     gp_Pnt& P = Poles(i);
 197     P.SetCoord(1,FP(j) / w); j++;
 198     P.SetCoord(2,FP(j) / w); j++;
 199     P.SetCoord(3,FP(j) / w); j++;
 200     j++;
 201   }
 202 }
 203
 204 // specialized allocator
 205 namespace
 206 {
 207
 208 class BinomAllocator
 209 {
 210 public:
 211
 212   //! Main constructor
 213   BinomAllocator (const Standard_Integer theMaxBinom)
 214   : myBinom (NULL),
 215     myMaxBinom (theMaxBinom)
 216   {
 217     Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
 218     Standard_Integer np1 = myMaxBinom + 1;
 219     myBinom = new Standard_Integer*[np1];
 220     myBinom[0] = new Standard_Integer[1];
 221     myBinom[0][0] = 1;
 222     for (i = 1; i < np1; ++i)
 223     {
 224       im1 = i - 1;
 225       ip1 = i + 1;
 226       id2 = i >> 1;
 227       md2 = im1 >> 1;
 228       md3 = ip1 >> 1;
 229       k   = 0;
 230       myBinom[i] = new Standard_Integer[ip1];
 231
 232       for (j = 0; j < id2; ++j)
 233       {
 234         myBinom[i][j] = k + myBinom[im1][j];
 235         k = myBinom[im1][j];
 236       }
 237       j = id2;
 238       if (j > md2) j = im1 - j;
 239       myBinom[i][id2] = k + myBinom[im1][j];
 240
 241       for (j = ip1 - md3; j < ip1; j++)
 242       {
 243         myBinom[i][j] = myBinom[i][i - j];
 244       }
 245     }
 246   }
 247
 248   //! Destructor
 249   ~BinomAllocator()
 250   {
 251     // free memory
 252     for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
 253     {
 254       delete[] myBinom[i];
 255     }
 256     delete[] myBinom;
 257   }
 258
 259   Standard_Real Value (const Standard_Integer N,
 260                        const Standard_Integer P) const
 261   {
 262     Standard_OutOfRange_Raise_if (N > myMaxBinom,
 263       "PLib, BinomAllocator: requested degree is greater than maximum supported");
 264     return Standard_Real (myBinom[N][P]);
 265   }
 266
 267 private:
 268   BinomAllocator (const BinomAllocator&);
 269   BinomAllocator& operator= (const BinomAllocator&);
 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   }
1405   case 15 : {
1406     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1407     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1408     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1409
1410     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1411     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1412     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1413
1414     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1415     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1416     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1417
1418     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1419     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1420     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1421
1422     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1423     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1424     *tmpRA = *tmpPA; tmpRA++;
1425     tmpPA -= 14;
1426     if (DerivativeRequest > 0 ) {
1427       Standard_Real *valRA;
1428       Standard_Integer ii, LocalRequest;
1429       Standard_Integer Index1, Index2;
1430       Standard_Integer MaxIndex1, MaxIndex2;
1431       if      (DerivativeRequest < Degree) {
1432         LocalRequest = DerivativeRequest;
1433         MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest;
1434       }
1435       else {
1436         LocalRequest = Degree;
1437         MaxIndex2 = MaxIndex1 = DegreeDimension;
1438       }
1439       MaxIndex2 -= 15;
1440
1441       for (ii = 1; ii <= LocalRequest; ii++) {
1442         *tmpRA = 0.0e0; tmpRA++;
1443         *tmpRA = 0.0e0; tmpRA++;
1444         *tmpRA = 0.0e0; tmpRA++;
1445
1446         *tmpRA = 0.0e0; tmpRA++;
1447         *tmpRA = 0.0e0; tmpRA++;
1448         *tmpRA = 0.0e0; tmpRA++;
1449
1450         *tmpRA = 0.0e0; tmpRA++;
1451         *tmpRA = 0.0e0; tmpRA++;
1452         *tmpRA = 0.0e0; tmpRA++;
1453
1454         *tmpRA = 0.0e0; tmpRA++;
1455         *tmpRA = 0.0e0; tmpRA++;
1456         *tmpRA = 0.0e0; tmpRA++;
1457
1458         *tmpRA = 0.0e0; tmpRA++;
1459         *tmpRA = 0.0e0; tmpRA++;
1460         *tmpRA = 0.0e0; tmpRA++;
1461       }
1462
1463       for (jj = Degree ; jj >  0 ; jj--) {
1464         tmpPA -= 15;
1465
1466         Index1 = MaxIndex1;
1467         Index2 = MaxIndex2;
1468
1469         for (ii = LocalRequest ; ii > 0 ; ii--) {
1470           valRA  = &RA[Index1];
1471           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1472           Index1 -= 15;
1473           Index2 -= 15;
1474         }
1475         valRA  = &RA[Index1];
1476         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1477
1478         Index1 = MaxIndex1 + 1;
1479         Index2 = MaxIndex2 + 1;
1480
1481         for (ii = LocalRequest ; ii > 0 ; ii--) {
1482           valRA  = &RA[Index1];
1483           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1484           Index1 -= 15;
1485           Index2 -= 15;
1486         }
1487         valRA  = &RA[Index1];
1488         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1489
1490         Index1 = MaxIndex1 + 2;
1491         Index2 = MaxIndex2 + 2;
1492
1493         for (ii = LocalRequest ; ii > 0 ; ii--) {
1494           valRA  = &RA[Index1];
1495           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1496           Index1 -= 15;
1497           Index2 -= 15;
1498         }
1499         valRA  = &RA[Index1];
1500         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1501
1502         Index1 = MaxIndex1 + 3;
1503         Index2 = MaxIndex2 + 3;
1504
1505         for (ii = LocalRequest ; ii > 0 ; ii--) {
1506           valRA  = &RA[Index1];
1507           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1508           Index1 -= 15;
1509           Index2 -= 15;
1510         }
1511         valRA  = &RA[Index1];
1512         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1513
1514         Index1 = MaxIndex1 + 4;
1515         Index2 = MaxIndex2 + 4;
1516
1517         for (ii = LocalRequest ; ii > 0 ; ii--) {
1518           valRA  = &RA[Index1];
1519           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1520           Index1 -= 15;
1521           Index2 -= 15;
1522         }
1523         valRA  = &RA[Index1];
1524         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1525
1526         Index1 = MaxIndex1 + 5;
1527         Index2 = MaxIndex2 + 5;
1528
1529         for (ii = LocalRequest ; ii > 0 ; ii--) {
1530           valRA  = &RA[Index1];
1531           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1532           Index1 -= 15;
1533           Index2 -= 15;
1534         }
1535         valRA  = &RA[Index1];
1536         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1537
1538         Index1 = MaxIndex1 + 6;
1539         Index2 = MaxIndex2 + 6;
1540
1541         for (ii = LocalRequest ; ii > 0 ; ii--) {
1542           valRA  = &RA[Index1];
1543           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1544           Index1 -= 15;
1545           Index2 -= 15;
1546         }
1547         valRA  = &RA[Index1];
1548         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1549
1550         Index1 = MaxIndex1 + 7;
1551         Index2 = MaxIndex2 + 7;
1552
1553         for (ii = LocalRequest ; ii > 0 ; ii--) {
1554           valRA  = &RA[Index1];
1555           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1556           Index1 -= 15;
1557           Index2 -= 15;
1558         }
1559         valRA  = &RA[Index1];
1560         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1561
1562         Index1 = MaxIndex1 + 8;
1563         Index2 = MaxIndex2 + 8;
1564
1565         for (ii = LocalRequest ; ii > 0 ; ii--) {
1566           valRA  = &RA[Index1];
1567           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1568           Index1 -= 15;
1569           Index2 -= 15;
1570         }
1571         valRA  = &RA[Index1];
1572         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1573
1574         Index1 = MaxIndex1 + 9;
1575         Index2 = MaxIndex2 + 9;
1576
1577         for (ii = LocalRequest ; ii > 0 ; ii--) {
1578           valRA  = &RA[Index1];
1579           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1580           Index1 -= 15;
1581           Index2 -= 15;
1582         }
1583         valRA  = &RA[Index1];
1584         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1585
1586         Index1 = MaxIndex1 + 10;
1587         Index2 = MaxIndex2 + 10;
1588
1589         for (ii = LocalRequest ; ii > 0 ; ii--) {
1590           valRA  = &RA[Index1];
1591           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1592           Index1 -= 15;
1593           Index2 -= 15;
1594         }
1595         valRA  = &RA[Index1];
1596         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1597
1598         Index1 = MaxIndex1 + 11;
1599         Index2 = MaxIndex2 + 11;
1600
1601         for (ii = LocalRequest ; ii > 0 ; ii--) {
1602           valRA  = &RA[Index1];
1603           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1604           Index1 -= 15;
1605           Index2 -= 15;
1606         }
1607         valRA  = &RA[Index1];
1608         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1609
1610         Index1 = MaxIndex1 + 12;
1611         Index2 = MaxIndex2 + 12;
1612
1613         for (ii = LocalRequest ; ii > 0 ; ii--) {
1614           valRA  = &RA[Index1];
1615           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1616           Index1 -= 15;
1617           Index2 -= 15;
1618         }
1619         valRA  = &RA[Index1];
1620         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1621
1622         Index1 = MaxIndex1 + 13;
1623         Index2 = MaxIndex2 + 13;
1624
1625         for (ii = LocalRequest ; ii > 0 ; ii--) {
1626           valRA  = &RA[Index1];
1627           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1628           Index1 -= 15;
1629           Index2 -= 15;
1630         }
1631         valRA  = &RA[Index1];
1632         *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1633
1634         Index1 = MaxIndex1 + 14;
1635         Index2 = MaxIndex2 + 14;
1636
1637         for (ii = LocalRequest ; ii > 0 ; ii--) {
1638           valRA  = &RA[Index1];
1639           *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1640           Index1 -= 15;
1641           Index2 -= 15;
1642         }
1643         valRA  = &RA[Index1];
1644         *valRA = Par * (*valRA) + (*tmpPA);
1645
1646         tmpPA -= 14;
1647       }
1648     }
1649     else {
1650
1651       for (jj = Degree  ; jj >  0 ; jj--) {
1652         tmpPA -= 15;
1653         tmpRA  = RA;
1654
1655         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1656         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1657         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1658
1659         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1660         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1661         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1662
1663         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1664         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1665         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1666
1667         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1668         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1669         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1670
1671         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1672         *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1673         *tmpRA = Par * (*tmpRA) + (*tmpPA);
1674         tmpPA -= 14;
1675       }
1676     }
1677     break;
1678   }
1679     default : {
1680       Standard_Integer kk ;
1681       for ( kk = 0 ; kk < Dimension ; kk++) {
1682         *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1683       }
1684       tmpPA -= Dimension;
1685       if (DerivativeRequest > 0 ) {
1686         Standard_Real *valRA;
1687         Standard_Integer ii, LocalRequest;
1688         Standard_Integer Index1, Index2;
1689         Standard_Integer MaxIndex1, MaxIndex2;
1690         if      (DerivativeRequest < Degree) {
1691           LocalRequest = DerivativeRequest;
1692           MaxIndex2 = MaxIndex1 = LocalRequest * Dimension;
1693         }
1694         else {
1695           LocalRequest = Degree;
1696           MaxIndex2 = MaxIndex1 = DegreeDimension;
1697         }
1698         MaxIndex2 -= Dimension;
1699
1700         for (ii = 1; ii <= MaxIndex1; ii++) {
1701           *tmpRA = 0.0e0; tmpRA++;
1702         }
1703
1704         for (jj = Degree ; jj >  0 ; jj--) {
1705           tmpPA -= Dimension;
1706
1707           for (kk = 0 ; kk < Dimension ; kk++) {
1708             Index1 = MaxIndex1 + kk;
1709             Index2 = MaxIndex2 + kk;
1710
1711             for (ii = LocalRequest ; ii > 0 ; ii--) {
1712               valRA  = &RA[Index1];
1713               *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1714               Index1 -= Dimension;
1715               Index2 -= Dimension;
1716             }
1717             valRA  = &RA[Index1];
1718             *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1719           }
1720           tmpPA -= Dimension;
1721         }
1722       }
1723       else {
1724
1725         for (jj = Degree  ; jj >  0 ; jj--) {
1726           tmpPA -= Dimension;
1727           tmpRA  = RA;
1728
1729           for (kk = 0 ; kk < Dimension ; kk++) {
1730             *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1731           }
1732           tmpPA -= Dimension;
1733         }
1734       }
1735     }
1736   }
1737 }
1738
1739 //=======================================================================
1740 //function : This evaluates a polynomial without derivative
1741 //purpose  :
1742 //=======================================================================
1743
1744 void  PLib::NoDerivativeEvalPolynomial(const Standard_Real    Par,
1745                                        const Standard_Integer Degree,
1746                                        const Standard_Integer Dimension,
1747                                        const Standard_Integer DegreeDimension,
1748                                        Standard_Real&         PolynomialCoeff,
1749                                        Standard_Real&         Results)
1750 {
1751   Standard_Integer jj;
1752   Standard_Real *RA = &Results ;
1753   Standard_Real *PA = &PolynomialCoeff ;
1754   Standard_Real *tmpRA = RA;
1755   Standard_Real *tmpPA = PA + DegreeDimension;
1756
1757   switch (Dimension) {
1758
1759   case 1 : {
1760     *tmpRA = *tmpPA;
1761
1762     for (jj = Degree  ; jj >  0 ; jj--) {
1763       tmpPA--;
1764
1765       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1766     }
1767     break;
1768   }
1769   case 2 : {
1770     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1771     *tmpRA = *tmpPA;
1772     tmpPA--;
1773
1774     for (jj = Degree  ; jj >  0 ; jj--) {
1775       tmpPA -= 2;
1776       tmpRA  = RA;
1777
1778       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1779       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1780       tmpPA--;
1781     }
1782     break;
1783   }
1784   case 3 : {
1785     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1786     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1787     *tmpRA = *tmpPA;
1788     tmpPA -= 2;
1789
1790     for (jj = Degree  ; jj >  0 ; jj--) {
1791       tmpPA -= 3;
1792       tmpRA  = RA;
1793
1794       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1795       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1796       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1797       tmpPA -= 2;
1798     }
1799     break;
1800   }
1801   case 6 : {
1802     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1803     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1804     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1805
1806     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1807     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1808     *tmpRA = *tmpPA;
1809     tmpPA -= 5;
1810
1811     for (jj = Degree  ; jj >  0 ; jj--) {
1812       tmpPA -= 6;
1813       tmpRA  = RA;
1814
1815       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1816       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1817       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1818
1819       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1820       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1821       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1822       tmpPA -= 5;
1823     }
1824     break;
1825   }
1826   case 9 : {
1827     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1828     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1829     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1830
1831     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1832     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1833     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1834
1835     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1836     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1837     *tmpRA = *tmpPA;
1838     tmpPA -= 8;
1839
1840     for (jj = Degree  ; jj >  0 ; jj--) {
1841       tmpPA -= 9;
1842       tmpRA  = RA;
1843
1844       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1845       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1846       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1847
1848       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1849       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1850       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1851
1852       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1853       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1854       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1855       tmpPA -= 8;
1856     }
1857     break;
1858   }
1859   case 12 : {
1860     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1861     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1862     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1863
1864     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1865     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1866     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1867
1868     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1869     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1870     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1871
1872     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1873     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1874     *tmpRA = *tmpPA;
1875     tmpPA -= 11;
1876
1877     for (jj = Degree  ; jj >  0 ; jj--) {
1878       tmpPA -= 12;
1879       tmpRA  = RA;
1880
1881       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1882       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1883       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1884
1885       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1886       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1887       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1888
1889       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1890       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1891       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1892
1893       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1894       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1895       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1896       tmpPA -= 11;
1897     }
1898     break;
1899   }
1900   case 15 : {
1901     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1902     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1903     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1904
1905     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1906     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1907     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1908
1909     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1910     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1911     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1912
1913     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1914     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1915     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1916
1917     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1918     *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1919     *tmpRA = *tmpPA;
1920     tmpPA -= 14;
1921
1922     for (jj = Degree  ; jj >  0 ; jj--) {
1923       tmpPA -= 15;
1924       tmpRA  = RA;
1925
1926       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1927       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1928       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1929
1930       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1931       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1932       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1933
1934       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1935       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1936       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1937
1938       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1939       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1940       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1941
1942       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1943       *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1944       *tmpRA = Par * (*tmpRA) + (*tmpPA);
1945       tmpPA -= 14;
1946     }
1947     break;
1948   }
1949     default : {
1950       Standard_Integer kk ;
1951       for ( kk = 0 ; kk < Dimension ; kk++) {
1952         *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1953       }
1954       tmpPA -= Dimension;
1955
1956       for (jj = Degree  ; jj >  0 ; jj--) {
1957         tmpPA -= Dimension;
1958         tmpRA  = RA;
1959
1960         for (kk = 0 ; kk < Dimension ; kk++) {
1961           *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1962         }
1963         tmpPA -= Dimension;
1964       }
1965     }
1966   }
1967 }
1968
1969 //=======================================================================
1970 //function : This evaluates a polynomial of 2 variables
1971 //purpose  : or its derivative at the requested orders
1972 //=======================================================================
1973
1974 void  PLib::EvalPoly2Var(const Standard_Real    UParameter,
1975                          const Standard_Real    VParameter,
1976                          const Standard_Integer UDerivativeRequest,
1977                          const Standard_Integer VDerivativeRequest,
1978                          const Standard_Integer UDegree,
1979                          const Standard_Integer VDegree,
1980                          const Standard_Integer Dimension,
1981                          Standard_Real&         PolynomialCoeff,
1982                          Standard_Real&         Results)
1983      //
1984      // the polynomial coefficients are assumed to be stored as follows :
1985      //                                                      0 0
1986      // [0]                  [Dimension -1]                 U V    coefficient
1987      //                                                      1 0
1988      // [Dimension]          [Dimension + Dimension -1]     U V    coefficient
1989      //                                                      2 0
1990      // [2 * Dimension]      [2 * Dimension + Dimension-1]  U V    coefficient
1991      //
1992      //   ...................................................
1993      //
1994      //
1995      //                                                      m 0
1996      // [m * Dimension]      [m * Dimension + Dimension-1]  U V    coefficient
1997      //
1998      //  where m = UDegree
1999      //
2000      //                                                           0 1
2001      // [(m+1) * Dimension]  [(m+1) * Dimension + Dimension-1]   U V   coefficient
2002      //
2003      //   ...................................................
2004      //
2005      //                                                          m 1
2006      // [2*m * Dimension]      [2*m * Dimension + Dimension-1]  U V    coefficient
2007      //
2008      //   ...................................................
2009      //
2010      //                                                          m n
2011      // [m*n * Dimension]      [m*n * Dimension + Dimension-1]  U V    coefficient
2012      //
2013      //  where n = VDegree
2014 {
2015   Standard_Integer Udim = (VDegree+1)*Dimension,
2016   index = Udim*UDerivativeRequest;
2017   TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
2018   TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
2019   Standard_Real * Result =  (Standard_Real *) &Curve.ChangeValue(1);
2020   Standard_Real * Digit  =  (Standard_Real *) &Point.ChangeValue(1);
2021   Standard_Real * ResultArray ;
2022   ResultArray = &Results ;
2023
2024   PLib::EvalPolynomial(UParameter,
2025                        UDerivativeRequest,
2026                        UDegree,
2027                        Udim,
2028                        PolynomialCoeff,
2029                        Result[0]);
2030
2031   PLib::EvalPolynomial(VParameter,
2032                        VDerivativeRequest,
2033                        VDegree,
2034                        Dimension,
2035                        Result[index],
2036                        Digit[0]);
2037
2038   index = Dimension*VDerivativeRequest;
2039
2040   for (Standard_Integer i=0;i<Dimension;i++) {
2041     ResultArray[i] = Digit[index+i];
2042   }
2043 }
2044
2045
2046
2047 //=======================================================================
2048 //function : This evaluates the lagrange polynomial and its derivatives
2049 //purpose  : up to the requested order that interpolates a series of
2050 //points of dimension <Dimension> with given assigned parameters
2051 //=======================================================================
2052
2053 Standard_Integer
2054 PLib::EvalLagrange(const Standard_Real                   Parameter,
2055                    const Standard_Integer                DerivativeRequest,
2056                    const Standard_Integer                Degree,
2057                    const Standard_Integer                Dimension,
2058                    Standard_Real&                        Values,
2059                    Standard_Real&                        Parameters,
2060                    Standard_Real&                        Results)
2061 {
2062   //
2063   // the points  are assumed to be stored as follows in the Values array :
2064   //
2065   // [0]              [Dimension -1]                first  point   coefficients
2066   //
2067   // [Dimension]      [Dimension + Dimension -1]    second point   coefficients
2068   //
2069   // [2 * Dimension]  [2 * Dimension + Dimension-1] third  point   coefficients
2070   //
2071   //   ...................................................
2072   //
2073   //
2074   //
2075   // [d * Dimension]  [d * Dimension + Dimension-1] d + 1 point   coefficients
2076   //
2077   //  where d is the Degree
2078   //
2079   //  The ParameterArray stores the parameter value assign to each point in
2080   //  order described above, that is
2081   //  [0]   is assign to first  point
2082   //  [1]   is assign to second point
2083   //
2084   Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
2085   Standard_Integer local_request = DerivativeRequest;
2086   Standard_Real  *ParameterArray;
2087   Standard_Real  difference;
2088   Standard_Real  *PointsArray;
2089   Standard_Real  *ResultArray ;
2090
2091   PointsArray    = &Values ;
2092   ParameterArray = &Parameters ;
2093   ResultArray = &Results ;
2094   if (local_request >= Degree) {
2095     local_request = Degree ;
2096   }
2097   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
2098   //
2099   //  Build the divided differences array
2100   //
2101
2102   for (ii = 0 ; ii <  (Degree + 1) * Dimension  ; ii++) {
2103     divided_differences_array[ii] = PointsArray[ii] ;
2104   }
2105
2106   for (ii = Degree ; ii >= 0   ; ii--) {
2107
2108     for (jj = Degree  ; jj > Degree - ii  ; jj--) {
2109       Index = jj * Dimension ;
2110       Index1 = Index - Dimension ;
2111
2112       for (kk = 0 ; kk < Dimension ; kk++) {
2113         divided_differences_array[Index + kk] -=
2114           divided_differences_array[Index1 + kk] ;
2115       }
2116       difference =
2117         ParameterArray[jj] - ParameterArray[jj - Degree -1 +  ii] ;
2118       if (Abs(difference) < RealSmall()) {
2119         ReturnCode = 1 ;
2120         goto FINISH ;
2121       }
2122       difference = 1.0e0 / difference ;
2123
2124       for (kk = 0 ; kk < Dimension ; kk++) {
2125         divided_differences_array[Index + kk] *= difference ;
2126       }
2127     }
2128   }
2129   //
2130   //
2131   // Evaluate the divided difference array polynomial which expresses as
2132   //
2133   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2134   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2135   //
2136   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2137   //
2138   //
2139   Index = Degree * Dimension ;
2140
2141   for (kk = 0 ; kk < Dimension ; kk++) {
2142     ResultArray[kk] = divided_differences_array[Index + kk] ;
2143   }
2144
2145   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
2146     ResultArray[ii] = 0.0e0 ;
2147   }
2148
2149   for (ii = Degree ; ii >= 1 ; ii--) {
2150     difference =  Parameter - ParameterArray[ii - 1] ;
2151
2152     for (jj = local_request ; jj > 0 ; jj--) {
2153       Index = jj * Dimension ;
2154       Index1 = Index - Dimension ;
2155
2156       for (kk = 0 ; kk < Dimension ; kk++) {
2157         ResultArray[Index + kk] *= difference ;
2158         ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
2159       }
2160     }
2161     Index = (ii -1) * Dimension ;
2162
2163     for (kk = 0 ; kk < Dimension ; kk++) {
2164       ResultArray[kk] *= difference ;
2165       ResultArray[kk] += divided_differences_array[Index+kk] ;
2166     }
2167   }
2168   FINISH :
2169     return (ReturnCode) ;
2170 }
2171
2172 //=======================================================================
2173 //function : This evaluates the hermite polynomial and its derivatives
2174 //purpose  : up to the requested order that interpolates a series of
2175 //points of dimension <Dimension> with given assigned parameters
2176 //=======================================================================
2177
2178 Standard_Integer PLib::EvalCubicHermite
2179 (const Standard_Real          Parameter,
2180  const Standard_Integer       DerivativeRequest,
2181  const Standard_Integer       Dimension,
2182  Standard_Real&               Values,
2183  Standard_Real&               Derivatives,
2184  Standard_Real&               theParameters,
2185  Standard_Real&               Results)
2186 {
2187   //
2188   // the points  are assumed to be stored as follows in the Values array :
2189   //
2190   // [0]            [Dimension -1]             first  point   coefficients
2191   //
2192   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
2193   //
2194   //
2195   // the derivatives  are assumed to be stored as follows in
2196   // the Derivatives array :
2197   //
2198   // [0]            [Dimension -1]             first  point   coefficients
2199   //
2200   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
2201   //
2202   //  The ParameterArray stores the parameter value assign to each point in
2203   //  order described above, that is
2204   //  [0]   is assign to first  point
2205   //  [1]   is assign to last   point
2206   //
2207   Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
2208   Standard_Integer local_request = DerivativeRequest ;
2209
2210   ReturnCode = 0 ;
2211   Degree = 3 ;
2212   Standard_Real  ParametersArray[4];
2213   Standard_Real  difference;
2214   Standard_Real  inverse;
2215   Standard_Real  *FirstLast;
2216   Standard_Real  *PointsArray;
2217   Standard_Real  *DerivativesArray;
2218   Standard_Real  *ResultArray ;
2219
2220   DerivativesArray = &Derivatives ;
2221   PointsArray    = &Values ;
2222   FirstLast = &theParameters ;
2223   ResultArray = &Results ;
2224   if (local_request >= Degree) {
2225     local_request = Degree ;
2226   }
2227   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
2228
2229   for (ii = 0, jj = 0  ; ii < 2 ; ii++, jj+= 2) {
2230     ParametersArray[jj] =
2231       ParametersArray[jj+1] = FirstLast[ii] ;
2232   }
2233   //
2234   //  Build the divided differences array
2235   //
2236   //
2237   // initialise it at the stage 2 of the building algorithm
2238   // for devided differences
2239   //
2240   inverse = FirstLast[1] - FirstLast[0] ;
2241   inverse = 1.0e0 / inverse ;
2242
2243   for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension  ;
2244        ii <  Dimension  ;
2245        ii++, jj++, kk++, pp++) {
2246     divided_differences_array[ii] = PointsArray[ii] ;
2247     divided_differences_array[kk] = inverse *
2248       (PointsArray[jj] - PointsArray[ii]) ;
2249     divided_differences_array[jj] = DerivativesArray[ii] ;
2250     divided_differences_array[pp] = DerivativesArray[jj] ;
2251   }
2252
2253   for (ii = 1 ; ii <= Degree   ; ii++) {
2254
2255     for (jj = Degree  ; jj >=  ii+1  ; jj--) {
2256       Index = jj * Dimension ;
2257       Index1 = Index - Dimension ;
2258
2259       for (kk = 0 ; kk < Dimension ; kk++) {
2260         divided_differences_array[Index + kk] -=
2261           divided_differences_array[Index1 + kk] ;
2262       }
2263
2264       for (kk = 0 ; kk < Dimension ; kk++) {
2265         divided_differences_array[Index + kk] *= inverse ;
2266       }
2267     }
2268   }
2269   //
2270   //
2271   // Evaluate the divided difference array polynomial which expresses as
2272   //
2273   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2274   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2275   //
2276   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2277   //
2278   //
2279   Index = Degree * Dimension ;
2280
2281   for (kk = 0 ; kk < Dimension ; kk++) {
2282     ResultArray[kk] = divided_differences_array[Index + kk] ;
2283   }
2284
2285   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
2286     ResultArray[ii] = 0.0e0 ;
2287   }
2288
2289   for (ii = Degree ; ii >= 1 ; ii--) {
2290     difference =  Parameter - ParametersArray[ii - 1] ;
2291
2292     for (jj = local_request ; jj > 0 ; jj--) {
2293       Index = jj * Dimension ;
2294       Index1 = Index - Dimension ;
2295
2296       for (kk = 0 ; kk < Dimension ; kk++) {
2297         ResultArray[Index + kk] *= difference ;
2298         ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
2299       }
2300     }
2301     Index = (ii -1) * Dimension ;
2302
2303     for (kk = 0 ; kk < Dimension ; kk++) {
2304       ResultArray[kk] *= difference ;
2305       ResultArray[kk] += divided_differences_array[Index+kk] ;
2306     }
2307   }
2308 //  FINISH :
2309     return (ReturnCode) ;
2310 }
2311
2312 //=======================================================================
2313 //function : HermiteCoefficients
2314 //purpose  : calcul des polynomes d'Hermite
2315 //=======================================================================
2316
2317 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
2318                                            const Standard_Real LastParameter,
2319                                            const Standard_Integer FirstOrder,
2320                                            const Standard_Integer LastOrder,
2321                                            math_Matrix& MatrixCoefs)
2322 {
2323   Standard_Integer NbCoeff =  FirstOrder +  LastOrder + 2, Ordre[2];
2324   Standard_Integer ii, jj, pp, cote, iof=0;
2325   Standard_Real Prod, TBorne = FirstParameter;
2326   math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
2327   math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
2328
2329   // Test de validites
2330
2331   if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
2332   Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
2333   if (D1 > 100 || D2 > 100) return Standard_False;
2334   D2 += D1;
2335   if (D2 < 0.01) return Standard_False;
2336   if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
2337
2338   // Calcul de la matrice a inverser (MAT)
2339
2340   Ordre[0] = FirstOrder+1;
2341   Ordre[1] = LastOrder+1;
2342
2343   for (cote=0; cote<=1; cote++) {
2344     Coeff.Init(1);
2345
2346     for (pp=1; pp<=Ordre[cote]; pp++) {
2347       ii = pp + iof;
2348       Prod = 1;
2349
2350       for (jj=pp; jj<=NbCoeff; jj++) {
2351         //        tout se passe dans les 3 lignes suivantes
2352         MAT(ii, jj) = Coeff(jj) * Prod;
2353         Coeff(jj) *= jj - pp;
2354         Prod      *= TBorne;
2355       }
2356     }
2357     TBorne = LastParameter;
2358     iof = Ordre[0];
2359   }
2360
2361   // resolution du systemes
2362   math_Gauss ResolCoeff(MAT, 1.0e-10);
2363   if (!ResolCoeff.IsDone()) return Standard_False;
2364
2365   for (ii=1; ii<=NbCoeff; ii++) {
2366     B(ii) = 1;
2367     ResolCoeff.Solve(B, Coeff);
2368     MatrixCoefs.SetRow( ii, Coeff);
2369     B(ii) = 0;
2370   }
2371   return Standard_True;
2372 }
2373
2374 //=======================================================================
2375 //function : CoefficientsPoles
2376 //purpose  :
2377 //=======================================================================
2378
2379 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt&   Coefs,
2380                               const TColStd_Array1OfReal& WCoefs,
2381                               TColgp_Array1OfPnt&         Poles,
2382                               TColStd_Array1OfReal&       Weights)
2383 {
2384   TColStd_Array1OfReal tempC(1,3*Coefs.Length());
2385   PLib::SetPoles(Coefs,tempC);
2386   TColStd_Array1OfReal tempP(1,3*Poles.Length());
2387   PLib::SetPoles(Coefs,tempP);
2388   PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
2389   PLib::GetPoles(tempP,Poles);
2390 }
2391
2392 //=======================================================================
2393 //function : CoefficientsPoles
2394 //purpose  :
2395 //=======================================================================
2396
2397 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
2398                               const TColStd_Array1OfReal& WCoefs,
2399                               TColgp_Array1OfPnt2d&       Poles,
2400                               TColStd_Array1OfReal&       Weights)
2401 {
2402   TColStd_Array1OfReal tempC(1,2*Coefs.Length());
2403   PLib::SetPoles(Coefs,tempC);
2404   TColStd_Array1OfReal tempP(1,2*Poles.Length());
2405   PLib::SetPoles(Coefs,tempP);
2406   PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
2407   PLib::GetPoles(tempP,Poles);
2408 }
2409
2410 //=======================================================================
2411 //function : CoefficientsPoles
2412 //purpose  :
2413 //=======================================================================
2414
2415 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
2416                               const TColStd_Array1OfReal& WCoefs,
2417                               TColStd_Array1OfReal&       Poles,
2418                               TColStd_Array1OfReal&       Weights)
2419 {
2420   PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
2421 }
2422
2423 //=======================================================================
2424 //function : CoefficientsPoles
2425 //purpose  :
2426 //=======================================================================
2427
2428 void PLib::CoefficientsPoles (const Standard_Integer      dim,
2429                               const TColStd_Array1OfReal& Coefs,
2430                               const TColStd_Array1OfReal& WCoefs,
2431                               TColStd_Array1OfReal&       Poles,
2432                               TColStd_Array1OfReal&       Weights)
2433 {
2434   Standard_Boolean rat = &WCoefs != NULL;
2435   Standard_Integer loc = Coefs.Lower();
2436   Standard_Integer lop = Poles.Lower();
2437   Standard_Integer lowc=0;
2438   Standard_Integer lowp=0;
2439   Standard_Integer upc = Coefs.Upper();
2440   Standard_Integer upp = Poles.Upper();
2441   Standard_Integer upwc=0;
2442   Standard_Integer upwp=0;
2443   Standard_Integer reflen = Coefs.Length()/dim;
2444   Standard_Integer i,j,k;
2445   //Les Extremites.
2446   if (rat) {
2447     lowc = WCoefs.Lower(); lowp = Weights.Lower();
2448     upwc = WCoefs.Upper(); upwp = Weights.Upper();
2449   }
2450
2451   for (i = 0; i < dim; i++){
2452     Poles (lop + i) = Coefs (loc + i);
2453     Poles (upp - i) = Coefs (upc - i);
2454   }
2455   if (rat) {
2456     Weights (lowp) = WCoefs (lowc);
2457     Weights (upwp) = WCoefs (upwc);
2458   }
2459
2460   Standard_Real Cnp;
2461   for (i = 2; i < reflen; i++ ) {
2462     Cnp = PLib::Bin(reflen - 1, i - 1);
2463     if (rat) Weights (lowp + i - 1) = WCoefs (lowc + i - 1) / Cnp;
2464
2465     for(j = 0; j < dim; j++){
2466       Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
2467     }
2468   }
2469
2470   for (i = 1; i <= reflen - 1; i++) {
2471
2472     for (j = reflen - 1; j >= i; j--) {
2473       if (rat) Weights (lowp + j) += Weights (lowp + j -1);
2474
2475       for(k = 0; k < dim; k++){
2476         Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
2477       }
2478     }
2479   }
2480   if (rat) {
2481
2482     for (i = 1; i <= reflen; i++) {
2483
2484       for(j = 0; j < dim; j++){
2485         Poles(lop + dim * (i-1) + j) /= Weights(lowp + i -1);
2486       }
2487     }
2488   }
2489 }
2490
2491 //=======================================================================
2492 //function : Trimming
2493 //purpose  :
2494 //=======================================================================
2495
2496 void PLib::Trimming(const Standard_Real   U1,
2497                     const Standard_Real   U2,
2498                     TColgp_Array1OfPnt&   Coefs,
2499                     TColStd_Array1OfReal& WCoefs)
2500 {
2501   TColStd_Array1OfReal temp(1,3*Coefs.Length());
2502   PLib::SetPoles(Coefs,temp);
2503   PLib::Trimming(U1,U2,3,temp,WCoefs);
2504   PLib::GetPoles(temp,Coefs);
2505 }
2506
2507 //=======================================================================
2508 //function : Trimming
2509 //purpose  :
2510 //=======================================================================
2511
2512 void PLib::Trimming(const Standard_Real   U1,
2513                     const Standard_Real   U2,
2514                     TColgp_Array1OfPnt2d& Coefs,
2515                     TColStd_Array1OfReal& WCoefs)
2516 {
2517   TColStd_Array1OfReal temp(1,2*Coefs.Length());
2518   PLib::SetPoles(Coefs,temp);
2519   PLib::Trimming(U1,U2,2,temp,WCoefs);
2520   PLib::GetPoles(temp,Coefs);
2521 }
2522
2523 //=======================================================================
2524 //function : Trimming
2525 //purpose  :
2526 //=======================================================================
2527
2528 void PLib::Trimming(const Standard_Real   U1,
2529                     const Standard_Real   U2,
2530                     TColStd_Array1OfReal& Coefs,
2531                     TColStd_Array1OfReal& WCoefs)
2532 {
2533   PLib::Trimming(U1,U2,1,Coefs,WCoefs);
2534 }
2535
2536 //=======================================================================
2537 //function : Trimming
2538 //purpose  :
2539 //=======================================================================
2540
2541 void PLib::Trimming(const Standard_Real U1,
2542                     const Standard_Real U2,
2543                     const Standard_Integer dim,
2544                     TColStd_Array1OfReal& Coefs,
2545                     TColStd_Array1OfReal& WCoefs)
2546 {
2547
2548   // principe :
2549   // on fait le changement de variable v = (u-U1) / (U2-U1)
2550   // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
2551   // decomposee sous la forme du schema iteratif de horner.
2552
2553   Standard_Real lsp = U2 - U1;
2554   Standard_Integer indc, indw=0;
2555   Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
2556   Standard_Integer len = Coefs.Length()/dim;
2557   Standard_Boolean rat = &WCoefs != NULL;
2558
2559   if (rat) {
2560     if(len != WCoefs.Length())
2561       Standard_Failure::Raise("PLib::Trimming : nbcoefs/dim != nbweights !!!");
2562     upw = WCoefs.Upper();
2563   }
2564   len --;
2565
2566   for (Standard_Integer i = 1; i <= len; i++) {
2567     Standard_Integer j ;
2568     indc = upc - dim*(i-1);
2569     if (rat) indw = upw - i + 1;
2570     //calcul du coefficient de degre le plus faible a l'iteration i
2571
2572     for( j = 0; j < dim; j++){
2573       Coefs(indc - dim + j) += U1 * Coefs(indc + j);
2574     }
2575     if (rat) WCoefs(indw - 1) += U1 * WCoefs(indw);
2576
2577     //calcul des coefficients intermediaires :
2578
2579     while (indc < upc){
2580       indc += dim;
2581
2582       for(Standard_Integer k = 0; k < dim; k++){
2583         Coefs(indc - dim + k) =
2584           U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
2585       }
2586       if (rat) {
2587         indw ++;
2588         WCoefs(indw - 1) = U1 * WCoefs(indw) + lsp * WCoefs(indw - 1);
2589       }
2590     }
2591
2592     //calcul du coefficient de degre le plus eleve :
2593
2594     for(j = 0; j < dim; j++){
2595       Coefs(upc + j) *= lsp;
2596     }
2597     if (rat) WCoefs(upw) *= lsp;
2598   }
2599 }
2600
2601 //=======================================================================
2602 //function : CoefficientsPoles
2603 //purpose  :
2604 // Modified: 21/10/1996 by PMN :  PolesCoefficient (PRO5852).
2605 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
2606 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
2607 //=======================================================================
2608
2609 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt&   Coefs,
2610                               const TColStd_Array2OfReal& WCoefs,
2611                               TColgp_Array2OfPnt&         Poles,
2612                               TColStd_Array2OfReal&       Weights)
2613 {
2614   Standard_Boolean rat = (&WCoefs != NULL);
2615   Standard_Integer LowerRow  = Poles.LowerRow();
2616   Standard_Integer UpperRow  = Poles.UpperRow();
2617   Standard_Integer LowerCol  = Poles.LowerCol();
2618   Standard_Integer UpperCol  = Poles.UpperCol();
2619   Standard_Integer ColLength = Poles.ColLength();
2620   Standard_Integer RowLength = Poles.RowLength();
2621
2622   // Bidouille pour retablir u et v pour les coefs calcules
2623   // par buildcache
2624 //  Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
2625
2626   Standard_Integer Row, Col;
2627   Standard_Real W, Cnp;
2628
2629   Standard_Integer I1, I2;
2630   Standard_Integer NPoleu , NPolev;
2631   gp_XYZ Temp;
2632
2633   for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
2634     Poles (NPoleu, LowerCol) =  Coefs (NPoleu, LowerCol);
2635     if (rat) {
2636       Weights (NPoleu, LowerCol) =  WCoefs (NPoleu, LowerCol);
2637     }
2638
2639     for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
2640       Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
2641       Temp = Coefs (NPoleu, Col).XYZ();
2642       Temp.Divide (Cnp);
2643       Poles (NPoleu, Col).SetXYZ (Temp);
2644       if (rat) {
2645         Weights (NPoleu, Col) = WCoefs (NPoleu, Col) /  Cnp;
2646       }
2647     }
2648     Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
2649     if (rat) {
2650       Weights (NPoleu, UpperCol) = WCoefs (NPoleu, UpperCol);
2651     }
2652
2653     for (I1 = 1; I1 <= RowLength - 1; I1++) {
2654
2655       for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
2656         Temp.SetLinearForm
2657           (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
2658         Poles (NPoleu, I2).SetXYZ (Temp);
2659         if (rat) Weights(NPoleu, I2) += Weights(NPoleu, I2-1);
2660       }
2661     }
2662   }
2663
2664   for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
2665
2666     for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
2667       Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
2668       Temp = Poles (Row, NPolev).XYZ();
2669       Temp.Divide (Cnp);
2670       Poles (Row, NPolev).SetXYZ (Temp);
2671       if (rat) Weights(Row, NPolev) /= Cnp;
2672     }
2673
2674     for (I1 = 1; I1 <= ColLength - 1; I1++) {
2675
2676       for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
2677         Temp.SetLinearForm
2678           (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
2679         Poles (I2, NPolev).SetXYZ (Temp);
2680         if (rat) Weights(I2, NPolev) += Weights(I2-1, NPolev);
2681       }
2682     }
2683   }
2684   if (rat) {
2685
2686     for (Row = LowerRow; Row <= UpperRow; Row++) {
2687
2688       for (Col = LowerCol; Col <= UpperCol; Col++) {
2689         W = Weights (Row, Col);
2690         Temp = Poles(Row, Col).XYZ();
2691         Temp.Divide (W);
2692         Poles(Row, Col).SetXYZ (Temp);
2693       }
2694     }
2695   }
2696 }
2697
2698 //=======================================================================
2699 //function : UTrimming
2700 //purpose  :
2701 //=======================================================================
2702
2703 void PLib::UTrimming(const Standard_Real U1,
2704                      const Standard_Real U2,
2705                      TColgp_Array2OfPnt& Coeffs,
2706                      TColStd_Array2OfReal& WCoeffs)
2707 {
2708   Standard_Boolean rat = &WCoeffs != NULL;
2709   Standard_Integer lr = Coeffs.LowerRow();
2710   Standard_Integer ur = Coeffs.UpperRow();
2711   Standard_Integer lc = Coeffs.LowerCol();
2712   Standard_Integer uc = Coeffs.UpperCol();
2713   TColgp_Array1OfPnt  Temp (lr,ur);
2714   TColStd_Array1OfReal Temw (lr,ur);
2715
2716   for (Standard_Integer icol = lc; icol <= uc; icol++) {
2717     Standard_Integer irow ;
2718     for ( irow = lr; irow <= ur; irow++) {
2719       Temp (irow) = Coeffs  (irow, icol);
2720       if (rat) Temw (irow) = WCoeffs (irow, icol);
2721     }
2722     if (rat) PLib::Trimming (U1, U2, Temp, Temw);
2723     else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
2724
2725     for (irow = lr; irow <= ur; irow++) {
2726       Coeffs  (irow, icol) = Temp (irow);
2727       if (rat) WCoeffs (irow, icol) = Temw (irow);
2728     }
2729   }
2730 }
2731
2732 //=======================================================================
2733 //function : VTrimming
2734 //purpose  :
2735 //=======================================================================
2736
2737 void PLib::VTrimming(const Standard_Real V1,
2738                      const Standard_Real V2,
2739                      TColgp_Array2OfPnt& Coeffs,
2740                      TColStd_Array2OfReal& WCoeffs)
2741 {
2742   Standard_Boolean rat = &WCoeffs != NULL;
2743   Standard_Integer lr = Coeffs.LowerRow();
2744   Standard_Integer ur = Coeffs.UpperRow();
2745   Standard_Integer lc = Coeffs.LowerCol();
2746   Standard_Integer uc = Coeffs.UpperCol();
2747   TColgp_Array1OfPnt  Temp (lc,uc);
2748   TColStd_Array1OfReal Temw (lc,uc);
2749
2750   for (Standard_Integer irow = lr; irow <= ur; irow++) {
2751     Standard_Integer icol ;
2752     for ( icol = lc; icol <= uc; icol++) {
2753       Temp (icol) = Coeffs  (irow, icol);
2754       if (rat) Temw (icol) = WCoeffs (irow, icol);
2755     }
2756     if (rat) PLib::Trimming (V1, V2, Temp, Temw);
2757     else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
2758
2759     for (icol = lc; icol <= uc; icol++) {
2760       Coeffs  (irow, icol) = Temp (icol);
2761       if (rat) WCoeffs (irow, icol) = Temw (icol);
2762     }
2763   }
2764 }
2765
2766 //=======================================================================
2767 //function : HermiteInterpolate
2768 //purpose  :
2769 //=======================================================================
2770
2771 Standard_Boolean PLib::HermiteInterpolate
2772 (const Standard_Integer Dimension,
2773  const Standard_Real FirstParameter,
2774  const Standard_Real LastParameter,
2775  const Standard_Integer FirstOrder,
2776  const Standard_Integer LastOrder,
2777  const TColStd_Array2OfReal& FirstConstr,
2778  const TColStd_Array2OfReal& LastConstr,
2779  TColStd_Array1OfReal& Coefficients)
2780 {
2781   Standard_Real Pattern[3][6];
2782
2783   // portage HP : il faut les initialiser 1 par 1
2784
2785   Pattern[0][0] = 1;
2786   Pattern[0][1] = 1;
2787   Pattern[0][2] = 1;
2788   Pattern[0][3] = 1;
2789   Pattern[0][4] = 1;
2790   Pattern[0][5] = 1;
2791   Pattern[1][0] = 0;
2792   Pattern[1][1] = 1;
2793   Pattern[1][2] = 2;
2794   Pattern[1][3] = 3;
2795   Pattern[1][4] = 4;
2796   Pattern[1][5] = 5;
2797   Pattern[2][0] = 0;
2798   Pattern[2][1] = 0;
2799   Pattern[2][2] = 2;
2800   Pattern[2][3] = 6;
2801   Pattern[2][4] = 12;
2802   Pattern[2][5] = 20;
2803
2804   math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
2805   //  The initialisation of the matrix A
2806   Standard_Integer irow ;
2807   for ( irow=0; irow<=FirstOrder; irow++) {
2808     Standard_Real FirstVal = 1.;
2809
2810     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2811       A(irow,icol) = Pattern[irow][icol]*FirstVal;
2812       if (irow <= icol) FirstVal *= FirstParameter;
2813     }
2814   }
2815
2816   for (irow=0; irow<=LastOrder; irow++) {
2817     Standard_Real LastVal  = 1.;
2818
2819     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2820       A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
2821       if (irow <= icol) LastVal *= LastParameter;
2822     }
2823   }
2824   //
2825   //  The filled matrix A for FirstOrder=LastOrder=2 is:
2826   //
2827   //      1   FP  FP**2   FP**3    FP**4     FP**5
2828   //      0   1   2*FP    3*FP**2  4*FP**3   5*FP**4        FP - FirstParameter
2829   //      0   0   2       6*FP     12*FP**2  20*FP**3
2830   //      1   LP  LP**2   LP**3    LP**4     LP**5
2831   //      0   1   2*LP    3*LP**2  4*LP**3   5*LP**4        LP - LastParameter
2832   //      0   0   2       6*LP     12*LP**2  20*LP**3
2833   //
2834   //  If FirstOrder or LastOrder <=2 then some rows and columns are missing.
2835   //  For example:
2836   //  If FirstOrder=1 then 3th row and 6th column are missing
2837   //  If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
2838
2839   math_Gauss Equations(A);
2840   //  cout << "A=" << A << endl;
2841
2842   for (Standard_Integer idim=1; idim<=Dimension; idim++) {
2843     //  cout << "idim=" << idim << endl;
2844
2845     math_Vector B(0,FirstOrder+LastOrder+1);
2846     Standard_Integer icol ;
2847     for ( icol=0; icol<=FirstOrder; icol++)
2848       B(icol) = FirstConstr(idim,icol);
2849
2850     for (icol=0; icol<=LastOrder; icol++)
2851       B(FirstOrder+1+icol) = LastConstr(idim,icol);
2852     //  cout << "B=" << B << endl;
2853
2854     //  The solving of equations system A * X = B. Then B = X
2855     Equations.Solve(B);
2856     //  cout << "After Solving" << endl << "B=" << B << endl;
2857
2858     if (Equations.IsDone()==Standard_False) return Standard_False;
2859
2860     //  the filling of the Coefficients
2861
2862     for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
2863       Coefficients(Dimension*icol+idim-1) = B(icol);
2864   }
2865   return Standard_True;
2866 }
2867
2868 //=======================================================================
2869 //function : JacobiParameters
2870 //purpose  :
2871 //=======================================================================
2872
2873 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
2874                             const Standard_Integer MaxDegree,
2875                             const Standard_Integer Code,
2876                             Standard_Integer& NbGaussPoints,
2877                             Standard_Integer& WorkDegree)
2878 {
2879   // ConstraintOrder: Ordre de contrainte aux extremites :
2880   //            C0 = contraintes de passage aux bornes;
2881   //            C1 = C0 + contraintes de derivees 1eres;
2882   //            C2 = C1 + contraintes de derivees 2ndes.
2883   // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
2884   //            d' approximation (doit etre superieur ou egal a
2885   //            2*NivConstr+2 et inferieur ou egal a 50).
2886   // Code:      Code d' init. des parametres de discretisation.
2887   //            (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
2888   //            = -5 Calcul tres rapide mais peu precis (8pts)
2889   //            = -4    '     '    '      '   '    '    (10pts)
2890   //            = -3    '     '    '      '   '    '    (15pts)
2891   //            = -2    '     '    '      '   '    '    (20pts)
2892   //            = -1    '     '    '      '   '    '    (25pts)
2893   //            = 1 calcul rapide avec precision moyenne (30pts).
2894   //            = 2 calcul rapide avec meilleure precision (40pts).
2895   //            = 3 calcul un peu plus lent avec bonne precision (50 pts).
2896   //            = 4 calcul lent avec la meilleure precision possible
2897   //             (61pts).
2898
2899   // The possible values of NbGaussPoints
2900
2901   const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
2902   NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
2903
2904   Standard_Integer NivConstr=0;
2905   switch (ConstraintOrder) {
2906   case GeomAbs_C0: NivConstr = 0; break;
2907   case GeomAbs_C1: NivConstr = 1; break;
2908   case GeomAbs_C2: NivConstr = 2; break;
2909   default:
2910     Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2911   }
2912   if (MaxDegree < 2*NivConstr+1)
2913     Standard_ConstructionError::Raise("Invalid MaxDegree");
2914
2915   if (Code >= 1)
2916     WorkDegree = MaxDegree + 9;
2917   else
2918     WorkDegree = MaxDegree + 6;
2919
2920   //---> Nbre mini de points necessaires.
2921   Standard_Integer IPMIN=0;
2922   if (WorkDegree < NDEG8)
2923     IPMIN=NDEG8;
2924   else if (WorkDegree < NDEG10)
2925     IPMIN=NDEG10;
2926   else if (WorkDegree < NDEG15)
2927     IPMIN=NDEG15;
2928   else if (WorkDegree < NDEG20)
2929     IPMIN=NDEG20;
2930   else if (WorkDegree < NDEG25)
2931     IPMIN=NDEG25;
2932   else if (WorkDegree < NDEG30)
2933     IPMIN=NDEG30;
2934   else if (WorkDegree < NDEG40)
2935     IPMIN=NDEG40;
2936   else if (WorkDegree < NDEG50)
2937     IPMIN=NDEG50;
2938   else if (WorkDegree < NDEG61)
2939     IPMIN=NDEG61;
2940   else
2941     Standard_ConstructionError::Raise("Invalid MaxDegree");
2942   // ---> Nbre de points voulus.
2943   Standard_Integer IWANT=0;
2944   switch (Code) {
2945   case -5: IWANT=NDEG8;  break;
2946   case -4: IWANT=NDEG10; break;
2947   case -3: IWANT=NDEG15; break;
2948   case -2: IWANT=NDEG20; break;
2949   case -1: IWANT=NDEG25; break;
2950   case  1: IWANT=NDEG30; break;
2951   case  2: IWANT=NDEG40; break;
2952   case  3: IWANT=NDEG50; break;
2953   case  4: IWANT=NDEG61; break;
2954   default:
2955     Standard_ConstructionError::Raise("Invalid Code");
2956   }
2957   //-->  NbGaussPoints est le nombre de points de discretisation de la fonction,
2958   //     il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
2959   //     NbGaussPoints doit etre superieur strictement a WorkDegree.
2960   NbGaussPoints = Max(IPMIN,IWANT);
2961   //  NbGaussPoints +=2;
2962 }
2963
2964 //=======================================================================
2965 //function : NivConstr
2966 //purpose  : translates from GeomAbs_Shape to Integer
2967 //=======================================================================
2968
2969  Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
2970 {
2971   Standard_Integer NivConstr=0;
2972   switch (ConstraintOrder) {
2973     case GeomAbs_C0: NivConstr = 0; break;
2974     case GeomAbs_C1: NivConstr = 1; break;
2975     case GeomAbs_C2: NivConstr = 2; break;
2976     default:
2977       Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2978   }
2979   return NivConstr;
2980 }
2981
2982 //=======================================================================
2983 //function : ConstraintOrder
2984 //purpose  : translates from Integer to GeomAbs_Shape
2985 //=======================================================================
2986
2987  GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
2988 {
2989   GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
2990   switch (NivConstr) {
2991     case 0: ConstraintOrder = GeomAbs_C0; break;
2992     case 1: ConstraintOrder = GeomAbs_C1; break;
2993     case 2: ConstraintOrder = GeomAbs_C2; break;
2994     default:
2995       Standard_ConstructionError::Raise("Invalid NivConstr");
2996   }
2997   return ConstraintOrder;
2998 }
2999
3000 //=======================================================================
3001 //function : EvalLength
3002 //purpose  :
3003 //=======================================================================
3004
3005  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
3006                       Standard_Real& PolynomialCoeff,
3007                       const Standard_Real U1, const Standard_Real U2,
3008                       Standard_Real& Length)
3009 {
3010   Standard_Integer i,j,idim, degdim;
3011   Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
3012
3013   Standard_Real *PolynomialArray = &PolynomialCoeff ;
3014
3015   Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
3016
3017   math_Vector GaussPoints(1,NbGaussPoints);
3018   math::GaussPoints(NbGaussPoints,GaussPoints);
3019
3020   math_Vector GaussWeights(1,NbGaussPoints);
3021   math::GaussWeights(NbGaussPoints,GaussWeights);
3022
3023   C1 = (U2 + U1) / 2.;
3024   C2 = (U2 - U1) / 2.;
3025
3026 //-----------------------------------------------------------
3027 //****** Integration - Boucle sur les intervalles de GAUSS **
3028 //-----------------------------------------------------------
3029
3030   Sum = 0;
3031
3032   for (j=1; j<=NbGaussPoints/2; j++) {
3033     // Integration en tenant compte de la symetrie
3034     Tran  = C2 * GaussPoints(j);
3035     X1 = C1 + Tran;
3036     X2 = C1 - Tran;
3037
3038     //****** Derivation sur la dimension de l'espace **
3039
3040     degdim =  Degree*Dimension;
3041     Der1 = Der2 = 0.;
3042     for (idim=0; idim<Dimension; idim++) {
3043       D1 = D2 = Degree * PolynomialArray [idim + degdim];
3044       for (i=Degree-1; i>=1; i--) {
3045         DD = i * PolynomialArray [idim + i*Dimension];
3046         D1 = D1 * X1 + DD;
3047         D2 = D2 * X2 + DD;
3048       }
3049       Der1 += D1 * D1;
3050       Der2 += D2 * D2;
3051     }
3052
3053     //****** Integration **
3054
3055     Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
3056
3057 //****** Fin de boucle dur les intervalles de GAUSS **
3058   }
3059   Length = Sum;
3060 }
3061
3062
3063 //=======================================================================
3064 //function : EvalLength
3065 //purpose  :
3066 //=======================================================================
3067
3068  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
3069                       Standard_Real& PolynomialCoeff,
3070                       const Standard_Real U1, const Standard_Real U2,
3071                       const Standard_Real Tol,
3072                       Standard_Real& Length, Standard_Real& Error)
3073 {
3074   Standard_Integer i;
3075   Standard_Integer NbSubInt = 1,    // Current number of subintervals
3076                    MaxNbIter = 13,  // Max number of iterations
3077                    NbIter    = 1;   // Current number of iterations
3078   Standard_Real    dU,OldLen,LenI;
3079
3080   PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
3081
3082   do {
3083     OldLen = Length;
3084     Length = 0.;
3085     NbSubInt *= 2;
3086     dU = (U2-U1)/NbSubInt;
3087     for (i=1; i<=NbSubInt; i++) {
3088       PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
3089       Length += LenI;
3090     }
3091     NbIter++;
3092     Error = Abs(OldLen-Length);
3093   }
3094   while (Error > Tol && NbIter <= MaxNbIter);
3095 }