// produced by Taka Kondo on 16 Feb. 2009, corrected 2011.11.25
// refer  Amaldi et al., Phys. Lett. 260(1991)447
#include <iostream>
#include <math.h>
void running(int isusy,float* a0,float q0, int nf, int nh, float q, float* a);

int main()
{
float a0[4] = { 0., 0.01684,  0.03353,  0.120 };
float a[4], ainv;
float mz = 91.1876, q0, q;
int isusy = 0;
float SUSY=1.E19,;   // case for Standard Model
int nf = 3, nh = 1 ;
for (int j=0 ; j < 96 ; j++){
   q = exp(j/5.*log(10.));
   if(isusy==0 && q > SUSY) {
      //std::cout << "----------SUSY is introduced here---------"<< std::endl;
      running(isusy, a0, q0, nf , nh, SUSY, a);
      for (int k=1 ; k <=3 ; k++) { 
         a0[k] = a[k];
         ainv = 1./a[k];
         std::cout<<k <<" "<< SUSY <<" "<< a[k] <<" "<< ainv << std::endl;
      }
     isusy = 1;
   }
   if( q <= SUSY) { nh = 1; q0 = mz;}
   if( q >  SUSY) { nh = 2; q0 = SUSY;}
   running(isusy, a0, q0, nf , nh, q, a);
   for (int k=1 ; k<4 ; k++) {
   ainv = 1./a[k];
   std::cout << k << " " << q << " " << a[k] << " " << ainv << std::endl;
   }
}
return 0;
}

//--------------function running------------
void running(int isusy, float* a0, float q0, int nf, int nh, float q, float* a)
{
float b[4],s[4],b2[4][4],s2[4][4];
float pi = 3.1415;
float x0, x,anew, ainv, da, rbeta, beta0, beta1;

b[1]=0.        +nf*4./3.    +nh*1./10.; 
b[2]=-22./3.   +nf*4./3.    +nh*1./6.; 
b[3]=-11.      +nf*4./3.    +0.;

s[1]=0.        +nf*2.    +nh*3./10.; 
s[2]=-6.       +nf*2.    +nh*1./2.; 
s[3]=-9.       +nf*2.    +0.;

b2[1][1]=0.     +nf*19./15.  +nh*9./50.; 
b2[1][2]=0.     +nf*3./5.    +nh*9./10.; 
b2[1][3]=0.     +nf*44./15.  +0.;
b2[2][1]=0.     +nf*1./5.    +nh*3./10.;
b2[2][2]=-136./3.+nf*49./3.   +nh*13./6.; 
b2[2][3]=0.     +nf*4.       +0.;
b2[3][1]=0.     +nf*11./30.  +0.;
b2[3][2]=0.     +nf*3./2.    +0.;
b2[3][3]=-102.  +nf*76./3.   +0.;

s2[1][1]=0.     +nf*38./15.  +nh*9./50.; 
s2[1][2]=0.     +nf*6./5.    +nh*9./10.; 
s2[1][3]=0.     +nf*88./15.  +0.;
s2[2][1]=0.     +nf*2./5.    +nh*3./10.;
s2[2][2]=-24.   +nf*14.      +nh*7./2.; 
s2[2][3]=0.     +nf*8.       +0.;
s2[3][1]=0.     +nf*11./15.  +0.;
s2[3][2]=0.     +nf*3.       +0.;
s2[3][3]=-54.   +nf*68./3.   +0.;
int m,n;
//------------------------------------------
   x = q*q;
   x0 = q0 * q0;
   for (int k=1 ; k<=3  ; k++) {
	a[k] = a0[k];
	if(k==1){m=2;n=3;}
	if(k==2){m=3,n=1;}
	if(k==3){m=1;n=2;}
        if(isusy == 0) {
  		beta0 = -1./(4*pi)*
		( b[k] + (b2[k][m]*a0[m]+b2[k][n]*a0[n])/(4*pi) );
		beta1 = -2.*b2[k][k]/(16*pi*pi);
	}
	else {
		beta0 = -1./(4*pi)*
		( s[k] + (s2[k][m]*a0[m]+s2[k][n]*a0[n])/(4*pi) );
		beta1 = -2.*s2[k][k]/(16*pi*pi);
	}
	rbeta = beta1/beta0;
	for (int i=0; i<10 ; i++) {
		anew = 1./(beta0*log(x/x0) + 1/a0[k] + rbeta*log( 
			(1./a[k] + rbeta)/(1/a0[k] + rbeta) ) );
		da = fabs(anew - a[k]);
		if(da <= 0.000001) continue;
//	std::cout << "k, a[k],anew,da="<<k<<" "<<a[k]<<" "
//	<<anew << " " << da <<std::endl;
		a[k] = anew;
	}
   }
   return ;
}