Zaawansowane CPP/Ćwiczenia 8: Metaprogramowanie

Ćwiczenie 1

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

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

Rozwiązanie

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

$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 $N^{M}$. Np.:

Pow<3,4>::val;

powinno mieć wartość 81.

Rozwiązanie

Implementujemy rekurencyjną definicję potęgi:

$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ź

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

Rozwiązanie

Do wzoru podanego w podpowiedzi dodajemy warunki brzegowe:

$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)$:

$sin<$N$>(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ą:

$sin<N>(x)=sin<N-1>(x)+(-1)^{N+1}{\frac {x^{2N-1}}{(2N-1)!}},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.

$\operatorname {inner}$ <N> $(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)

${\begin{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}\\\end{aligned}}$

Tablica $A$ jest reprezentowana w pamięci zgodnie z konwencją $C$, tzn. wiersz po wierszu: elementowi $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.