projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022550: Fixing data races
[occt.git] / src / PLib / PLib.cxx
CommitLineData
7fd59977 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
7fd59977 9#include <PLib.ixx>
41194117 10#include <PLib_LocalArray.hxx>
7fd59977 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
41194117 42class BinomAllocator
7fd59977 43{
41194117
K		 44public:
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 {
7fd59977 58 im1 = i - 1;
59 ip1 = i + 1;
60 id2 = i >> 1;
61 md2 = im1 >> 1;
62 md3 = ip1 >> 1;
63 k = 0;
41194117 64 myBinom[i] = new Standard_Integer[ip1];
7fd59977 65
41194117
K		 66 for (j = 0; j < id2; ++j)
67 {
68 myBinom[i][j] = k + myBinom[im1][j];
69 k = myBinom[im1][j];
7fd59977 70 }
71 j = id2;
72 if (j > md2) j = im1 - j;
41194117 73 myBinom[i][id2] = k + myBinom[im1][j];
7fd59977 74
41194117
K		 75 for (j = ip1 - md3; j < ip1; j++)
76 {
77 myBinom[i][j] = myBinom[i][i - j];
7fd59977 78 }
79 }
7fd59977 80 }
7fd59977 81
41194117
K		 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 }
7fd59977 92
41194117
K		 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]);
7fd59977 99 }
41194117
K		 100
101private:
102 Standard_Integer** myBinom;
103 Standard_Integer myMaxBinom;
104
105};
106
107namespace
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
119Standard_Real PLib::Bin(const Standard_Integer N,
120 const Standard_Integer P)
121{
122 return THE_BINOM.Value (N, P);
7fd59977 123}
124
125//=======================================================================
126//function : RationalDerivative
127//purpose :
128//=======================================================================
129
130void 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;
41194117
K		 179 PLib_LocalArray binomial_array;
180 PLib_LocalArray derivative_storage;
7fd59977 181 if (Dimension == 3) {
182 Standard_Integer DeRequest1 = DerivativeRequest + 1;
183 Standard_Integer MinDegRequ = DerivativeRequest;
184 if (MinDegRequ > Degree) MinDegRequ = Degree;
41194117 185 binomial_array.Allocate (DeRequest1);
7fd59977 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;
41194117 192 derivative_storage.Allocate (DimDeRequ1);
7fd59977 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;
41194117 267 binomial_array.Allocate (DeRequest1);
7fd59977 268
269 for (ii = 0 ; ii < DeRequest1 ; ii++) {
270 binomial_array[ii] = 1.0e0 ;
271 }
272 if (!All) {
273 Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
41194117 274 derivative_storage.Allocate (DimDeRequ1);
7fd59977 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
356void 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
41194117
K		 399 PLib_LocalArray binomial_array (DeRequest1);
400 PLib_LocalArray derivative_storage;
7fd59977 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
476void 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
1577void 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
1807void 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
1879static Standard_Integer storage_divided = 0 ;
1880static 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
1888Standard_Integer
1889PLib::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 ;
41194117
K		 1931 }
1932 PLib_LocalArray divided_differences_array ((Degree + 1) * Dimension);
7fd59977 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
2013Standard_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 }
41194117 2062 PLib_LocalArray divided_differences_array ((Degree + 1) * Dimension);
7fd59977 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
2152Standard_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
2214void 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
2232void 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
2250void 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
2263void 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;
7fd59977 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
2331void 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
2347void 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
2363void 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
2376void 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
2444void 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;
7fd59977 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 }
7fd59977 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
2538void 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
2572void 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
2606Standard_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
2708void 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}
Open CASCADE Technology repository