#include <stdio.h>
#include <math.h>
typedef double R;
void pm(int n, char * m, R a[n][n]){
   int i; for(i=0; i<n; ++i) {int j; for(j=0; j<n; ++j) printf("%12.8f ", a[i][j]);}
   printf(" %s\n", m);}

void expm(int n, R x[n][n], R a[n][n]){
  R c = 1;
  R z[n][n]; {int j=n; while(j--) {int k=n; while(k--) {
     z[j][k] = j==k; a[j][k] = 0;}}}
  while(1){
    {int j=n; while(j--) {int k=n; while(k--)
      if(fabs(z[j][k]) > 1.e-10) goto more;}} return; more:
    {int j=n; while(j--) {int k=n; while(k--) a[j][k] += z[j][k];}}
    {R p[n][n];
      {int j=n; while(j--) {int k=n; while(k--)
        {R s=0; int i=n; while(i--) s += z[j][i]*x[i][k]; p[j][k] = s/c;}}}
    {int j=n; while(j--) {int k=n; while(k--) z[j][k] = p[j][k];}}}
  c += 1.;}}

int main(){if(0) {R t[3][3] = {{2, 0, 0}, {0, 3, 0}, {0, 0, 4}};
  R a[3][3]; expm(3, t, a); pm(3, "ans", a);
  printf("%e %e %e\n", exp(2), exp(3), exp(4));}
  if(1) {R P=2, Q=1; const int n = 2;
    R ti[] = {8, 14}; R C[2] = {10, 3};
    R c[][2] = {{-P, P}, {P, -P-Q}};
    R M[2][2];
    {int j=n; while(j--) {int k=n; while(k--) M[j][k] = c[j][k]/C[j];}}
    {R t; for(t=0; t<2.05; t += 0.1) {R y[n][n], a[n][n];
      {int j=n; while(j--) {int k=n; while(k--) y[j][k] = t*M[j][k];}}
      expm(n, y, a);
      printf("t=%5.1f ", t); {int j; for(j=0; j<n; ++j) {R s = 0;
          int k=n; while(k--) s += a[j][k]*ti[k]; printf("%12.8f ", s);}
        printf("\n");}}}}
  return 0;}