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 #include <PLib.hxx>
  18
  19 #include <GeomAbs_Shape.hxx>
  20 #include <math.hxx>
  21 #include <math_Gauss.hxx>
  22 #include <math_Matrix.hxx>
  23 #include <NCollection_LocalArray.hxx>
  24 #include <Standard_ConstructionError.hxx>
  25
  26 // To convert points array into Real ..
  27 // *********************************
  28 //=======================================================================
  29 //function : SetPoles
  30 //purpose  :
  31 //=======================================================================
  32 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
  33                     TColStd_Array1OfReal& FP)
  34 {
  35   Standard_Integer j      = FP   .Lower();
  36   Standard_Integer PLower = Poles.Lower();
  37   Standard_Integer PUpper = Poles.Upper();
  38
  39   for (Standard_Integer i = PLower; i <= PUpper; i++) {
  40     const gp_Pnt2d& P = Poles(i);
  41     FP(j) = P.Coord(1); j++;
  42     FP(j) = P.Coord(2); j++;
  43   }
  44 }
  45
  46 //=======================================================================
  47 //function : SetPoles
  48 //purpose  :
  49 //=======================================================================
  50
  51 void PLib::SetPoles(const TColgp_Array1OfPnt2d&       Poles,
  52                     const TColStd_Array1OfReal& Weights,
  53                     TColStd_Array1OfReal&       FP)
  54 {
  55   Standard_Integer j      = FP   .Lower();
  56   Standard_Integer PLower = Poles.Lower();
  57   Standard_Integer PUpper = Poles.Upper();
  58
  59   for (Standard_Integer i = PLower; i <= PUpper; i++) {
  60     Standard_Real w = Weights(i);
  61     const gp_Pnt2d& P = Poles(i);
  62     FP(j) = P.Coord(1) * w; j++;
  63     FP(j) = P.Coord(2) * w; j++;
  64     FP(j) =              w; j++;
  65   }
  66 }
  67
  68 //=======================================================================
  69 //function : GetPoles
  70 //purpose  :
  71 //=======================================================================
  72
  73 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
  74                     TColgp_Array1OfPnt2d& Poles)
  75 {
  76   Standard_Integer j      = FP   .Lower();
  77   Standard_Integer PLower = Poles.Lower();
  78   Standard_Integer PUpper = Poles.Upper();
  79
  80   for (Standard_Integer i = PLower; i <= PUpper; i++) {
  81     gp_Pnt2d& P = Poles(i);
  82     P.SetCoord(1,FP(j)); j++;
  83     P.SetCoord(2,FP(j)); j++;
  84   }
  85 }
  86
  87 //=======================================================================
  88 //function : GetPoles
  89 //purpose  :
  90 //=======================================================================
  91
  92 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
  93                     TColgp_Array1OfPnt2d&       Poles,
  94                     TColStd_Array1OfReal&       Weights)
  95 {
  96   Standard_Integer j      = FP   .Lower();
  97   Standard_Integer PLower = Poles.Lower();
  98   Standard_Integer PUpper = Poles.Upper();
  99
 100   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 101     Standard_Real w = FP(j + 2);
 102     Weights(i) = w;
 103     gp_Pnt2d& P = Poles(i);
 104     P.SetCoord(1,FP(j) / w); j++;
 105     P.SetCoord(2,FP(j) / w); j++;
 106     j++;
 107   }
 108 }
 109
 110 //=======================================================================
 111 //function : SetPoles
 112 //purpose  :
 113 //=======================================================================
 114
 115 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
 116                     TColStd_Array1OfReal& FP)
 117 {
 118   Standard_Integer j      = FP   .Lower();
 119   Standard_Integer PLower = Poles.Lower();
 120   Standard_Integer PUpper = Poles.Upper();
 121
 122   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 123     const gp_Pnt& P = Poles(i);
 124     FP(j) = P.Coord(1); j++;
 125     FP(j) = P.Coord(2); j++;
 126     FP(j) = P.Coord(3); j++;
 127   }
 128 }
 129
 130 //=======================================================================
 131 //function : SetPoles
 132 //purpose  :
 133 //=======================================================================
 134
 135 void PLib::SetPoles(const TColgp_Array1OfPnt&   Poles,
 136                     const TColStd_Array1OfReal& Weights,
 137                     TColStd_Array1OfReal&       FP)
 138 {
 139   Standard_Integer j      = FP   .Lower();
 140   Standard_Integer PLower = Poles.Lower();
 141   Standard_Integer PUpper = Poles.Upper();
 142
 143   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 144     Standard_Real w = Weights(i);
 145     const gp_Pnt& P = Poles(i);
 146     FP(j) = P.Coord(1) * w; j++;
 147     FP(j) = P.Coord(2) * w; j++;
 148     FP(j) = P.Coord(3) * w; j++;
 149     FP(j) =              w; j++;
 150   }
 151 }
 152
 153 //=======================================================================
 154 //function : GetPoles
 155 //purpose  :
 156 //=======================================================================
 157
 158 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
 159                     TColgp_Array1OfPnt&         Poles)
 160 {
 161   Standard_Integer j      = FP   .Lower();
 162   Standard_Integer PLower = Poles.Lower();
 163   Standard_Integer PUpper = Poles.Upper();
 164
 165   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 166     gp_Pnt& P = Poles(i);
 167     P.SetCoord(1,FP(j)); j++;
 168     P.SetCoord(2,FP(j)); j++;
 169     P.SetCoord(3,FP(j)); j++;
 170   }
 171 }
 172
 173 //=======================================================================
 174 //function : GetPoles
 175 //purpose  :
 176 //=======================================================================
 177
 178 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
 179                     TColgp_Array1OfPnt&         Poles,
 180                     TColStd_Array1OfReal&       Weights)
 181 {
 182   Standard_Integer j      = FP   .Lower();
 183   Standard_Integer PLower = Poles.Lower();
 184   Standard_Integer PUpper = Poles.Upper();
 185
 186   for (Standard_Integer i = PLower; i <= PUpper; i++) {
 187     Standard_Real w = FP(j + 3);
 188     Weights(i) = w;
 189     gp_Pnt& P = Poles(i);
 190     P.SetCoord(1,FP(j) / w); j++;
 191     P.SetCoord(2,FP(j) / w); j++;
 192     P.SetCoord(3,FP(j) / w); j++;
 193     j++;
 194   }
 195 }
 196
 197 // specialized allocator
 198 namespace
 199 {
 200
 201 class BinomAllocator
 202 {
 203 public:
 204
 205   //! Main constructor
 206   BinomAllocator (const Standard_Integer theMaxBinom)
 207   : myBinom (NULL),
 208     myMaxBinom (theMaxBinom)
 209   {
 210     Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
 211     Standard_Integer np1 = myMaxBinom + 1;
 212     myBinom = new Standard_Integer*[np1];
 213     myBinom[0] = new Standard_Integer[1];
 214     myBinom[0][0] = 1;
 215     for (i = 1; i < np1; ++i)
 216     {
 217       im1 = i - 1;
 218       ip1 = i + 1;
 219       id2 = i >> 1;
 220       md2 = im1 >> 1;
 221       md3 = ip1 >> 1;
 222       k   = 0;
 223       myBinom[i] = new Standard_Integer[ip1];
 224
 225       for (j = 0; j < id2; ++j)
 226       {
 227         myBinom[i][j] = k + myBinom[im1][j];
 228         k = myBinom[im1][j];
 229       }
 230       j = id2;
 231       if (j > md2) j = im1 - j;
 232       myBinom[i][id2] = k + myBinom[im1][j];
 233
 234       for (j = ip1 - md3; j < ip1; j++)
 235       {
 236         myBinom[i][j] = myBinom[i][i - j];
 237       }
 238     }
 239   }
 240
 241   //! Destructor
 242   ~BinomAllocator()
 243   {
 244     // free memory
 245     for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
 246     {
 247       delete[] myBinom[i];
 248     }
 249     delete[] myBinom;
 250   }
 251
 252   Standard_Real Value (const Standard_Integer N,
 253                        const Standard_Integer P) const
 254   {
 255     Standard_OutOfRange_Raise_if (N > myMaxBinom,
 256       "PLib, BinomAllocator: requested degree is greater than maximum supported");
 257     return Standard_Real (myBinom[N][P]);
 258   }
 259
 260 private:
 261   BinomAllocator (const BinomAllocator&);
 262   BinomAllocator& operator= (const BinomAllocator&);
 263
 264 private:
 265   Standard_Integer**  myBinom;
 266   Standard_Integer myMaxBinom;
 267
 268 };
 269
 270   // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
 271   //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
 272   static BinomAllocator THE_BINOM (25 + 1);
 273 }
 274
 275 //=======================================================================
 276 //function : Bin
 277 //purpose  :
 278 //=======================================================================
 279
 280 Standard_Real PLib::Bin(const Standard_Integer N,
 281                         const Standard_Integer P)
 282 {
 283   return THE_BINOM.Value (N, P);
 284 }
 285
 286 //=======================================================================
 287 //function : RationalDerivative
 288 //purpose  :
 289 //=======================================================================
 290
 291 void  PLib::RationalDerivative(const Standard_Integer Degree,
 292                                const Standard_Integer DerivativeRequest,
 293                                const Standard_Integer Dimension,
 294                                Standard_Real& Ders,
 295                                Standard_Real& RDers,
 296                                const Standard_Boolean All)
 297 {
 298   //
 299   // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
 300   //
 301   //  We Write  u = fv
 302   //  Let C(N,P) be the binomial
 303   //
 304   //        then we have
 305   //
 306   //         (q)                             (p)   (q-p)
 307   //        u    =     SUM          C (q,p) f     v
 308   //                p = 0 to q
 309   //
 310   //
 311   //        Therefore
 312   //
 313   //
 314   //         (q)         (   (q)                               (p)   (q-p)   )
 315   //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
 316   //                     (          p = 0 to q-1                             )
 317   //
 318   //
 319   // make arrays for the binomial since computing it each time could raise a performance
 320   // issue
 321   // As oppose to the method below the <Der> array is organized in the following
 322   // fashion :
 323   //
 324   //  u (1)     u (2) ....   u  (Dimension)     v (1)
 325   //
 326   //   (1)       (1)          (1)                (1)
 327   //  u (1)     u (2) ....   u  (Dimension)     v (1)
 328   //
 329   //   ............................................
 330   //
 331   //   (Degree)  (Degree)     (Degree)           (Degree)
 332   //  u (1)     u (2) ....   u  (Dimension)     v (1)
 333   //
 334   //
 335   Standard_Real Inverse;
 336   Standard_Real *PolesArray = &Ders;
 337   Standard_Real *RationalArray = &RDers;
 338   Standard_Real Factor ;
 339   Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
 340   NCollection_LocalArray<Standard_Real> binomial_array;
 341   NCollection_LocalArray<Standard_Real> derivative_storage;
 342   if (Dimension == 3) {
 343     Standard_Integer DeRequest1 = DerivativeRequest + 1;
 344     Standard_Integer MinDegRequ = DerivativeRequest;
 345     if (MinDegRequ > Degree) MinDegRequ = Degree;
 346     binomial_array.Allocate (DeRequest1);
 347
 348     for (ii = 0 ; ii < DeRequest1 ; ii++) {
 349       binomial_array[ii] = 1.0e0 ;
 350     }
 351     if (!All) {
 352       Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
 353       derivative_storage.Allocate (DimDeRequ1);
 354       RationalArray = derivative_storage ;
 355     }
 356
 357     Inverse = 1.0e0 / PolesArray[3]  ;
 358     Index = 0 ;
 359     Index2 = - 6;
 360     OtherIndex = 0 ;
 361
 362     for (ii = 0 ; ii <= MinDegRequ ; ii++) {
 363       Index2 += 3;
 364       Index1  = Index2;
 365       RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
 366       RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
 367       RationalArray[Index] = PolesArray[OtherIndex];
 368       Index      -= 2;
 369       OtherIndex += 2;
 370
 371       for (jj = ii - 1 ; jj >= 0 ; jj--) {
 372         Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
 373         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 374         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 375         RationalArray[Index] -=  Factor * RationalArray[Index1];
 376         Index  -= 2;
 377         Index1 -= 5;
 378       }
 379
 380       for (jj = ii ; jj >=  1 ; jj--) {
 381         binomial_array[jj] += binomial_array[jj - 1] ;
 382       }
 383       RationalArray[Index] *= Inverse ; Index++;
 384       RationalArray[Index] *= Inverse ; Index++;
 385       RationalArray[Index] *= Inverse ; Index++;
 386     }
 387
 388     for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
 389       Index2 += 3;
 390       Index1 = Index2;
 391       RationalArray[Index] = 0.0e0; Index++;
 392       RationalArray[Index] = 0.0e0; Index++;
 393       RationalArray[Index] = 0.0e0;
 394       Index -= 2;
 395
 396       for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
 397         Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
 398         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 399         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 400         RationalArray[Index] -=  Factor * RationalArray[Index1];
 401         Index  -= 2;
 402         Index1 -= 5;
 403       }
 404
 405       for (jj = ii ; jj >=  1 ; jj--) {
 406         binomial_array[jj] += binomial_array[jj - 1] ;
 407       }
 408       RationalArray[Index] *= Inverse; Index++;
 409       RationalArray[Index] *= Inverse; Index++;
 410       RationalArray[Index] *= Inverse; Index++;
 411     }
 412
 413     if (!All) {
 414       RationalArray = &RDers ;
 415       Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
 416       RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
 417       RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
 418       RationalArray[2] = derivative_storage[DimDeRequ];
 419     }
 420   }
 421   else {
 422     Standard_Integer kk;
 423     Standard_Integer Dimension1 = Dimension + 1;
 424     Standard_Integer Dimension2 = Dimension << 1;
 425     Standard_Integer DeRequest1 = DerivativeRequest + 1;
 426     Standard_Integer MinDegRequ = DerivativeRequest;
 427     if (MinDegRequ > Degree) MinDegRequ = Degree;
 428     binomial_array.Allocate (DeRequest1);
 429
 430     for (ii = 0 ; ii < DeRequest1 ; ii++) {
 431       binomial_array[ii] = 1.0e0 ;
 432     }
 433     if (!All) {
 434       Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
 435       derivative_storage.Allocate (DimDeRequ1);
 436       RationalArray = derivative_storage ;
 437     }
 438
 439     Inverse = 1.0e0 / PolesArray[Dimension]  ;
 440     Index = 0 ;
 441     Index2 = - Dimension2;
 442     OtherIndex = 0 ;
 443
 444     for (ii = 0 ; ii <= MinDegRequ ; ii++) {
 445       Index2 += Dimension;
 446       Index1  = Index2;
 447
 448       for (kk = 0 ; kk < Dimension ; kk++) {
 449         RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
 450       }
 451       Index -= Dimension;
 452       ++OtherIndex;
 453
 454       for (jj = ii - 1 ; jj >= 0 ; jj--) {
 455         Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
 456
 457         for (kk = 0 ; kk < Dimension ; kk++) {
 458           RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 459         }
 460         Index  -= Dimension ;
 461         Index1 -= Dimension2 ;
 462       }
 463
 464       for (jj = ii ; jj >=  1 ; jj--) {
 465         binomial_array[jj] += binomial_array[jj - 1] ;
 466       }
 467
 468       for (kk = 0 ; kk < Dimension ; kk++) {
 469         RationalArray[Index] *= Inverse ; Index++;
 470       }
 471     }
 472
 473     for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
 474       Index2 += Dimension;
 475       Index1  = Index2;
 476
 477       for (kk = 0 ; kk < Dimension ; kk++) {
 478         RationalArray[Index] = 0.0e0 ; Index++;
 479       }
 480       Index -= Dimension;
 481
 482       for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
 483         Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
 484
 485         for (kk = 0 ; kk < Dimension ; kk++) {
 486           RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 487         }
 488         Index  -= Dimension ;
 489         Index1 -= Dimension2 ;
 490       }
 491
 492       for (jj = ii ; jj >=  1 ; jj--) {
 493         binomial_array[jj] += binomial_array[jj - 1] ;
 494       }
 495
 496       for (kk = 0 ; kk < Dimension ; kk++) {
 497         RationalArray[Index] *= Inverse; Index++;
 498       }
 499     }
 500
 501     if (!All) {
 502       RationalArray = &RDers ;
 503       Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
 504
 505       for (kk = 0 ; kk < Dimension ; kk++) {
 506         RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
 507       }
 508     }
 509   }
 510 }
 511
 512 //=======================================================================
 513 //function : RationalDerivatives
 514 //purpose  : Uses Homogeneous Poles Derivatives and Deivatives of Weights
 515 //=======================================================================
 516
 517 void  PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
 518                                 const Standard_Integer  Dimension,
 519                                 Standard_Real& PolesDerivates,
 520                                 // must be an array with
 521                                 // (DerivativeRequest + 1) * Dimension slots
 522                                 Standard_Real& WeightsDerivates,
 523                                 // must be an array with
 524                                 // (DerivativeRequest + 1) slots
 525                                 Standard_Real& RationalDerivates)
 526 {
 527   //
 528   // Our purpose is to compute f = (u/v) derivated N times
 529   //
 530   //  We Write  u = fv
 531   //  Let C(N,P) be the binomial
 532   //
 533   //        then we have
 534   //
 535   //         (q)                             (p)   (q-p)
 536   //        u    =     SUM          C (q,p) f     v
 537   //                p = 0 to q
 538   //
 539   //
 540   //        Therefore
 541   //
 542   //
 543   //         (q)         (   (q)                               (p)   (q-p)   )
 544   //        f    = (1/v) (  u    -     SUM            C (q,p) f     v        )
 545   //                     (          p = 0 to q-1                             )
 546   //
 547   //
 548   // make arrays for the binomial since computing it each time could
 549   // raize a performance issue
 550   //
 551   Standard_Real Inverse;
 552   Standard_Real *PolesArray = &PolesDerivates;
 553   Standard_Real *WeightsArray = &WeightsDerivates;
 554   Standard_Real *RationalArray = &RationalDerivates;
 555   Standard_Real Factor ;
 556
 557   Standard_Integer ii, Index, Index1, Index2, jj;
 558   Standard_Integer DeRequest1 = DerivativeRequest + 1;
 559
 560   NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
 561   NCollection_LocalArray<Standard_Real> derivative_storage;
 562
 563   for (ii = 0 ; ii < DeRequest1 ; ii++) {
 564     binomial_array[ii] = 1.0e0 ;
 565   }
 566   Inverse = 1.0e0 / WeightsArray[0]  ;
 567   if (Dimension == 3) {
 568     Index  = 0 ;
 569     Index2 = - 6 ;
 570
 571     for (ii = 0 ; ii < DeRequest1 ; ii++) {
 572       Index2 += 3;
 573       Index1 = Index2;
 574       RationalArray[Index] = PolesArray[Index] ; Index++;
 575       RationalArray[Index] = PolesArray[Index] ; Index++;
 576       RationalArray[Index] = PolesArray[Index] ;
 577       Index -= 2;
 578
 579       for (jj = ii - 1 ; jj >= 0 ; jj--) {
 580         Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
 581         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 582         RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 583         RationalArray[Index] -=  Factor * RationalArray[Index1];
 584         Index  -= 2;
 585         Index1 -= 5;
 586       }
 587
 588       for (jj = ii ; jj >=  1 ; jj--) {
 589         binomial_array[jj] += binomial_array[jj - 1] ;
 590       }
 591       RationalArray[Index] *= Inverse ; Index++;
 592       RationalArray[Index] *= Inverse ; Index++;
 593       RationalArray[Index] *= Inverse ; Index++;
 594     }
 595   }
 596   else {
 597     Standard_Integer kk;
 598     Standard_Integer Dimension2 = Dimension << 1;
 599     Index  = 0 ;
 600     Index2 = - Dimension2;
 601
 602     for (ii = 0 ; ii < DeRequest1 ; ii++) {
 603       Index2 += Dimension;
 604       Index1  = Index2;
 605
 606       for (kk = 0 ; kk < Dimension ; kk++) {
 607         RationalArray[Index] = PolesArray[Index]; Index++;
 608       }
 609       Index  -= Dimension;
 610
 611       for (jj = ii - 1 ; jj >= 0 ; jj--) {
 612         Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
 613
 614         for (kk = 0 ; kk < Dimension ; kk++) {
 615           RationalArray[Index] -=  Factor * RationalArray[Index1]; Index++; Index1++;
 616         }
 617         Index  -= Dimension;
 618         Index1 -= Dimension2;
 619       }
 620
 621       for (jj = ii ; jj >=  1 ; jj--) {
 622         binomial_array[jj] += binomial_array[jj - 1] ;
 623       }
 624
 625       for (kk = 0 ; kk < Dimension ; kk++) {
 626         RationalArray[Index] *= Inverse ; Index++;
 627       }
 628     }
 629   }
 630 }
 631
 632 //=======================================================================
 633 // Auxiliary template functions used for optimized evaluation of polynome
 634 // and its derivatives for smaller dimensions of the polynome
 635 //=======================================================================
 636
 637 namespace {
 638   // recursive template for evaluating value or first derivative
 639   template<int dim>
 640   inline void eval_step1 (double* poly, double par, double* coef)
 641   {
 642     eval_step1<dim - 1> (poly, par, coef);
 643     poly[dim] = poly[dim] * par + coef[dim];
 644   }
 645
 646   // recursion end
 647   template<>
 648   inline void eval_step1<-1> (double*, double, double*)
 649   {
 650   }
 651
 652   // recursive template for evaluating second derivative
 653   template<int dim>
 654   inline void eval_step2 (double* poly, double par, double* coef)
 655   {
 656     eval_step2<dim - 1> (poly, par, coef);
 657     poly[dim] = poly[dim] * par + coef[dim] * 2.;
 658   }
 659
 660   // recursion end
 661   template<>
 662   inline void eval_step2<-1> (double*, double, double*)
 663   {
 664   }
 665
 666   // evaluation of only value
 667   template<int dim>
 668   inline void eval_poly0 (double* aRes, double* aCoeffs, int Degree, double Par)
 669   {
 670     Standard_Real* aRes0 = aRes;
 671     memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
 672
 673     for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
 674     {
 675       aCoeffs -= dim;
 676       // Calculating the value of the polynomial
 677       eval_step1<dim-1> (aRes0, Par, aCoeffs);
 678     }
 679   }
 680
 681   // evaluation of value and first derivative
 682   template<int dim>
 683   inline void eval_poly1 (double* aRes, double* aCoeffs, int Degree, double Par)
 684   {
 685     Standard_Real* aRes0 = aRes;
 686     Standard_Real* aRes1 = aRes + dim;
 687
 688     memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
 689     memset(aRes1, 0, sizeof(Standard_Real) * dim);
 690
 691     for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
 692     {
 693       aCoeffs -= dim;
 694       // Calculating derivatives of the polynomial
 695       eval_step1<dim-1> (aRes1, Par, aRes0);
 696       // Calculating the value of the polynomial
 697       eval_step1<dim-1> (aRes0, Par, aCoeffs);
 698     }
 699   }
 700
 701   // evaluation of value and first and second derivatives
 702   template<int dim>
 703   inline void eval_poly2 (double* aRes, double* aCoeffs, int Degree, double Par)
 704   {
 705     Standard_Real* aRes0 = aRes;
 706     Standard_Real* aRes1 = aRes + dim;
 707     Standard_Real* aRes2 = aRes + 2 * dim;
 708
 709     memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
 710     memset(aRes1, 0, sizeof(Standard_Real) * dim);
 711     memset(aRes2, 0, sizeof(Standard_Real) * dim);
 712
 713     for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
 714     {
 715       aCoeffs -= dim;
 716       // Calculating second derivatives of the polynomial
 717       eval_step2<dim-1> (aRes2, Par, aRes1);
 718       // Calculating derivatives of the polynomial
 719       eval_step1<dim-1> (aRes1, Par, aRes0);
 720       // Calculating the value of the polynomial
 721       eval_step1<dim-1> (aRes0, Par, aCoeffs);
 722     }
 723   }
 724 }
 725
 726 //=======================================================================
 727 //function : This evaluates a polynomial and its derivatives
 728 //purpose  : up to the requested order
 729 //=======================================================================
 730
 731 void  PLib::EvalPolynomial(const Standard_Real    Par,
 732                            const Standard_Integer DerivativeRequest,
 733                            const Standard_Integer Degree,
 734                            const Standard_Integer Dimension,
 735                                  Standard_Real&   PolynomialCoeff,
 736                                  Standard_Real&   Results)
 737      //
 738      // the polynomial coefficients are assumed to be stored as follows :
 739      //                                                      0
 740      // [0]                  [Dimension -1]                 X     coefficient
 741      //                                                      1
 742      // [Dimension]          [Dimension + Dimension -1]     X     coefficient
 743      //                                                      2
 744      // [2 * Dimension]      [2 * Dimension + Dimension-1]  X     coefficient
 745      //
 746      //   ...................................................
 747      //
 748      //
 749      //                                                      d
 750      // [d * Dimension]      [d * Dimension + Dimension-1]  X     coefficient
 751      //
 752      //  where d is the Degree
 753      //
 754 {
 755   Standard_Real* aCoeffs = &PolynomialCoeff + Degree * Dimension;
 756   Standard_Real* aRes = &Results;
 757   Standard_Real* anOriginal;
 758   Standard_Integer ind = 0;
 759   switch (DerivativeRequest)
 760   {
 761   case 1:
 762   {
 763     switch (Dimension)
 764     {
 765     case 1:  eval_poly1<1>  (aRes, aCoeffs, Degree, Par); break;
 766     case 2:  eval_poly1<2>  (aRes, aCoeffs, Degree, Par); break;
 767     case 3:  eval_poly1<3>  (aRes, aCoeffs, Degree, Par); break;
 768     case 4:  eval_poly1<4>  (aRes, aCoeffs, Degree, Par); break;
 769     case 5:  eval_poly1<5>  (aRes, aCoeffs, Degree, Par); break;
 770     case 6:  eval_poly1<6>  (aRes, aCoeffs, Degree, Par); break;
 771     case 7:  eval_poly1<7>  (aRes, aCoeffs, Degree, Par); break;
 772     case 8:  eval_poly1<8>  (aRes, aCoeffs, Degree, Par); break;
 773     case 9:  eval_poly1<9>  (aRes, aCoeffs, Degree, Par); break;
 774     case 10: eval_poly1<10> (aRes, aCoeffs, Degree, Par); break;
 775     case 11: eval_poly1<11> (aRes, aCoeffs, Degree, Par); break;
 776     case 12: eval_poly1<12> (aRes, aCoeffs, Degree, Par); break;
 777     case 13: eval_poly1<13> (aRes, aCoeffs, Degree, Par); break;
 778     case 14: eval_poly1<14> (aRes, aCoeffs, Degree, Par); break;
 779     case 15: eval_poly1<15> (aRes, aCoeffs, Degree, Par); break;
 780     default:
 781       {
 782         Standard_Real* aRes0 = aRes;
 783         Standard_Real* aRes1 = aRes + Dimension;
 784
 785         memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * Dimension);
 786         memset(aRes1, 0, sizeof(Standard_Real) * Dimension);
 787
 788         for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
 789         {
 790           aCoeffs -= Dimension;
 791           // Calculating derivatives of the polynomial
 792           for (ind = 0; ind < Dimension; ind++)
 793             aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
 794           // Calculating the value of the polynomial
 795           for (ind = 0; ind < Dimension; ind++)
 796             aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
 797         }
 798       }
 799     }
 800     break;
 801   }
 802   case 2:
 803   {
 804     switch (Dimension)
 805     {
 806     case 1:  eval_poly2<1>  (aRes, aCoeffs, Degree, Par); break;
 807     case 2:  eval_poly2<2>  (aRes, aCoeffs, Degree, Par); break;
 808     case 3:  eval_poly2<3>  (aRes, aCoeffs, Degree, Par); break;
 809     case 4:  eval_poly2<4>  (aRes, aCoeffs, Degree, Par); break;
 810     case 5:  eval_poly2<5>  (aRes, aCoeffs, Degree, Par); break;
 811     case 6:  eval_poly2<6>  (aRes, aCoeffs, Degree, Par); break;
 812     case 7:  eval_poly2<7>  (aRes, aCoeffs, Degree, Par); break;
 813     case 8:  eval_poly2<8>  (aRes, aCoeffs, Degree, Par); break;
 814     case 9:  eval_poly2<9>  (aRes, aCoeffs, Degree, Par); break;
 815     case 10: eval_poly2<10> (aRes, aCoeffs, Degree, Par); break;
 816     case 11: eval_poly2<11> (aRes, aCoeffs, Degree, Par); break;
 817     case 12: eval_poly2<12> (aRes, aCoeffs, Degree, Par); break;
 818     case 13: eval_poly2<13> (aRes, aCoeffs, Degree, Par); break;
 819     case 14: eval_poly2<14> (aRes, aCoeffs, Degree, Par); break;
 820     case 15: eval_poly2<15> (aRes, aCoeffs, Degree, Par); break;
 821     default:
 822       {
 823         Standard_Real* aRes0 = aRes;
 824         Standard_Real* aRes1 = aRes + Dimension;
 825         Standard_Real* aRes2 = aRes1 + Dimension;
 826
 827         // Nullify the results
 828         Standard_Integer aSize = 2 * Dimension;
 829         memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
 830         memset(aRes1, 0, sizeof(Standard_Real) * aSize);
 831
 832         for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
 833         {
 834           aCoeffs -= Dimension;
 835           // Calculating derivatives of the polynomial
 836           for (ind = 0; ind < Dimension; ind++)
 837             aRes2[ind] = aRes2[ind] * Par + aRes1[ind] * 2.0;
 838           for (ind = 0; ind < Dimension; ind++)
 839             aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
 840           // Calculating the value of the polynomial
 841           for (ind = 0; ind < Dimension; ind++)
 842             aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
 843         }
 844         break;
 845       }
 846     }
 847     break;
 848   }
 849   default:
 850     {
 851       // Nullify the results
 852       Standard_Integer aResSize = (1 + DerivativeRequest) * Dimension;
 853       memset(aRes, 0, sizeof(Standard_Real) * aResSize);
 854
 855       for (Standard_Integer aDeg = 0; aDeg <= Degree; aDeg++)
 856       {
 857         aRes = &Results + aResSize - Dimension;
 858         // Calculating derivatives of the polynomial
 859         for (Standard_Integer aDeriv = DerivativeRequest; aDeriv > 0; aDeriv--)
 860         {
 861           anOriginal = aRes - Dimension; // pointer to the derivative minus 1
 862           for (ind = 0; ind < Dimension; ind++)
 863             aRes[ind] = aRes[ind] * Par + anOriginal[ind] * aDeriv;
 864           aRes = anOriginal;
 865         }
 866         // Calculating the value of the polynomial
 867         for (ind = 0; ind < Dimension; ind++)
 868           aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
 869         aCoeffs -= Dimension;
 870       }
 871     }
 872   }
 873 }
 874
 875 //=======================================================================
 876 //function : This evaluates a polynomial without derivative
 877 //purpose  :
 878 //=======================================================================
 879
 880 void  PLib::NoDerivativeEvalPolynomial(const Standard_Real    Par,
 881                                        const Standard_Integer Degree,
 882                                        const Standard_Integer Dimension,
 883                                        const Standard_Integer DegreeDimension,
 884                                        Standard_Real&         PolynomialCoeff,
 885                                        Standard_Real&         Results)
 886 {
 887   Standard_Real* aCoeffs = &PolynomialCoeff + DegreeDimension;
 888   Standard_Real* aRes = &Results;
 889
 890   switch (Dimension)
 891   {
 892   case 1:  eval_poly0<1>  (aRes, aCoeffs, Degree, Par); break;
 893   case 2:  eval_poly0<2>  (aRes, aCoeffs, Degree, Par); break;
 894   case 3:  eval_poly0<3>  (aRes, aCoeffs, Degree, Par); break;
 895   case 4:  eval_poly0<4>  (aRes, aCoeffs, Degree, Par); break;
 896   case 5:  eval_poly0<5>  (aRes, aCoeffs, Degree, Par); break;
 897   case 6:  eval_poly0<6>  (aRes, aCoeffs, Degree, Par); break;
 898   case 7:  eval_poly0<7>  (aRes, aCoeffs, Degree, Par); break;
 899   case 8:  eval_poly0<8>  (aRes, aCoeffs, Degree, Par); break;
 900   case 9:  eval_poly0<9>  (aRes, aCoeffs, Degree, Par); break;
 901   case 10: eval_poly0<10> (aRes, aCoeffs, Degree, Par); break;
 902   case 11: eval_poly0<11> (aRes, aCoeffs, Degree, Par); break;
 903   case 12: eval_poly0<12> (aRes, aCoeffs, Degree, Par); break;
 904   case 13: eval_poly0<13> (aRes, aCoeffs, Degree, Par); break;
 905   case 14: eval_poly0<14> (aRes, aCoeffs, Degree, Par); break;
 906   case 15: eval_poly0<15> (aRes, aCoeffs, Degree, Par); break;
 907   default:
 908     {
 909       memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
 910       for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
 911       {
 912         aCoeffs -= Dimension;
 913         for (Standard_Integer ind = 0; ind < Dimension; ind++)
 914           aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
 915       }
 916     }
 917   }
 918 }
 919
 920 //=======================================================================
 921 //function : This evaluates a polynomial of 2 variables
 922 //purpose  : or its derivative at the requested orders
 923 //=======================================================================
 924
 925 void  PLib::EvalPoly2Var(const Standard_Real    UParameter,
 926                          const Standard_Real    VParameter,
 927                          const Standard_Integer UDerivativeRequest,
 928                          const Standard_Integer VDerivativeRequest,
 929                          const Standard_Integer UDegree,
 930                          const Standard_Integer VDegree,
 931                          const Standard_Integer Dimension,
 932                          Standard_Real&         PolynomialCoeff,
 933                          Standard_Real&         Results)
 934      //
 935      // the polynomial coefficients are assumed to be stored as follows :
 936      //                                                      0 0
 937      // [0]                  [Dimension -1]                 U V    coefficient
 938      //                                                      1 0
 939      // [Dimension]          [Dimension + Dimension -1]     U V    coefficient
 940      //                                                      2 0
 941      // [2 * Dimension]      [2 * Dimension + Dimension-1]  U V    coefficient
 942      //
 943      //   ...................................................
 944      //
 945      //
 946      //                                                      m 0
 947      // [m * Dimension]      [m * Dimension + Dimension-1]  U V    coefficient
 948      //
 949      //  where m = UDegree
 950      //
 951      //                                                           0 1
 952      // [(m+1) * Dimension]  [(m+1) * Dimension + Dimension-1]   U V   coefficient
 953      //
 954      //   ...................................................
 955      //
 956      //                                                          m 1
 957      // [2*m * Dimension]      [2*m * Dimension + Dimension-1]  U V    coefficient
 958      //
 959      //   ...................................................
 960      //
 961      //                                                          m n
 962      // [m*n * Dimension]      [m*n * Dimension + Dimension-1]  U V    coefficient
 963      //
 964      //  where n = VDegree
 965 {
 966   Standard_Integer Udim = (VDegree+1)*Dimension,
 967   index = Udim*UDerivativeRequest;
 968   TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
 969   TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
 970   Standard_Real * Result =  (Standard_Real *) &Curve.ChangeValue(1);
 971   Standard_Real * Digit  =  (Standard_Real *) &Point.ChangeValue(1);
 972   Standard_Real * ResultArray ;
 973   ResultArray = &Results ;
 974
 975   PLib::EvalPolynomial(UParameter,
 976                        UDerivativeRequest,
 977                        UDegree,
 978                        Udim,
 979                        PolynomialCoeff,
 980                        Result[0]);
 981
 982   PLib::EvalPolynomial(VParameter,
 983                        VDerivativeRequest,
 984                        VDegree,
 985                        Dimension,
 986                        Result[index],
 987                        Digit[0]);
 988
 989   index = Dimension*VDerivativeRequest;
 990
 991   for (Standard_Integer i=0;i<Dimension;i++) {
 992     ResultArray[i] = Digit[index+i];
 993   }
 994 }
 995
 996
 997
 998 //=======================================================================
 999 //function : This evaluates the lagrange polynomial and its derivatives
