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