NGSolve
4.9
|
00001 #ifndef FILE_RECURSIVE_POL 00002 #define FILE_RECURSIVE_POL 00003 00004 /*********************************************************************/ 00005 /* File: recursive_pol.hpp */ 00006 /* Author: Start */ 00007 /* Date: 6. Feb. 2003 */ 00008 /*********************************************************************/ 00009 00010 namespace ngfem 00011 { 00012 00013 /* 00014 Recursive Polynomials 00015 */ 00016 00017 00018 00019 00020 00021 template <class REC, int N> 00022 class CEvalFO 00023 { 00024 public: 00025 template <class S, class T> 00026 ALWAYS_INLINE static void Eval (S x, T & values, S & p1, S & p2) 00027 { 00028 S p3; 00029 CEvalFO<REC,N-1>::Eval (x, values, p2, p3); 00030 values[N] = p1 = ( REC::A(N) * x + REC::B(N)) * p2 + REC::C(N) * p3; 00031 } 00032 00033 00034 template <class S, class Sc, class T> 00035 ALWAYS_INLINE static void EvalMult (S x, Sc c, T & values, S & p1, S & p2) 00036 { 00037 S p3; 00038 CEvalFO<REC,N-1>::EvalMult (x, c, values, p2, p3); 00039 values[N] = p1 = ( REC::A(N) * x + REC::B(N)) * p2 + REC::C(N) * p3; 00040 } 00041 00042 00043 template <class S, class Sy, class T> 00044 ALWAYS_INLINE static void EvalScaled (S x, Sy y, T & values, S & p1, S & p2) 00045 { 00046 S p3; 00047 CEvalFO<REC,N-1>::EvalScaled (x, y, values, p2, p3); 00048 values[N] = p1 = ( REC::A(N) * x + REC::B(N) * y) * p2 + REC::C(N)*y*y * p3; 00049 } 00050 00051 00052 template <class S, class Sy, class Sc, class T> 00053 ALWAYS_INLINE static void EvalScaledMult (S x, Sy y, Sc c, T & values, S & p1, S & p2) 00054 { 00055 S p3; 00056 CEvalFO<REC,N-1>::EvalScaledMult (x, y, c, values, p2, p3); 00057 values[N] = p1 = ( REC::A(N) * x + REC::B(N) * y) * p2 + REC::C(N)*y*y * p3; 00058 } 00059 00060 00061 }; 00062 00063 00064 template <class REC> 00065 class CEvalFO<REC, -1> 00066 { 00067 public: 00068 template <class S, class T> 00069 ALWAYS_INLINE static void Eval (S x, T & values, S & /* p1 */, S & /* p2 */) 00070 { ; } 00071 00072 template <class S, class Sc, class T> 00073 ALWAYS_INLINE static void EvalMult (S x, Sc c, T & values, S & /* p1 */, S & /* p2 */) 00074 { ; } 00075 00076 00077 template <class S, class Sy, class T> 00078 ALWAYS_INLINE static void EvalScaled (S x, Sy y, T & values, S & /* p1 */, S & /* p2 */) 00079 { ; } 00080 00081 template <class S, class Sy, class Sc, class T> 00082 ALWAYS_INLINE static void EvalScaledMult (S x, Sy y, Sc c, T & values, S & /* p1 */, S & /* p2 */) 00083 { ; } 00084 00085 }; 00086 00087 00088 template <class REC> 00089 class CEvalFO<REC, 0> 00090 { 00091 public: 00092 template <class S, class T> 00093 ALWAYS_INLINE static void Eval (S x, T & values, S & p1, S & /* p2 */) 00094 { 00095 values[0] = p1 = REC::P0(x); 00096 } 00097 00098 template <class S, class Sc, class T> 00099 ALWAYS_INLINE static void EvalMult (S x, Sc c, T & values, S & p1, S & /* p2 */) 00100 { 00101 values[0] = p1 = c * REC::P0(x); 00102 } 00103 00104 00105 template <class S, class Sy, class T> 00106 ALWAYS_INLINE static void EvalScaled (S x, Sy y, T & values, S & p1, S & /* p2 */) 00107 { 00108 values[0] = p1 = REC::P0(x); 00109 } 00110 00111 template <class S, class Sy, class Sc, class T> 00112 ALWAYS_INLINE static void EvalScaledMult (S x, Sy y, Sc c, T & values, S & p1, S & /* p2 */) 00113 { 00114 values[0] = p1 = c * REC::P0(x); 00115 } 00116 00117 }; 00118 00119 template <class REC> 00120 class CEvalFO<REC, 1> 00121 { 00122 public: 00123 template <class S, class T> 00124 ALWAYS_INLINE static void Eval (S x, T & values, S & p1, S & p2) 00125 { 00126 values[0] = p2 = REC::P0(x); 00127 values[1] = p1 = REC::P1(x); 00128 } 00129 00130 template <class S, class Sc, class T> 00131 ALWAYS_INLINE static void EvalMult (S x, Sc c, T & values, S & p1, S & p2) 00132 { 00133 values[0] = p2 = c * REC::P0(x); 00134 values[1] = p1 = c * REC::P1(x); 00135 } 00136 00137 template <class S, class Sy, class T> 00138 ALWAYS_INLINE static void EvalScaled (S x, Sy y, T & values, S & p1, S & p2) 00139 { 00140 values[0] = p2 = REC::P0(x); 00141 values[1] = p1 = REC::P1(x); 00142 } 00143 00144 template <class S, class Sy, class Sc, class T> 00145 ALWAYS_INLINE static void EvalScaledMult (S x, Sy y, Sc c, T & values, S & p1, S & p2) 00146 { 00147 values[0] = p2 = c * REC::P0(x); 00148 values[1] = p1 = c * REC::P1(x); 00149 } 00150 }; 00151 00152 00153 00154 00155 00156 // P_i = (a_i x + b_i) P_{i-1} + c_i P_{i-2} 00157 template<class REC> 00158 class RecursivePolynomial 00159 { 00160 public: 00161 template <int N, class S, class T> 00162 ALWAYS_INLINE static void EvalFO (S x, T & values) 00163 { 00164 S p1, p2; 00165 CEvalFO<REC, N>::Eval (x, values, p1, p2); 00166 } 00167 00168 00169 template <class S> 00170 ALWAYS_INLINE static void EvalNext (int i, S x, S & p1, S & p2) 00171 { 00172 S pnew = (REC::A(i) * x + REC::B(i)) * p1 + REC::C(i) * p2; 00173 p2 = p1; 00174 p1 = pnew; 00175 } 00176 00177 template <class S, class Sy> 00178 ALWAYS_INLINE static void EvalScaledNext (int i, S x, Sy y, S & p1, S & p2) 00179 { 00180 S pnew = (REC::A(i) * x + REC::B(i) * y) * p1 + REC::C(i) * y*y*p2; 00181 p2 = p1; 00182 p1 = pnew; 00183 } 00184 00185 template <class S, class T> 00186 static void Eval (int n, S x, T & values) 00187 { 00188 EvalMult (n, x, 1.0, values); 00189 } 00190 00191 00192 template <class S, class Sc, class T> 00193 static void EvalMult (int n, S x, Sc c, T & values) 00194 { 00195 S p1, p2; 00196 00197 if (n < 0) return; 00198 00199 values[0] = p2 = c * REC::P0(x); 00200 if (n < 1) return; 00201 00202 values[1] = p1 = c * REC::P1(x); 00203 if (n < 2) return; 00204 00205 EvalNext(2, x, p1, p2); 00206 values[2] = p1; 00207 if (n < 3) return; 00208 00209 EvalNext(3, x, p1, p2); 00210 values[3] = p1; 00211 if (n < 4) return; 00212 00213 EvalNext(4, x, p1, p2); 00214 values[4] = p1; 00215 if (n < 5) return; 00216 00217 EvalNext(5, x, p1, p2); 00218 values[5] = p1; 00219 if (n < 6) return; 00220 00221 EvalNext(6, x, p1, p2); 00222 values[6] = p1; 00223 if (n < 7) return; 00224 00225 EvalNext(7, x, p1, p2); 00226 values[7] = p1; 00227 if (n < 8) return; 00228 00229 EvalNext(8, x, p1, p2); 00230 values[8] = p1; 00231 if (n < 9) return; 00232 00233 EvalNext(9, x, p1, p2); 00234 values[9] = p1; 00235 if (n < 10) return; 00236 00237 EvalNext(10, x, p1, p2); 00238 values[10] = p1; 00239 if (n < 11) return; 00240 00241 EvalNext(11, x, p1, p2); 00242 values[11] = p1; 00243 if (n < 12) return; 00244 00245 00246 for (int i = 12; i <= n; i++) 00247 { 00248 EvalNext (i, x, p1, p2); 00249 values[i] = p1; 00250 } 00251 00252 /* 00253 const REC & pol = static_cast<const REC&> (*this); 00254 00255 if (n < 0) return; 00256 00257 switch (n) 00258 { 00259 case 0: EvalMultFO<0> (x, c, values); return; 00260 case 1: EvalMultFO<1> (x, c, values); return; 00261 case 2: EvalMultFO<2> (x, c, values); return; 00262 case 3: EvalMultFO<3> (x, c, values); return; 00263 case 4: EvalMultFO<4> (x, c, values); return; 00264 case 5: EvalMultFO<5> (x, c, values); return; 00265 case 6: EvalMultFO<6> (x, c, values); return; 00266 case 7: EvalMultFO<7> (x, c, values); return; 00267 default: 00268 S p1, p2, p3; 00269 CEvalFO<REC, 8>::EvalMult (pol, x, c, values, p1, p2); 00270 00271 for (int i = 9; i <= n; i++) 00272 { 00273 p3 = p2; p2 = p1; 00274 p1 = (REC::A(i) * x + REC::B(i)) * p2 + REC::C(i) * p3; 00275 values[i] = p1; 00276 } 00277 } 00278 */ 00279 } 00280 00281 template <int N, class S, class Sc, class T> 00282 ALWAYS_INLINE static void EvalMultFO (S x, Sc c, T & values) 00283 { 00284 S p1, p2; 00285 CEvalFO<REC, N>::EvalMult (x, c, values, p1, p2); 00286 } 00287 00288 00289 00290 00291 00292 template <class S, class Sy, class T> 00293 static void EvalScaled (int n, S x, Sy y, T & values) 00294 { 00295 EvalScaledMult (n, x, y, 1.0, values); 00296 } 00297 00298 template <int N, class S, class Sy, class T> 00299 ALWAYS_INLINE static void EvalScaledFO (S x, Sy y, T & values) 00300 { 00301 S p1, p2; 00302 CEvalFO<REC, N>::EvalScaled (x, y, values, p1, p2); 00303 } 00304 00305 00306 00307 template <class S, class Sy, class Sc, class T> 00308 static void EvalScaledMult (int n, S x, Sy y, Sc c, T & values) 00309 { 00310 S p1, p2; 00311 00312 if (n < 0) return; 00313 00314 values[0] = p2 = c * REC::P0(x); 00315 if (n < 1) return; 00316 00317 values[1] = p1 = c * REC::P1(x); 00318 if (n < 2) return; 00319 00320 EvalScaledNext(2, x, y, p1, p2); 00321 values[2] = p1; 00322 if (n < 3) return; 00323 00324 EvalScaledNext(3, x, y, p1, p2); 00325 values[3] = p1; 00326 if (n < 4) return; 00327 00328 EvalScaledNext(4, x, y, p1, p2); 00329 values[4] = p1; 00330 if (n < 5) return; 00331 00332 EvalScaledNext(5, x, y, p1, p2); 00333 values[5] = p1; 00334 if (n < 6) return; 00335 00336 EvalScaledNext(6, x, y, p1, p2); 00337 values[6] = p1; 00338 if (n < 7) return; 00339 00340 EvalScaledNext(7, x, y, p1, p2); 00341 values[7] = p1; 00342 if (n < 8) return; 00343 00344 EvalScaledNext(8, x, y, p1, p2); 00345 values[8] = p1; 00346 00347 for (int i = 9; i <= n; i++) 00348 { 00349 EvalScaledNext (i, x, y, p1, p2); 00350 values[i] = p1; 00351 } 00352 /* 00353 const REC & pol = static_cast<const REC&> (*this); 00354 00355 if (n < 0) return; 00356 00357 switch (n) 00358 { 00359 case 0: EvalScaledMultFO<0> (x, y, c, values); return; 00360 case 1: EvalScaledMultFO<1> (x, y, c, values); return; 00361 case 2: EvalScaledMultFO<2> (x, y, c, values); return; 00362 case 3: EvalScaledMultFO<3> (x, y, c, values); return; 00363 case 4: EvalScaledMultFO<4> (x, y, c, values); return; 00364 case 5: EvalScaledMultFO<5> (x, y, c, values); return; 00365 case 6: EvalScaledMultFO<6> (x, y, c, values); return; 00366 case 7: EvalScaledMultFO<7> (x, y, c, values); return; 00367 default: 00368 S p1, p2, p3; 00369 CEvalFO<REC, 8>::EvalScaledMult (pol, x, y, c, values, p1, p2); 00370 00371 for (int i = 9; i <= n; i++) 00372 { 00373 p3 = p2; p2 = p1; 00374 p1 = (REC::A(i) * x + REC::B(i) * y) * p2 + REC::C(i)*y*y * p3; 00375 values[i] = p1; 00376 } 00377 } 00378 */ 00379 } 00380 00381 template <int N, class S, class Sy, class Sc, class T> 00382 ALWAYS_INLINE static void EvalScaledMultFO (S x, Sy y, Sc c,T & values) 00383 { 00384 S p1, p2; 00385 CEvalFO<REC, N>::EvalScaledMult (x, y, c, values, p1, p2); 00386 } 00387 }; 00388 00389 00390 00391 class LegendrePolynomial : public RecursivePolynomial<LegendrePolynomial> 00392 { 00393 public: 00394 LegendrePolynomial () { ; } 00395 00396 template <class S, class T> 00397 inline LegendrePolynomial (int n, S x, T & values) 00398 { 00399 Eval (n, x, values); 00400 } 00401 00402 template <class S> 00403 static S P0(S x) { return S(1.0); } 00404 template <class S> 00405 static S P1(S x) { return x; } 00406 00407 static double A (int i) { return 2.0-1.0/i; } 00408 static double B (int i) { return 0; } 00409 static double C (int i) { return 1.0/i-1.0; } 00410 }; 00411 00412 00413 00414 00415 00416 00417 template <int al, int be> 00418 class JacobiPolynomialFix : public RecursivePolynomial<JacobiPolynomialFix<al,be> > 00419 { 00420 public: 00421 JacobiPolynomialFix () { ; } 00422 00423 template <class S, class T> 00424 inline JacobiPolynomialFix (int n, S x, T & values) 00425 { 00426 Eval (n, x, values); 00427 } 00428 00429 00430 00431 template <class S> 00432 static ALWAYS_INLINE S P0(S x) { return S(1.0); } 00433 template <class S> 00434 static ALWAYS_INLINE S P1(S x) { return 0.5 * (2*(al+1)+(al+be+2)*(x-1)); } 00435 00436 static ALWAYS_INLINE double A (int i) 00437 { i--; return (2.0*i+al+be)*(2*i+al+be+1)*(2*i+al+be+2) / ( 2 * (i+1) * (i+al+be+1) * (2*i+al+be)); } 00438 static ALWAYS_INLINE double B (int i) 00439 { i--; return (2.0*i+al+be+1)*(al*al-be*be) / ( 2 * (i+1) * (i+al+be+1) * (2*i+al+be)); } 00440 static ALWAYS_INLINE double C (int i) 00441 { i--; return -2.0*(i+al)*(i+be) * (2*i+al+be+2) / ( 2 * (i+1) * (i+al+be+1) * (2*i+al+be)); } 00442 }; 00443 00444 00445 00446 00447 00448 /* 00449 class JacobiPolynomial2 : public RecursivePolynomial<JacobiPolynomial2> 00450 { 00451 double al, be; 00452 public: 00453 JacobiPolynomial2 (double aal, double abe) : al(aal), be(abe) { ; } 00454 00455 template <class S> 00456 S P0(S x) const { return S(1.0); } 00457 template <class S> 00458 S P1(S x) const { return 0.5 * (2*(al+1)+(al+be+2)*(x-1)); } 00459 00460 double A (int i) const 00461 { i--; return (2.0*i+al+be)*(2*i+al+be+1)*(2*i+al+be+2) / ( 2 * (i+1) * (i+al+be+1) * (2*i+al+be)); } 00462 double B (int i) const 00463 { i--; return (2.0*i+al+be+1)*(al*al-be*be) / ( 2 * (i+1) * (i+al+be+1) * (2*i+al+be)); } 00464 double C (int i) const 00465 { i--; return -2.0*(i+al)*(i+be) * (2*i+al+be+2) / ( 2 * (i+1) * (i+al+be+1) * (2*i+al+be)); } 00466 }; 00467 */ 00468 00469 00470 00471 00472 00473 00474 00475 00476 00477 00478 00479 00490 template <class S, class T> 00491 inline void LegendrePolynomial1 (int n, S x, T & values) 00492 { 00493 S p1, p2; 00494 00495 if (n < 0) return; 00496 values[0] = 1.0; 00497 if (n < 1) return; 00498 values[1] = x; 00499 if (n < 2) return; 00500 values[2] = p2 = 1.5 * sqr(x) - 0.5; 00501 if (n < 3) return; 00502 values[3] = p1 = ( (5.0/3.0) * p2 - (2.0/3.0)) * x; 00503 if (n < 4) return; 00504 00505 for (int j=4; j < n; j+=2) 00506 { 00507 double invj = 1.0 / j; 00508 p2 *= (invj-1); 00509 p2 += (2-invj) * x * p1; 00510 values[j] = p2; 00511 00512 invj = 1.0 / (j+1); 00513 p1 *= (invj-1); 00514 p1 += (2-invj) * x * p2; 00515 values[j+1] = p1; 00516 } 00517 00518 if (n % 2 == 0) 00519 { 00520 double invn = 1.0 / n; 00521 values[n] = (2-invn)*x*p1 - (1-invn)*p2; 00522 } 00523 00524 /* 00525 S p1 = 1.0, p2 = 0.0, p3; 00526 00527 if (n >= 0) 00528 values[0] = 1.0; 00529 00530 for (int j=1; j<=n; j++) 00531 { 00532 p3 = p2; p2 = p1; 00533 p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j; 00534 values[j] = p1; 00535 } 00536 */ 00537 } 00538 00539 00540 template <class S, class Sc, class T> 00541 inline void LegendrePolynomialMult (int n, S x, Sc c , T & values) 00542 { 00543 LegendrePolynomial leg; 00544 leg.EvalMult (n, x, c, values); 00545 } 00546 00547 00548 00549 00550 00551 template <class S, class T> 00552 inline void JacobiPolynomial (int n, S x, double alpha, double beta, T & values) 00553 { 00554 S p1 = 1.0, p2 = 0.0, p3; 00555 00556 if (n >= 0) 00557 p2 = values[0] = 1.0; 00558 if (n >= 1) 00559 p1 = values[1] = 0.5 * (2*(alpha+1)+(alpha+beta+2)*(x-1)); 00560 00561 for (int i = 1; i < n; i++) 00562 { 00563 p3 = p2; p2=p1; 00564 p1 = 00565 1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) * 00566 ( 00567 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) + 00568 (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 00569 * p2 00570 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) * p3 00571 ); 00572 values[i+1] = p1; 00573 } 00574 } 00575 00576 00577 00578 00579 00580 template <class S, class Sc, class T> 00581 inline void JacobiPolynomialMult (int n, S x, double alpha, double beta, Sc c, T & values) 00582 { 00583 S p1 = c, p2 = 0.0, p3; 00584 00585 if (n >= 0) 00586 p2 = values[0] = c; 00587 if (n >= 1) 00588 p1 = values[1] = 0.5 * c * (2*(alpha+1)+(alpha+beta+2)*(x-1)); 00589 00590 for (int i = 1; i < n; i++) 00591 { 00592 p3 = p2; p2=p1; 00593 p1 = 00594 1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) * 00595 ( 00596 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) + 00597 (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 00598 * p2 00599 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) * p3 00600 ); 00601 values[i+1] = p1; 00602 } 00603 } 00604 00605 00606 00607 00608 00609 00610 00611 template <class S, class St, class T> 00612 inline void ScaledJacobiPolynomial (int n, S x, St t, double alpha, double beta, T & values) 00613 { 00614 /* 00615 S p1 = 1.0, p2 = 0.0, p3; 00616 00617 if (n >= 0) values[0] = 1.0; 00618 */ 00619 00620 S p1 = 1.0, p2 = 0.0, p3; 00621 00622 if (n >= 0) 00623 p2 = values[0] = 1.0; 00624 if (n >= 1) 00625 p1 = values[1] = 0.5 * (2*(alpha+1)*t+(alpha+beta+2)*(x-t)); 00626 00627 for (int i=1; i < n; i++) 00628 { 00629 p3 = p2; p2=p1; 00630 p1 = 00631 1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) * 00632 ( 00633 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) * t + 00634 (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 00635 * p2 00636 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) * t * t * p3 00637 ); 00638 values[i+1] = p1; 00639 } 00640 } 00641 00642 00643 00644 00645 00646 00647 template <class S, class St, class Sc, class T> 00648 inline void ScaledJacobiPolynomialMult (int n, S x, St t, double alpha, double beta, Sc c, T & values) 00649 { 00650 /* 00651 S p1 = 1.0, p2 = 0.0, p3; 00652 if (n >= 0) values[0] = 1.0; 00653 */ 00654 00655 S p1 = c, p2 = 0.0, p3; 00656 00657 if (n >= 0) 00658 p2 = values[0] = c; 00659 if (n >= 1) 00660 p1 = values[1] = 0.5 * c * (2*(alpha+1)*t+(alpha+beta+2)*(x-t)); 00661 00662 for (int i=1; i < n; i++) 00663 { 00664 p3 = p2; p2=p1; 00665 p1 = 00666 1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) * 00667 ( 00668 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) * t + 00669 (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 00670 * p2 00671 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) * t * t * p3 00672 ); 00673 values[i+1] = p1; 00674 } 00675 } 00676 00677 00678 00679 00680 00681 00682 00683 00684 00685 00686 00687 #ifdef V1 00688 00691 template <class S, class Sc, class T> 00692 inline void LegendrePolynomialMult (int n, S x, Sc c , T & values) 00693 { 00694 S p1, p2; 00695 00696 if (n < 0) return; 00697 values[0] = c; 00698 if (n < 1) return; 00699 values[1] = c * x; 00700 if (n < 2) return; 00701 values[2] = p2 = c * (1.5 * sqr(x) - 0.5); 00702 if (n < 3) return; 00703 values[3] = p1 = ( (5.0/3.0) * p2 - (2.0/3.0) * c) * x; 00704 if (n < 4) return; 00705 00706 for (int j=4; j < n; j+=2) 00707 { 00708 double invj = 1.0 / j; 00709 p2 *= (invj-1); 00710 p2 += (2-invj) * x * p1; 00711 values[j] = p2; 00712 00713 invj = 1.0 / (j+1); 00714 p1 *= (invj-1); 00715 p1 += (2-invj) * x * p2; 00716 values[j+1] = p1; 00717 } 00718 00719 if (n % 2 == 0) 00720 { 00721 double invn = 1.0 / n; 00722 values[n] = (2-invn)*x*p1 - (1-invn)*p2; 00723 } 00724 00725 /* 00726 S p1 = 1.0, p2 = 0.0, p3; 00727 00728 if (n >= 0) 00729 values[0] = 1.0; 00730 00731 for (int j=1; j<=n; j++) 00732 { 00733 p3 = p2; p2 = p1; 00734 p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j; 00735 values[j] = p1; 00736 } 00737 */ 00738 } 00739 00740 00741 00742 00743 00744 00745 00746 template <int n> 00747 class LegendrePolynomialFO 00748 { 00749 public: 00750 template <class S, class T> 00751 static void Eval (S x, T & values) 00752 { 00753 LegendrePolynomialFO<n-1>::Eval (x, values); 00754 values[n] = (2.0*n-1)/n * x * values[n-1] - (n-1.0)/n * values[n-2]; 00755 } 00756 00757 template <class S, class Sc, class T> 00758 static void EvalMult (S x, Sc c, T & values) 00759 { 00760 LegendrePolynomialFO<n-1>::EvalMult (x, c, values); 00761 values[n] = (2.0*n-1)/n * x * values[n-1] - (n-1.0)/n * values[n-2]; 00762 } 00763 00764 template <class S, class Sy, class T> 00765 static void EvalScaled (S x, Sy y, T & values) 00766 { 00767 LegendrePolynomialFO<n-1>::EvalScaled (x, y, values); 00768 values[n] = (2.0*n-1)/n * x * values[n-1] - (n-1.0)/n * y*y * values[n-2]; 00769 } 00770 00771 template <class S, class Sy, class Sc, class T> 00772 static void EvalScaledMult (S x, Sy y, Sc c, T & values) 00773 { 00774 LegendrePolynomialFO<n-1>::EvalScaledMult (x, y, c, values); 00775 values[n] = (2.0*n-1)/n * x * values[n-1] - (n-1.0)/n * y*y * values[n-2]; 00776 } 00777 00778 }; 00779 00780 00781 template <> class LegendrePolynomialFO<-1> 00782 { 00783 public: 00784 template <class S, class T> 00785 static void Eval (S x, T & values) 00786 { ; } 00787 00788 template <class S, class Sc, class T> 00789 static void EvalMult (S x, Sc c, T & values) 00790 { ; } 00791 00792 template <class S, class Sy, class T> 00793 static void EvalScaled (S x, S y, T & values) 00794 { ; } 00795 00796 template <class S, class Sy, class Sc, class T> 00797 static void EvalScaledMult (S x, Sy y, Sc c, T & values) 00798 { ; } 00799 }; 00800 00801 template <> class LegendrePolynomialFO<0> 00802 { 00803 public: 00804 template <class S, class T> 00805 static void Eval (S x, T & values) 00806 { 00807 values[0] = 1; 00808 } 00809 00810 template <class S, class Sc, class T> 00811 static void EvalMult (S x, Sc c, T & values) 00812 { 00813 values[0] = c; 00814 } 00815 00816 template <class S, class Sy, class T> 00817 static void EvalScaled (S x, Sy y, T & values) 00818 { 00819 values[0] = 1; 00820 } 00821 00822 template <class S, class Sy, class Sc, class T> 00823 static void EvalScaledMult (S x, Sy y, Sc c, T & values) 00824 { 00825 values[0] = c; 00826 } 00827 }; 00828 00829 template <> class LegendrePolynomialFO<1> 00830 { 00831 public: 00832 template <class S, class T> 00833 static void Eval (S x, T & values) 00834 { 00835 values[0] = 1; 00836 values[1] = x; 00837 } 00838 00839 template <class S, class Sc, class T> 00840 static void EvalMult (S x, Sc c, T & values) 00841 { 00842 values[0] = c; 00843 values[1] = c*x; 00844 } 00845 00846 template <class S, class Sy, class T> 00847 static void EvalScaled (S x, Sy y, T & values) 00848 { 00849 values[0] = 1; 00850 values[1] = x; 00851 } 00852 00853 template <class S, class Sy, class Sc, class T> 00854 static void EvalScaledMult (S x, Sy y, Sc c, T & values) 00855 { 00856 values[0] = c; 00857 values[1] = c*x; 00858 } 00859 }; 00860 #endif 00861 00862 00863 00864 00865 00866 #ifdef V1 00867 template <int n, int alpha, int beta> 00868 class JacobiPolynomialFO 00869 { 00870 public: 00871 template <class S, class T> 00872 static void Eval (S x, T & values) 00873 { 00874 JacobiPolynomialFO<n-1, alpha, beta>::Eval (x, values); 00875 int i = n-1; 00876 values[n] = 00877 1.0 / ( 2 * (i+1) * (i+alpha+beta+1) * (2*i+alpha+beta) ) * 00878 ( 00879 ( (2*i+alpha+beta+1)*(alpha*alpha-beta*beta) + 00880 (2*i+alpha+beta)*(2*i+alpha+beta+1)*(2*i+alpha+beta+2) * x) 00881 * values[n-1] 00882 - 2*(i+alpha)*(i+beta) * (2*i+alpha+beta+2) 00883 * values[n-2] 00884 ); 00885 // (2.0*n-1)/n * x * values[n-1] - (n-1.0)/n * values[n-2]; 00886 } 00887 /* 00888 template <class S, class Sc, class T> 00889 static void EvalMult (S x, Sc c, T & values) 00890 { 00891 JacobiPolynomialFO<n-1>::EvalMult (x, c, values); 00892 values[n] = (2.0*n-1)/n * x * values[n-1] - (n-1.0)/n * values[n-2]; 00893 } 00894 */ 00895 }; 00896 00897 00898 template <int alpha, int beta> class JacobiPolynomialFO<-1, alpha, beta> 00899 { 00900 public: 00901 template <class S, class T> 00902 static void Eval (S x, T & values) 00903 { ; } 00904 /* 00905 template <class S, class Sc, class T> 00906 static void EvalMult (S x, Sc c, T & values) 00907 { ; } 00908 */ 00909 }; 00910 00911 template <int alpha, int beta> class JacobiPolynomialFO<0, alpha, beta> 00912 { 00913 public: 00914 template <class S, class T> 00915 static void Eval (S x, T & values) 00916 { 00917 values[0] = 1; 00918 } 00919 /* 00920 template <class S, class Sc, class T> 00921 static void EvalMult (S x, Sc c, T & values) 00922 { 00923 values[0] = c; 00924 } 00925 */ 00926 }; 00927 00928 template <int alpha, int beta> class JacobiPolynomialFO<1, alpha, beta> 00929 { 00930 public: 00931 template <class S, class T> 00932 static void Eval (S x, T & values) 00933 { 00934 values[0] = 1; 00935 values[1] = 0.5 * (2*(alpha+1)+(alpha+beta+2)*(x-1)); 00936 } 00937 /* 00938 template <class S, class Sc, class T> 00939 static void EvalMult (S x, Sc c, T & values) 00940 { 00941 values[0] = c; 00942 values[1] = c*x; 00943 } 00944 */ 00945 }; 00946 #endif 00947 00948 00949 00950 00951 00952 // compute J_j^{2i+alpha0, beta} (x), for i+j <= n 00953 00954 template <class S, class T> 00955 void DubinerJacobiPolynomials1 (int n, S x, int alpha0, int beta, T & values) 00956 { 00957 for (int i = 0; i <= n; i++) 00958 JacobiPolynomial (n-i, x, alpha0+2*i, beta, values.Row(i)); 00959 } 00960 00961 00962 template <int n, int i, int alpha0, int beta> 00963 class DubinerJacobiPolynomialsFO 00964 { 00965 public: 00966 template <class S, class T> 00967 static void Eval (S x, T & values) 00968 { 00969 DubinerJacobiPolynomialsFO<n, i-1, alpha0, beta>::Eval (x, values); 00970 // JacobiPolynomialFO<n-i, alpha0+2*i, beta>::Eval (x, values.Row(i)); 00971 00972 JacobiPolynomialFix<alpha0+2*i, beta> jac; 00973 S p1, p2; 00974 CEvalFO<JacobiPolynomialFix<alpha0+2*i, beta>, n-i>::Eval (x, values.Row(i), p1, p2); 00975 } 00976 00977 00978 template <class S, class St, class T> 00979 static void EvalScaled (S x, St t, T & values) 00980 { 00981 DubinerJacobiPolynomialsFO<n, i-1, alpha0, beta>::EvalScaled (x, t, values); 00982 // JacobiPolynomialFO<n-i, alpha0+2*i, beta>::Eval (x, values.Row(i)); 00983 00984 JacobiPolynomialFix<alpha0+2*i, beta> jac; 00985 S p1, p2; 00986 CEvalFO<JacobiPolynomialFix<alpha0+2*i, beta>, n-i>::EvalScaled (x, t, values.Row(i), p1, p2); 00987 } 00988 00989 00990 }; 00991 00992 template <int n, int alpha0, int beta> 00993 class DubinerJacobiPolynomialsFO<n, -1, alpha0, beta> 00994 { 00995 public: 00996 template <class S, class T> 00997 static void Eval (S x, T & values) 00998 { ; } 00999 template <class S, class St, class T> 01000 static void EvalScaled (S x, St t, T & values) 01001 { ; } 01002 }; 01003 01004 01005 01006 01007 01008 01009 01010 template <int n, int i, int alpha0, int beta> 01011 class DubinerJacobiPolynomialsPowFO 01012 { 01013 public: 01014 template <class S, class T> 01015 static S Eval (S x, T & values) 01016 { 01017 S power = DubinerJacobiPolynomialsPowFO<n, i-1, alpha0, beta>::Eval (x, values); 01018 S p1, p2; 01019 CEvalFO<JacobiPolynomialFix<alpha0+2*i, beta>, n-i>::EvalMult (x, power, values.Row(i), p1, p2); 01020 return power * (1-x)/2; 01021 } 01022 01023 01024 template <class S, class St, class T> 01025 static S EvalScaled (S x, St t, T & values) 01026 { 01027 S power = DubinerJacobiPolynomialsPowFO<n, i-1, alpha0, beta>::EvalScaled (x, t, values); 01028 S p1, p2; 01029 CEvalFO<JacobiPolynomialFix<alpha0+2*i, beta>, n-i>::EvalScaledMult (x, t, power, values.Row(i), p1, p2); 01030 return power * (1-x)/2; 01031 } 01032 }; 01033 01034 template <int n, int alpha0, int beta> 01035 class DubinerJacobiPolynomialsPowFO<n, -1, alpha0, beta> 01036 { 01037 public: 01038 template <class S, class T> 01039 static S Eval (S x, T & values) 01040 { return 1.0; } 01041 template <class S, class St, class T> 01042 static S EvalScaled (S x, St t, T & values) 01043 { return 1.0; } 01044 }; 01045 01046 01047 01048 01049 01050 01051 01052 01053 01054 01055 01056 01057 01058 01059 01060 01061 01062 template <int ALPHA0, int BETA, class S, class T> 01063 void DubinerJacobiPolynomials2 (int n, S x, T & values) 01064 { 01065 switch (n) 01066 { 01067 case 0: DubinerJacobiPolynomialsFO<0, 0, ALPHA0, BETA>::Eval (x, values); break; 01068 case 1: DubinerJacobiPolynomialsFO<1, 1, ALPHA0, BETA>::Eval (x, values); break; 01069 case 2: DubinerJacobiPolynomialsFO<2, 2, ALPHA0, BETA>::Eval (x, values); break; 01070 case 3: DubinerJacobiPolynomialsFO<3, 3, ALPHA0, BETA>::Eval (x, values); break; 01071 case 4: DubinerJacobiPolynomialsFO<4, 4, ALPHA0, BETA>::Eval (x, values); break; 01072 case 5: DubinerJacobiPolynomialsFO<5, 5, ALPHA0, BETA>::Eval (x, values); break; 01073 case 6: DubinerJacobiPolynomialsFO<6, 6, ALPHA0, BETA>::Eval (x, values); break; 01074 case 7: DubinerJacobiPolynomialsFO<7, 7, ALPHA0, BETA>::Eval (x, values); break; 01075 case 8: DubinerJacobiPolynomialsFO<8, 8, ALPHA0, BETA>::Eval (x, values); break; 01076 case 9: DubinerJacobiPolynomialsFO<9, 9, ALPHA0, BETA>::Eval (x, values); break; 01077 case 10: DubinerJacobiPolynomialsFO<10, 10, ALPHA0, BETA>::Eval (x, values); break; 01078 default: DubinerJacobiPolynomials1 (n, x, ALPHA0, BETA, values); 01079 } 01080 } 01081 01082 01083 01084 template <int ALPHA0, int BETA, int DIAG, int ORDER = DIAG> 01085 class DubinerJacobiPolynomialsDiag 01086 { 01087 public: 01088 template<class S, class Thelp, class T> 01089 ALWAYS_INLINE static void Step (S x, Thelp & help, T & values) 01090 { 01091 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, DIAG, ORDER-1>::Step (x, help, values); 01092 typedef JacobiPolynomialFix<ALPHA0+2*(DIAG-ORDER), BETA> REC; 01093 01094 if (ORDER == 0) 01095 help[DIAG-ORDER][0] = REC::P0(x); 01096 else if (ORDER == 1) 01097 { 01098 help[DIAG-ORDER][0] = REC::P1(x); 01099 help[DIAG-ORDER][1] = REC::P0(x); 01100 } 01101 else 01102 { 01103 REC::EvalNext (ORDER, x, help[DIAG-ORDER][0], help[DIAG-ORDER][1]); 01104 } 01105 values(DIAG-ORDER, ORDER) = help[DIAG-ORDER][0]; 01106 } 01107 }; 01108 01109 template <int ALPHA0, int BETA, int DIAG> 01110 class DubinerJacobiPolynomialsDiag<ALPHA0, BETA, DIAG, -1> 01111 { 01112 public: 01113 template<class S, class Thelp, class T> 01114 ALWAYS_INLINE static void Step (S x, Thelp & help, T & values) {;} 01115 }; 01116 01117 01118 template <int ALPHA0, int BETA, class S, class T> 01119 void DubinerJacobiPolynomials (int n, S x, T & values) 01120 { 01121 S help[20][2]; 01122 if (n < 0) return; 01123 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 0>::Step (x, help, values); 01124 01125 if (n < 1) return; 01126 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 1>::Step (x, help, values); 01127 01128 if (n < 2) return; 01129 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 2>::Step (x, help, values); 01130 01131 if (n < 3) return; 01132 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 3>::Step (x, help, values); 01133 01134 if (n < 4) return; 01135 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 4>::Step (x, help, values); 01136 01137 if (n < 5) return; 01138 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 5>::Step (x, help, values); 01139 01140 if (n < 6) return; 01141 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 6>::Step (x, help, values); 01142 01143 if (n < 7) return; 01144 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 7>::Step (x, help, values); 01145 01146 if (n < 8) return; 01147 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 8>::Step (x, help, values); 01148 01149 if (n < 9) return; 01150 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 9>::Step (x, help, values); 01151 01152 if (n < 10) return; 01153 DubinerJacobiPolynomialsDiag<ALPHA0, BETA, 10>::Step (x, help, values); 01154 01155 if (n < 11) return; 01156 01157 DubinerJacobiPolynomials1 (n, x, ALPHA0, BETA, values); 01158 } 01159 01160 01161 01162 01163 01164 01165 template <class S, class St, class T> 01166 void DubinerJacobiPolynomialsScaled1 (int n, S x, St t, int alpha0, int beta, T & values) 01167 { 01168 for (int i = 0; i <= n; i++) 01169 ScaledJacobiPolynomial (n-i, x, t, alpha0+2*i, beta, values.Row(i)); 01170 } 01171 01172 01173 01174 01175 template <int ALPHA0, int BETA, class S, class St, class T> 01176 void DubinerJacobiPolynomialsScaled (int n, S x, St t, T & values) 01177 { 01178 switch (n) 01179 { 01180 case 0: DubinerJacobiPolynomialsFO<0, 0, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01181 case 1: DubinerJacobiPolynomialsFO<1, 1, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01182 case 2: DubinerJacobiPolynomialsFO<2, 2, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01183 case 3: DubinerJacobiPolynomialsFO<3, 3, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01184 case 4: DubinerJacobiPolynomialsFO<4, 4, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01185 case 5: DubinerJacobiPolynomialsFO<5, 5, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01186 case 6: DubinerJacobiPolynomialsFO<6, 6, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01187 case 7: DubinerJacobiPolynomialsFO<7, 7, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01188 case 8: DubinerJacobiPolynomialsFO<8, 8, ALPHA0, BETA>::EvalScaled (x, t, values); break; 01189 default: DubinerJacobiPolynomialsScaled1 (n, x, t, ALPHA0, BETA, values); 01190 } 01191 } 01192 01193 01194 01195 01196 01197 01198 01199 01200 01201 01202 01203 01204 01205 01206 01207 01208 01209 class DubinerBasis 01210 { 01211 public: 01212 template <class S, class T> 01213 static void Eval (int n, S x, S y, T & values) 01214 { 01215 EvalMult (n, x, y, 1.0, values); 01216 } 01217 01218 template <class S, class Sc, class T> 01219 static void EvalMult (int n, S x, S y, Sc c, T & values) 01220 { 01221 ArrayMem<S, 20> poly(n+1); 01222 ArrayMem<S, 400> polx_mem( sqr(n+1) ); 01223 FlatMatrix<S> polx(n+1,n+1, &polx_mem[0]); 01224 01225 LegendrePolynomial leg; 01226 leg.EvalScaledMult (n, 2*y+x-1, 1-x, c, poly); 01227 01228 DubinerJacobiPolynomials<1,0> (n, 2*x-1, polx); 01229 01230 for (int i = 0, ii = 0; i <= n; i++) 01231 for (int j = 0; j <= n-i; j++) 01232 values[ii++] = poly[j] * polx(j, i); 01233 } 01234 01235 01236 template <class S, class St, class Sc, class T> 01237 static void EvalScaledMult (int n, S x, S y, St t, Sc c, T & values) 01238 { 01239 ArrayMem<S, 20> poly(n+1); 01240 ArrayMem<S, 400> polx_mem( sqr(n+1) ); 01241 FlatMatrix<S> polx(n+1,n+1, &polx_mem[0]); 01242 01243 LegendrePolynomial leg; 01244 leg.EvalScaledMult (n, 2*y+x-1, t-x, c, poly); 01245 01246 DubinerJacobiPolynomialsScaled<1,0> (n, 2*x-1, t, polx); 01247 01248 for (int i = 0, ii = 0; i <= n; i++) 01249 for (int j = 0; j <= n-i; j++) 01250 values[ii++] = poly[j] * polx(j, i); 01251 } 01252 01253 01254 }; 01255 01256 01257 01258 01259 01260 01261 01262 01263 01264 01265 01266 template <class S, class T> 01267 inline void GegenbauerPolynomial (int n, S x, double lam, T & values) 01268 { 01269 S p1 = 1.0, p2 = 0.0, p3; 01270 01271 if (n >= 0) 01272 values[0] = 1.0; 01273 01274 for (int j=1; j<=n; j++) 01275 { 01276 p3=p2; p2=p1; 01277 p1=( 2.0*(j+lam-1.0) * x * p2 - (j+2*lam-2.0) * p3) / j; 01278 values[j] = p1; 01279 } 01280 } 01281 01282 01292 template <class S, class T> 01293 inline void IntegratedLegendrePolynomial (int n, S x, T & values) 01294 { 01295 S p1 = -1.0; 01296 S p2 = 0.0; 01297 S p3; 01298 01299 if (n >= 0) 01300 values[0] = (S)-1.0; 01301 01302 for (int j=1; j<=n; j++) 01303 { 01304 p3=p2; p2=p1; 01305 p1=( (2*j-3) * x * p2 - (j-3) * p3) / j; 01306 values[j] = p1; 01307 01308 } 01309 } 01310 01311 01312 #ifdef USEDXXX 01313 01317 template <class S, class T> 01318 inline void DerivedLegendrePolynomial (int n, S x, T & values) 01319 { 01320 GegenbauerPolynomial<S,T> (n, x, 1.5, values); 01321 } 01322 #endif 01323 01324 01325 01326 01327 01328 01341 template <class S, class T> 01342 inline void HermitePolynomial (int n, S x, T & values) 01343 { 01344 S p1, p2, p3; 01345 01346 p2 = 0; 01347 if (n >= 0) 01348 p1 = values[0] = 1.0; 01349 for (int j=1; j<=n; j++) 01350 { 01351 p3 = p2; p2 = p1; 01352 p1 = 2*x*p2 - 2*(j-1)*p3; 01353 values[j] = p1; 01354 } 01355 } 01356 01357 01358 01359 01360 01361 01362 01363 01364 01365 01366 01367 01382 template <class Sx, class Sy, class T> 01383 inline void TriangleExtensionMonomial (int n, Sx x, Sy y, T & values) 01384 { 01385 Sx p1 = -1.0, p2 = 0.0, p3; 01386 values[0] = -1.0; 01387 Sy fy = (1-y)*(1-y); 01388 for (int j=1; j<=n; j++) 01389 { 01390 p3=p2; p2=p1; 01391 p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j; 01392 values[j] = p1; 01393 } 01394 } 01395 01396 template <class Sx, class Sy, class T> 01397 inline void DiffTriangleExtensionMonomial (int n, Sx x, Sy y, T & values) 01398 { 01399 Array<AutoDiff<2> > ad_values(n+1); 01400 AutoDiff<2> ad_x(x, 0); 01401 AutoDiff<2> ad_y(y, 1); 01402 01403 TriangleExtensionMonomial (n, ad_x, ad_y, ad_values); 01404 01405 for (int i = 0; i <= n; i++) 01406 for (int j = 0; j < 2; j++) 01407 values(i,j) = ad_values[i].DValue(j); 01408 } 01409 01410 01411 01415 template <class Sx, class Sy, class T> 01416 inline void TriangleExtensionJacobi (int n, Sx x, Sy y, T & values) 01417 { 01418 if ( (1-y) != 0.0) 01419 { 01420 int j; 01421 01422 JacobiPolynomial (n-2, x / (1-y), 2, 2, values); 01423 Sy fac = (1.-x-y) * (1.+x-y); 01424 for (j = 0; j <= n-2; j++) 01425 { 01426 values[j] *= fac; 01427 fac *= 1-y; 01428 } 01429 for (j = n; j >= 2; j--) 01430 values[j] = values[j-2]; 01431 if (n >= 0) values[0] = 0; 01432 if (n >= 1) values[1] = 0; 01433 } 01434 else 01435 { 01436 for (int j = 0; j <= n; j++) 01437 values[j] = 0; 01438 } 01439 } 01440 01441 template <class Sx, class Sy, class T> 01442 inline void DiffTriangleExtensionJacobi (int n, Sx x, Sy y, T & values) 01443 { 01444 Array<AutoDiff<2> > ad_values(n+1); 01445 AutoDiff<2> ad_x(x, 0); 01446 AutoDiff<2> ad_y(y, 1); 01447 01448 TriangleExtensionJacobi (n, ad_x, ad_y, ad_values); 01449 for (int i = 0; i <= n; i++) 01450 for (int j = 0; j < 2; j++) 01451 values(i,j) = ad_values[i].DValue(j); 01452 } 01453 01454 01455 01456 01457 01458 01459 01463 template <class Sx, class Sy, class T> 01464 inline void TriangleExtensionOpt (int n, Sx x, Sy y, T & values) 01465 { 01466 if (y < 1e-10) 01467 { 01468 IntegratedLegendrePolynomial (n, x, values); 01469 } 01470 else 01471 { 01472 Array<Sx> ge1(n+2); 01473 Array<Sx> ge2(n+2); 01474 Array<Sx> ge3(n+2); 01475 Array<Sx> ge4(n+2); 01476 01477 GegenbauerPolynomial (n+1, Sx(-1.0), -1.5, ge1); 01478 GegenbauerPolynomial (n+1, x-y, -1.5, ge2); 01479 GegenbauerPolynomial (n+1, x+y, -1.5, ge3); 01480 GegenbauerPolynomial (n+1, Sx(1.0), -1.5, ge4); 01481 01482 for (int i = 0; i <= n; i++) 01483 values[i] = 1.0/3.0 * 01484 ( (2*y/(1+x+y)/(1+x+y) - y/2) * ge1[i+1] + 01485 (-1/(2*y) + 2*y/(1-x+y)/(1-x+y)) * ge2[i+1] + 01486 (1/(2*y) - 2*y/(1+x+y)/(1+x+y) ) * ge3[i+1] + 01487 (-2*y/(1-x+y)/(1-x+y) + y/2 ) * ge4[i+1] ); 01488 } 01489 } 01490 01491 template <class Sx, class Sy, class T> 01492 inline void DiffTriangleExtensionOpt (int n, Sx x, Sy y, T & values) 01493 { 01494 Array<AutoDiff<2> > ad_values(n+1); 01495 AutoDiff<2> ad_x(x, 0); 01496 AutoDiff<2> ad_y(y, 1); 01497 01498 if (y < 1e-10) 01499 { 01500 values = 0.; 01501 } 01502 else 01503 { 01504 TriangleExtensionOpt (n, ad_x, ad_y, ad_values); 01505 01506 for (int i = 0; i <= n; i++) 01507 for (int j = 0; j < 2; j++) 01508 values(i,j) = ad_values[i].DValue(j); 01509 } 01510 } 01511 01512 01513 01514 template <class S1, class S2, class S3> 01515 inline void StdOp (S1 & v1, const S2 & tt, const S3 & v2, double fac) 01516 { 01517 v1 = fac * (v1*tt - v2) + v2; 01518 // v1 = fac * (v1*tt) + (1-fac) * v2; 01519 } 01520 01521 template <int D> 01522 inline void StdOp (AutoDiff<D> & v1, const AutoDiff<D> & tt, const AutoDiff<D> & v2, double fac) 01523 { 01524 for (int j = 0; j < D; j++) 01525 v1.DValue(j) = fac * (v1.DValue(j) * tt.Value() + v1.Value() * tt.DValue(j) - v2.DValue(j)) + v2.DValue(j); 01526 v1.Value() = fac * (v1.Value()*tt.Value()-v2.Value()) + v2.Value(); 01527 } 01528 01529 01530 /* 01531 E_i(x,y) = P_i(x/t) * t^i 01532 */ 01533 template <class Sx, class St, class T> 01534 inline void ScaledLegendrePolynomial (int n, Sx x, St t, T & values) 01535 { 01536 // St tt = t*t; 01537 St tt = sqr(t); 01538 01539 Sx p1, p2; 01540 01541 if (n < 0) return; 01542 values[0] = p2 = 1.0; 01543 if (n < 1) return; 01544 values[1] = p1 = x; 01545 if (n < 2) return; 01546 01547 for (int j=2; j < n; j+=2) 01548 { 01549 /* 01550 double invj = 1.0/j; 01551 p2 *= (invj-1) * tt; 01552 p2 += (2-invj) * x * p1; 01553 values[j] = p2; 01554 01555 double invj2 = 1.0/(j+1); 01556 p1 *= (invj2-1) * tt; 01557 p1 += (2-invj2) * x * p2; 01558 values[j+1] = p1; 01559 */ 01560 01561 StdOp (p2, tt, x*p1, 1.0/j-1); 01562 values[j] = p2; 01563 StdOp (p1, tt, x*p2, 1.0/(j+1)-1); 01564 values[j+1] = p1; 01565 } 01566 01567 if (n % 2 == 0) 01568 { 01569 double invn = 1.0/n; 01570 values[n] = (2-invn)*x*p1 - (1-invn) * tt*p2; 01571 } 01572 01573 01574 /* 01575 if (n < 0) return; 01576 values[0] = 1.0; 01577 if (n < 1) return; 01578 values[1] = x; 01579 if (n < 2) return; 01580 values[2] = p2 = 1.5 * x * x - 0.5 * tt; 01581 if (n < 3) return; 01582 values[3] = p1 = (5.0/3.0) * p2 * x - (2.0/3.0) * tt * x; 01583 if (n < 4) return; 01584 01585 for (int j=4; j < n; j+=2) 01586 { 01587 double invj = 1.0/j; 01588 p2 *= (invj-1) * tt; 01589 p2 += (2-invj) * x * p1; 01590 values[j] = p2; 01591 01592 invj = 1.0/(j+1); 01593 p1 *= (invj-1) * tt; 01594 p1 += (2-invj) * x * p2; 01595 values[j+1] = p1; 01596 } 01597 01598 if (n % 2 == 0) 01599 { 01600 double invn = 1.0/n; 01601 values[n] = (2-invn)*x*p1 - (invn-1) * tt*p2; 01602 } 01603 */ 01604 01605 01606 01607 01608 01609 /* 01610 Sx p1 = 1.0, p2 = 0.0, p3; 01611 01612 if (n>=0) values[0] = 1.0; 01613 01614 for (int j=1; j<=n; j++) 01615 { 01616 p3=p2; p2=p1; 01617 p1=((2.0*j-1.0) * x*p2 - tt*(j-1.0)*p3)/j; 01618 values[j] = p1; 01619 } 01620 */ 01621 } 01622 01623 01624 01625 /* 01626 E_i(x,y) = c * P_i(x/t) * t^i 01627 */ 01628 template <class Sx, class St, class Sc, class T> 01629 inline void ScaledLegendrePolynomialMult (int n, Sx x, St t, Sc c, T & values) 01630 { 01631 St tt = sqr(t); 01632 Sx p1, p2; 01633 01634 if (n < 0) return; 01635 values[0] = p2 = c; 01636 if (n < 1) return; 01637 values[1] = p1 = c * x; 01638 if (n < 2) return; 01639 01640 for (int j=2; j < n; j+=2) 01641 { 01642 StdOp (p2, tt, x*p1, 1.0/j-1); 01643 values[j] = p2; 01644 StdOp (p1, tt, x*p2, 1.0/(j+1)-1); 01645 values[j+1] = p1; 01646 /* 01647 double invj = 1.0/j; 01648 p2 *= (invj-1) * tt; 01649 p2 += (2-invj) * x * p1; 01650 values[j] = p2; 01651 01652 invj = 1.0/(j+1); 01653 p1 *= (invj-1) * tt; 01654 p1 += (2-invj) * x * p2; 01655 values[j+1] = p1; 01656 */ 01657 } 01658 01659 if (n % 2 == 0) 01660 { 01661 StdOp (p2, tt, x*p1, 1.0/n-1); 01662 values[n] = p2; 01663 01664 // double invn = 1.0/n; 01665 // values[n] = (2-invn)*x*p1 - (1-invn) * tt*p2; 01666 } 01667 } 01668 01669 01670 01671 01672 01673 01674 01675 01676 01677 01678 01679 01680 01681 01682 01683 01684 01685 01686 01687 template <class T> 01688 inline void DiffScaledLegendrePolynomial (int n, double x, double t, T & values) 01689 { 01690 ArrayMem<AutoDiff<2>,10> ad_values(n+1); 01691 AutoDiff<2> ad_x(x, 0); 01692 AutoDiff<2> ad_t(t, 1); 01693 01694 ScaledLegendrePolynomial(n, ad_x, ad_t, ad_values); 01695 01696 for (int i = 0; i <= n; i++) 01697 for (int j = 0; j < 2; j++) 01698 values(i,j) = ad_values[i].DValue(j); 01699 } 01700 01701 01702 template <class Sx, class St, class T> 01703 inline void ScaledIntegratedLegendrePolynomial (int n, Sx x, St t, T & values) 01704 { 01705 Sx p1 = -1.0; 01706 Sx p2 = 0.0; 01707 Sx p3; 01708 St tt = t*t; 01709 if (n >= 0) 01710 values[0] = -1.0; 01711 01712 for (int j=1; j<=n; j++) 01713 { 01714 p3=p2; p2=p1; 01715 p1=((2.0*j-3.0) * x*p2 - tt*(j-3.0)*p3)/j; 01716 values[j] = p1; 01717 } 01718 } 01719 01720 01721 01722 01723 01724 01725 template <class T> 01726 inline void ScaledLegendrePolynomialandDiff(int n, double x, double t, T & P, T & Px, T & Pt) 01727 { 01728 /* 01729 ArrayMem<AutoDiff<2>,10> ad_values(n+1); 01730 AutoDiff<2> ad_x(x, 0); 01731 AutoDiff<2> ad_t(t, 1); 01732 01733 ScaledLegendrePolynomial(n, ad_x, ad_t, ad_values); 01734 01735 for (int i = 0; i <= n; i++) 01736 { 01737 P[i] = ad_values[i].Value(); 01738 Px[i] = ad_values[i].DValue(0); 01739 Pt[i] = ad_values[i].DValue(1); 01740 } 01741 */ 01742 if(n>=0) 01743 { 01744 P[0] = 1.; 01745 Px[0] = 0.; 01746 Pt[0] = 0.; 01747 if(n>=1) 01748 { 01749 P[1] = x; 01750 Px[1] = 1.; 01751 Pt[1] = 0.; 01752 } 01753 01754 double px0 = 0., px1 = 0., px2 =1.; 01755 double sqt = t*t; 01756 for(int l = 2; l<=n; l++) 01757 { 01758 px0=px1; 01759 px1=px2; 01760 px2= ( (2*l-1)*x*px1 - l*sqt*px0)/(l-1); 01761 Px[l] = px2; 01762 Pt[l] = -t*px1; 01763 P[l] = (x*px2-sqt*px1)/l; 01764 } 01765 } 01766 01767 } 01768 01769 template <class T> 01770 inline void LegendrePolynomialandDiff(int n, double x, T & P, T & Px) 01771 { 01772 /* 01773 ArrayMem<AutoDiff<1>,10> ad_values(n+1); 01774 AutoDiff<1> ad_x(x, 0); 01775 LegendrePolynomial(n, ad_x, ad_values); 01776 01777 for (int i = 0; i <= n; i++) 01778 { 01779 P[i] = ad_values[i].Value(); 01780 Px[i] = ad_values[i].DValue(0); 01781 } 01782 01783 (*testout) << "P = " << endl << P << endl 01784 << "Px = " << endl << Px << endl; 01785 */ 01786 if(n>=0) 01787 { 01788 P[0] = 1.; 01789 Px[0] = 0.; 01790 if(n>=1) 01791 { 01792 P[1] = x; 01793 Px[1] = 1.; 01794 } 01795 double px0 = 0., px1 = 0., px2 =1.; 01796 for(int l = 2; l<=n; l++) 01797 { 01798 px0=px1; 01799 px1=px2; 01800 01801 px2= ( (2*l-1)*x*px1 - l*px0)/(l-1); 01802 Px[l] = px2; 01803 P[l] = (x*px2 - px1)/l; 01804 } 01805 } 01806 } 01807 01808 01809 01810 01811 #ifdef OLD 01812 01813 01814 /* 01815 u(0) = 0, u(1) = 1, 01816 min \int_0^1 (1-x)^{DIM-1} (u')^2 dx 01817 01818 representation as 01819 \sum c_i P_i(2x-1) 01820 */ 01821 01822 template <int DIM> 01823 class VertexExtensionOptimal 01824 { 01825 enum { SIZE = 50 }; 01826 static double coefs[SIZE][SIZE]; 01827 static bool initialized; 01828 public: 01829 01830 VertexExtensionOptimal (); 01831 01832 template <typename Tx> 01833 inline static Tx Calc (int p, Tx x) 01834 { 01835 Tx p1 = 1.0, p2 = 0.0, p3; 01836 Tx sum = 0; 01837 01838 x = 2*x-1; 01839 01840 if (p >= 0) 01841 sum += coefs[0][p]; 01842 01843 for (int j=1; j<=p; j++) 01844 { 01845 p3 = p2; p2 = p1; 01846 p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j; 01847 sum += coefs[j][p] * p1; 01848 } 01849 01850 return sum; 01851 } 01852 01853 inline static double CalcDeriv (int p, double x) 01854 { 01855 AutoDiff<1> p1 = 1.0, p2 = 0.0, p3; 01856 AutoDiff<1> sum = 0; 01857 AutoDiff<1> adx (x, 0); // \nabla adx = e_0 01858 01859 adx = 2.0*adx-1; 01860 01861 if (p >= 0) 01862 sum += coefs[0][p]; 01863 01864 for (int j=1; j<=p; j++) 01865 { 01866 p3 = p2; p2 = p1; 01867 p1 = ((2.0*j-1.0)*adx*p2 - (j-1.0)*p3) / j; 01868 sum += coefs[j][p] * p1; 01869 } 01870 01871 return sum.DValue(0); 01872 } 01873 01874 01875 /* 01876 // Based on Jacobi pols, 3D only 01877 template <typename Tx> 01878 inline static Tx Calc (int p, Tx x) 01879 { 01880 ArrayMem<Tx,20> jacpol(p+1); 01881 01882 JacobiPolynomial (p, 2*x-1, 1, -1, jacpol); 01883 01884 Tx sum = 0; 01885 for (int j = 0; j <= p; j++) 01886 sum += coefs[j][p] * jacpol[j]; 01887 01888 return sum; 01889 } 01890 01891 inline static double CalcDeriv (int p, double x) 01892 { 01893 ArrayMem<double,20> jacpol(p+1); 01894 01895 JacobiPolynomial (p, 2*x-1, 2, 0, jacpol); 01896 01897 double sum = 0; 01898 for (int j = 1; j <= p; j++) 01899 sum += coefs[j][p] * 0.5 * (j+1) * jacpol[j-1]; 01900 01901 return sum; 01902 } 01903 */ 01904 }; 01905 01906 01907 01908 01909 /* 01910 u(-1) = 0, u(1) = 1, 01911 min \int_-1^1 (1-x)^{DIM-1} (u')^2 dx 01912 */ 01913 01914 template <class Tx, class T> 01915 inline void LowEnergyVertexPolynomials2D (int n, Tx x, T & values) 01916 { 01917 JacobiPolynomial (n, x, 0, -1, values); 01918 Tx sum1 = 0.0, sum2 = 0.0; 01919 for (int i = 1; i <= n; i++) 01920 { 01921 sum1 += 1.0/i; 01922 sum2 += values[i] / i; 01923 values[i] = sum2/sum1; 01924 } 01925 values[0] = 1; 01926 } 01927 01928 template <class Tx, class T> 01929 inline void LowEnergyVertexPolynomials3D (int n, Tx x, T & values) 01930 { 01931 JacobiPolynomial (n, x, 1, -1, values); 01932 Tx sum = 0.0; 01933 for (int i = 1; i <= n; i++) 01934 { 01935 sum += (2.0*i+1)/(i+1) * values[i]; 01936 values[i] = 1.0/(i*(i+2)) * sum; 01937 } 01938 values[0] = 1; 01939 } 01940 01941 01942 01943 01944 01945 class VertexStandard 01946 { 01947 public: 01948 template <typename Tx> 01949 inline static Tx Calc (int p, Tx x) 01950 { 01951 return x; 01952 } 01953 01954 inline static double CalcDeriv (int p, double x) 01955 { 01956 return 1; 01957 } 01958 }; 01959 01960 #endif 01961 01962 01977 class IntegratedLegendreMonomialExt 01978 { 01979 enum { SIZE = 1000 }; 01980 static double coefs[SIZE][2]; 01981 01982 public: 01983 01984 static void CalcCoeffs () 01985 { 01986 for (int j = 1; j < SIZE; j++) 01987 { 01988 coefs[j][0] = double(2*j-3)/j; 01989 coefs[j][1] = double(j-3)/j; 01990 } 01991 } 01992 01993 template <class Sx, class Sy, class T> 01994 inline static int CalcScaled (int n, Sx x, Sy y, T & values) 01995 { 01996 Sy fy = y*y; 01997 Sx p3 = 0; 01998 Sx p2 = -1; 01999 Sx p1 = x; 02000 02001 for (int j=2; j<=n; j++) 02002 { 02003 p3=p2; p2=p1; 02004 p1 = double(2*j-3)/j * x * p2 - double(j-3)/j * fy * p3; 02005 values[j-2] = p1; 02006 } 02007 02008 return n-1; 02009 } 02010 02011 02012 template <int n, class Sx, class Sy, class T> 02013 inline static int CalcScaled (Sx x, Sy y, T & values) 02014 { 02015 Sy fy = y*y; 02016 Sx p3 = 0; 02017 Sx p2 = -1; 02018 Sx p1 = x; 02019 02020 for (int j=2; j<=n; j++) 02021 { 02022 p3=p2; p2=p1; 02023 p1 = double(2*j-3)/j * x * p2 - double(j-3)/j * fy * p3; 02024 values[j-2] = p1; 02025 } 02026 02027 return n-1; 02028 } 02029 02030 02031 02032 02033 02034 template <class Sx, class Sy, class T> 02035 inline static int CalcTrigExt (int n, Sx x, Sy y, T & values) 02036 { 02037 Sy fy = (1-y)*(1-y); 02038 Sx p3 = 0; 02039 Sx p2 = -1; 02040 Sx p1 = x; 02041 02042 for (int j=2; j<=n; j++) 02043 { 02044 p3=p2; p2=p1; 02045 p1 = double(2*j-3)/j * x * p2 - double(j-3)/j * fy * p3; 02046 values[j-2] = p1; 02047 } 02048 02049 return n-1; 02050 } 02051 02052 template <class Sx, class Sy, class Sf, class T> 02053 inline static int CalcTrigExtMult (int n, Sx x, Sy y, Sf fac, T & values) 02054 { 02055 Sy fy = (1-y)*(1-y); 02056 Sx p3 = 0; 02057 Sx p2 = -fac; 02058 Sx p1 = x * fac; 02059 02060 for (int j=2; j<=n; j++) 02061 { 02062 p3=p2; p2=p1; 02063 // p1=( (2*j-3) * x * p2 - (j-3) * fy * p3) / j; 02064 p1 = double(2*j-3)/j * x * p2 - double(j-3)/j * fy * p3; 02065 // p1= coefs[j][0] * x * p2 - coefs[j][1] * fy * p3; 02066 values[j-2] = p1; 02067 } 02068 02069 return n-1; 02070 } 02071 02072 02073 02074 02075 template <class T> 02076 inline static int CalcTrigExtDeriv (int n, double x, double y, T & values) 02077 { 02078 double fy = (1-y)*(1-y); 02079 double p3 = 0, p3x = 0, p3y = 0; 02080 double p2 = -1, p2x = 0, p2y = 0; 02081 double p1 = x, p1x = 1, p1y = 0; 02082 02083 for (int j=2; j<=n; j++) 02084 { 02085 p3=p2; p3x = p2x; p3y = p2y; 02086 p2=p1; p2x = p1x; p2y = p1y; 02087 double c1 = (2.0*j-3) / j; 02088 double c2 = (j-3.0) / j; 02089 02090 p1 = c1 * x * p2 - c2 * fy * p3; 02091 p1x = c1 * p2 + c1 * x * p2x - c2 * fy * p3x; 02092 p1y = c1 * x * p2y - (c2 * 2 * (y-1) * p3 + c2 * fy * p3y); 02093 values (j-2, 0) = p1x; 02094 values (j-2, 1) = p1y; 02095 } 02096 return n-1; 02097 } 02098 02099 02100 template <class Sx, class T> 02101 inline static int Calc (int n, Sx x, T & values) 02102 { 02103 Sx p3 = 0; 02104 Sx p2 = -1; 02105 Sx p1 = x; 02106 02107 for (int j=2; j<=n; j++) 02108 { 02109 p3=p2; p2=p1; 02110 p1=( (2*j-3) * x * p2 - (j-3) * p3) / j; 02111 values[j-2] = p1; 02112 } 02113 return n-1; 02114 } 02115 02116 template <class Sx, class Sf, class T> 02117 inline static int CalcMult (int n, Sx x, Sf fac, T & values) 02118 { 02119 Sx p3 = 0; 02120 Sx p2 = -fac; 02121 Sx p1 = x*fac; 02122 02123 for (int j=2; j<=n; j++) 02124 { 02125 p3=p2; p2=p1; 02126 p1=( (2*j-3) * x * p2 - (j-3) * p3) / j; 02127 values[j-2] = p1; 02128 } 02129 return n-1; 02130 } 02131 02132 02133 02134 template <class T> 02135 inline static int CalcDeriv (int n, double x, T & values) 02136 { 02137 double p1 = 1.0, p2 = 0.0, p3; 02138 02139 for (int j=1; j<=n-1; j++) 02140 { 02141 p3 = p2; p2 = p1; 02142 p1 = ((2.0*j-1.0)*x*p2 - (j-1.0)*p3) / j; 02143 values[j-1] = p1; 02144 } 02145 return n-1; 02146 } 02147 02148 02149 }; 02150 02151 02152 02153 02154 02155 02156 /* Conversion of orthogonal polynomials */ 02157 02158 // given: \sum u_i P^{al,0} 02159 // find: \sum v_i P^{al-1, 0} 02160 template <class T> 02161 void ConvertJacobiReduceAlpha (int n, int alpha, T & inout) 02162 { 02163 for (int i = n; i >= 1; i--) 02164 { 02165 double val = inout(i) / (i+alpha); 02166 inout(i) = (2*i+alpha) * val; 02167 inout(i-1) += i * val; 02168 } 02169 } 02170 02171 // given: \sum u_i (1-x) P^{al+1,0} 0 <= i < n 02172 // find: \sum v_i P^{al, 0} 0 <= i <= n 02173 template <class T> 02174 void ConvertJacobiReduceAlphaFactor (int n, double alpha, T & inout) 02175 { 02176 inout(n) = 0; 02177 for (int i = n; i > 0; i--) 02178 { 02179 double val = inout(i-1) / (i+alpha/2); 02180 inout(i-1) = (i+alpha) * val; 02181 inout(i) -= i * val; 02182 } 02183 } 02184 02185 02186 02187 02188 02189 // Differentiate Jacobi 02190 // (P_i^alpha)' 02191 02192 template <class T> 02193 void DifferentiateJacobi (int n, double alpha, T & inout) 02194 { 02195 for (int i = n; i >= 1; i--) 02196 { 02197 double val = inout(i); 02198 inout(i-1) -= alpha*(2*i+alpha-1) / ( (i+alpha)*(2*i+alpha-2)) * val; 02199 if (i > 1) 02200 inout(i-2) += (i-1)*(2*i+alpha) / ( (i+alpha)*(2*i+alpha-2)) * val; 02201 02202 inout(i) *= (2*i+alpha)*(2*i+alpha-1) / ( 2 * (i+alpha) ); 02203 } 02204 for (int i = 0; i < n; i++) 02205 inout(i) = inout(i+1); 02206 inout(n) = 0; 02207 } 02208 02209 02210 template <class T> 02211 void DifferentiateJacobiTrans (int n, double alpha, T & inout) 02212 { 02213 for (int i = n-1; i >= 0; i--) 02214 inout(i+1) = inout(i); 02215 inout(0) = 0; 02216 02217 for (int i = 1; i <= n; i++) 02218 { 02219 inout(i) *= (2*i+alpha)*(2*i+alpha-1) / ( 2 * (i+alpha) ); 02220 02221 inout(i) -= alpha*(2*i+alpha-1) / ( (i+alpha)*(2*i+alpha-2)) * inout(i-1); 02222 if (i > 1) 02223 inout(i) += (i-1)*(2*i+alpha) / ( (i+alpha)*(2*i+alpha-2)) * inout(i-2); 02224 } 02225 } 02226 02227 02228 template <class T> 02229 void DifferentiateLegendre (int n, T & inout) 02230 { 02231 for (int i = n; i >= 1; i--) 02232 { 02233 if (i > 1) inout(i-2) += inout(i); 02234 inout(i) *= 2*i-1; 02235 } 02236 for (int i = 0; i < n; i++) 02237 inout(i) = inout(i+1); 02238 inout(n) = 0; 02239 } 02240 02241 template <class T> 02242 void DifferentiateLegendreTrans (int n, T & inout) 02243 { 02244 for (int i = n-1; i >= 0; i--) 02245 inout(i+1) = inout(i); 02246 inout(0) = 0; 02247 02248 for (int i = 1; i <= n; i++) 02249 { 02250 inout(i) *= (2*i-1); 02251 if (i > 1) 02252 inout(i) += inout(i-2); 02253 } 02254 } 02255 02256 02257 02258 02259 02260 02261 class ConvertJacobi 02262 { 02263 typedef double d2[2]; 02264 static Array<d2*> coefs_increasealpha; 02265 static Array<d2*> coefs_reducealpha; 02266 static Array<d2*> coefs_reducealphafac; 02267 02268 public: 02269 ConvertJacobi (); 02270 ~ConvertJacobi (); 02271 02272 template <class T> 02273 static void ReduceAlpha (int n, int alpha, T & inout) // alpha of input 02274 { 02275 d2 * c = coefs_reducealpha[alpha]; 02276 02277 double val = inout[n]; 02278 for (int i = n; i >= 1; i--) 02279 { 02280 inout[i] = c[i][1] * val; 02281 val = c[i][0] * val + inout[i-1]; 02282 } 02283 inout[0] = val; 02284 02285 /* 02286 for (int i = n; i >= 1; i--) 02287 { 02288 inout[i-1] += c[i][0] * inout[i]; 02289 inout[i] *= c[i][1]; 02290 } 02291 */ 02292 } 02293 02294 template <class T> 02295 static void ReduceAlphaFactor (int n, int alpha, T & inout) // alpha of output 02296 { 02297 d2 * c = coefs_reducealphafac[alpha]; 02298 02299 inout(n) = c[n][0] * inout[n-1]; 02300 for (int i = n-1; i >= 1; i--) 02301 inout[i] = c[i+1][1] * inout[i] + c[i][0] * inout[i-1]; 02302 inout[0] = c[1][1] * inout[0]; 02303 02304 /* 02305 for (int i = n; i >= 1; i--) 02306 { 02307 inout[i] += c[i][0] * inout[i-1]; 02308 inout[i-1] *= c[i][1]; 02309 } 02310 */ 02311 } 02312 02313 02314 template <class T> 02315 static void ReduceAlphaTrans (int n, int alpha, T & inout) // alpha of input 02316 { 02317 d2 * c = coefs_reducealpha[alpha]; 02318 02319 /* 02320 for (int i = n; i >= 1; i--) 02321 { 02322 inout[i-1] += c[i][0] * inout[i]; 02323 inout[i] *= c[i][1]; 02324 } 02325 */ 02326 for (int i = 1; i <= n; i++) 02327 { 02328 inout[i] *= c[i][1]; 02329 inout[i] += c[i][0] * inout[i-1]; 02330 } 02331 02332 } 02333 02334 template <class T> 02335 static void ReduceAlphaFactorTrans (int n, int alpha, T & inout) // alpha of output 02336 { 02337 d2 * c = coefs_reducealphafac[alpha]; 02338 /* 02339 for (int i = n; i >= 1; i--) 02340 { 02341 inout[i] += c[i][0] * inout[i-1]; 02342 inout[i-1] *= c[i][1]; 02343 } 02344 */ 02345 for (int i = 1; i <= n; i++) 02346 { 02347 inout[i-1] *= c[i][1]; 02348 inout[i-1] += c[i][0] * inout[i]; 02349 } 02350 inout[n] = 0; 02351 } 02352 02353 02354 02355 02356 02357 02358 02359 02360 02361 02362 02363 02364 02365 02366 02367 02368 02369 // reduce, fac 02370 // P_i^alpha(x) (1-x)/2 = c_i^{alpha} P_i^{\alpha-1} + hatc_i^{alpha} P_{i+1}^{\alpha-1} 02371 static double c (int i, int alpha) { return double(i+alpha)/double(2*i+alpha+1); } 02372 static double hatc (int i, int alpha) { return -double(i+1)/double(2*i+alpha+1); } 02373 02374 // increase alpha 02375 // P_i^alpha(x) = d_i^{alpha} P_i^{\alpha+1} + hatd_i^{alpha} P_{i-1}^{\alpha+1} 02376 static double d (int i, int alpha) { return double(i+alpha+1)/double(2*i+alpha+1); } 02377 static double hatd (int i, int alpha) { return -double(i)/double(2*i+alpha+1); } 02378 02379 // decrease alpha 02380 // P_i^alpha(x) = e_i^{alpha} P_i^{\alpha-1} + hate_i^{alpha} P_{i-1}^{\alpha} 02381 static double e (int i, int alpha) { return double(2*i+alpha)/double(i+alpha); } 02382 static double hate (int i, int alpha) { return double(i)/double(i+alpha); } 02383 02384 02385 static Array<d2*> coefs_c, coefs_d, coefs_e; 02386 02387 // alpha,beta of input, 02388 // order of input 02389 // reduce alpha, reduce beta, reduce factor (1-x)/2 (1-y)/2 02390 template <class T> 02391 static void TriangularReduceFactor (int order, int alpha, int beta, T & inout) 02392 { 02393 for (int i = 0; i <= order+1; i++) 02394 inout(i, order+1-i) = 0; 02395 02396 for (int j = order; j >= 0; j--) 02397 { 02398 d2 * calpha = coefs_c[alpha+2*j]; 02399 d2 * cbeta = coefs_c[beta]; 02400 d2 * dalpha = coefs_d[alpha+2*j]; 02401 02402 for (int i = order-j; i >= 0; i--) 02403 { 02404 double val = inout(i,j); 02405 inout(i,j) = val * calpha[i][0] * cbeta[j][0]; 02406 inout(i+1,j) += val * calpha[i][1] * cbeta[j][0]; 02407 inout(i,j+1) += val * dalpha[i][0] * cbeta[j][1]; 02408 if (i > 0) 02409 inout(i-1,j+1) += val * dalpha[i][1] * cbeta[j][1]; 02410 } 02411 02412 /* 02413 for (int i = order-j; i >= 0; i--) 02414 { 02415 double val = inout(i,j); 02416 inout(i,j) = val * c(i,alpha+2*j) * c(j,beta); 02417 inout(i+1,j) += val * hatc(i, alpha+2*j) * c(j, beta); 02418 inout(i,j+1) += val * d(i, alpha+2*j) * hatc (j,beta); 02419 if (i > 0) 02420 inout(i-1,j+1) += val * hatd(i, alpha+2*j) * hatc(j,beta); 02421 } 02422 */ 02423 } 02424 } 02425 02426 // alpha,beta of input, 02427 // const alpha, dec beta 02428 // order is constant 02429 template <class T, class S> 02430 static void TriangularReduceBeta (int order, int alpha, int beta, T & inout, S & hv) 02431 { 02432 d2 * ebeta = coefs_e[beta]; 02433 02434 for (int j = order; j > 0; j--) 02435 { 02436 d2 * calpha = coefs_c[alpha+2*j]; 02437 02438 hv = 0.0; 02439 for (int i = 0; i <= order-j; i++) 02440 { 02441 double val = inout(i,j); 02442 inout(i,j) = val * ebeta[j][0]; 02443 02444 hv(i) += val * calpha[i][0] * ebeta[j][1]; 02445 hv(i+1) += val * calpha[i][1] * ebeta[j][1]; 02446 } 02447 02448 ReduceAlpha (order-j+1, alpha+2*j-1, hv); 02449 inout.Col(j-1) += hv; 02450 } 02451 } 02452 }; 02453 02454 02455 02456 02457 02458 02459 02460 02461 02462 02463 02464 02465 // template meta-programming 02466 02467 template <int N, int AL> 02468 class TReduceAlpha 02469 { 02470 public: 02471 enum { FLOP = TReduceAlpha<N-1,AL>::FLOP + 2 }; 02472 02473 template <class T> 02474 static ALWAYS_INLINE void Do (T & inout) 02475 { 02476 inout[N-1] += double(N)/double(N+AL) * inout[N]; 02477 inout[N] *= double(2*N+AL)/double(N+AL); 02478 TReduceAlpha<N-1,AL>::Do(inout); 02479 } 02480 02481 template <class T> 02482 static ALWAYS_INLINE void Trans (T & inout) 02483 { 02484 TReduceAlpha<N-1,AL>::Trans(inout); 02485 inout[N] *= double(2*N+AL)/double(N+AL); 02486 inout[N] += double(N)/double(N+AL) * inout[N-1]; 02487 } 02488 }; 02489 02490 template <int AL> 02491 class TReduceAlpha<0,AL> 02492 { 02493 public: 02494 enum { FLOP = 0 }; 02495 template <class T> 02496 static void Do (T & inout) { ; } 02497 02498 template <class T> 02499 static void Trans (T & inout) { ; } 02500 }; 02501 02502 02503 02504 02505 02506 template <int N, int AL, int HIGHEST = 1> 02507 class TReduceAlphaFactor 02508 { 02509 public: 02510 enum { FLOP = TReduceAlphaFactor<N-1,AL>::FLOP + 2 }; 02511 02512 template <class T> 02513 static ALWAYS_INLINE void Do (T & inout) 02514 { 02515 if (HIGHEST) 02516 inout[N] = double(-N)/double(2*N+AL) * inout[N-1]; 02517 else 02518 inout[N] += double(-N)/double(2*N+AL) * inout[N-1]; 02519 inout[N-1] *= double(N+AL)/double(2*N+AL); 02520 TReduceAlphaFactor<N-1,AL,0>::Do(inout); 02521 } 02522 02523 template <class T> 02524 static ALWAYS_INLINE void Trans (T & inout) 02525 { 02526 TReduceAlphaFactor<N-1,AL,0>::Trans(inout); 02527 inout[N-1] *= double(N+AL)/double(2*N+AL); 02528 inout[N-1] += double(-N)/double(2*N+AL) * inout[N]; 02529 } 02530 }; 02531 02532 template <int AL, int HIGHEST> 02533 class TReduceAlphaFactor<0,AL,HIGHEST> 02534 { 02535 public: 02536 enum { FLOP = 0 }; 02537 02538 template <class T> 02539 static inline void Do (T & inout) 02540 { 02541 if (HIGHEST) inout[0] = 0.0; 02542 } 02543 02544 template <class T> 02545 static inline void Trans (T & inout) { ; } 02546 }; 02547 02548 02549 02550 02551 02552 02553 /* 02554 template <class T> 02555 void DifferentiateJacobi (int n, double alpha, T & inout) 02556 { 02557 for (int i = n; i >= 1; i--) 02558 { 02559 double val = inout(i); 02560 inout(i-1) -= alpha*(2*i+alpha-1) / ( (i+alpha)*(2*i+alpha-2)) * val; 02561 if (i > 1) 02562 inout(i-2) += (i-1)*(2*i+alpha) / ( (i+alpha)*(2*i+alpha-2)) * val; 02563 02564 inout(i) *= (2*i+alpha)*(2*i+alpha-1) / ( 2 * (i+alpha) ); 02565 } 02566 for (int i = 0; i < n; i++) 02567 inout(i) = inout(i+1); 02568 inout(n) = 0; 02569 } 02570 */ 02571 02572 02573 02574 template <int N, int AL> 02575 class TDifferentiateJacobi 02576 { 02577 public: 02578 enum { FLOP = TDifferentiateJacobi<N-1,AL>::FLOP + 3 }; 02579 02580 // (P_i^Al)' = c1 P_i^AL + c2 (P_{i-1}^AL)' + c3 (P_{i-2}^AL)' 02581 02582 static double c1() { return double ( (2*N+AL)*(2*N+AL-1) ) / double( 2 * (N+AL) ); } 02583 static double c2() { return -double (AL*(2*N+AL-1)) / double( (N+AL)*(2*N+AL-2)); } 02584 static double c3() { return double((N-1)*(2*N+AL)) / double( (N+AL)*(2*N+AL-2)); } 02585 /* 02586 enum { c2 = -double (AL*(2*N+AL-1)) / double( (N+AL)*(2*N+AL-2)) }; 02587 enum { c3 = double((N-1)*(2*N+AL)) / double( (N+AL)*(2*N+AL-2)) }; 02588 */ 02589 template <class T> 02590 static ALWAYS_INLINE void Do (T & inout) 02591 { 02592 double val = inout(N); 02593 inout(N-1) += c2() * val; 02594 if (N > 1) 02595 inout(N-2) += c3() * val; 02596 02597 inout(N) *= c1(); 02598 /* 02599 double val = inout(N); 02600 inout(N-1) -= double (AL*(2*N+AL-1)) / double( (N+AL)*(2*N+AL-2)) * val; 02601 if (N > 1) 02602 inout(N-2) += double((N-1)*(2*N+AL)) / double( (N+AL)*(2*N+AL-2)) * val; 02603 02604 inout(N) *= double ( (2*N+AL)*(2*N+AL-1) ) / double( 2 * (N+AL) ); 02605 */ 02606 TDifferentiateJacobi<N-1,AL>::Do(inout); 02607 02608 inout(N-1) = inout(N); 02609 inout(N) = 0; 02610 } 02611 02612 template <class T> 02613 static ALWAYS_INLINE void Trans (T & inout) 02614 { 02615 inout(N) = inout(N-1); 02616 inout(N-1) = 0; 02617 02618 TDifferentiateJacobi<N-1,AL>::Trans(inout); 02619 02620 inout(N) *= double ((2*N+AL)*(2*N+AL-1)) / ( 2 * (N+AL) ); 02621 02622 inout(N) -= double (AL*(2*N+AL-1)) / ( (N+AL)*(2*N+AL-2)) * inout(N-1); 02623 if (N > 1) 02624 inout(N) += double ((N-1)*(2*N+AL)) / ( (N+AL)*(2*N+AL-2)) * inout(N-2); 02625 02626 /* 02627 double val = inout(N); 02628 inout(N-1) -= double (AL*(2*N+AL-1)) / double( (N+AL)*(2*N+AL-2)) * val; 02629 if (N > 1) 02630 inout(N-2) += double((N-1)*(2*N+AL)) / double( (N+AL)*(2*N+AL-2)) * val; 02631 02632 inout(N) *= double ( (2*N+AL)*(2*N+AL-1) ) / double( 2 * (N+AL) ); 02633 */ 02634 02635 } 02636 02637 02638 }; 02639 02640 template <int AL> 02641 class TDifferentiateJacobi<0,AL> 02642 { 02643 public: 02644 enum { FLOP = 0 }; 02645 02646 template <class T> 02647 static inline void Do (T & inout) { ; } 02648 02649 template <class T> 02650 static inline void Trans (T & inout) { ; } 02651 }; 02652 02653 02654 02655 02656 02657 02658 02659 02660 02661 template <int N> 02662 class TDifferentiateLegendre 02663 { 02664 public: 02665 enum { FLOP = TDifferentiateLegendre<N-1>::FLOP + 3 }; 02666 02667 // (P_i^Al)' = c1 P_i^AL + c2 (P_{i-1}^AL)' + c3 (P_{i-2}^AL)' 02668 02669 template <class T> 02670 static ALWAYS_INLINE void Do (T & inout) 02671 { 02672 double val = inout(N); 02673 if (N > 1) inout(N-2) += val; 02674 02675 inout(N) *= double (2*N-1); 02676 TDifferentiateLegendre<N-1>::Do(inout); 02677 02678 inout(N-1) = inout(N); 02679 inout(N) = 0; 02680 } 02681 02682 template <class T> 02683 static ALWAYS_INLINE void Trans (T & inout) 02684 { 02685 inout(N) = inout(N-1); 02686 inout(N-1) = 0; 02687 02688 TDifferentiateLegendre<N-1>::Trans(inout); 02689 02690 inout(N) *= double (2*N-1); 02691 if (N > 1) inout(N) += inout(N-2); 02692 } 02693 02694 }; 02695 02696 template <> 02697 class TDifferentiateLegendre<0> 02698 { 02699 public: 02700 enum { FLOP = 0 }; 02701 02702 template <class T> 02703 static inline void Do (T & inout) { ; } 02704 02705 template <class T> 02706 static inline void Trans (T & inout) { ; } 02707 }; 02708 02709 02710 02711 02712 02713 02714 02715 02716 02717 02718 02719 02720 02721 02722 02723 02724 02725 02726 02727 02728 template <int N, int J, int I, int AL, int BE> 02729 class TTriangleReduceFactorCol 02730 { 02731 // reduce, fac 02732 // P_i^alpha(x) (1-x)/2 = c_i^{alpha} P_i^{\alpha-1} + hatc_i^{alpha} P_{i+1}^{\alpha-1} 02733 static double c (int i, int alpha) { return double(i+alpha)/double(2*i+alpha+1); } 02734 static double hatc (int i, int alpha) { return -double(i+1)/double(2*i+alpha+1); } 02735 02736 // increase alpha 02737 // P_i^alpha(x) = d_i^{alpha} P_i^{\alpha+1} + hatd_i^{alpha} P_{i-1}^{\alpha+1} 02738 static double d (int i, int alpha) { return double(i+alpha+1)/double(2*i+alpha+1); } 02739 static double hatd (int i, int alpha) { return -double(i)/double(2*i+alpha+1); } 02740 02741 public: 02742 enum { FLOP = TTriangleReduceFactorCol<N,J,I-1,AL,BE>::FLOP + 4 }; 02743 02744 template <class T> 02745 static ALWAYS_INLINE void Do (T & inout) 02746 { 02747 double val = inout(I,J); 02748 02749 inout(I,J) = val * c(I,AL+2*J) * c(J,BE); 02750 inout(I+1,J) += val * hatc(I, AL+2*J) * c(J, BE); 02751 inout(I,J+1) += val * d(I, AL+2*J) * hatc (J, BE); 02752 if (I > 0) 02753 inout(I-1,J+1) += val * hatd(I, AL+2*J) * hatc(J,BE); 02754 02755 TTriangleReduceFactorCol<N,J,I-1,AL,BE> ::Do(inout); 02756 } 02757 02758 02759 template <class T> 02760 static ALWAYS_INLINE void Trans (T & inout) 02761 { 02762 TTriangleReduceFactorCol<N,J,I-1,AL,BE> ::Trans(inout); 02763 02764 double val = 02765 c(I,AL+2*J) * c(J,BE) * inout(I,J) 02766 + hatc(I, AL+2*J) * c(J, BE) * inout(I+1,J) 02767 + d(I, AL+2*J) * hatc (J, BE) * inout(I,J+1); 02768 02769 if (I > 0) 02770 val += hatd(I, AL+2*J) * hatc(J,BE) * inout(I-1,J+1); 02771 02772 inout(I,J) = val; 02773 } 02774 02775 }; 02776 02777 02778 template <int N, int J, int AL, int BE> 02779 class TTriangleReduceFactorCol<N,J,-1,AL,BE> 02780 { 02781 public: 02782 enum { FLOP = 0 }; 02783 template <class T> 02784 static void Do (T & inout) { ; } 02785 template <class T> 02786 static void Trans (T & inout) { ; } 02787 }; 02788 02789 02790 template <int N, int J, int AL, int BE> 02791 class TTriangleReduceFactor 02792 { 02793 public: 02794 enum { FLOP = TTriangleReduceFactor<N,J-1,AL,BE>::FLOP + 02795 TTriangleReduceFactorCol<N,J,N-J,AL,BE>::FLOP }; 02796 02797 template <class T> 02798 static ALWAYS_INLINE void Do (T & inout) 02799 { 02800 if (J == N) 02801 for (int i = 0; i <= N+1; i++) 02802 inout(i, N+1-i) = 0; 02803 02804 TTriangleReduceFactorCol<N,J,N-J,AL,BE> ::Do(inout); 02805 TTriangleReduceFactor<N,J-1,AL,BE> ::Do(inout); 02806 } 02807 02808 template <class T> 02809 static ALWAYS_INLINE void Trans (T & inout) 02810 { 02811 /* 02812 if (J == N) 02813 for (int i = 0; i <= N+1; i++) 02814 inout(i, N+1-i) = 0; 02815 */ 02816 TTriangleReduceFactor<N,J-1,AL,BE> ::Trans(inout); 02817 TTriangleReduceFactorCol<N,J,N-J,AL,BE> ::Trans(inout); 02818 } 02819 }; 02820 02821 02822 template <int N, int AL, int BE> 02823 class TTriangleReduceFactor<N,-1,AL,BE> 02824 { 02825 public: 02826 enum { FLOP = 0 }; 02827 template <class T> 02828 static void Do (T & inout) { ; } 02829 template <class T> 02830 static void Trans (T & inout) { ; } 02831 }; 02832 02833 02834 02835 02836 02837 02838 02839 02840 /* 02841 02842 02843 template <int N, int J, int AL, int BE> 02844 class TTriangleReduce 02845 { 02846 // reduce, fac 02847 // P_i^alpha(x) (1-x)/2 = c_i^{alpha} P_i^{\alpha-1} + hatc_i^{alpha} P_{i+1}^{\alpha-1} 02848 static double c (int i, int alpha) { return double(i+alpha)/double(2*i+alpha+1); } 02849 static double hatc (int i, int alpha) { return -double(i+1)/double(2*i+alpha+1); } 02850 02851 // decrease alpha 02852 // P_i^alpha(x) = e_i^{alpha} P_i^{\alpha-1} + hate_i^{alpha} P_{i-1}^{\alpha} 02853 static double e (int i, int alpha) { return double(2*i+alpha)/double(i+alpha); } 02854 static double hate (int i, int alpha) { return double(i)/double(i+alpha); } 02855 02856 02857 public: 02858 enum { FLOP = TTriangleReduce<N,J-1,AL,BE>::FLOP + 3*(N-J+1) + 02859 TReduceAlpha<N-J+1, AL+2*J-1>::FLOP 02860 }; 02861 02862 02863 template <class T> 02864 static ALWAYS_INLINE void Do (T & inout) 02865 { 02866 Vec<N-J+2> hv; 02867 hv(0) = 0.0; 02868 02869 for (int i = 0; i <= N-J; i++) 02870 { 02871 double val = inout(i,J); 02872 inout(i,J) = val * e(J,BE); 02873 hv(i) += val * c(i, AL+2*J) * hate(J, BE); 02874 hv(i+1) = val * hatc(i, AL+2*J) * hate(J, BE); 02875 } 02876 02877 TReduceAlpha<N-J+1, AL+2*J-1>::Do(hv); 02878 02879 for (int i = 0; i <= N-J+1; i++) 02880 inout(i,J-1) += hv(i); 02881 02882 TTriangleReduce<N,J-1,AL,BE> ::Do(inout); 02883 02884 TReduceAlpha<N-J, AL+2*J>::Do(inout.Col(J)); 02885 } 02886 02887 02888 02889 02890 template <class T> 02891 static ALWAYS_INLINE void Trans (T & inout) 02892 { 02893 Vec<N-J+2> hv; 02894 02895 TReduceAlpha<N-J, AL+2*J>::Trans(inout.Col(J)); 02896 02897 TTriangleReduce<N,J-1,AL,BE>::Trans(inout); 02898 02899 for (int i = 0; i <= N-J+1; i++) 02900 hv(i) = inout(i,J-1); 02901 02902 TReduceAlpha<N-J+1, AL+2*J-1>::Trans(hv); 02903 02904 for (int i = 0; i <= N-J; i++) 02905 { 02906 inout(i,J) = 02907 e(J,BE) * inout(i,J) 02908 + c(i, AL+2*J) * hate(J, BE) * hv(i) 02909 + hatc(i, AL+2*J) * hate(J, BE) * hv(i+1); 02910 } 02911 } 02912 }; 02913 02914 02915 template <int N, int AL, int BE> 02916 class TTriangleReduce<N,0,AL,BE> 02917 { 02918 public: 02919 enum { FLOP = 0 }; 02920 template <class T> 02921 static void Do (T & inout) 02922 { 02923 TReduceAlpha<N, AL>::Do(inout.Col(0)); 02924 } 02925 02926 template <class T> 02927 static void Trans (T & inout) 02928 { 02929 TReduceAlpha<N, AL>::Trans(inout.Col(0)); 02930 } 02931 }; 02932 */ 02933 02934 02935 02936 02937 02938 02939 02940 02941 template <int N, int I, int J, int AL, int BE> 02942 class TTriangleReduceLoop2New 02943 { 02944 public: 02945 template <class T> 02946 static ALWAYS_INLINE void Do (T & inout) 02947 { 02948 if (BE+I == 0) cout << "is 0" << endl; 02949 double fac = 1.0 / ( (BE + I)*(4 + AL-1 + 2*I-2 + J-1) ); 02950 double val = inout(J,I); 02951 inout(J,I) = fac * (BE + 2*I)*(5 + AL-1 + 2*I-2 + 2*J-2) * val; 02952 if (I >= 1) 02953 { 02954 inout(J,I-1) += fac * I*(3 + AL-1 + 2*I-2 + J-1) * val; 02955 inout(1+J,I-1) -= fac * I*(2 + J-1) * val; 02956 } 02957 if (J >= 1) 02958 inout(J-1, I) += fac * (BE + I)*(1 + J-1) * val; 02959 02960 TTriangleReduceLoop2New<N,I,J-1,AL,BE>::Do(inout); 02961 } 02962 02963 02964 template <class T> 02965 static ALWAYS_INLINE void Trans (T & inout) 02966 { 02967 TTriangleReduceLoop2New<N,I,J-1,AL,BE>::Trans(inout); 02968 02969 double fac = 1.0 / ( (BE + I)*(4 + AL-1 + 2*I-2 + J-1) ); 02970 double val = inout(J,I) * ( fac * (BE + 2*I)*(5 + AL-1 + 2*I-2 + 2*J-2) ); 02971 if (I >= 1) 02972 { 02973 val += inout(J,I-1) * ( fac * I*(3 + AL-1 + 2*I-2 + J-1) ); 02974 val -= inout(1+J,I-1) * ( fac * I*(2 + J-1)); 02975 } 02976 if (J >= 1) 02977 val += inout(J-1, I) * (fac * (BE + I)*(1 + J-1)); 02978 02979 inout(J,I) = val; 02980 } 02981 02982 }; 02983 02984 template <int N, int AL, int I, int BE> 02985 class TTriangleReduceLoop2New<N,I,-1,AL,BE> 02986 { 02987 public: 02988 enum { FLOP = 0 }; 02989 template <class T> 02990 static void Do (T & inout) { ; } 02991 02992 template <class T> 02993 static void Trans (T & inout) { ; } 02994 }; 02995 02996 02997 02998 template <int N, int I, int AL, int BE> 02999 class TTriangleReduceNew 03000 { 03001 public: 03002 enum { FLOP = TTriangleReduceNew<N,I-1,AL,BE>::FLOP + 4*(N-I+1) }; 03003 03004 03005 template <class T> 03006 static ALWAYS_INLINE void Do (T & inout) 03007 { 03008 TTriangleReduceLoop2New<N,I,N-I,AL,BE>::Do(inout); 03009 TTriangleReduceNew<N,I-1,AL,BE>::Do(inout); 03010 } 03011 03012 template <class T> 03013 static ALWAYS_INLINE void Trans (T & inout) 03014 { 03015 TTriangleReduceNew<N,I-1,AL,BE>::Trans(inout); 03016 TTriangleReduceLoop2New<N,I,N-I,AL,BE>::Trans(inout); 03017 } 03018 }; 03019 03020 template <int N, int AL, int BE> 03021 class TTriangleReduceNew<N,-1,AL,BE> 03022 { 03023 public: 03024 enum { FLOP = 0 }; 03025 03026 template <class T> 03027 static void Do (T & inout) { ; } 03028 03029 template <class T> 03030 static void Trans (T & inout) { ; } 03031 }; 03032 03033 03034 03035 03036 03037 /* 03038 03039 template <int N, int I, int AL, int BE> 03040 class TTriangleReduceNew 03041 { 03042 public: 03043 03044 template <class T> 03045 static ALWAYS_INLINE void Do (T & inout) 03046 { 03047 for (int J = N-I; J >= 0; J--) 03048 { 03049 double fac = 1.0 / ( (2 + BE-1 + I-1)*(4 + AL-1 + 2*I-2 + J-1) ); 03050 double val = inout(J,I); 03051 inout(J,I) = fac * (3 + BE-1 + 2*I-2)*(5 + AL-1 + 2*I-2 + 2*J-2) * val; 03052 if (I >= 1) 03053 { 03054 inout(1+J,I-1) += -fac * I*(2 + J-1) * val; 03055 inout(J,I-1) += fac * I*(3 + AL-1 + 2*I-2 + J-1) * val; 03056 } 03057 if (J >= 1) 03058 inout(J-1, I) += fac * (2 + BE-1 + I-1)*(1 + J-1) * val; 03059 } 03060 03061 TTriangleReduceNew<N,I-1,AL,BE>::Do(inout); 03062 } 03063 03064 03065 template <class T> 03066 static ALWAYS_INLINE void Trans (T & inout) 03067 { 03068 TTriangleReduceNew<N,I-1,AL,BE>::Trans(inout); 03069 03070 for (int J = 0; J <= N-I; J++) 03071 { 03072 double fac = 1.0 / ( (2 + BE-1 + I-1)*(4 + AL-1 + 2*I-2 + J-1) ); 03073 double val = inout(J,I) * ( fac * (3 + BE-1 + 2*I-2)*(5 + AL-1 + 2*I-2 + 2*J-2) ); 03074 if (I >= 1) 03075 { 03076 val += inout(J,I-1) * ( fac * I*(3 + AL-1 + 2*I-2 + J-1) ); 03077 val -= inout(1+J,I-1) * ( fac * I*(2 + J-1)); 03078 } 03079 03080 if (J >= 1) 03081 val += inout(J-1, I) * (fac * (2 + BE-1 + I-1)*(1 + J-1)); 03082 03083 inout(J,I) = val; 03084 } 03085 } 03086 03087 }; 03088 03089 template <int N, int AL, int BE> 03090 class TTriangleReduceNew<N,-1,AL,BE> 03091 { 03092 public: 03093 enum { FLOP = 0 }; 03094 template <class T> 03095 static void Do (T & inout) 03096 { ; } 03097 03098 template <class T> 03099 static void Trans (T & inout) 03100 { ; } 03101 }; 03102 */ 03103 03104 03105 03106 03107 03108 03109 03110 03111 03112 03113 03114 03115 } 03116 03117 03118 #endif