computes exp(t a), a a 3x3 double matrix

expm3

int expm3(double t, double *a, double *aa, double X[3], double dis)
/*
	computes exp(t a)

	input
	t	double
	a 	3x3 double matrix
	aa	a*a	
	dis	negative, if matrix has complex conjugate pair, zero if two eigenvalues are equal and positive, if all eigenvalues are real and different 
	X	double[3], containing either:
			the three real eigenvalues in ascending order
			or
			first the real eigenvalue, then the real component and the absolute value of the complex component of the complex conjugate eigenvalues

	output
	aa	the exponential

	return	unspecified integer

*/
{
	
	double la1; // and la2 la3: eigenvalues, la2 and la3 might be complex
	double r1, r3, R; // and r2:  putzer coefficients, r2 might be complex
	int ca = 0;
	
	
	/* aa = a.a */
	int i, j, k;
	for (i = 0; i < 9; i++) aa[i] = 0.0;
	
	for (i = 0; i < 3; i++)
	  for (k = 0; k < 3; k++)
	    for (j = 0; j < 3; j++)
		 aa[i*3 + j] += a[i*3 + k]*a[k*3 + j];


        if (dis == 0.) // real roots, at least x[1]=x[2]
        {   
	    double la2, r2; 
	    la1 = X[0];
	    la2 = X[1];
            if (cabs(la1 - la2) < DBL_EPSILON ) // la1 = la2 = la3
            {
    		ca = 1;
	        r1 = exp(la1 * t);
		r2 = t * exp(la1 * t);
		r3 = 0.5 * t * t * exp(la1 * t);
  		R = r2 - la2*r3; 

	    } else
	    {
       		ca = 2;
	    	    r1 = exp(la1 * t);
        	    r2 = (expm1(la1 * t) - expm1(la2 * t)) / (la1 - la2); // la1 != la2
        	    r3 = (exp(la1 * t) - exp(la2 * t)*(1.+t*(la1-la2))) / square(la1-la2);
       		    R = r2 - la2*r3; 
       	    }
        
        } else if (dis > 0.0) // three distinct real roots
        {
	    double la2, la3; 
	    double  r2;
		ca = 3;
	    la1 = X[0];
	    la2 = X[1];
            la3 = X[2];
            r1 = exp(la1 * t);
            r2 = (expm1(la1 * t) - expm1(la2 * t))/ (la1 - la2); // expm1(a*t) - expm1(b*t) is more stable then  exp(a*t) - exp(b*t)
            r3 = -((la2 - la3)*exp(la1*t) + (la3 - la1)*exp(la2*t) + (la1 - la2)*exp(la3*t))
            		/ ((la1 - la2)*(la2 - la3)*(la3 - la1));
            R = r2 - la2*r3;
 
	} else //   One real and a complex conjugate pair of roots, all different
	{
	    double complex la2, la3; // only these values are possibly complex
	    double complex r2;
		ca = 4;
	    
	    la1 =  X[0];
            la2 =  X[1] + I*X[2];
            la3 =  X[1] - I*X[2];
        
    	    r1 = exp(la1 * t);
            r2 = (cexp(la1 * t) - cexp(la2 * t))/ ((complex) la1 - la2); // possibly complex, cexpm1 still missing, so use cexp(a*t) - cexp(b*t)
            r3 = -((la2 - la3)*exp(la1*t) + (la3 - la1)*cexp(la2*t) + (la1 - la2)*cexp(la3*t)) // r3 is real, see short comm. by r.r. huilgol
            		/ ((la1 - la2)*(la2 - la3)*(la3 - la1));
            R = r2 - la2*r3; // also real
     }


	for (i = 0;i < 3; i++) for (j = 0;j < 3;j++)
	{
		aa[i*3+j] = r3 * aa[i*3+j] + (R - la1*r3)* a[i*3+j] ;
		if (i == j) aa[i*3+j] += r1 - la1 * R; 
					//= r1 I + (A - la1 I)*(r2 I  + r3*(A - la2 I))
	}
	return ca;

}