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)!,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},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,x^{(2n+1)}=x(x^2{)}^n$

Rozwiązanie

Do wzoru podanego w podpowiedzi dodajemy warunki brzegowe:

$x^0=1,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.

$\mathrm{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{array}{rl}{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}\\ & ⋮& \\ 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{array}$

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.