1000 //purpose  : up to the requested order that interpolates a series of
1001 //points of dimension <Dimension> with given assigned parameters
1002 //=======================================================================
1003
1004 Standard_Integer
1005 PLib::EvalLagrange(const Standard_Real                   Parameter,
1006                    const Standard_Integer                DerivativeRequest,
1007                    const Standard_Integer                Degree,
1008                    const Standard_Integer                Dimension,
1009                    Standard_Real&                        Values,
1010                    Standard_Real&                        Parameters,
1011                    Standard_Real&                        Results)
1012 {
1013   //
1014   // the points  are assumed to be stored as follows in the Values array :
1015   //
1016   // [0]              [Dimension -1]                first  point   coefficients
1017   //
1018   // [Dimension]      [Dimension + Dimension -1]    second point   coefficients
1019   //
1020   // [2 * Dimension]  [2 * Dimension + Dimension-1] third  point   coefficients
1021   //
1022   //   ...................................................
1023   //
1024   //
1025   //
1026   // [d * Dimension]  [d * Dimension + Dimension-1] d + 1 point   coefficients
1027   //
1028   //  where d is the Degree
1029   //
1030   //  The ParameterArray stores the parameter value assign to each point in
1031   //  order described above, that is
1032   //  [0]   is assign to first  point
1033   //  [1]   is assign to second point
1034   //
1035   Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
1036   Standard_Integer local_request = DerivativeRequest;
1037   Standard_Real  *ParameterArray;
1038   Standard_Real  difference;
1039   Standard_Real  *PointsArray;
1040   Standard_Real  *ResultArray ;
1041
1042   PointsArray    = &Values ;
1043   ParameterArray = &Parameters ;
1044   ResultArray = &Results ;
1045   if (local_request >= Degree) {
1046     local_request = Degree ;
1047   }
1048   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1049   //
1050   //  Build the divided differences array
1051   //
1052
1053   for (ii = 0 ; ii <  (Degree + 1) * Dimension  ; ii++) {
1054     divided_differences_array[ii] = PointsArray[ii] ;
1055   }
1056
1057   for (ii = Degree ; ii >= 0   ; ii--) {
1058
1059     for (jj = Degree  ; jj > Degree - ii  ; jj--) {
1060       Index = jj * Dimension ;
1061       Index1 = Index - Dimension ;
1062
1063       for (kk = 0 ; kk < Dimension ; kk++) {
1064         divided_differences_array[Index + kk] -=
1065           divided_differences_array[Index1 + kk] ;
1066       }
1067       difference =
1068         ParameterArray[jj] - ParameterArray[jj - Degree -1 +  ii] ;
1069       if (Abs(difference) < RealSmall()) {
1070         ReturnCode = 1 ;
1071         goto FINISH ;
1072       }
1073       difference = 1.0e0 / difference ;
1074
1075       for (kk = 0 ; kk < Dimension ; kk++) {
1076         divided_differences_array[Index + kk] *= difference ;
1077       }
1078     }
1079   }
1080   //
1081   //
1082   // Evaluate the divided difference array polynomial which expresses as
1083   //
1084   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1085   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1086   //
1087   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1088   //
1089   //
1090   Index = Degree * Dimension ;
1091
1092   for (kk = 0 ; kk < Dimension ; kk++) {
1093     ResultArray[kk] = divided_differences_array[Index + kk] ;
1094   }
1095
1096   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
1097     ResultArray[ii] = 0.0e0 ;
1098   }
1099
1100   for (ii = Degree ; ii >= 1 ; ii--) {
1101     difference =  Parameter - ParameterArray[ii - 1] ;
1102
1103     for (jj = local_request ; jj > 0 ; jj--) {
1104       Index = jj * Dimension ;
1105       Index1 = Index - Dimension ;
1106
1107       for (kk = 0 ; kk < Dimension ; kk++) {
1108         ResultArray[Index + kk] *= difference ;
1109         ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
1110       }
1111     }
1112     Index = (ii -1) * Dimension ;
1113
1114     for (kk = 0 ; kk < Dimension ; kk++) {
1115       ResultArray[kk] *= difference ;
1116       ResultArray[kk] += divided_differences_array[Index+kk] ;
1117     }
1118   }
1119   FINISH :
1120     return (ReturnCode) ;
1121 }
1122
1123 //=======================================================================
1124 //function : This evaluates the hermite polynomial and its derivatives
1125 //purpose  : up to the requested order that interpolates a series of
1126 //points of dimension <Dimension> with given assigned parameters
1127 //=======================================================================
1128
1129 Standard_Integer PLib::EvalCubicHermite
1130 (const Standard_Real          Parameter,
1131  const Standard_Integer       DerivativeRequest,
1132  const Standard_Integer       Dimension,
1133  Standard_Real&               Values,
1134  Standard_Real&               Derivatives,
1135  Standard_Real&               theParameters,
1136  Standard_Real&               Results)
1137 {
1138   //
1139   // the points  are assumed to be stored as follows in the Values array :
1140   //
1141   // [0]            [Dimension -1]             first  point   coefficients
1142   //
1143   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
1144   //
1145   //
1146   // the derivatives  are assumed to be stored as follows in
1147   // the Derivatives array :
1148   //
1149   // [0]            [Dimension -1]             first  point   coefficients
1150   //
1151   // [Dimension]    [Dimension + Dimension -1] last point   coefficients
1152   //
1153   //  The ParameterArray stores the parameter value assign to each point in
1154   //  order described above, that is
1155   //  [0]   is assign to first  point
1156   //  [1]   is assign to last   point
1157   //
1158   Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
1159   Standard_Integer local_request = DerivativeRequest ;
1160
1161   ReturnCode = 0 ;
1162   Degree = 3 ;
1163   Standard_Real  ParametersArray[4];
1164   Standard_Real  difference;
1165   Standard_Real  inverse;
1166   Standard_Real  *FirstLast;
1167   Standard_Real  *PointsArray;
1168   Standard_Real  *DerivativesArray;
1169   Standard_Real  *ResultArray ;
1170
1171   DerivativesArray = &Derivatives ;
1172   PointsArray    = &Values ;
1173   FirstLast = &theParameters ;
1174   ResultArray = &Results ;
1175   if (local_request >= Degree) {
1176     local_request = Degree ;
1177   }
1178   NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1179
1180   for (ii = 0, jj = 0  ; ii < 2 ; ii++, jj+= 2) {
1181     ParametersArray[jj] =
1182       ParametersArray[jj+1] = FirstLast[ii] ;
1183   }
1184   //
1185   //  Build the divided differences array
1186   //
1187   //
1188   // initialise it at the stage 2 of the building algorithm
1189   // for divided differences
1190   //
1191   inverse = FirstLast[1] - FirstLast[0] ;
1192   inverse = 1.0e0 / inverse ;
1193
1194   for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension  ;
1195        ii <  Dimension  ;
1196        ii++, jj++, kk++, pp++) {
1197     divided_differences_array[ii] = PointsArray[ii] ;
1198     divided_differences_array[kk] = inverse *
1199       (PointsArray[jj] - PointsArray[ii]) ;
1200     divided_differences_array[jj] = DerivativesArray[ii] ;
1201     divided_differences_array[pp] = DerivativesArray[jj] ;
1202   }
1203
1204   for (ii = 1 ; ii <= Degree   ; ii++) {
1205
1206     for (jj = Degree  ; jj >=  ii+1  ; jj--) {
1207       Index = jj * Dimension ;
1208       Index1 = Index - Dimension ;
1209
1210       for (kk = 0 ; kk < Dimension ; kk++) {
1211         divided_differences_array[Index + kk] -=
1212           divided_differences_array[Index1 + kk] ;
1213       }
1214
1215       for (kk = 0 ; kk < Dimension ; kk++) {
1216         divided_differences_array[Index + kk] *= inverse ;
1217       }
1218     }
1219   }
1220   //
1221   //
1222   // Evaluate the divided difference array polynomial which expresses as
1223   //
1224   //  P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1225   //         + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1226   //
1227   // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1228   //
1229   //
1230   Index = Degree * Dimension ;
1231
1232   for (kk = 0 ; kk < Dimension ; kk++) {
1233     ResultArray[kk] = divided_differences_array[Index + kk] ;
1234   }
1235
1236   for (ii = Dimension ; ii < (local_request + 1)  * Dimension ; ii++) {
1237     ResultArray[ii] = 0.0e0 ;
1238   }
1239
1240   for (ii = Degree ; ii >= 1 ; ii--) {
1241     difference =  Parameter - ParametersArray[ii - 1] ;
1242
1243     for (jj = local_request ; jj > 0 ; jj--) {
1244       Index = jj * Dimension ;
1245       Index1 = Index - Dimension ;
1246
1247       for (kk = 0 ; kk < Dimension ; kk++) {
1248         ResultArray[Index + kk] *= difference ;
1249         ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
1250       }
1251     }
1252     Index = (ii -1) * Dimension ;
1253
1254     for (kk = 0 ; kk < Dimension ; kk++) {
1255       ResultArray[kk] *= difference ;
1256       ResultArray[kk] += divided_differences_array[Index+kk] ;
1257     }
1258   }
1259 //  FINISH :
1260     return (ReturnCode) ;
1261 }
1262
1263 //=======================================================================
1264 //function : HermiteCoefficients
1265 //purpose  : calcul des polynomes d'Hermite
1266 //=======================================================================
1267
1268 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
1269                                            const Standard_Real LastParameter,
1270                                            const Standard_Integer FirstOrder,
1271                                            const Standard_Integer LastOrder,
1272                                            math_Matrix& MatrixCoefs)
1273 {
1274   Standard_Integer NbCoeff =  FirstOrder +  LastOrder + 2, Ordre[2];
1275   Standard_Integer ii, jj, pp, cote, iof=0;
1276   Standard_Real Prod, TBorne = FirstParameter;
1277   math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
1278   math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
1279
1280   // Test de validites
1281
1282   if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
1283   Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
1284   if (D1 > 100 || D2 > 100) return Standard_False;
1285   D2 += D1;
1286   if (D2 < 0.01) return Standard_False;
1287   if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
1288
1289   // Calcul de la matrice a inverser (MAT)
1290
1291   Ordre[0] = FirstOrder+1;
1292   Ordre[1] = LastOrder+1;
1293
1294   for (cote=0; cote<=1; cote++) {
1295     Coeff.Init(1);
1296
1297     for (pp=1; pp<=Ordre[cote]; pp++) {
1298       ii = pp + iof;
1299       Prod = 1;
1300
1301       for (jj=pp; jj<=NbCoeff; jj++) {
1302         //        tout se passe dans les 3 lignes suivantes
1303         MAT(ii, jj) = Coeff(jj) * Prod;
1304         Coeff(jj) *= jj - pp;
1305         Prod      *= TBorne;
1306       }
1307     }
1308     TBorne = LastParameter;
1309     iof = Ordre[0];
1310   }
1311
1312   // resolution du systemes
1313   math_Gauss ResolCoeff(MAT, 1.0e-10);
1314   if (!ResolCoeff.IsDone()) return Standard_False;
1315
1316   for (ii=1; ii<=NbCoeff; ii++) {
1317     B(ii) = 1;
1318     ResolCoeff.Solve(B, Coeff);
1319     MatrixCoefs.SetRow( ii, Coeff);
1320     B(ii) = 0;
1321   }
1322   return Standard_True;
1323 }
1324
1325 //=======================================================================
1326 //function : CoefficientsPoles
1327 //purpose  :
1328 //=======================================================================
1329
1330 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt&   Coefs,
1331                               const TColStd_Array1OfReal* WCoefs,
1332                               TColgp_Array1OfPnt&         Poles,
1333                               TColStd_Array1OfReal*       Weights)
1334 {
1335   TColStd_Array1OfReal tempC(1,3*Coefs.Length());
1336   PLib::SetPoles(Coefs,tempC);
1337   TColStd_Array1OfReal tempP(1,3*Poles.Length());
1338   PLib::SetPoles(Coefs,tempP);
1339   PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
1340   PLib::GetPoles(tempP,Poles);
1341 }
1342
1343 //=======================================================================
1344 //function : CoefficientsPoles
1345 //purpose  :
1346 //=======================================================================
1347
1348 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
1349                               const TColStd_Array1OfReal* WCoefs,
1350                               TColgp_Array1OfPnt2d&       Poles,
1351                               TColStd_Array1OfReal*       Weights)
1352 {
1353   TColStd_Array1OfReal tempC(1,2*Coefs.Length());
1354   PLib::SetPoles(Coefs,tempC);
1355   TColStd_Array1OfReal tempP(1,2*Poles.Length());
1356   PLib::SetPoles(Coefs,tempP);
1357   PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
1358   PLib::GetPoles(tempP,Poles);
1359 }
1360
1361 //=======================================================================
1362 //function : CoefficientsPoles
1363 //purpose  :
1364 //=======================================================================
1365
1366 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
1367                               const TColStd_Array1OfReal* WCoefs,
1368                               TColStd_Array1OfReal&       Poles,
1369                               TColStd_Array1OfReal*       Weights)
1370 {
1371   PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
1372 }
1373
1374 //=======================================================================
1375 //function : CoefficientsPoles
1376 //purpose  :
1377 //=======================================================================
1378
1379 void PLib::CoefficientsPoles (const Standard_Integer      dim,
1380                               const TColStd_Array1OfReal& Coefs,
1381                               const TColStd_Array1OfReal* WCoefs,
1382                               TColStd_Array1OfReal&       Poles,
1383                               TColStd_Array1OfReal*       Weights)
1384 {
1385   Standard_Boolean rat = WCoefs != NULL;
1386   Standard_Integer loc = Coefs.Lower();
1387   Standard_Integer lop = Poles.Lower();
1388   Standard_Integer lowc=0;
1389   Standard_Integer lowp=0;
1390   Standard_Integer upc = Coefs.Upper();
1391   Standard_Integer upp = Poles.Upper();
1392   Standard_Integer upwc=0;
1393   Standard_Integer upwp=0;
1394   Standard_Integer reflen = Coefs.Length()/dim;
1395   Standard_Integer i,j,k;
1396   //Les Extremites.
1397   if (rat) {
1398     lowc = WCoefs->Lower(); lowp = Weights->Lower();
1399     upwc = WCoefs->Upper(); upwp = Weights->Upper();
1400   }
1401
1402   for (i = 0; i < dim; i++){
1403     Poles (lop + i) = Coefs (loc + i);
1404     Poles (upp - i) = Coefs (upc - i);
1405   }
1406   if (rat) {
1407     (*Weights) (lowp) = (*WCoefs) (lowc);
1408     (*Weights) (upwp) = (*WCoefs) (upwc);
1409   }
1410
1411   Standard_Real Cnp;
1412   for (i = 2; i < reflen; i++ ) {
1413     Cnp = PLib::Bin(reflen - 1, i - 1);
1414     if (rat) (*Weights)(lowp + i - 1) = (*WCoefs)(lowc + i - 1) / Cnp;
1415
1416     for(j = 0; j < dim; j++){
1417       Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
1418     }
1419   }
1420
1421   for (i = 1; i <= reflen - 1; i++) {
1422
1423     for (j = reflen - 1; j >= i; j--) {
1424       if (rat) (*Weights)(lowp + j) += (*Weights)(lowp + j - 1);
1425
1426       for(k = 0; k < dim; k++){
1427         Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
1428       }
1429     }
1430   }
1431   if (rat) {
1432
1433     for (i = 1; i <= reflen; i++) {
1434
1435       for(j = 0; j < dim; j++){
1436         Poles(lop + dim * (i-1) + j) /= (*Weights)(lowp + i -1);
1437       }
1438     }
1439   }
1440 }
1441
1442 //=======================================================================
1443 //function : Trimming
1444 //purpose  :
1445 //=======================================================================
1446
1447 void PLib::Trimming(const Standard_Real   U1,
1448                     const Standard_Real   U2,
1449                     TColgp_Array1OfPnt&   Coefs,
1450                     TColStd_Array1OfReal* WCoefs)
1451 {
1452   TColStd_Array1OfReal temp(1,3*Coefs.Length());
1453   PLib::SetPoles(Coefs,temp);
1454   PLib::Trimming(U1,U2,3,temp,WCoefs);
1455   PLib::GetPoles(temp,Coefs);
1456 }
1457
1458 //=======================================================================
1459 //function : Trimming
1460 //purpose  :
1461 //=======================================================================
1462
1463 void PLib::Trimming(const Standard_Real   U1,
1464                     const Standard_Real   U2,
1465                     TColgp_Array1OfPnt2d& Coefs,
1466                     TColStd_Array1OfReal* WCoefs)
1467 {
1468   TColStd_Array1OfReal temp(1,2*Coefs.Length());
1469   PLib::SetPoles(Coefs,temp);
1470   PLib::Trimming(U1,U2,2,temp,WCoefs);
1471   PLib::GetPoles(temp,Coefs);
1472 }
1473
1474 //=======================================================================
1475 //function : Trimming
1476 //purpose  :
1477 //=======================================================================
1478
1479 void PLib::Trimming(const Standard_Real   U1,
1480                     const Standard_Real   U2,
1481                     TColStd_Array1OfReal& Coefs,
1482                     TColStd_Array1OfReal* WCoefs)
1483 {
1484   PLib::Trimming(U1,U2,1,Coefs,WCoefs);
1485 }
1486
1487 //=======================================================================
1488 //function : Trimming
1489 //purpose  :
1490 //=======================================================================
1491
1492 void PLib::Trimming(const Standard_Real U1,
1493                     const Standard_Real U2,
1494                     const Standard_Integer dim,
1495                     TColStd_Array1OfReal& Coefs,
1496                     TColStd_Array1OfReal* WCoefs)
1497 {
1498
1499   // principe :
1500   // on fait le changement de variable v = (u-U1) / (U2-U1)
1501   // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
1502   // decomposee sous la forme du schema iteratif de horner.
1503
1504   Standard_Real lsp = U2 - U1;
1505   Standard_Integer indc, indw=0;
1506   Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
1507   Standard_Integer len = Coefs.Length()/dim;
1508   Standard_Boolean rat = WCoefs != NULL;
1509
1510   if (rat) {
1511     if(len != WCoefs->Length())
1512       throw Standard_Failure("PLib::Trimming : nbcoefs/dim != nbweights !!!");
1513     upw = WCoefs->Upper();
1514   }
1515   len --;
1516
1517   for (Standard_Integer i = 1; i <= len; i++) {
1518     Standard_Integer j ;
1519     indc = upc - dim*(i-1);
1520     if (rat) indw = upw - i + 1;
1521     //calcul du coefficient de degre le plus faible a l'iteration i
1522
1523     for( j = 0; j < dim; j++){
1524       Coefs(indc - dim + j) += U1 * Coefs(indc + j);
1525     }
1526     if (rat) (*WCoefs)(indw - 1) += U1 * (*WCoefs)(indw);
1527
1528     //calcul des coefficients intermediaires :
1529
1530     while (indc < upc){
1531       indc += dim;
1532
1533       for(Standard_Integer k = 0; k < dim; k++){
1534         Coefs(indc - dim + k) =
1535           U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
1536       }
1537       if (rat) {
1538         indw ++;
1539         (*WCoefs)(indw - 1) = U1 * (*WCoefs)(indw) + lsp * (*WCoefs)(indw - 1);
1540       }
1541     }
1542
1543     //calcul du coefficient de degre le plus eleve :
1544
1545     for(j = 0; j < dim; j++){
1546       Coefs(upc + j) *= lsp;
1547     }
1548     if (rat) (*WCoefs)(upw) *= lsp;
1549   }
1550 }
1551
1552 //=======================================================================
1553 //function : CoefficientsPoles
1554 //purpose  :
1555 // Modified: 21/10/1996 by PMN :  PolesCoefficient (PRO5852).
1556 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
1557 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
1558 //=======================================================================
1559
1560 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt&   Coefs,
1561                               const TColStd_Array2OfReal* WCoefs,
1562                               TColgp_Array2OfPnt&         Poles,
1563                               TColStd_Array2OfReal*       Weights)
1564 {
1565   Standard_Boolean rat = (WCoefs != NULL);
1566   Standard_Integer LowerRow  = Poles.LowerRow();
1567   Standard_Integer UpperRow  = Poles.UpperRow();
1568   Standard_Integer LowerCol  = Poles.LowerCol();
1569   Standard_Integer UpperCol  = Poles.UpperCol();
1570   Standard_Integer ColLength = Poles.ColLength();
1571   Standard_Integer RowLength = Poles.RowLength();
1572
1573   // Bidouille pour retablir u et v pour les coefs calcules
1574   // par buildcache
1575 //  Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
1576
1577   Standard_Integer Row, Col;
1578   Standard_Real W, Cnp;
1579
1580   Standard_Integer I1, I2;
1581   Standard_Integer NPoleu , NPolev;
1582   gp_XYZ Temp;
1583
1584   for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
1585     Poles (NPoleu, LowerCol) =  Coefs (NPoleu, LowerCol);
1586     if (rat) {
1587       (*Weights) (NPoleu, LowerCol) =  (*WCoefs) (NPoleu, LowerCol);
1588     }
1589
1590     for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
1591       Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
1592       Temp = Coefs (NPoleu, Col).XYZ();
1593       Temp.Divide (Cnp);
1594       Poles (NPoleu, Col).SetXYZ (Temp);
1595       if (rat) {
1596         (*Weights) (NPoleu, Col) = (*WCoefs) (NPoleu, Col) /  Cnp;
1597       }
1598     }
1599     Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
1600     if (rat) {
1601       (*Weights) (NPoleu, UpperCol) = (*WCoefs) (NPoleu, UpperCol);
1602     }
1603
1604     for (I1 = 1; I1 <= RowLength - 1; I1++) {
1605
1606       for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
1607         Temp.SetLinearForm
1608           (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
1609         Poles (NPoleu, I2).SetXYZ (Temp);
1610         if (rat) (*Weights)(NPoleu, I2) += (*Weights)(NPoleu, I2-1);
1611       }
1612     }
1613   }
1614
1615   for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
1616
1617     for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
1618       Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
1619       Temp = Poles (Row, NPolev).XYZ();
1620       Temp.Divide (Cnp);
1621       Poles (Row, NPolev).SetXYZ (Temp);
1622       if (rat) (*Weights)(Row, NPolev) /= Cnp;
1623     }
1624
1625     for (I1 = 1; I1 <= ColLength - 1; I1++) {
1626
1627       for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
1628         Temp.SetLinearForm
1629           (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
1630         Poles (I2, NPolev).SetXYZ (Temp);
1631         if (rat) (*Weights)(I2, NPolev) += (*Weights)(I2-1, NPolev);
1632       }
1633     }
1634   }
1635   if (rat) {
1636
1637     for (Row = LowerRow; Row <= UpperRow; Row++) {
1638
1639       for (Col = LowerCol; Col <= UpperCol; Col++) {
1640         W = (*Weights) (Row, Col);
1641         Temp = Poles(Row, Col).XYZ();
1642         Temp.Divide (W);
1643         Poles(Row, Col).SetXYZ (Temp);
1644       }
1645     }
1646   }
1647 }
1648
1649 //=======================================================================
1650 //function : UTrimming
1651 //purpose  :
1652 //=======================================================================
1653
1654 void PLib::UTrimming(const Standard_Real U1,
1655                      const Standard_Real U2,
1656                      TColgp_Array2OfPnt& Coeffs,
1657                      TColStd_Array2OfReal* WCoeffs)
1658 {
1659   Standard_Boolean rat = WCoeffs != NULL;
1660   Standard_Integer lr = Coeffs.LowerRow();
1661   Standard_Integer ur = Coeffs.UpperRow();
1662   Standard_Integer lc = Coeffs.LowerCol();
1663   Standard_Integer uc = Coeffs.UpperCol();
1664   TColgp_Array1OfPnt  Temp (lr,ur);
1665   TColStd_Array1OfReal Temw (lr,ur);
1666
1667   for (Standard_Integer icol = lc; icol <= uc; icol++) {
1668     Standard_Integer irow ;
1669     for ( irow = lr; irow <= ur; irow++) {
1670       Temp (irow) = Coeffs  (irow, icol);
1671       if (rat) Temw (irow) = (*WCoeffs) (irow, icol);
1672     }
1673     if (rat) PLib::Trimming (U1, U2, Temp, &Temw);
1674     else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
1675
1676     for (irow = lr; irow <= ur; irow++) {
1677       Coeffs  (irow, icol) = Temp (irow);
1678       if (rat) (*WCoeffs) (irow, icol) = Temw (irow);
1679     }
1680   }
1681 }
1682
1683 //=======================================================================
1684 //function : VTrimming
1685 //purpose  :
1686 //=======================================================================
1687
1688 void PLib::VTrimming(const Standard_Real V1,
1689                      const Standard_Real V2,
1690                      TColgp_Array2OfPnt& Coeffs,
1691                      TColStd_Array2OfReal* WCoeffs)
1692 {
1693   Standard_Boolean rat = WCoeffs != NULL;
1694   Standard_Integer lr = Coeffs.LowerRow();
1695   Standard_Integer ur = Coeffs.UpperRow();
1696   Standard_Integer lc = Coeffs.LowerCol();
1697   Standard_Integer uc = Coeffs.UpperCol();
1698   TColgp_Array1OfPnt  Temp (lc,uc);
1699   TColStd_Array1OfReal Temw (lc,uc);
1700
1701   for (Standard_Integer irow = lr; irow <= ur; irow++) {
1702     Standard_Integer icol ;
1703     for ( icol = lc; icol <= uc; icol++) {
1704       Temp (icol) = Coeffs  (irow, icol);
1705       if (rat) Temw (icol) = (*WCoeffs) (irow, icol);
1706     }
1707     if (rat) PLib::Trimming (V1, V2, Temp, &Temw);
1708     else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
1709
1710     for (icol = lc; icol <= uc; icol++) {
1711       Coeffs  (irow, icol) = Temp (icol);
1712       if (rat) (*WCoeffs) (irow, icol) = Temw (icol);
1713     }
1714   }
1715 }
1716
1717 //=======================================================================
1718 //function : HermiteInterpolate
1719 //purpose  :
1720 //=======================================================================
1721
1722 Standard_Boolean PLib::HermiteInterpolate
1723 (const Standard_Integer Dimension,
1724  const Standard_Real FirstParameter,
1725  const Standard_Real LastParameter,
1726  const Standard_Integer FirstOrder,
1727  const Standard_Integer LastOrder,
1728  const TColStd_Array2OfReal& FirstConstr,
1729  const TColStd_Array2OfReal& LastConstr,
1730  TColStd_Array1OfReal& Coefficients)
1731 {
1732   Standard_Real Pattern[3][6];
1733
1734   // portage HP : il faut les initialiser 1 par 1
1735
1736   Pattern[0][0] = 1;
1737   Pattern[0][1] = 1;
1738   Pattern[0][2] = 1;
1739   Pattern[0][3] = 1;
1740   Pattern[0][4] = 1;
1741   Pattern[0][5] = 1;
1742   Pattern[1][0] = 0;
1743   Pattern[1][1] = 1;
1744   Pattern[1][2] = 2;
1745   Pattern[1][3] = 3;
1746   Pattern[1][4] = 4;
1747   Pattern[1][5] = 5;
1748   Pattern[2][0] = 0;
1749   Pattern[2][1] = 0;
1750   Pattern[2][2] = 2;
1751   Pattern[2][3] = 6;
1752   Pattern[2][4] = 12;
1753   Pattern[2][5] = 20;
1754
1755   math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
1756   //  The initialisation of the matrix A
1757   Standard_Integer irow ;
1758   for ( irow=0; irow<=FirstOrder; irow++) {
1759     Standard_Real FirstVal = 1.;
1760
1761     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1762       A(irow,icol) = Pattern[irow][icol]*FirstVal;
1763       if (irow <= icol) FirstVal *= FirstParameter;
1764     }
1765   }
1766
1767   for (irow=0; irow<=LastOrder; irow++) {
1768     Standard_Real LastVal  = 1.;
1769
1770     for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1771       A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
1772       if (irow <= icol) LastVal *= LastParameter;
1773     }
1774   }
1775   //
1776   //  The filled matrix A for FirstOrder=LastOrder=2 is:
1777   //
1778   //      1   FP  FP**2   FP**3    FP**4     FP**5
1779   //      0   1   2*FP    3*FP**2  4*FP**3   5*FP**4        FP - FirstParameter
1780   //      0   0   2       6*FP     12*FP**2  20*FP**3
1781   //      1   LP  LP**2   LP**3    LP**4     LP**5
1782   //      0   1   2*LP    3*LP**2  4*LP**3   5*LP**4        LP - LastParameter
1783   //      0   0   2       6*LP     12*LP**2  20*LP**3
1784   //
1785   //  If FirstOrder or LastOrder <=2 then some rows and columns are missing.
1786   //  For example:
1787   //  If FirstOrder=1 then 3th row and 6th column are missing
1788   //  If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
1789
1790   math_Gauss Equations(A);
1791   //  std::cout << "A=" << A << std::endl;
1792
1793   for (Standard_Integer idim=1; idim<=Dimension; idim++) {
1794     //  std::cout << "idim=" << idim << std::endl;
1795
1796     math_Vector B(0,FirstOrder+LastOrder+1);
1797     Standard_Integer icol ;
1798     for ( icol=0; icol<=FirstOrder; icol++)
1799       B(icol) = FirstConstr(idim,icol);
1800
1801     for (icol=0; icol<=LastOrder; icol++)
1802       B(FirstOrder+1+icol) = LastConstr(idim,icol);
1803     //  std::cout << "B=" << B << std::endl;
1804
1805     //  The solving of equations system A * X = B. Then B = X
1806     Equations.Solve(B);
1807     //  std::cout << "After Solving" << std::endl << "B=" << B << std::endl;
1808
1809     if (Equations.IsDone()==Standard_False) return Standard_False;
1810
1811     //  the filling of the Coefficients
1812
1813     for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
1814       Coefficients(Dimension*icol+idim-1) = B(icol);
1815   }
1816   return Standard_True;
1817 }
1818
1819 //=======================================================================
1820 //function : JacobiParameters
1821 //purpose  :
1822 //=======================================================================
1823
1824 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
1825                             const Standard_Integer MaxDegree,
1826                             const Standard_Integer Code,
1827                             Standard_Integer& NbGaussPoints,
1828                             Standard_Integer& WorkDegree)
1829 {
1830   // ConstraintOrder: Ordre de contrainte aux extremites :
1831   //            C0 = contraintes de passage aux bornes;
1832   //            C1 = C0 + contraintes de derivees 1eres;
1833   //            C2 = C1 + contraintes de derivees 2ndes.
1834   // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
1835   //            d' approximation (doit etre superieur ou egal a
1836   //            2*NivConstr+2 et inferieur ou egal a 50).
1837   // Code:      Code d' init. des parametres de discretisation.
1838   //            (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
1839   //            = -5 Calcul tres rapide mais peu precis (8pts)
1840   //            = -4    '     '    '      '   '    '    (10pts)
1841   //            = -3    '     '    '      '   '    '    (15pts)
1842   //            = -2    '     '    '      '   '    '    (20pts)
1843   //            = -1    '     '    '      '   '    '    (25pts)
1844   //            = 1 calcul rapide avec precision moyenne (30pts).
1845   //            = 2 calcul rapide avec meilleure precision (40pts).
1846   //            = 3 calcul un peu plus lent avec bonne precision (50 pts).
1847   //            = 4 calcul lent avec la meilleure precision possible
1848   //             (61pts).
1849
1850   // The possible values of NbGaussPoints
1851
1852   const Standard_Integer NDEG8=8,   NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
1853   NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
1854
1855   Standard_Integer NivConstr=0;
1856   switch (ConstraintOrder) {
1857   case GeomAbs_C0: NivConstr = 0; break;
1858   case GeomAbs_C1: NivConstr = 1; break;
1859   case GeomAbs_C2: NivConstr = 2; break;
1860   default:
1861     throw Standard_ConstructionError("Invalid ConstraintOrder");
1862   }
1863   if (MaxDegree < 2*NivConstr+1)
1864     throw Standard_ConstructionError("Invalid MaxDegree");
1865
1866   if (Code >= 1)
1867     WorkDegree = MaxDegree + 9;
1868   else
1869     WorkDegree = MaxDegree + 6;
1870
1871   //---> Nbre mini de points necessaires.
1872   Standard_Integer IPMIN=0;
1873   if (WorkDegree < NDEG8)
1874     IPMIN=NDEG8;
1875   else if (WorkDegree < NDEG10)
1876     IPMIN=NDEG10;
1877   else if (WorkDegree < NDEG15)
1878     IPMIN=NDEG15;
1879   else if (WorkDegree < NDEG20)
1880     IPMIN=NDEG20;
1881   else if (WorkDegree < NDEG25)
1882     IPMIN=NDEG25;
1883   else if (WorkDegree < NDEG30)
1884     IPMIN=NDEG30;
1885   else if (WorkDegree < NDEG40)
1886     IPMIN=NDEG40;
1887   else if (WorkDegree < NDEG50)
1888     IPMIN=NDEG50;
1889   else if (WorkDegree < NDEG61)
1890     IPMIN=NDEG61;
1891   else
1892     throw Standard_ConstructionError("Invalid MaxDegree");
1893   // ---> Nbre de points voulus.
1894   Standard_Integer IWANT=0;
1895   switch (Code) {
1896   case -5: IWANT=NDEG8;  break;
1897   case -4: IWANT=NDEG10; break;
1898   case -3: IWANT=NDEG15; break;
1899   case -2: IWANT=NDEG20; break;
1900   case -1: IWANT=NDEG25; break;
1901   case  1: IWANT=NDEG30; break;
1902   case  2: IWANT=NDEG40; break;
1903   case  3: IWANT=NDEG50; break;
1904   case  4: IWANT=NDEG61; break;
1905   default:
1906     throw Standard_ConstructionError("Invalid Code");
1907   }
1908   //-->  NbGaussPoints est le nombre de points de discretisation de la fonction,
1909   //     il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
1910   //     NbGaussPoints doit etre superieur strictement a WorkDegree.
1911   NbGaussPoints = Max(IPMIN,IWANT);
1912   //  NbGaussPoints +=2;
1913 }
1914
1915 //=======================================================================
1916 //function : NivConstr
1917 //purpose  : translates from GeomAbs_Shape to Integer
1918 //=======================================================================
1919
1920  Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
1921 {
1922   Standard_Integer NivConstr=0;
1923   switch (ConstraintOrder) {
1924     case GeomAbs_C0: NivConstr = 0; break;
1925     case GeomAbs_C1: NivConstr = 1; break;
1926     case GeomAbs_C2: NivConstr = 2; break;
1927     default:
1928       throw Standard_ConstructionError("Invalid ConstraintOrder");
1929   }
1930   return NivConstr;
1931 }
1932
1933 //=======================================================================
1934 //function : ConstraintOrder
1935 //purpose  : translates from Integer to GeomAbs_Shape
1936 //=======================================================================
1937
1938  GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
1939 {
1940   GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
1941   switch (NivConstr) {
1942     case 0: ConstraintOrder = GeomAbs_C0; break;
1943     case 1: ConstraintOrder = GeomAbs_C1; break;
1944     case 2: ConstraintOrder = GeomAbs_C2; break;
1945     default:
1946       throw Standard_ConstructionError("Invalid NivConstr");
1947   }
1948   return ConstraintOrder;
1949 }
1950
1951 //=======================================================================
1952 //function : EvalLength
1953 //purpose  :
1954 //=======================================================================
1955
1956  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
1957                       Standard_Real& PolynomialCoeff,
1958                       const Standard_Real U1, const Standard_Real U2,
1959                       Standard_Real& Length)
1960 {
1961   Standard_Integer i,j,idim, degdim;
1962   Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
1963
1964   Standard_Real *PolynomialArray = &PolynomialCoeff ;
1965
1966   Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
1967
1968   math_Vector GaussPoints(1,NbGaussPoints);
1969   math::GaussPoints(NbGaussPoints,GaussPoints);
1970
1971   math_Vector GaussWeights(1,NbGaussPoints);
1972   math::GaussWeights(NbGaussPoints,GaussWeights);
1973
1974   C1 = (U2 + U1) / 2.;
1975   C2 = (U2 - U1) / 2.;
1976
1977 //-----------------------------------------------------------
1978 //****** Integration - Boucle sur les intervalles de GAUSS **
1979 //-----------------------------------------------------------
1980
1981   Sum = 0;
1982
1983   for (j=1; j<=NbGaussPoints/2; j++) {
1984     // Integration en tenant compte de la symetrie
1985     Tran  = C2 * GaussPoints(j);
1986     X1 = C1 + Tran;
1987     X2 = C1 - Tran;
1988
1989     //****** Derivation sur la dimension de l'espace **
1990
1991     degdim =  Degree*Dimension;
1992     Der1 = Der2 = 0.;
1993     for (idim=0; idim<Dimension; idim++) {
1994       D1 = D2 = Degree * PolynomialArray [idim + degdim];
1995       for (i=Degree-1; i>=1; i--) {
1996         DD = i * PolynomialArray [idim + i*Dimension];
1997         D1 = D1 * X1 + DD;
1998         D2 = D2 * X2 + DD;
1999       }
2000       Der1 += D1 * D1;
2001       Der2 += D2 * D2;
2002     }
2003
2004     //****** Integration **
2005
2006     Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
2007
2008 //****** Fin de boucle dur les intervalles de GAUSS **
2009   }
2010   Length = Sum;
2011 }
2012
2013
2014 //=======================================================================
2015 //function : EvalLength
2016 //purpose  :
2017 //=======================================================================
2018
2019  void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2020                       Standard_Real& PolynomialCoeff,
2021                       const Standard_Real U1, const Standard_Real U2,
2022                       const Standard_Real Tol,
2023                       Standard_Real& Length, Standard_Real& Error)
2024 {
2025   Standard_Integer i;
2026   Standard_Integer NbSubInt = 1,    // Current number of subintervals
2027                    MaxNbIter = 13,  // Max number of iterations
2028                    NbIter    = 1;   // Current number of iterations
2029   Standard_Real    dU,OldLen,LenI;
2030
2031   PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
2032
2033   do {
2034     OldLen = Length;
2035     Length = 0.;
2036     NbSubInt *= 2;
2037     dU = (U2-U1)/NbSubInt;
2038     for (i=1; i<=NbSubInt; i++) {
2039       PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
2040       Length += LenI;
2041     }
2042     NbIter++;
2043     Error = Abs(OldLen-Length);
2044   }
2045   while (Error > Tol && NbIter <= MaxNbIter);
2046 }