src/PLib/PLib.cxx

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