/* 
 * "rustic" math library, for muLinux
 * (C) by M. Andreoli & E. Cavalli

 * This library emulate the math library. Linking the math.o, 
   avoid the use of -lm in compiling. This orrendous piece of
   code use Taylor series expansions and other dirty hack.
	USE AT OWN RISK!
 */

#define PI 3.14159265359
#define LN_10 2.30258509  /* ln(10) */
#define SQRT_E 1.6487212707

#define double double 

#define ABS(x)  ( (x)>=0?(x):(-x) )


/* exp() */

double exp2( double x)  
{
double x2,x212;
if(x<0.5) {
x2=x*x;
x212=12.0+x2;
return (x212+6.0*x)/(x212-6.0*x);
}
else {
x-=0.5;
x2=x*x;
x212=12.0+x2;
return SQRT_E*(x212+6.0*x)/(x212-6.0*x);
}

}

double exp(double x)
{
double r,rest;
int i,m,sign;

sign= (int)(x>=0?1.0:-1.0);
x=sign*x;

m=(int)x;
rest=x-1.0*(int)x;
r=1.0;
for (i=0; i<m;i++) r=r*2.71828182; 
return sign>0? r*exp2(rest) : exp2(-rest)/r;
}


/* cos(x) */

double cos(double x)
{
double y,x2,x4;
int n;
x=ABS(x);

n=(int) (x/(2.0*PI));
y=x - n*PI*2.0-PI;

x2=y*y;
x4=x2*x2;
return -(1. - 0.4563492063492063*x2 + 
     0.0207010582010582*x4)/
   (1. + 0.04365079365079364*x2 + 
     0.00085978835978836*x4);

}

/* sin(x) */                                                                    
                                                                                
double sin(double x)                                                         
{                                                                               
return cos(PI/2.0-x); 
}                                                                               

/* ceil(x) */

double ceil(double x)
{
int i;
i=(int) x;
return i;
}

/* log(x) if x is near x */

double log_serie(double y)
{
double x,p,sum;
int i,sign, iter=80;

if (  (y-1.0 < 0.00001) && ( y-1.0 > -0.00001 ) ) return 0.0; 

if ( y > 1.0 ) {
	x=y-1;
	sum=x; sign=-1; p=x*x;
	for(i=2; i<=iter; i++) {
	sum=sum +sign*p/i;
	p*=x; sign*=-1;
	}

} else {
	x=1-y;
        sum=x; p=x*x;
        for(i=2; i<=iter; i++) {
        sum=sum +p/i;
        p*=x; 
        }
	sum*=-1;
	}

return sum;
}

/* log(x) : a dirty hack! */

double log( double x)
{
char buf[64];
float man, ex;
sprintf(buf,"%2.2E",x);
buf[4]='\0'; buf[8]='\0';
man=atof(buf); ex=atof(buf+5);

/* reduce mantissa < 1 */
man/=10.0; ex++;
return log_serie(man) + ex*LN_10;
}


/* sqrt(x) */

double sqrt(double x)
{
if ( x < 0.0 ) { printf("sqrt: x < 0\n"); exit(1);}
return (double) exp(0.5*log(x));
}

/*  pow(x,y) */

double pow (double b, double x)
{
if ( b < 0.0 ) { printf("pow: base < 0\n"); exit(1);}
return (double) exp(x*log(b));
}