NGSolve  4.9
fem/recursive_pol.hpp
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