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