Zaawansowane CPP/Ćwiczenia 8: Metaprogramowanie

From Studia Informatyczne

Ćwiczenie 1

Napisz szablon funkcji lub klasy wyliczający funkcję silnia:

\displaystyle  n!=n(n-1)(n-2)\cdots 1

Rozwiązanie

Rozwiązanie jest bezpośrednim zastosowaniem rekurencyjnej definicji funkcji silnia:

\displaystyle  n!=n*(n-1)!,\quad 0!=1
template<size_t N> struct factorial {
   enum {val=N*factorial<N-1>::val};
};
template<> struct factorial <0>{ enum {val=1}; };

Patrz plik factorial.h.

Ćwiczenie 2

Zaimplementuj szablon Pow<N,M> obliczający \displaystyle N^M. Np.:

Pow<3,4>::val;

powinno mieć wartość 81.

Rozwiązanie

Implementujemy rekurencyjną definicję potęgi:

\displaystyle  N^M=N*N^{M-1},\quad N^0=1
template<size_t N,size_t M> struct pow {
   enum {val=N*pow<N,M-1>::val};
};
template<size_t N> struct pow<N,0> { enum {val=1}; };

Dla przyspieszenia dodajemy dwie specjalizacje:

template<size_t M> struct pow<1,M>{
   enum {val=1};
};
template<size_t M> struct pow<0,M>{ enum {val=0}; };

Powyższa specjalizacja spowoduje, że konkretyzacja Pow<0,0> będzie niejednoznaczna. Możemy to tak zostawić, bo 0^0 jest nieokreślone. Jeśli jednak zdecydujemy się na jakąś wartość, np. zero, to należy dodać jeszcze jedną specjalizację:

template<> struct pow<0,0>{
   enum {val=0};
};

Patrz plik pown.h.

Ćwiczenie 3

Wymyśl i zaimplementuj jako metaprogram szybszy algorytm funkcji pow(x).

Podpowiedź

\displaystyle  x^{(2 n)} = (x^2)^n,\quad x^{(2 n+1)} = x (x^2)^n

Rozwiązanie

Do wzoru podanego w podpowiedzi dodajemy warunki brzegowe:

\displaystyle  x^0=1,\quad x^1=x

i implementujemy bezpośrednio:

template<size_t N> double pow(double x) {
   return ((N%2)?x:1)*pow<N/2>(x*x);
}
template<> double pow<1>(double x) {return x;};
template<> double pow<0>(double x) {return 1.0;};

Patrz plik powx.h.

Ćwiczenie 4

Napisz szablon generujący pierwsze N wyrazów rozwinięcia funkcji sin(x):

\displaystyle  sin<N> \displaystyle  (x)=x-\frac{x^3}{3!}+\cdots +(-1)^{N+1}\frac{(x)^{2N-1}}{(2N-1)!}

Możesz skorzystać z rozwiązań wcześniejszych zadań.

Rozwiązanie

Jak zwykle najpierw piszemy definicję rekurencyjną:

\displaystyle  sin<N> (x)=sin<N-1> (x)+(-1)^{N+1}\frac{x^{2N-1}}{(2N-1)!},  \displaystyle sin<0> (x)=0

i implementujemy ją korzystając z wcześniej zdefiniowanych szablonów:

template<int N> double sin(double x) {
   return sin<N-1>(x) + (N%2?1:-1)*pow<(2*N-1)>(x)/factorial<2*N-1>::val;
}
template<> double sin<0>(double x) { return 0; }
template<int N> double inner(double *a, double *b) { return (*a)*(*b)+inner<N-1>(++a,++b); }
template<> double inner<1>(double *a, double *b) { return (*a)*(*b); }

Patrz plik inner.h.

Ćwiczenie 5

Napisz szablon generujący funkcję implementującą iloczyn skalarny dwu wektorów:

template<size_t N> double inner(double *x, double *y);

Parametrem szablonu ma być dlugość mnożonych wektorów.

\displaystyle  \operatorname{inner} <N> \displaystyle  (x,y)=x_1 y_1+\cdots y_N x_N

Ćwiczenie 6

Rozszerz powyższy szablon tak, aby również typ elementów wektora był parametrem szablonu:

template<size_t N, typename T> T dot(T *x, T *y);

Podpowiedź

Konieczne będzie użycie specjalizacji częściowej. W tym celu użyj pomocniczej klasy (specjalizacja częściowa jest niedozwolona dla funkcji).

Rozwiązanie

Ponieważ warunek brzegowy N=1 będzie wymagał specjalizacji częściowej, implementujemy iloczyn skalarny jako składową klasy:

template<size_t N,typename T = double > struct Inner {
   static T dot(T *a,T *b) {
       return (*a)*(*b)+Inner<N-1,T>::dot(++a,++b);
   }
}
template<typename T > struct Inner<1,T> { static T dot(T *a,T *b) { return (*a)*(*b); } };

a następnie dodajemy funkcję wywołującą tę metodę:

template<size_t N,typename T> T dot( T *a,T *b) {
   return Inner<N,T>::dot(a,b);
};

Patrz plik inner.h.

Ćwiczenie 7

Napisz szablon generujący funkcję implementującą iloczyn macierzy NxM i wektora o M elementach:

 void matrix_v<N>(double *A,double *v,double *u)
\displaystyle \aligned u_0&= A_{0,0} v_0+A_{0,1} v_1+\cdots+A_{0,M-1} v_{M-1}\\ u_1&= A_{1,0} v_0+A_{1,1} v_1+\cdots+A_{1,M-1} v_{M-1}\\ &\vdots&\\ u_{N-1}&= A_{N-1,0} v_0+A_{N-1,1} v_1+\cdots+A_{N-1,M-1} v_{M-1}\\ \endaligned

Tablica \displaystyle A jest reprezentowana w pamięci zgodnie z konwencją C, tzn. wiersz po wierszu: elementowi \displaystyle A_{i,j} odpowiada A[M*i+j].

Rozwiązanie

Najpierw definiujemy szablon, który wykonuje mnożenie I-tego wiersza macierzy A przez wektor v:

template<int M,size_t I> double row_vec(double *A, double *v) {
   return inner<M>(A+I*M,v);
}

Następnie implementujemy rekurencyjnie mnożenie kolejnych wierszy macierzy:

template<int N,int M> struct matrix_vec_c {
   static void matrix_vec(double *A,double *v,double *u) {
       u[N-1]=row_vec<M,N-1>(A,v);
       matrix_vec_c<N-1,M>::matrix_vec(A,v,u);
   }
};
template<int M> struct matrix_vec_c<0,M> { static void matrix_vec(double *A,double *v,double *u) {} };

Używamy szablonu klasy aby móc wykorzystać specjalizację częściową. Tak jak w poprzednim zadaniou opakowujemy wywołanie metody w szablon funkcji:

template<size_t N,size_t M>
inline void matrix_vec(double *u,double *A,double *v) {
    matrix_vec_c<N,M>::matrix_vec(A,v,u);
}

Całość kodu znajduje się w pliku matrix.h.