diff --git a/src/bsp/lk/lib/libm/e_acos.c b/src/bsp/lk/lib/libm/e_acos.c
new file mode 100644
index 0000000..8ba672a
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_acos.c
@@ -0,0 +1,104 @@
+
+/* @(#)e_acos.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_acos.c,v 1.10 2005/02/04 18:26:05 das Exp $";
+#endif
+
+/* __ieee754_acos(x)
+ * Method :                  
+ *	acos(x)  = pi/2 - asin(x)
+ *	acos(-x) = pi/2 + asin(x)
+ * For |x|<=0.5
+ *	acos(x) = pi/2 - (x + x*x^2*R(x^2))	(see asin.c)
+ * For x>0.5
+ * 	acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
+ *		= 2asin(sqrt((1-x)/2))  
+ *		= 2s + 2s*z*R(z) 	...z=(1-x)/2, s=sqrt(z)
+ *		= 2f + (2c + 2s*z*R(z))
+ *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
+ *     for f so that f+c ~ sqrt(z).
+ * For x<-0.5
+ *	acos(x) = pi - 2asin(sqrt((1-|x|)/2))
+ *		= pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
+ *
+ * Special cases:
+ *	if x is NaN, return x itself;
+ *	if |x|>1, return NaN with invalid signal.
+ *
+ * Function needed: sqrt
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+one=  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+pi =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+pio2_hi =  1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
+pio2_lo =  6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
+pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+double
+__ieee754_acos(double x)
+{
+	double z,p,q,r,w,s,c,df;
+	int32_t hx,ix;
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x3ff00000) {	/* |x| >= 1 */
+	    u_int32_t lx;
+	    GET_LOW_WORD(lx,x);
+	    if(((ix-0x3ff00000)|lx)==0) {	/* |x|==1 */
+		if(hx>0) return 0.0;		/* acos(1) = 0  */
+		else return pi+2.0*pio2_lo;	/* acos(-1)= pi */
+	    }
+	    return (x-x)/(x-x);		/* acos(|x|>1) is NaN */
+	}
+	if(ix<0x3fe00000) {	/* |x| < 0.5 */
+	    if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
+	    z = x*x;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    return pio2_hi - (x - (pio2_lo-x*r));
+	} else  if (hx<0) {		/* x < -0.5 */
+	    z = (one+x)*0.5;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    s = sqrt(z);
+	    r = p/q;
+	    w = r*s-pio2_lo;
+	    return pi - 2.0*(s+w);
+	} else {			/* x > 0.5 */
+	    z = (one-x)*0.5;
+	    s = sqrt(z);
+	    df = s;
+	    SET_LOW_WORD(df,0);
+	    c  = (z-df*df)/(s+df);
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    w = r*s+c;
+	    return 2.0*(df+w);
+	}
+}
diff --git a/src/bsp/lk/lib/libm/e_acosf.c b/src/bsp/lk/lib/libm/e_acosf.c
new file mode 100644
index 0000000..a11f48e
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_acosf.c
@@ -0,0 +1,81 @@
+/* e_acosf.c -- float version of e_acos.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_acosf.c,v 1.7 2002/05/28 17:03:12 alfred Exp $";
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+static const float
+one =  1.0000000000e+00, /* 0x3F800000 */
+pi =  3.1415925026e+00, /* 0x40490fda */
+pio2_hi =  1.5707962513e+00, /* 0x3fc90fda */
+pio2_lo =  7.5497894159e-08, /* 0x33a22168 */
+pS0 =  1.6666667163e-01, /* 0x3e2aaaab */
+pS1 = -3.2556581497e-01, /* 0xbea6b090 */
+pS2 =  2.0121252537e-01, /* 0x3e4e0aa8 */
+pS3 = -4.0055535734e-02, /* 0xbd241146 */
+pS4 =  7.9153501429e-04, /* 0x3a4f7f04 */
+pS5 =  3.4793309169e-05, /* 0x3811ef08 */
+qS1 = -2.4033949375e+00, /* 0xc019d139 */
+qS2 =  2.0209457874e+00, /* 0x4001572d */
+qS3 = -6.8828397989e-01, /* 0xbf303361 */
+qS4 =  7.7038154006e-02; /* 0x3d9dc62e */
+
+float
+__ieee754_acosf(float x)
+{
+	float z,p,q,r,w,s,c,df;
+	int32_t hx,ix;
+	GET_FLOAT_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix==0x3f800000) {		/* |x|==1 */
+	    if(hx>0) return 0.0;	/* acos(1) = 0  */
+	    else return pi+(float)2.0*pio2_lo;	/* acos(-1)= pi */
+	} else if(ix>0x3f800000) {	/* |x| >= 1 */
+	    return (x-x)/(x-x);		/* acos(|x|>1) is NaN */
+	}
+	if(ix<0x3f000000) {	/* |x| < 0.5 */
+	    if(ix<=0x23000000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
+	    z = x*x;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    return pio2_hi - (x - (pio2_lo-x*r));
+	} else  if (hx<0) {		/* x < -0.5 */
+	    z = (one+x)*(float)0.5;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    s = __ieee754_sqrtf(z);
+	    r = p/q;
+	    w = r*s-pio2_lo;
+	    return pi - (float)2.0*(s+w);
+	} else {			/* x > 0.5 */
+	    int32_t idf;
+	    z = (one-x)*(float)0.5;
+	    s = __ieee754_sqrtf(z);
+	    df = s;
+	    GET_FLOAT_WORD(idf,df);
+	    SET_FLOAT_WORD(df,idf&0xfffff000);
+	    c  = (z-df*df)/(s+df);
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    w = r*s+c;
+	    return (float)2.0*(df+w);
+	}
+}
diff --git a/src/bsp/lk/lib/libm/e_asin.c b/src/bsp/lk/lib/libm/e_asin.c
new file mode 100644
index 0000000..1ba7026
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_asin.c
@@ -0,0 +1,113 @@
+
+/* @(#)e_asin.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_asin.c,v 1.11 2005/02/04 18:26:05 das Exp $";
+#endif
+
+/* __ieee754_asin(x)
+ * Method :                  
+ *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ *	we approximate asin(x) on [0,0.5] by
+ *		asin(x) = x + x*x^2*R(x^2)
+ *	where
+ *		R(x^2) is a rational approximation of (asin(x)-x)/x^3 
+ *	and its remez error is bounded by
+ *		|(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
+ *
+ *	For x in [0.5,1]
+ *		asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ *	Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
+ *	then for x>0.98
+ *		asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *			= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
+ *	For x<=0.98, let pio4_hi = pio2_hi/2, then
+ *		f = hi part of s;
+ *		c = sqrt(z) - f = (z-f*f)/(s+f) 	...f+c=sqrt(z)
+ *	and
+ *		asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *			= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
+ *			= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
+ *
+ * Special cases:
+ *	if x is NaN, return x itself;
+ *	if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+huge =  1.000e+300,
+pio2_hi =  1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
+pio2_lo =  6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
+pio4_hi =  7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
+	/* coefficient for R(x^2) */
+pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+double
+__ieee754_asin(double x)
+{
+	double t=0.0,w,p,q,c,r,s;
+	int32_t hx,ix;
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>= 0x3ff00000) {		/* |x|>= 1 */
+	    u_int32_t lx;
+	    GET_LOW_WORD(lx,x);
+	    if(((ix-0x3ff00000)|lx)==0)
+		    /* asin(1)=+-pi/2 with inexact */
+		return x*pio2_hi+x*pio2_lo;	
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
+	} else if (ix<0x3fe00000) {	/* |x|<0.5 */
+	    if(ix<0x3e400000) {		/* if |x| < 2**-27 */
+		if(huge+x>one) return x;/* return x with inexact if x!=0*/
+	    } else 
+		t = x*x;
+		p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+		q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+		w = p/q;
+		return x+x*w;
+	}
+	/* 1> |x|>= 0.5 */
+	w = one-fabs(x);
+	t = w*0.5;
+	p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+	q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+	s = sqrt(t);
+	if(ix>=0x3FEF3333) { 	/* if |x| > 0.975 */
+	    w = p/q;
+	    t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
+	} else {
+	    w  = s;
+	    SET_LOW_WORD(w,0);
+	    c  = (t-w*w)/(s+w);
+	    r  = p/q;
+	    p  = 2.0*s*r-(pio2_lo-2.0*c);
+	    q  = pio4_hi-2.0*w;
+	    t  = pio4_hi-(p-q);
+	}    
+	if(hx>0) return t; else return -t;    
+}
diff --git a/src/bsp/lk/lib/libm/e_asinf.c b/src/bsp/lk/lib/libm/e_asinf.c
new file mode 100644
index 0000000..1405faf
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_asinf.c
@@ -0,0 +1,84 @@
+/* e_asinf.c -- float version of e_asin.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_asinf.c,v 1.9 2005/12/04 13:52:46 bde Exp $";
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+static const float
+one =  1.0000000000e+00, /* 0x3F800000 */
+huge =  1.000e+30,
+pio2_hi =  1.5707962513e+00, /* 0x3fc90fda */
+pio2_lo =  7.5497894159e-08, /* 0x33a22168 */
+pio4_hi =  7.8539812565e-01, /* 0x3f490fda */
+	/* coefficient for R(x^2) */
+pS0 =  1.6666667163e-01, /* 0x3e2aaaab */
+pS1 = -3.2556581497e-01, /* 0xbea6b090 */
+pS2 =  2.0121252537e-01, /* 0x3e4e0aa8 */
+pS3 = -4.0055535734e-02, /* 0xbd241146 */
+pS4 =  7.9153501429e-04, /* 0x3a4f7f04 */
+pS5 =  3.4793309169e-05, /* 0x3811ef08 */
+qS1 = -2.4033949375e+00, /* 0xc019d139 */
+qS2 =  2.0209457874e+00, /* 0x4001572d */
+qS3 = -6.8828397989e-01, /* 0xbf303361 */
+qS4 =  7.7038154006e-02; /* 0x3d9dc62e */
+
+float
+__ieee754_asinf(float x)
+{
+	float t=0.0,w,p,q,c,r,s;
+	int32_t hx,ix;
+	GET_FLOAT_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix==0x3f800000) {
+		/* asin(1)=+-pi/2 with inexact */
+	    return x*pio2_hi+x*pio2_lo;
+	} else if(ix> 0x3f800000) {	/* |x|>= 1 */
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */
+	} else if (ix<0x3f000000) {	/* |x|<0.5 */
+	    if(ix<0x32000000) {		/* if |x| < 2**-27 */
+		if(huge+x>one) return x;/* return x with inexact if x!=0*/
+	    } else
+		t = x*x;
+		p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+		q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+		w = p/q;
+		return x+x*w;
+	}
+	/* 1> |x|>= 0.5 */
+	w = one-fabsf(x);
+	t = w*(float)0.5;
+	p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+	q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+	s = __ieee754_sqrtf(t);
+	if(ix>=0x3F79999A) { 	/* if |x| > 0.975 */
+	    w = p/q;
+	    t = pio2_hi-((float)2.0*(s+s*w)-pio2_lo);
+	} else {
+	    int32_t iw;
+	    w  = s;
+	    GET_FLOAT_WORD(iw,w);
+	    SET_FLOAT_WORD(w,iw&0xfffff000);
+	    c  = (t-w*w)/(s+w);
+	    r  = p/q;
+	    p  = (float)2.0*s*r-(pio2_lo-(float)2.0*c);
+	    q  = pio4_hi-(float)2.0*w;
+	    t  = pio4_hi-(p-q);
+	}
+	if(hx>0) return t; else return -t;
+}
diff --git a/src/bsp/lk/lib/libm/e_atan2.c b/src/bsp/lk/lib/libm/e_atan2.c
new file mode 100644
index 0000000..e91d309
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_atan2.c
@@ -0,0 +1,131 @@
+
+/* @(#)e_atan2.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __ieee754_atan2(y,x)
+ * Method :
+ *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
+ *	2. Reduce x to positive by (if x and y are unexceptional): 
+ *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
+ *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
+ *
+ * Special cases:
+ *
+ *	ATAN2((anything), NaN ) is NaN;
+ *	ATAN2(NAN , (anything) ) is NaN;
+ *	ATAN2(+-0, +(anything but NaN)) is +-0  ;
+ *	ATAN2(+-0, -(anything but NaN)) is +-pi ;
+ *	ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+ *	ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+ *	ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+ *	ATAN2(+-INF,+INF ) is +-pi/4 ;
+ *	ATAN2(+-INF,-INF ) is +-3pi/4;
+ *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
+ * to produce the hexadecimal values shown.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static volatile double
+tiny  = 1.0e-300;
+static const double
+zero  = 0.0,
+pi_o_4  = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
+pi_o_2  = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
+pi      = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
+static volatile double
+pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
+
+double
+__ieee754_atan2(double y, double x)
+{
+	double z;
+	int32_t k,m,hx,hy,ix,iy;
+	u_int32_t lx,ly;
+
+	EXTRACT_WORDS(hx,lx,x);
+	ix = hx&0x7fffffff;
+	EXTRACT_WORDS(hy,ly,y);
+	iy = hy&0x7fffffff;
+	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
+	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
+	   return x+y;
+	if((hx-0x3ff00000|lx)==0) return atan(y);   /* x=1.0 */
+	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
+
+    /* when y = 0 */
+	if((iy|ly)==0) {
+	    switch(m) {
+		case 0: 
+		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
+		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
+		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
+	    }
+	}
+    /* when x = 0 */
+	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
+	    
+    /* when x is INF */
+	if(ix==0x7ff00000) {
+	    if(iy==0x7ff00000) {
+		switch(m) {
+		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
+		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
+		    case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
+		    case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
+		}
+	    } else {
+		switch(m) {
+		    case 0: return  zero  ;	/* atan(+...,+INF) */
+		    case 1: return -zero  ;	/* atan(-...,+INF) */
+		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
+		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
+		}
+	    }
+	}
+    /* when y is INF */
+	if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
+
+    /* compute y/x */
+	k = (iy-ix)>>20;
+	if(k > 60) {		 	/* |y/x| >  2**60 */
+	    z=pi_o_2+0.5*pi_lo;
+	    m&=1;
+	}
+	else if(hx<0&&k<-60) z=0.0; 	/* 0 > |y|/x > -2**-60 */
+	else z=atan(fabs(y/x));		/* safe to do y/x */
+	switch (m) {
+	    case 0: return       z  ;	/* atan(+,+) */
+	    case 1: return      -z  ;	/* atan(-,+) */
+	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
+	    default: /* case 3 */
+	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+	}
+}
+
+#if SUPPORT_LONG_DOUBLE
+#if LDBL_MANT_DIG == 53
+__weak_reference(atan2, atan2l);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/e_exp.c b/src/bsp/lk/lib/libm/e_exp.c
new file mode 100644
index 0000000..c03c9de
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_exp.c
@@ -0,0 +1,166 @@
+
+/* @(#)e_exp.c 1.6 04/04/22 */
+/*
+ * ====================================================
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __ieee754_exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ *   1. Argument reduction:
+ *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ *	Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.  
+ *
+ *      Here r will be represented as r = hi-lo for better 
+ *	accuracy.
+ *
+ *   2. Approximation of exp(r) by a special rational function on
+ *	the interval [0,0.34658]:
+ *	Write
+ *	    R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ *      We use a special Remes algorithm on [0,0.34658] to generate 
+ * 	a polynomial of degree 5 to approximate R. The maximum error 
+ *	of this polynomial approximation is bounded by 2**-59. In
+ *	other words,
+ *	    R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ *  	(where z=r*r, and the values of P1 to P5 are listed below)
+ *	and
+ *	    |                  5          |     -59
+ *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2 
+ *	    |                             |
+ *	The computation of exp(r) thus becomes
+ *                             2*r
+ *		exp(r) = 1 + -------
+ *		              R - r
+ *                                 r*R1(r)	
+ *		       = 1 + r + ----------- (for better accuracy)
+ *		                  2 - R1(r)
+ *	where
+ *			         2       4             10
+ *		R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
+ *	
+ *   3. Scale back to obtain exp(x):
+ *	From step 1, we have
+ *	   exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ *	exp(INF) is INF, exp(NaN) is NaN;
+ *	exp(-INF) is 0, and
+ *	for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *	For IEEE double 
+ *	    if x >  7.09782712893383973096e+02 then exp(x) overflow
+ *	    if x < -7.45133219101941108420e+02 then exp(x) underflow
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+one	= 1.0,
+halF[2]	= {0.5,-0.5,},
+o_threshold=  7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
+u_threshold= -7.45133219101941108420e+02,  /* 0xc0874910, 0xD52D3051 */
+ln2HI[2]   ={ 6.93147180369123816490e-01,  /* 0x3fe62e42, 0xfee00000 */
+	     -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
+ln2LO[2]   ={ 1.90821492927058770002e-10,  /* 0x3dea39ef, 0x35793c76 */
+	     -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
+invln2 =  1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
+P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
+
+static volatile double
+huge	= 1.0e+300,
+twom1000= 9.33263618503218878990e-302;     /* 2**-1000=0x01700000,0*/
+
+double
+__ieee754_exp(double x)	/* default IEEE double exp */
+{
+	double y,hi=0.0,lo=0.0,c,t,twopk;
+	int32_t k=0,xsb;
+	u_int32_t hx;
+
+	GET_HIGH_WORD(hx,x);
+	xsb = (hx>>31)&1;		/* sign bit of x */
+	hx &= 0x7fffffff;		/* high word of |x| */
+
+    /* filter out non-finite argument */
+	if(hx >= 0x40862E42) {			/* if |x|>=709.78... */
+            if(hx>=0x7ff00000) {
+	        u_int32_t lx;
+		GET_LOW_WORD(lx,x);
+		if(((hx&0xfffff)|lx)!=0)
+		     return x+x; 		/* NaN */
+		else return (xsb==0)? x:0.0;	/* exp(+-inf)={inf,0} */
+	    }
+	    if(x > o_threshold) return huge*huge; /* overflow */
+	    if(x < u_threshold) return twom1000*twom1000; /* underflow */
+	}
+
+    /* argument reduction */
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */ 
+	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
+		hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
+	    } else {
+		k  = (int)(invln2*x+halF[xsb]);
+		t  = k;
+		hi = x - t*ln2HI[0];	/* t*ln2HI is exact here */
+		lo = t*ln2LO[0];
+	    }
+	    STRICT_ASSIGN(double, x, hi - lo);
+	} 
+	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
+	    if(huge+x>one) return one+x;/* trigger inexact */
+	}
+	else k = 0;
+
+    /* x is now in primary range */
+	t  = x*x;
+	if(k >= -1021)
+	    INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0);
+	else
+	    INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0);
+	c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+	if(k==0) 	return one-((x*c)/(c-2.0)-x); 
+	else 		y = one-((lo-(x*c)/(2.0-c))-hi);
+	if(k >= -1021) {
+	    if (k==1024) return y*2.0*0x1p1023;
+	    return y*twopk;
+	} else {
+	    return y*twopk*twom1000;
+	}
+}
+
+#ifdef SUPPORT_LONG_DOUBLE
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(exp, expl);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/e_fmod.c b/src/bsp/lk/lib/libm/e_fmod.c
new file mode 100644
index 0000000..720aa03
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_fmod.c
@@ -0,0 +1,132 @@
+
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* 
+ * __ieee754_fmod(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double one = 1.0, Zero[] = {0.0, -0.0,};
+
+double
+__ieee754_fmod(double x, double y)
+{
+	int32_t n,hx,hy,hz,ix,iy,sx,i;
+	u_int32_t lx,ly,lz;
+
+	EXTRACT_WORDS(hx,lx,x);
+	EXTRACT_WORDS(hy,ly,y);
+	sx = hx&0x80000000;		/* sign of x */
+	hx ^=sx;		/* |x| */
+	hy &= 0x7fffffff;	/* |y| */
+
+    /* purge off exception values */
+	if((hy|ly)==0||(hx>=0x7ff00000)||	/* y=0,or x not finite */
+	  ((hy|((ly|-ly)>>31))>0x7ff00000))	/* or y is NaN */
+	    return (x*y)/(x*y);
+	if(hx<=hy) {
+	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
+	    if(lx==ly) 
+		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
+	}
+
+    /* determine ix = ilogb(x) */
+	if(hx<0x00100000) {	/* subnormal x */
+	    if(hx==0) {
+		for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+	    } else {
+		for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
+	    }
+	} else ix = (hx>>20)-1023;
+
+    /* determine iy = ilogb(y) */
+	if(hy<0x00100000) {	/* subnormal y */
+	    if(hy==0) {
+		for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+	    } else {
+		for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
+	    }
+	} else iy = (hy>>20)-1023;
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	if(ix >= -1022) 
+	    hx = 0x00100000|(0x000fffff&hx);
+	else {		/* subnormal x, shift x to normal */
+	    n = -1022-ix;
+	    if(n<=31) {
+	        hx = (hx<<n)|(lx>>(32-n));
+	        lx <<= n;
+	    } else {
+		hx = lx<<(n-32);
+		lx = 0;
+	    }
+	}
+	if(iy >= -1022) 
+	    hy = 0x00100000|(0x000fffff&hy);
+	else {		/* subnormal y, shift y to normal */
+	    n = -1022-iy;
+	    if(n<=31) {
+	        hy = (hy<<n)|(ly>>(32-n));
+	        ly <<= n;
+	    } else {
+		hy = ly<<(n-32);
+		ly = 0;
+	    }
+	}
+
+    /* fix point fmod */
+	n = ix - iy;
+	while(n--) {
+	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	    if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
+	    else {
+	    	if((hz|lz)==0) 		/* return sign(x)*0 */
+		    return Zero[(u_int32_t)sx>>31];
+	    	hx = hz+hz+(lz>>31); lx = lz+lz;
+	    }
+	}
+	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	if(hz>=0) {hx=hz;lx=lz;}
+
+    /* convert back to floating value and restore the sign */
+	if((hx|lx)==0) 			/* return sign(x)*0 */
+	    return Zero[(u_int32_t)sx>>31];
+	while(hx<0x00100000) {		/* normalize x */
+	    hx = hx+hx+(lx>>31); lx = lx+lx;
+	    iy -= 1;
+	}
+	if(iy>= -1022) {	/* normalize output */
+	    hx = ((hx-0x00100000)|((iy+1023)<<20));
+	    INSERT_WORDS(x,hx|sx,lx);
+	} else {		/* subnormal output */
+	    n = -1022 - iy;
+	    if(n<=20) {
+		lx = (lx>>n)|((u_int32_t)hx<<(32-n));
+		hx >>= n;
+	    } else if (n<=31) {
+		lx = (hx<<(32-n))|(lx>>n); hx = sx;
+	    } else {
+		lx = hx>>(n-32); hx = sx;
+	    }
+	    INSERT_WORDS(x,hx|sx,lx);
+	    x *= one;		/* create necessary signal */
+	}
+	return x;		/* exact output */
+}
diff --git a/src/bsp/lk/lib/libm/e_log.c b/src/bsp/lk/lib/libm/e_log.c
new file mode 100644
index 0000000..9425b85
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_log.c
@@ -0,0 +1,149 @@
+
+/* @(#)e_log.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __ieee754_log(x)
+ * Return the logrithm of x
+ *
+ * Method :                  
+ *   1. Argument Reduction: find k and f such that 
+ *			x = 2^k * (1+f), 
+ *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *   2. Approximation of log(1+f).
+ *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *	     	 = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate 
+ * 	a polynomial of degree 14 to approximate R The maximum error 
+ *	of this polynomial approximation is bounded by 2**-58.45. In
+ *	other words,
+ *		        2      4      6      8      10      12      14
+ *	    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+ *  	(the values of Lg1 to Lg7 are listed in the program)
+ *	and
+ *	    |      2          14          |     -58.45
+ *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
+ *	    |                             |
+ *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *	In order to guarantee error in log below 1ulp, we compute log
+ *	by
+ *		log(1+f) = f - s*(f - R)	(if f is not too large)
+ *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
+ *	
+ *	3. Finally,  log(x) = k*ln2 + log(1+f).  
+ *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *	   Here ln2 is split into two floating point number: 
+ *			ln2_hi + ln2_lo,
+ *	   where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *	log(x) is NaN with signal if x < 0 (including -INF) ; 
+ *	log(+INF) is +INF; log(0) is -INF with signal;
+ *	log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
+ * to produce the hexadecimal values shown.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
+ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
+two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
+Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
+Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
+Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
+Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
+Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
+Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
+Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+
+static const double zero   =  0.0;
+static volatile double vzero = 0.0;
+
+double
+__ieee754_log(double x)
+{
+	double hfsq,f,s,z,R,w,t1,t2,dk;
+	int32_t k,hx,i,j;
+	u_int32_t lx;
+
+	EXTRACT_WORDS(hx,lx,x);
+
+	k=0;
+	if (hx < 0x00100000) {			/* x < 2**-1022  */
+	    if (((hx&0x7fffffff)|lx)==0) 
+		return -two54/vzero;		/* log(+-0)=-inf */
+	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
+	    k -= 54; x *= two54; /* subnormal number, scale up x */
+	    GET_HIGH_WORD(hx,x);
+	} 
+	if (hx >= 0x7ff00000) return x+x;
+	k += (hx>>20)-1023;
+	hx &= 0x000fffff;
+	i = (hx+0x95f64)&0x100000;
+	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
+	k += (i>>20);
+	f = x-1.0;
+	if((0x000fffff&(2+hx))<3) {	/* -2**-20 <= f < 2**-20 */
+	    if(f==zero) {
+		if(k==0) {
+		    return zero;
+		} else {
+		    dk=(double)k;
+		    return dk*ln2_hi+dk*ln2_lo;
+		}
+	    }
+	    R = f*f*(0.5-0.33333333333333333*f);
+	    if(k==0) return f-R; else {dk=(double)k;
+	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
+	}
+ 	s = f/(2.0+f); 
+	dk = (double)k;
+	z = s*s;
+	i = hx-0x6147a;
+	w = z*z;
+	j = 0x6b851-hx;
+	t1= w*(Lg2+w*(Lg4+w*Lg6)); 
+	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
+	i |= j;
+	R = t2+t1;
+	if(i>0) {
+	    hfsq=0.5*f*f;
+	    if(k==0) return f-(hfsq-s*(hfsq+R)); else
+		     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
+	} else {
+	    if(k==0) return f-s*(f-R); else
+		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
+	}
+}
+
+#ifdef SUPPORT_LONG_DOUBLE
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log, logl);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/e_pow.c b/src/bsp/lk/lib/libm/e_pow.c
new file mode 100644
index 0000000..d213132
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_pow.c
@@ -0,0 +1,304 @@
+/* @(#)e_pow.c 1.5 04/04/22 SMI */
+/*
+ * ====================================================
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_pow.c,v 1.11 2005/02/04 18:26:06 das Exp $";
+#endif
+
+/* __ieee754_pow(x,y) return x**y
+ *
+ *		      n
+ * Method:  Let x =  2   * (1+f)
+ *	1. Compute and return log2(x) in two pieces:
+ *		log2(x) = w1 + w2,
+ *	   where w1 has 53-24 = 29 bit trailing zeros.
+ *	2. Perform y*log2(x) = n+y' by simulating muti-precision 
+ *	   arithmetic, where |y'|<=0.5.
+ *	3. Return x**y = 2**n*exp(y'*log2)
+ *
+ * Special cases:
+ *	1.  (anything) ** 0  is 1
+ *	2.  (anything) ** 1  is itself
+ *	3.  (anything) ** NAN is NAN
+ *	4.  NAN ** (anything except 0) is NAN
+ *	5.  +-(|x| > 1) **  +INF is +INF
+ *	6.  +-(|x| > 1) **  -INF is +0
+ *	7.  +-(|x| < 1) **  +INF is +0
+ *	8.  +-(|x| < 1) **  -INF is +INF
+ *	9.  +-1         ** +-INF is NAN
+ *	10. +0 ** (+anything except 0, NAN)               is +0
+ *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
+ *	12. +0 ** (-anything except 0, NAN)               is +INF
+ *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
+ *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
+ *	15. +INF ** (+anything except 0,NAN) is +INF
+ *	16. +INF ** (-anything except 0,NAN) is +0
+ *	17. -INF ** (anything)  = -0 ** (-anything)
+ *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+ *	19. (-anything except 0 and inf) ** (non-integer) is NAN
+ *
+ * Accuracy:
+ *	pow(x,y) returns x**y nearly rounded. In particular
+ *			pow(integer,integer)
+ *	always returns the correct integer provided it is 
+ *	representable.
+ *
+ * Constants :
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
+ * to produce the hexadecimal values shown.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+bp[] = {1.0, 1.5,},
+dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
+dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
+zero    =  0.0,
+one	=  1.0,
+two	=  2.0,
+two53	=  9007199254740992.0,	/* 0x43400000, 0x00000000 */
+huge	=  1.0e300,
+tiny    =  1.0e-300,
+	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
+L1  =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
+L2  =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
+L3  =  3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
+L4  =  2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
+L5  =  2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
+L6  =  2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
+P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5   =  4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
+lg2  =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
+lg2_h  =  6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
+lg2_l  = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
+ovt =  8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
+cp    =  9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
+cp_h  =  9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
+cp_l  = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
+ivln2    =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
+ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
+ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
+
+double
+__ieee754_pow(double x, double y)
+{
+	double z,ax,z_h,z_l,p_h,p_l;
+	double y1,t1,t2,r,s,t,u,v,w;
+	int32_t i,j,k,yisint,n;
+	int32_t hx,hy,ix,iy;
+	u_int32_t lx,ly;
+
+	EXTRACT_WORDS(hx,lx,x);
+	EXTRACT_WORDS(hy,ly,y);
+	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
+
+    /* y==zero: x**0 = 1 */
+	if((iy|ly)==0) return one; 	
+
+    /* +-NaN return x+y */
+	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
+	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 
+		return x+y;	
+
+    /* determine if y is an odd int when x < 0
+     * yisint = 0	... y is not an integer
+     * yisint = 1	... y is an odd int
+     * yisint = 2	... y is an even int
+     */
+	yisint  = 0;
+	if(hx<0) {	
+	    if(iy>=0x43400000) yisint = 2; /* even integer y */
+	    else if(iy>=0x3ff00000) {
+		k = (iy>>20)-0x3ff;	   /* exponent */
+		if(k>20) {
+		    j = ly>>(52-k);
+		    if((j<<(52-k))==ly) yisint = 2-(j&1);
+		} else if(ly==0) {
+		    j = iy>>(20-k);
+		    if((j<<(20-k))==iy) yisint = 2-(j&1);
+		}
+	    }		
+	} 
+
+    /* special value of y */
+	if(ly==0) { 	
+	    if (iy==0x7ff00000) {	/* y is +-inf */
+	        if(((ix-0x3ff00000)|lx)==0)
+		    return  y - y;	/* inf**+-1 is NaN */
+	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
+		    return (hy>=0)? y: zero;
+	        else			/* (|x|<1)**-,+inf = inf,0 */
+		    return (hy<0)?-y: zero;
+	    } 
+	    if(iy==0x3ff00000) {	/* y is  +-1 */
+		if(hy<0) return one/x; else return x;
+	    }
+	    if(hy==0x40000000) return x*x; /* y is  2 */
+	    if(hy==0x3fe00000) {	/* y is  0.5 */
+		if(hx>=0)	/* x >= +0 */
+		return sqrt(x);	
+	    }
+	}
+
+	ax   = fabs(x);
+    /* special value of x */
+	if(lx==0) {
+	    if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
+		z = ax;			/*x is +-0,+-inf,+-1*/
+		if(hy<0) z = one/z;	/* z = (1/|x|) */
+		if(hx<0) {
+		    if(((ix-0x3ff00000)|yisint)==0) {
+			z = (z-z)/(z-z); /* (-1)**non-int is NaN */
+		    } else if(yisint==1) 
+			z = -z;		/* (x<0)**odd = -(|x|**odd) */
+		}
+		return z;
+	    }
+	}
+    
+    /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be
+	n = (hx>>31)+1;
+       but ANSI C says a right shift of a signed negative quantity is
+       implementation defined.  */
+	n = ((u_int32_t)hx>>31)-1;
+
+    /* (x<0)**(non-int) is NaN */
+	if((n|yisint)==0) return (x-x)/(x-x);
+
+	s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+	if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
+
+    /* |y| is huge */
+	if(iy>0x41e00000) { /* if |y| > 2**31 */
+	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow */
+		if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
+		if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
+	    }
+	/* over/underflow if x is not close to one */
+	    if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
+	    if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;
+	/* now |1-x| is tiny <= 2**-20, suffice to compute 
+	   log(x) by x-x^2/2+x^3/3-x^4/4 */
+	    t = ax-one;		/* t has 20 trailing zeros */
+	    w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
+	    u = ivln2_h*t;	/* ivln2_h has 21 sig. bits */
+	    v = t*ivln2_l-w*ivln2;
+	    t1 = u+v;
+	    SET_LOW_WORD(t1,0);
+	    t2 = v-(t1-u);
+	} else {
+	    double ss,s2,s_h,s_l,t_h,t_l;
+	    n = 0;
+	/* take care subnormal number */
+	    if(ix<0x00100000)
+		{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
+	    n  += ((ix)>>20)-0x3ff;
+	    j  = ix&0x000fffff;
+	/* determine interval */
+	    ix = j|0x3ff00000;		/* normalize ix */
+	    if(j<=0x3988E) k=0;		/* |x|<sqrt(3/2) */
+	    else if(j<0xBB67A) k=1;	/* |x|<sqrt(3)   */
+	    else {k=0;n+=1;ix -= 0x00100000;}
+	    SET_HIGH_WORD(ax,ix);
+
+	/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
+	    u = ax-bp[k];		/* bp[0]=1.0, bp[1]=1.5 */
+	    v = one/(ax+bp[k]);
+	    ss = u*v;
+	    s_h = ss;
+	    SET_LOW_WORD(s_h,0);
+	/* t_h=ax+bp[k] High */
+	    t_h = zero;
+	    SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
+	    t_l = ax - (t_h-bp[k]);
+	    s_l = v*((u-s_h*t_h)-s_h*t_l);
+	/* compute log(ax) */
+	    s2 = ss*ss;
+	    r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
+	    r += s_l*(s_h+ss);
+	    s2  = s_h*s_h;
+	    t_h = 3.0+s2+r;
+	    SET_LOW_WORD(t_h,0);
+	    t_l = r-((t_h-3.0)-s2);
+	/* u+v = ss*(1+...) */
+	    u = s_h*t_h;
+	    v = s_l*t_h+t_l*ss;
+	/* 2/(3log2)*(ss+...) */
+	    p_h = u+v;
+	    SET_LOW_WORD(p_h,0);
+	    p_l = v-(p_h-u);
+	    z_h = cp_h*p_h;		/* cp_h+cp_l = 2/(3*log2) */
+	    z_l = cp_l*p_h+p_l*cp+dp_l[k];
+	/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
+	    t = (double)n;
+	    t1 = (((z_h+z_l)+dp_h[k])+t);
+	    SET_LOW_WORD(t1,0);
+	    t2 = z_l-(((t1-t)-dp_h[k])-z_h);
+	}
+
+    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
+	y1  = y;
+	SET_LOW_WORD(y1,0);
+	p_l = (y-y1)*t1+y*t2;
+	p_h = y1*t1;
+	z = p_l+p_h;
+	EXTRACT_WORDS(j,i,z);
+	if (j>=0x40900000) {				/* z >= 1024 */
+	    if(((j-0x40900000)|i)!=0)			/* if z > 1024 */
+		return s*huge*huge;			/* overflow */
+	    else {
+		if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
+	    }
+	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
+	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
+		return s*tiny*tiny;		/* underflow */
+	    else {
+		if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
+	    }
+	}
+    /*
+     * compute 2**(p_h+p_l)
+     */
+	i = j&0x7fffffff;
+	k = (i>>20)-0x3ff;
+	n = 0;
+	if(i>0x3fe00000) {		/* if |z| > 0.5, set n = [z+0.5] */
+	    n = j+(0x00100000>>(k+1));
+	    k = ((n&0x7fffffff)>>20)-0x3ff;	/* new k for n */
+	    t = zero;
+	    SET_HIGH_WORD(t,n&~(0x000fffff>>k));
+	    n = ((n&0x000fffff)|0x00100000)>>(20-k);
+	    if(j<0) n = -n;
+	    p_h -= t;
+	} 
+	t = p_l+p_h;
+	SET_LOW_WORD(t,0);
+	u = t*lg2_h;
+	v = (p_l-(t-p_h))*lg2+t*lg2_l;
+	z = u+v;
+	w = v-(z-u);
+	t  = z*z;
+	t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+	r  = (z*t1)/(t1-two)-(w+z*w);
+	z  = one-(r-z);
+	GET_HIGH_WORD(j,z);
+	j += (n<<20);
+	if((j>>20)<=0) z = scalbn(z,n);	/* subnormal output */
+	else SET_HIGH_WORD(z,j);
+	return s*z;
+}
diff --git a/src/bsp/lk/lib/libm/e_powf.c b/src/bsp/lk/lib/libm/e_powf.c
new file mode 100644
index 0000000..41f08dd
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_powf.c
@@ -0,0 +1,247 @@
+/* e_powf.c -- float version of e_pow.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_powf.c,v 1.12 2004/06/01 19:33:30 bde Exp $";
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+static const float
+bp[] = {1.0, 1.5,},
+dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
+dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
+zero    =  0.0,
+one	=  1.0,
+two	=  2.0,
+two24	=  16777216.0,	/* 0x4b800000 */
+huge	=  1.0e30,
+tiny    =  1.0e-30,
+	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
+L1  =  6.0000002384e-01, /* 0x3f19999a */
+L2  =  4.2857143283e-01, /* 0x3edb6db7 */
+L3  =  3.3333334327e-01, /* 0x3eaaaaab */
+L4  =  2.7272811532e-01, /* 0x3e8ba305 */
+L5  =  2.3066075146e-01, /* 0x3e6c3255 */
+L6  =  2.0697501302e-01, /* 0x3e53f142 */
+P1   =  1.6666667163e-01, /* 0x3e2aaaab */
+P2   = -2.7777778450e-03, /* 0xbb360b61 */
+P3   =  6.6137559770e-05, /* 0x388ab355 */
+P4   = -1.6533901999e-06, /* 0xb5ddea0e */
+P5   =  4.1381369442e-08, /* 0x3331bb4c */
+lg2  =  6.9314718246e-01, /* 0x3f317218 */
+lg2_h  =  6.93145752e-01, /* 0x3f317200 */
+lg2_l  =  1.42860654e-06, /* 0x35bfbe8c */
+ovt =  4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
+cp    =  9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
+cp_h  =  9.6179199219e-01, /* 0x3f763800 =head of cp */
+cp_l  =  4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
+ivln2    =  1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
+ivln2_h  =  1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
+ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
+
+float
+__ieee754_powf(float x, float y)
+{
+	float z,ax,z_h,z_l,p_h,p_l;
+	float y1,t1,t2,r,s,sn,t,u,v,w;
+	int32_t i,j,k,yisint,n;
+	int32_t hx,hy,ix,iy,is;
+
+	GET_FLOAT_WORD(hx,x);
+	GET_FLOAT_WORD(hy,y);
+	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
+
+    /* y==zero: x**0 = 1 */
+	if(iy==0) return one;
+
+    /* +-NaN return x+y */
+	if(ix > 0x7f800000 ||
+	   iy > 0x7f800000)
+		return x+y;
+
+    /* determine if y is an odd int when x < 0
+     * yisint = 0	... y is not an integer
+     * yisint = 1	... y is an odd int
+     * yisint = 2	... y is an even int
+     */
+	yisint  = 0;
+	if(hx<0) {
+	    if(iy>=0x4b800000) yisint = 2; /* even integer y */
+	    else if(iy>=0x3f800000) {
+		k = (iy>>23)-0x7f;	   /* exponent */
+		j = iy>>(23-k);
+		if((j<<(23-k))==iy) yisint = 2-(j&1);
+	    }
+	}
+
+    /* special value of y */
+	if (iy==0x7f800000) {	/* y is +-inf */
+	    if (ix==0x3f800000)
+	        return  y - y;	/* inf**+-1 is NaN */
+	    else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
+	        return (hy>=0)? y: zero;
+	    else			/* (|x|<1)**-,+inf = inf,0 */
+	        return (hy<0)?-y: zero;
+	}
+	if(iy==0x3f800000) {	/* y is  +-1 */
+	    if(hy<0) return one/x; else return x;
+	}
+	if(hy==0x40000000) return x*x; /* y is  2 */
+	if(hy==0x3f000000) {	/* y is  0.5 */
+	    if(hx>=0)	/* x >= +0 */
+	    return __ieee754_sqrtf(x);
+	}
+
+	ax   = fabsf(x);
+    /* special value of x */
+	if(ix==0x7f800000||ix==0||ix==0x3f800000){
+	    z = ax;			/*x is +-0,+-inf,+-1*/
+	    if(hy<0) z = one/z;	/* z = (1/|x|) */
+	    if(hx<0) {
+		if(((ix-0x3f800000)|yisint)==0) {
+		    z = (z-z)/(z-z); /* (-1)**non-int is NaN */
+		} else if(yisint==1)
+		    z = -z;		/* (x<0)**odd = -(|x|**odd) */
+	    }
+	    return z;
+	}
+
+	n = ((u_int32_t)hx>>31)-1;
+
+    /* (x<0)**(non-int) is NaN */
+	if((n|yisint)==0) return (x-x)/(x-x);
+
+	sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+	if((n|(yisint-1))==0) sn = -one;/* (-ve)**(odd int) */
+
+    /* |y| is huge */
+	if(iy>0x4d000000) { /* if |y| > 2**27 */
+	/* over/underflow if x is not close to one */
+	    if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+	    if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
+	/* now |1-x| is tiny <= 2**-20, suffice to compute
+	   log(x) by x-x^2/2+x^3/3-x^4/4 */
+	    t = ax-1;		/* t has 20 trailing zeros */
+	    w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
+	    u = ivln2_h*t;	/* ivln2_h has 16 sig. bits */
+	    v = t*ivln2_l-w*ivln2;
+	    t1 = u+v;
+	    GET_FLOAT_WORD(is,t1);
+	    SET_FLOAT_WORD(t1,is&0xfffff000);
+	    t2 = v-(t1-u);
+	} else {
+	    float s2,s_h,s_l,t_h,t_l;
+	    n = 0;
+	/* take care subnormal number */
+	    if(ix<0x00800000)
+		{ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
+	    n  += ((ix)>>23)-0x7f;
+	    j  = ix&0x007fffff;
+	/* determine interval */
+	    ix = j|0x3f800000;		/* normalize ix */
+	    if(j<=0x1cc471) k=0;	/* |x|<sqrt(3/2) */
+	    else if(j<0x5db3d7) k=1;	/* |x|<sqrt(3)   */
+	    else {k=0;n+=1;ix -= 0x00800000;}
+	    SET_FLOAT_WORD(ax,ix);
+
+	/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
+	    u = ax-bp[k];		/* bp[0]=1.0, bp[1]=1.5 */
+	    v = one/(ax+bp[k]);
+	    s = u*v;
+	    s_h = s;
+	    GET_FLOAT_WORD(is,s_h);
+	    SET_FLOAT_WORD(s_h,is&0xfffff000);
+	/* t_h=ax+bp[k] High */
+	    is = ((ix>>1)&0xfffff000)|0x20000000;
+	    SET_FLOAT_WORD(t_h,is+0x00400000+(k<<21));
+	    t_l = ax - (t_h-bp[k]);
+	    s_l = v*((u-s_h*t_h)-s_h*t_l);
+	/* compute log(ax) */
+	    s2 = s*s;
+	    r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
+	    r += s_l*(s_h+s);
+	    s2  = s_h*s_h;
+	    t_h = (float)3.0+s2+r;
+	    GET_FLOAT_WORD(is,t_h);
+	    SET_FLOAT_WORD(t_h,is&0xfffff000);
+	    t_l = r-((t_h-(float)3.0)-s2);
+	/* u+v = s*(1+...) */
+	    u = s_h*t_h;
+	    v = s_l*t_h+t_l*s;
+	/* 2/(3log2)*(s+...) */
+	    p_h = u+v;
+	    GET_FLOAT_WORD(is,p_h);
+	    SET_FLOAT_WORD(p_h,is&0xfffff000);
+	    p_l = v-(p_h-u);
+	    z_h = cp_h*p_h;		/* cp_h+cp_l = 2/(3*log2) */
+	    z_l = cp_l*p_h+p_l*cp+dp_l[k];
+	/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
+	    t = (float)n;
+	    t1 = (((z_h+z_l)+dp_h[k])+t);
+	    GET_FLOAT_WORD(is,t1);
+	    SET_FLOAT_WORD(t1,is&0xfffff000);
+	    t2 = z_l-(((t1-t)-dp_h[k])-z_h);
+	}
+
+    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
+	GET_FLOAT_WORD(is,y);
+	SET_FLOAT_WORD(y1,is&0xfffff000);
+	p_l = (y-y1)*t1+y*t2;
+	p_h = y1*t1;
+	z = p_l+p_h;
+	GET_FLOAT_WORD(j,z);
+	if (j>0x43000000)				/* if z > 128 */
+	    return sn*huge*huge;			/* overflow */
+	else if (j==0x43000000) {			/* if z == 128 */
+	    if(p_l+ovt>z-p_h) return sn*huge*huge;	/* overflow */
+	}
+	else if ((j&0x7fffffff)>0x43160000)		/* z <= -150 */
+	    return sn*tiny*tiny;			/* underflow */
+	else if (j==0xc3160000){			/* z == -150 */
+	    if(p_l<=z-p_h) return sn*tiny*tiny;		/* underflow */
+	}
+    /*
+     * compute 2**(p_h+p_l)
+     */
+	i = j&0x7fffffff;
+	k = (i>>23)-0x7f;
+	n = 0;
+	if(i>0x3f000000) {		/* if |z| > 0.5, set n = [z+0.5] */
+	    n = j+(0x00800000>>(k+1));
+	    k = ((n&0x7fffffff)>>23)-0x7f;	/* new k for n */
+	    SET_FLOAT_WORD(t,n&~(0x007fffff>>k));
+	    n = ((n&0x007fffff)|0x00800000)>>(23-k);
+	    if(j<0) n = -n;
+	    p_h -= t;
+	}
+	t = p_l+p_h;
+	GET_FLOAT_WORD(is,t);
+	SET_FLOAT_WORD(t,is&0xffff8000);
+	u = t*lg2_h;
+	v = (p_l-(t-p_h))*lg2+t*lg2_l;
+	z = u+v;
+	w = v-(z-u);
+	t  = z*z;
+	t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+	r  = (z*t1)/(t1-two)-(w+z*w);
+	z  = one-(r-z);
+	GET_FLOAT_WORD(j,z);
+	j += (n<<23);
+	if((j>>23)<=0) z = scalbnf(z,n);	/* subnormal output */
+	else SET_FLOAT_WORD(z,j);
+	return sn*z;
+}
diff --git a/src/bsp/lk/lib/libm/e_rem_pio2.c b/src/bsp/lk/lib/libm/e_rem_pio2.c
new file mode 100644
index 0000000..86fde2d
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_rem_pio2.c
@@ -0,0 +1,189 @@
+
+/* @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __ieee754_rem_pio2(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * invpio2:  53 bits of 2/pi
+ * pio2_1:   first  33 bit of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 33 bit of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  33 bit of pi/2
+ * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
+pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
+pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
+pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
+pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
+pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
+
+#ifdef INLINE_REM_PIO2
+static __inline __always_inline
+#endif
+int
+__ieee754_rem_pio2(double x, double *y)
+{
+    double z,w,t,r,fn;
+    double tx[3],ty[2];
+    int32_t e0,i,j,nx,n,ix,hx;
+    u_int32_t low;
+
+    GET_HIGH_WORD(hx,x);        /* high word of x */
+    ix = hx&0x7fffffff;
+#if 0 /* Must be handled in caller. */
+    if (ix<=0x3fe921fb)  /* |x| ~<= pi/4 , no need for reduction */
+        {y[0] = x; y[1] = 0; return 0;}
+#endif
+    if (ix <= 0x400f6a7a) {     /* |x| ~<= 5pi/4 */
+        if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
+            goto medium;        /* cancellation -- use medium case */
+        if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
+            if (hx > 0) {
+                z = x - pio2_1; /* one round good to 85 bits */
+                y[0] = z - pio2_1t;
+                y[1] = (z-y[0])-pio2_1t;
+                return 1;
+            } else {
+                z = x + pio2_1;
+                y[0] = z + pio2_1t;
+                y[1] = (z-y[0])+pio2_1t;
+                return -1;
+            }
+        } else {
+            if (hx > 0) {
+                z = x - 2*pio2_1;
+                y[0] = z - 2*pio2_1t;
+                y[1] = (z-y[0])-2*pio2_1t;
+                return 2;
+            } else {
+                z = x + 2*pio2_1;
+                y[0] = z + 2*pio2_1t;
+                y[1] = (z-y[0])+2*pio2_1t;
+                return -2;
+            }
+        }
+    }
+    if (ix <= 0x401c463b) {     /* |x| ~<= 9pi/4 */
+        if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
+            if (ix == 0x4012d97c)   /* |x| ~= 3pi/2 */
+                goto medium;
+            if (hx > 0) {
+                z = x - 3*pio2_1;
+                y[0] = z - 3*pio2_1t;
+                y[1] = (z-y[0])-3*pio2_1t;
+                return 3;
+            } else {
+                z = x + 3*pio2_1;
+                y[0] = z + 3*pio2_1t;
+                y[1] = (z-y[0])+3*pio2_1t;
+                return -3;
+            }
+        } else {
+            if (ix == 0x401921fb)   /* |x| ~= 4pi/2 */
+                goto medium;
+            if (hx > 0) {
+                z = x - 4*pio2_1;
+                y[0] = z - 4*pio2_1t;
+                y[1] = (z-y[0])-4*pio2_1t;
+                return 4;
+            } else {
+                z = x + 4*pio2_1;
+                y[0] = z + 4*pio2_1t;
+                y[1] = (z-y[0])+4*pio2_1t;
+                return -4;
+            }
+        }
+    }
+    if (ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
+medium:
+        /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+        STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52);
+        fn = fn-0x1.8p52;
+#ifdef HAVE_EFFICIENT_IRINT
+        n  = irint(fn);
+#else
+        n  = (int32_t)fn;
+#endif
+        r  = x-fn*pio2_1;
+        w  = fn*pio2_1t;    /* 1st round good to 85 bit */
+        {
+            u_int32_t high;
+            j  = ix>>20;
+            y[0] = r-w;
+            GET_HIGH_WORD(high,y[0]);
+            i = j-((high>>20)&0x7ff);
+            if (i>16) { /* 2nd iteration needed, good to 118 */
+                t  = r;
+                w  = fn*pio2_2;
+                r  = t-w;
+                w  = fn*pio2_2t-((t-r)-w);
+                y[0] = r-w;
+                GET_HIGH_WORD(high,y[0]);
+                i = j-((high>>20)&0x7ff);
+                if (i>49)  { /* 3rd iteration need, 151 bits acc */
+                    t  = r; /* will cover all possible cases */
+                    w  = fn*pio2_3;
+                    r  = t-w;
+                    w  = fn*pio2_3t-((t-r)-w);
+                    y[0] = r-w;
+                }
+            }
+        }
+        y[1] = (r-y[0])-w;
+        return n;
+    }
+    /*
+     * all other (large) arguments
+     */
+    if (ix>=0x7ff00000) {       /* x is inf or NaN */
+        y[0]=y[1]=x-x;
+        return 0;
+    }
+    /* set z = scalbn(|x|,ilogb(x)-23) */
+    GET_LOW_WORD(low,x);
+    e0  = (ix>>20)-1046;    /* e0 = ilogb(z)-23; */
+    INSERT_WORDS(z, ix - ((int32_t)(e0<<20)), low);
+    for (i=0; i<2; i++) {
+        tx[i] = (double)((int32_t)(z));
+        z     = (z-tx[i])*two24;
+    }
+    tx[2] = z;
+    nx = 3;
+    while (tx[nx-1]==zero) nx--; /* skip zero term */
+    n  =  __kernel_rem_pio2(tx,ty,e0,nx,1);
+    if (hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;}
+    y[0] = ty[0];
+    y[1] = ty[1];
+    return n;
+}
diff --git a/src/bsp/lk/lib/libm/e_rem_pio2f.c b/src/bsp/lk/lib/libm/e_rem_pio2f.c
new file mode 100644
index 0000000..b068466
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_rem_pio2f.c
@@ -0,0 +1,86 @@
+/* e_rem_pio2f.c -- float version of e_rem_pio2.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Debugged and optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __ieee754_rem_pio2f(x,y)
+ *
+ * return the remainder of x rem pi/2 in *y
+ * use double precision for everything except passing x
+ * use __kernel_rem_pio2() for large x
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * invpio2:  53 bits of 2/pi
+ * pio2_1:   first 25 bits of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ */
+
+static const double
+invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+pio2_1  =  1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
+pio2_1t =  1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
+
+#ifdef INLINE_REM_PIO2F
+static __inline __always_inline
+#endif
+int
+__ieee754_rem_pio2f(float x, double *y)
+{
+    double w,r,fn;
+    double tx[1],ty[1];
+    float z;
+    int32_t e0,n,ix,hx;
+
+    GET_FLOAT_WORD(hx,x);
+    ix = hx&0x7fffffff;
+    /* 33+53 bit pi is good enough for medium size */
+    if (ix<0x4dc90fdb) {    /* |x| ~< 2^28*(pi/2), medium size */
+        /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+        STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52);
+        fn = fn-0x1.8p52;
+#ifdef HAVE_EFFICIENT_IRINT
+        n  = irint(fn);
+#else
+        n  = (int32_t)fn;
+#endif
+        r  = x-fn*pio2_1;
+        w  = fn*pio2_1t;
+        *y = r-w;
+        return n;
+    }
+    /*
+     * all other (large) arguments
+     */
+    if (ix>=0x7f800000) {       /* x is inf or NaN */
+        *y=x-x;
+        return 0;
+    }
+    /* set z = scalbn(|x|,ilogb(|x|)-23) */
+    e0 = (ix>>23)-150;      /* e0 = ilogb(|x|)-23; */
+    SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
+    tx[0] = z;
+    n  =  __kernel_rem_pio2(tx,ty,e0,1,0);
+    if (hx<0) {*y = -ty[0]; return -n;}
+    *y = ty[0];
+    return n;
+}
diff --git a/src/bsp/lk/lib/libm/e_sqrt.c b/src/bsp/lk/lib/libm/e_sqrt.c
new file mode 100644
index 0000000..1abcb01
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_sqrt.c
@@ -0,0 +1,468 @@
+
+/* @(#)e_sqrt.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __ieee754_sqrt(x)
+ * Return correctly rounded sqrt.
+ *           ------------------------------------------
+ *       |  Use the hardware sqrt if you have one |
+ *           ------------------------------------------
+ * Method:
+ *   Bit by bit method using integer arithmetic. (Slow, but portable)
+ *   1. Normalization
+ *  Scale x to y in [1,4) with even powers of 2:
+ *  find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
+ *      sqrt(x) = 2^k * sqrt(y)
+ *   2. Bit by bit computation
+ *  Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+ *       i                           0
+ *                                     i+1         2
+ *      s  = 2*q , and  y  =  2   * ( y - q  ).     (1)
+ *       i      i            i                 i
+ *
+ *  To compute q    from q , one checks whether
+ *          i+1       i
+ *
+ *                -(i+1) 2
+ *          (q + 2      ) <= y.         (2)
+ *                i
+ *                                -(i+1)
+ *  If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+ *                 i+1   i             i+1   i
+ *
+ *  With some algebric manipulation, it is not difficult to see
+ *  that (2) is equivalent to
+ *                             -(i+1)
+ *          s  +  2       <= y          (3)
+ *           i                i
+ *
+ *  The advantage of (3) is that s  and y  can be computed by
+ *                    i      i
+ *  the following recurrence formula:
+ *      if (3) is false
+ *
+ *      s     =  s  ,   y    = y   ;            (4)
+ *       i+1      i      i+1    i
+ *
+ *      otherwise,
+ *                         -i                     -(i+1)
+ *      s     =  s  + 2  ,  y    = y  -  s  - 2         (5)
+ *           i+1      i          i+1    i     i
+ *
+ *  One may easily use induction to prove (4) and (5).
+ *  Note. Since the left hand side of (3) contain only i+2 bits,
+ *        it does not necessary to do a full (53-bit) comparison
+ *        in (3).
+ *   3. Final rounding
+ *  After generating the 53 bits result, we compute one more bit.
+ *  Together with the remainder, we can decide whether the
+ *  result is exact, bigger than 1/2ulp, or less than 1/2ulp
+ *  (it will never equal to 1/2ulp).
+ *  The rounding mode can be detected by checking whether
+ *  huge + tiny is equal to huge, and whether huge - tiny is
+ *  equal to huge for some floating point number "huge" and "tiny".
+ *
+ * Special cases:
+ *  sqrt(+-0) = +-0     ... exact
+ *  sqrt(inf) = inf
+ *  sqrt(-ve) = NaN     ... with invalid signal
+ *  sqrt(NaN) = NaN     ... with invalid signal for signaling NaN
+ *
+ * Other methods : see the appended file at the end of the program below.
+ *---------------
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#if defined(LK) && ARCH_ARM && ARM_WITH_VFP
+/* use ARM w/VFP sqrt instruction */
+double
+__ieee754_sqrt(double x)
+{
+    double res;
+
+    __asm__("vsqrt.f64 %0, %1" : "=w"(res) : "w"(x));
+
+    return res;
+}
+
+#else
+
+static  const double    one = 1.0, tiny=1.0e-300;
+
+double
+__ieee754_sqrt(double x)
+{
+    double z;
+    int32_t sign = (int)0x80000000;
+    int32_t ix0,s0,q,m,t,i;
+    u_int32_t r,t1,s1,ix1,q1;
+
+    EXTRACT_WORDS(ix0,ix1,x);
+
+    /* take care of Inf and NaN */
+    if ((ix0&0x7ff00000)==0x7ff00000) {
+        return x*x+x;       /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+                       sqrt(-inf)=sNaN */
+    }
+    /* take care of zero */
+    if (ix0<=0) {
+        if (((ix0&(~sign))|ix1)==0) return x; /* sqrt(+-0) = +-0 */
+        else if (ix0<0)
+            return (x-x)/(x-x);     /* sqrt(-ve) = sNaN */
+    }
+    /* normalize x */
+    m = (ix0>>20);
+    if (m==0) {             /* subnormal x */
+        while (ix0==0) {
+            m -= 21;
+            ix0 |= (ix1>>11);
+            ix1 <<= 21;
+        }
+        for (i=0; (ix0&0x00100000)==0; i++) ix0<<=1;
+        m -= i-1;
+        ix0 |= (ix1>>(32-i));
+        ix1 <<= i;
+    }
+    m -= 1023;  /* unbias exponent */
+    ix0 = (ix0&0x000fffff)|0x00100000;
+    if (m&1) {  /* odd m, double x to make it even */
+        ix0 += ix0 + ((ix1&sign)>>31);
+        ix1 += ix1;
+    }
+    m >>= 1;    /* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+    ix0 += ix0 + ((ix1&sign)>>31);
+    ix1 += ix1;
+    q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
+    r = 0x00200000;     /* r = moving bit from right to left */
+
+    while (r!=0) {
+        t = s0+r;
+        if (t<=ix0) {
+            s0   = t+r;
+            ix0 -= t;
+            q   += r;
+        }
+        ix0 += ix0 + ((ix1&sign)>>31);
+        ix1 += ix1;
+        r>>=1;
+    }
+
+    r = sign;
+    while (r!=0) {
+        t1 = s1+r;
+        t  = s0;
+        if ((t<ix0)||((t==ix0)&&(t1<=ix1))) {
+            s1  = t1+r;
+            if (((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
+            ix0 -= t;
+            if (ix1 < t1) ix0 -= 1;
+            ix1 -= t1;
+            q1  += r;
+        }
+        ix0 += ix0 + ((ix1&sign)>>31);
+        ix1 += ix1;
+        r>>=1;
+    }
+
+    /* use floating add to find out rounding direction */
+    if ((ix0|ix1)!=0) {
+        z = one-tiny; /* trigger inexact flag */
+        if (z>=one) {
+            z = one+tiny;
+            if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
+            else if (z>one) {
+                if (q1==(u_int32_t)0xfffffffe) q+=1;
+                q1+=2;
+            } else
+                q1 += (q1&1);
+        }
+    }
+    ix0 = (q>>1)+0x3fe00000;
+    ix1 =  q1>>1;
+    if ((q&1)==1) ix1 |= sign;
+    ix0 += (m <<20);
+    INSERT_WORDS(z,ix0,ix1);
+    return z;
+}
+#endif
+
+#if SUPPORT_LONG_DOUBLE
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sqrt, sqrtl);
+#endif
+#endif
+/*
+Other methods  (use floating-point arithmetic)
+-------------
+(This is a copy of a drafted paper by Prof W. Kahan
+and K.C. Ng, written in May, 1986)
+
+    Two algorithms are given here to implement sqrt(x)
+    (IEEE double precision arithmetic) in software.
+    Both supply sqrt(x) correctly rounded. The first algorithm (in
+    Section A) uses newton iterations and involves four divisions.
+    The second one uses reciproot iterations to avoid division, but
+    requires more multiplications. Both algorithms need the ability
+    to chop results of arithmetic operations instead of round them,
+    and the INEXACT flag to indicate when an arithmetic operation
+    is executed exactly with no roundoff error, all part of the
+    standard (IEEE 754-1985). The ability to perform shift, add,
+    subtract and logical AND operations upon 32-bit words is needed
+    too, though not part of the standard.
+
+A.  sqrt(x) by Newton Iteration
+
+   (1)  Initial approximation
+
+    Let x0 and x1 be the leading and the trailing 32-bit words of
+    a floating point number x (in IEEE double format) respectively
+
+        1    11          52               ...widths
+       ------------------------------------------------------
+    x: |s|    e     |         f             |
+       ------------------------------------------------------
+          msb    lsb  msb                     lsb ...order
+
+
+         ------------------------        ------------------------
+    x0:  |s|   e    |    f1     |    x1: |          f2           |
+         ------------------------        ------------------------
+
+    By performing shifts and subtracts on x0 and x1 (both regarded
+    as integers), we obtain an 8-bit approximation of sqrt(x) as
+    follows.
+
+        k  := (x0>>1) + 0x1ff80000;
+        y0 := k - T1[31&(k>>15)].   ... y ~ sqrt(x) to 8 bits
+    Here k is a 32-bit integer and T1[] is an integer array containing
+    correction terms. Now magically the floating value of y (y's
+    leading 32-bit word is y0, the value of its trailing word is 0)
+    approximates sqrt(x) to almost 8-bit.
+
+    Value of T1:
+    static int T1[32]= {
+    0,  1024,   3062,   5746,   9193,   13348,  18162,  23592,
+    29598,  36145,  43202,  50740,  58733,  67158,  75992,  85215,
+    83599,  71378,  60428,  50647,  41945,  34246,  27478,  21581,
+    16499,  12183,  8588,   5674,   3403,   1742,   661,    130,};
+
+    (2) Iterative refinement
+
+    Apply Heron's rule three times to y, we have y approximates
+    sqrt(x) to within 1 ulp (Unit in the Last Place):
+
+        y := (y+x/y)/2      ... almost 17 sig. bits
+        y := (y+x/y)/2      ... almost 35 sig. bits
+        y := y-(y-x/y)/2    ... within 1 ulp
+
+
+    Remark 1.
+        Another way to improve y to within 1 ulp is:
+
+        y := (y+x/y)        ... almost 17 sig. bits to 2*sqrt(x)
+        y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
+
+                2
+                (x-y )*y
+        y := y + 2* ----------  ...within 1 ulp
+                   2
+                 3y  + x
+
+
+    This formula has one division fewer than the one above; however,
+    it requires more multiplications and additions. Also x must be
+    scaled in advance to avoid spurious overflow in evaluating the
+    expression 3y*y+x. Hence it is not recommended uless division
+    is slow. If division is very slow, then one should use the
+    reciproot algorithm given in section B.
+
+    (3) Final adjustment
+
+    By twiddling y's last bit it is possible to force y to be
+    correctly rounded according to the prevailing rounding mode
+    as follows. Let r and i be copies of the rounding mode and
+    inexact flag before entering the square root program. Also we
+    use the expression y+-ulp for the next representable floating
+    numbers (up and down) of y. Note that y+-ulp = either fixed
+    point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+    mode.
+
+        I := FALSE; ... reset INEXACT flag I
+        R := RZ;    ... set rounding mode to round-toward-zero
+        z := x/y;   ... chopped quotient, possibly inexact
+        If(not I) then {    ... if the quotient is exact
+            if(z=y) {
+                I := i;  ... restore inexact flag
+                R := r;  ... restore rounded mode
+                return sqrt(x):=y.
+            } else {
+            z := z - ulp;   ... special rounding
+            }
+        }
+        i := TRUE;      ... sqrt(x) is inexact
+        If (r=RN) then z=z+ulp  ... rounded-to-nearest
+        If (r=RP) then {    ... round-toward-+inf
+            y = y+ulp; z=z+ulp;
+        }
+        y := y+z;       ... chopped sum
+        y0:=y0-0x00100000;  ... y := y/2 is correctly rounded.
+            I := i;         ... restore inexact flag
+            R := r;         ... restore rounded mode
+            return sqrt(x):=y.
+
+    (4) Special cases
+
+    Square root of +inf, +-0, or NaN is itself;
+    Square root of a negative number is NaN with invalid signal.
+
+
+B.  sqrt(x) by Reciproot Iteration
+
+   (1)  Initial approximation
+
+    Let x0 and x1 be the leading and the trailing 32-bit words of
+    a floating point number x (in IEEE double format) respectively
+    (see section A). By performing shifs and subtracts on x0 and y0,
+    we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
+
+        k := 0x5fe80000 - (x0>>1);
+        y0:= k - T2[63&(k>>14)].    ... y ~ 1/sqrt(x) to 7.8 bits
+
+    Here k is a 32-bit integer and T2[] is an integer array
+    containing correction terms. Now magically the floating
+    value of y (y's leading 32-bit word is y0, the value of
+    its trailing word y1 is set to zero) approximates 1/sqrt(x)
+    to almost 7.8-bit.
+
+    Value of T2:
+    static int T2[64]= {
+    0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
+    0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
+    0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
+    0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
+    0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
+    0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
+    0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
+    0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
+
+    (2) Iterative refinement
+
+    Apply Reciproot iteration three times to y and multiply the
+    result by x to get an approximation z that matches sqrt(x)
+    to about 1 ulp. To be exact, we will have
+        -1ulp < sqrt(x)-z<1.0625ulp.
+
+    ... set rounding mode to Round-to-nearest
+       y := y*(1.5-0.5*x*y*y)   ... almost 15 sig. bits to 1/sqrt(x)
+       y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
+    ... special arrangement for better accuracy
+       z := x*y         ... 29 bits to sqrt(x), with z*y<1
+       z := z + 0.5*z*(1-z*y)   ... about 1 ulp to sqrt(x)
+
+    Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
+    (a) the term z*y in the final iteration is always less than 1;
+    (b) the error in the final result is biased upward so that
+        -1 ulp < sqrt(x) - z < 1.0625 ulp
+        instead of |sqrt(x)-z|<1.03125ulp.
+
+    (3) Final adjustment
+
+    By twiddling y's last bit it is possible to force y to be
+    correctly rounded according to the prevailing rounding mode
+    as follows. Let r and i be copies of the rounding mode and
+    inexact flag before entering the square root program. Also we
+    use the expression y+-ulp for the next representable floating
+    numbers (up and down) of y. Note that y+-ulp = either fixed
+    point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+    mode.
+
+    R := RZ;        ... set rounding mode to round-toward-zero
+    switch(r) {
+        case RN:        ... round-to-nearest
+           if(x<= z*(z-ulp)...chopped) z = z - ulp; else
+           if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
+           break;
+        case RZ:case RM:    ... round-to-zero or round-to--inf
+           R:=RP;       ... reset rounding mod to round-to-+inf
+           if(x<z*z ... rounded up) z = z - ulp; else
+           if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
+           break;
+        case RP:        ... round-to-+inf
+           if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
+           if(x>z*z ...chopped) z = z+ulp;
+           break;
+    }
+
+    Remark 3. The above comparisons can be done in fixed point. For
+    example, to compare x and w=z*z chopped, it suffices to compare
+    x1 and w1 (the trailing parts of x and w), regarding them as
+    two's complement integers.
+
+    ...Is z an exact square root?
+    To determine whether z is an exact square root of x, let z1 be the
+    trailing part of z, and also let x0 and x1 be the leading and
+    trailing parts of x.
+
+    If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
+        I := 1;     ... Raise Inexact flag: z is not exact
+    else {
+        j := 1 - [(x0>>20)&1]   ... j = logb(x) mod 2
+        k := z1 >> 26;      ... get z's 25-th and 26-th
+                        fraction bits
+        I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
+    }
+    R:= r       ... restore rounded mode
+    return sqrt(x):=z.
+
+    If multiplication is cheaper then the foregoing red tape, the
+    Inexact flag can be evaluated by
+
+        I := i;
+        I := (z*z!=x) or I.
+
+    Note that z*z can overwrite I; this value must be sensed if it is
+    True.
+
+    Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
+    zero.
+
+            --------------------
+        z1: |        f2        |
+            --------------------
+        bit 31         bit 0
+
+    Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
+    or even of logb(x) have the following relations:
+
+    -------------------------------------------------
+    bit 27,26 of z1     bit 1,0 of x1   logb(x)
+    -------------------------------------------------
+    00          00      odd and even
+    01          01      even
+    10          10      odd
+    10          00      even
+    11          01      even
+    -------------------------------------------------
+
+    (4) Special cases (see (4) of Section A).
+
+ */
+
diff --git a/src/bsp/lk/lib/libm/e_sqrtf.c b/src/bsp/lk/lib/libm/e_sqrtf.c
new file mode 100644
index 0000000..5d5412a
--- /dev/null
+++ b/src/bsp/lk/lib/libm/e_sqrtf.c
@@ -0,0 +1,106 @@
+/* e_sqrtf.c -- float version of e_sqrt.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD$";
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+static  const float one = 1.0, tiny=1.0e-30;
+
+#if defined(LK) && ARCH_ARM && ARM_WITH_VFP
+/* use ARM w/VFP sqrt instruction */
+float
+__ieee754_sqrtf(float x)
+{
+    float res;
+
+    __asm__("vsqrt.f32 %0, %1" : "=t"(res) : "t"(x));
+
+    return res;
+}
+
+#else
+
+float
+__ieee754_sqrtf(float x)
+{
+    float z;
+    int32_t sign = (int)0x80000000;
+    int32_t ix,s,q,m,t,i;
+    u_int32_t r;
+
+    GET_FLOAT_WORD(ix,x);
+
+    /* take care of Inf and NaN */
+    if ((ix&0x7f800000)==0x7f800000) {
+        return x*x+x;       /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+                       sqrt(-inf)=sNaN */
+    }
+    /* take care of zero */
+    if (ix<=0) {
+        if ((ix&(~sign))==0) return x; /* sqrt(+-0) = +-0 */
+        else if (ix<0)
+            return (x-x)/(x-x);     /* sqrt(-ve) = sNaN */
+    }
+    /* normalize x */
+    m = (ix>>23);
+    if (m==0) {             /* subnormal x */
+        for (i=0; (ix&0x00800000)==0; i++) ix<<=1;
+        m -= i-1;
+    }
+    m -= 127;   /* unbias exponent */
+    ix = (ix&0x007fffff)|0x00800000;
+    if (m&1) /* odd m, double x to make it even */
+        ix += ix;
+    m >>= 1;    /* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+    ix += ix;
+    q = s = 0;      /* q = sqrt(x) */
+    r = 0x01000000;     /* r = moving bit from right to left */
+
+    while (r!=0) {
+        t = s+r;
+        if (t<=ix) {
+            s    = t+r;
+            ix  -= t;
+            q   += r;
+        }
+        ix += ix;
+        r>>=1;
+    }
+
+    /* use floating add to find out rounding direction */
+    if (ix!=0) {
+        z = one-tiny; /* trigger inexact flag */
+        if (z>=one) {
+            z = one+tiny;
+            if (z>one)
+                q += 2;
+            else
+                q += (q&1);
+        }
+    }
+    ix = (q>>1)+0x3f000000;
+    ix += (m <<23);
+    SET_FLOAT_WORD(z,ix);
+    return z;
+}
+
+#endif
+
diff --git a/src/bsp/lk/lib/libm/include/math.h b/src/bsp/lk/lib/libm/include/math.h
new file mode 100644
index 0000000..64c166f
--- /dev/null
+++ b/src/bsp/lk/lib/libm/include/math.h
@@ -0,0 +1,474 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ * $FreeBSD$
+ */
+
+#ifndef _MATH_H_
+#define _MATH_H_
+
+#include <sys/cdefs.h>
+#include <limits.h>
+
+__BEGIN_DECLS
+#pragma GCC visibility push(default)
+
+/*
+ * ANSI/POSIX
+ */
+extern const union __infinity_un {
+    unsigned char   __uc[8];
+    double      __ud;
+} __infinity;
+
+extern const union __nan_un {
+    unsigned char   __uc[sizeof(float)];
+    float       __uf;
+} __nan;
+
+#if __GNUC_PREREQ__(3, 3) || (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 800)
+#define __MATH_BUILTIN_CONSTANTS
+#endif
+
+#if __GNUC_PREREQ__(3, 0) && !defined(__INTEL_COMPILER)
+#define __MATH_BUILTIN_RELOPS
+#endif
+
+#ifdef __MATH_BUILTIN_CONSTANTS
+#define HUGE_VAL    __builtin_huge_val()
+#else
+#define HUGE_VAL    (__infinity.__ud)
+#endif
+
+#if __ISO_C_VISIBLE >= 1999
+#define FP_ILOGB0   (-INT_MAX) /* Android-changed */
+#define FP_ILOGBNAN INT_MAX /* Android-changed */
+
+#ifdef __MATH_BUILTIN_CONSTANTS
+#define HUGE_VALF   __builtin_huge_valf()
+#define HUGE_VALL   __builtin_huge_vall()
+#define INFINITY    __builtin_inff()
+#define NAN     __builtin_nanf("")
+#else
+#define HUGE_VALF   (float)HUGE_VAL
+#define HUGE_VALL   (long double)HUGE_VAL
+#define INFINITY    HUGE_VALF
+#define NAN     (__nan.__uf)
+#endif /* __MATH_BUILTIN_CONSTANTS */
+
+#define MATH_ERRNO  1
+#define MATH_ERREXCEPT  2
+#define math_errhandling    MATH_ERREXCEPT
+
+#define FP_FAST_FMAF    1
+#ifdef __ia64__
+#define FP_FAST_FMA 1
+#define FP_FAST_FMAL    1
+#endif
+
+/* Symbolic constants to classify floating point numbers. */
+#define FP_INFINITE 0x01
+#define FP_NAN      0x02
+#define FP_NORMAL   0x04
+#define FP_SUBNORMAL    0x08
+#define FP_ZERO     0x10
+#define fpclassify(x) \
+    ((sizeof (x) == sizeof (float)) ? __fpclassifyf(x) \
+    : (sizeof (x) == sizeof (double)) ? __fpclassifyd(x) \
+    : __fpclassifyl(x))
+
+#define isfinite(x)                 \
+    ((sizeof (x) == sizeof (float)) ? __isfinitef(x)    \
+    : (sizeof (x) == sizeof (double)) ? __isfinite(x)   \
+    : __isfinitel(x))
+#define isinf(x)                    \
+    ((sizeof (x) == sizeof (float)) ? __isinff(x)   \
+    : (sizeof (x) == sizeof (double)) ? isinf(x)    \
+    : __isinfl(x))
+#define isnan(x)                    \
+    ((sizeof (x) == sizeof (float)) ? __isnanf(x)   \
+    : (sizeof (x) == sizeof (double)) ? isnan(x)    \
+    : __isnanl(x))
+#define isnormal(x)                 \
+    ((sizeof (x) == sizeof (float)) ? __isnormalf(x)    \
+    : (sizeof (x) == sizeof (double)) ? __isnormal(x)   \
+    : __isnormall(x))
+
+#ifdef __MATH_BUILTIN_RELOPS
+#define isgreater(x, y)     __builtin_isgreater((x), (y))
+#define isgreaterequal(x, y)    __builtin_isgreaterequal((x), (y))
+#define isless(x, y)        __builtin_isless((x), (y))
+#define islessequal(x, y)   __builtin_islessequal((x), (y))
+#define islessgreater(x, y) __builtin_islessgreater((x), (y))
+#define isunordered(x, y)   __builtin_isunordered((x), (y))
+#else
+#define isgreater(x, y)     (!isunordered((x), (y)) && (x) > (y))
+#define isgreaterequal(x, y)    (!isunordered((x), (y)) && (x) >= (y))
+#define isless(x, y)        (!isunordered((x), (y)) && (x) < (y))
+#define islessequal(x, y)   (!isunordered((x), (y)) && (x) <= (y))
+#define islessgreater(x, y) (!isunordered((x), (y)) && \
+                    ((x) > (y) || (y) > (x)))
+#define isunordered(x, y)   (isnan(x) || isnan(y))
+#endif /* __MATH_BUILTIN_RELOPS */
+
+#define signbit(x)                  \
+    ((sizeof (x) == sizeof (float)) ? __signbitf(x) \
+    : (sizeof (x) == sizeof (double)) ? __signbit(x)    \
+    : __signbitl(x))
+
+typedef double __double_t;
+typedef __double_t double_t;
+typedef float __float_t;
+typedef __float_t float_t;
+#endif /* __ISO_C_VISIBLE >= 1999 */
+
+/*
+ * XOPEN/SVID
+ */
+#if __BSD_VISIBLE || __XSI_VISIBLE
+#define M_E     2.7182818284590452354   /* e */
+#define M_LOG2E     1.4426950408889634074   /* log 2e */
+#define M_LOG10E    0.43429448190325182765  /* log 10e */
+#define M_LN2       0.69314718055994530942  /* log e2 */
+#define M_LN10      2.30258509299404568402  /* log e10 */
+#define M_PI        3.14159265358979323846  /* pi */
+#define M_PI_2      1.57079632679489661923  /* pi/2 */
+#define M_PI_4      0.78539816339744830962  /* pi/4 */
+#define M_1_PI      0.31830988618379067154  /* 1/pi */
+#define M_2_PI      0.63661977236758134308  /* 2/pi */
+#define M_2_SQRTPI  1.12837916709551257390  /* 2/sqrt(pi) */
+#define M_SQRT2     1.41421356237309504880  /* sqrt(2) */
+#define M_SQRT1_2   0.70710678118654752440  /* 1/sqrt(2) */
+
+#define MAXFLOAT    ((float)3.40282346638528860e+38)
+extern int signgam;
+#endif /* __BSD_VISIBLE || __XSI_VISIBLE */
+
+#if __BSD_VISIBLE
+#if 0
+/* Old value from 4.4BSD-Lite math.h; this is probably better. */
+#define HUGE        HUGE_VAL
+#else
+#define HUGE        MAXFLOAT
+#endif
+#endif /* __BSD_VISIBLE */
+
+/*
+ * Most of these functions depend on the rounding mode and have the side
+ * effect of raising floating-point exceptions, so they are not declared
+ * as __pure2.  In C99, FENV_ACCESS affects the purity of these functions.
+ */
+
+/*
+ * ANSI/POSIX
+ */
+int __fpclassifyd(double) __pure2;
+int __fpclassifyf(float) __pure2;
+int __fpclassifyl(long double) __pure2;
+int __isfinitef(float) __pure2;
+int __isfinite(double) __pure2;
+int __isfinitel(long double) __pure2;
+int __isinff(float) __pure2;
+int __isinfl(long double) __pure2;
+int __isnanf(float) __pure2;
+int __isnanl(long double) __pure2;
+int __isnormalf(float) __pure2;
+int __isnormal(double) __pure2;
+int __isnormall(long double) __pure2;
+int __signbit(double) __pure2;
+int __signbitf(float) __pure2;
+int __signbitl(long double) __pure2;
+
+double  acos(double);
+double  asin(double);
+double  atan(double);
+double  atan2(double, double);
+double  cos(double);
+double  sin(double);
+double  tan(double);
+
+double  cosh(double);
+double  sinh(double);
+double  tanh(double);
+
+double  exp(double);
+double  frexp(double, int *);   /* fundamentally !__pure2 */
+double  ldexp(double, int);
+double  log(double);
+double  log10(double);
+double  modf(double, double *); /* fundamentally !__pure2 */
+
+double  pow(double, double);
+double  sqrt(double);
+
+double  ceil(double);
+double  fabs(double) __pure2;
+double  floor(double);
+double  fmod(double, double);
+
+/*
+ * These functions are not in C90.
+ */
+#if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE
+double  acosh(double);
+double  asinh(double);
+double  atanh(double);
+double  cbrt(double);
+double  erf(double);
+double  erfc(double);
+double  exp2(double);
+double  expm1(double);
+double  fma(double, double, double);
+double  hypot(double, double);
+int ilogb(double) __pure2;
+int (isinf)(double) __pure2;
+int (isnan)(double) __pure2;
+double  lgamma(double);
+long long llrint(double);
+long long llround(double);
+double  log1p(double);
+double  log2(double);
+double  logb(double);
+long    lrint(double);
+long    lround(double);
+double  nan(const char *) __pure2;
+double  nextafter(double, double);
+double  remainder(double, double);
+double  remquo(double, double, int *);
+double  rint(double);
+#endif /* __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999 || __XSI_VISIBLE */
+
+#if __BSD_VISIBLE || __XSI_VISIBLE
+double  j0(double);
+double  j1(double);
+double  jn(int, double);
+double  y0(double);
+double  y1(double);
+double  yn(int, double);
+
+#if __XSI_VISIBLE <= 500 || __BSD_VISIBLE
+double  gamma(double);
+#endif
+
+#if __XSI_VISIBLE <= 600 || __BSD_VISIBLE
+double  scalb(double, double);
+#endif
+#endif /* __BSD_VISIBLE || __XSI_VISIBLE */
+
+#if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999
+double  copysign(double, double) __pure2;
+double  fdim(double, double);
+double  fmax(double, double) __pure2;
+double  fmin(double, double) __pure2;
+double  nearbyint(double);
+double  round(double);
+double  scalbln(double, long);
+double  scalbn(double, int);
+double  tgamma(double);
+double  trunc(double);
+#endif
+
+/*
+ * BSD math library entry points
+ */
+#if __BSD_VISIBLE
+double  drem(double, double);
+int finite(double) __pure2;
+int isnanf(float) __pure2;
+long double significandl(long double);
+
+/*
+ * Reentrant version of gamma & lgamma; passes signgam back by reference
+ * as the second argument; user must allocate space for signgam.
+ */
+double  gamma_r(double, int *);
+double  lgamma_r(double, int *);
+
+/*
+ * IEEE Test Vector
+ */
+double  significand(double);
+#endif /* __BSD_VISIBLE */
+
+/* float versions of ANSI/POSIX functions */
+#if __ISO_C_VISIBLE >= 1999
+float   acosf(float);
+float   asinf(float);
+float   atanf(float);
+float   atan2f(float, float);
+float   cosf(float);
+float   sinf(float);
+float   tanf(float);
+
+float   coshf(float);
+float   sinhf(float);
+float   tanhf(float);
+
+float   exp2f(float);
+float   expf(float);
+float   expm1f(float);
+float   frexpf(float, int *);   /* fundamentally !__pure2 */
+int ilogbf(float) __pure2;
+float   ldexpf(float, int);
+float   log10f(float);
+float   log1pf(float);
+float   log2f(float);
+float   logf(float);
+float   modff(float, float *);  /* fundamentally !__pure2 */
+
+float   powf(float, float);
+float   sqrtf(float);
+
+float   ceilf(float);
+float   fabsf(float) __pure2;
+float   floorf(float);
+float   fmodf(float, float);
+float   roundf(float);
+
+float   erff(float);
+float   erfcf(float);
+float   hypotf(float, float);
+float   lgammaf(float);
+float   tgammaf(float);
+
+float   acoshf(float);
+float   asinhf(float);
+float   atanhf(float);
+float   cbrtf(float);
+float   logbf(float);
+float   copysignf(float, float) __pure2;
+long long llrintf(float);
+long long llroundf(float);
+long    lrintf(float);
+long    lroundf(float);
+float   nanf(const char *) __pure2;
+float   nearbyintf(float);
+float   nextafterf(float, float);
+float   remainderf(float, float);
+float   remquof(float, float, int *);
+float   rintf(float);
+float   scalblnf(float, long);
+float   scalbnf(float, int);
+float   truncf(float);
+
+float   fdimf(float, float);
+float   fmaf(float, float, float);
+float   fmaxf(float, float) __pure2;
+float   fminf(float, float) __pure2;
+#endif
+
+/*
+ * float versions of BSD math library entry points
+ */
+#if __BSD_VISIBLE
+float   dremf(float, float);
+int finitef(float) __pure2;
+float   gammaf(float);
+float   j0f(float);
+float   j1f(float);
+float   jnf(int, float);
+float   scalbf(float, float);
+float   y0f(float);
+float   y1f(float);
+float   ynf(int, float);
+
+/*
+ * Float versions of reentrant version of gamma & lgamma; passes
+ * signgam back by reference as the second argument; user must
+ * allocate space for signgam.
+ */
+float   gammaf_r(float, int *);
+float   lgammaf_r(float, int *);
+
+/*
+ * float version of IEEE Test Vector
+ */
+float   significandf(float);
+#endif  /* __BSD_VISIBLE */
+
+/*
+ * long double versions of ISO/POSIX math functions
+ */
+#if __ISO_C_VISIBLE >= 1999
+long double acoshl(long double);
+long double acosl(long double);
+long double asinhl(long double);
+long double asinl(long double);
+long double atan2l(long double, long double);
+long double atanhl(long double);
+long double atanl(long double);
+long double cbrtl(long double);
+long double ceill(long double);
+long double copysignl(long double, long double) __pure2;
+long double coshl(long double);
+long double cosl(long double);
+long double erfcl(long double);
+long double erfl(long double);
+long double exp2l(long double);
+long double expl(long double);
+long double expm1l(long double);
+long double fabsl(long double) __pure2;
+long double fdiml(long double, long double);
+long double floorl(long double);
+long double fmal(long double, long double, long double);
+long double fmaxl(long double, long double) __pure2;
+long double fminl(long double, long double) __pure2;
+long double fmodl(long double, long double);
+long double frexpl(long double value, int *); /* fundamentally !__pure2 */
+long double hypotl(long double, long double);
+int     ilogbl(long double) __pure2;
+long double ldexpl(long double, int);
+long double lgammal(long double);
+long long   llrintl(long double);
+long long   llroundl(long double);
+long double log10l(long double);
+long double log1pl(long double);
+long double log2l(long double);
+long double logbl(long double);
+long double logl(long double);
+long        lrintl(long double);
+long        lroundl(long double);
+long double modfl(long double, long double *); /* fundamentally !__pure2 */
+long double nanl(const char *) __pure2;
+long double nearbyintl(long double);
+long double nextafterl(long double, long double);
+double      nexttoward(double, long double);
+float       nexttowardf(float, long double);
+long double nexttowardl(long double, long double);
+long double powl(long double, long double);
+long double remainderl(long double, long double);
+long double remquol(long double, long double, int *);
+long double rintl(long double);
+long double roundl(long double);
+long double scalblnl(long double, long);
+long double scalbnl(long double, int);
+long double sinhl(long double);
+long double sinl(long double);
+long double sqrtl(long double);
+long double tanhl(long double);
+long double tanl(long double);
+long double tgammal(long double);
+long double truncl(long double);
+
+#endif /* __ISO_C_VISIBLE >= 1999 */
+
+#if defined(_GNU_SOURCE)
+void sincos(double, double*, double*);
+void sincosf(float, float*, float*);
+void sincosl(long double, long double*, long double*);
+#endif /* _GNU_SOURCE */
+
+#pragma GCC visibility pop
+__END_DECLS
+
+#endif /* !_MATH_H_ */
diff --git a/src/bsp/lk/lib/libm/include/sys/cdefs.h b/src/bsp/lk/lib/libm/include/sys/cdefs.h
new file mode 100644
index 0000000..4681e7c
--- /dev/null
+++ b/src/bsp/lk/lib/libm/include/sys/cdefs.h
@@ -0,0 +1,144 @@
+/*  $NetBSD: cdefs.h,v 1.58 2004/12/11 05:59:00 christos Exp $  */
+
+/*
+ * Copyright (c) 1991, 1993
+ *  The Regents of the University of California.  All rights reserved.
+ *
+ * This code is derived from software contributed to Berkeley by
+ * Berkeley Software Design, Inc.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ *  @(#)cdefs.h 8.8 (Berkeley) 1/9/95
+ */
+
+/*
+ * This is a minimum version of cdefs.h for LK, this might be
+ * modified if necessary
+ */
+#ifndef _SYS_CDEFS_H_
+#define _SYS_CDEFS_H_
+
+/*-
+ * Deal with _ANSI_SOURCE:
+ * If it is defined, and no other compilation environment is explicitly
+ * requested, then define our internal feature-test macros to zero.  This
+ * makes no difference to the preprocessor (undefined symbols in preprocessing
+ * expressions are defined to have value zero), but makes it more convenient for
+ * a test program to print out the values.
+ *
+ * If a program mistakenly defines _ANSI_SOURCE and some other macro such as
+ * _POSIX_C_SOURCE, we will assume that it wants the broader compilation
+ * environment (and in fact we will never get here).
+ */
+#if defined(_ANSI_SOURCE)   /* Hide almost everything. */
+#define __POSIX_VISIBLE     0
+#define __XSI_VISIBLE       0
+#define __BSD_VISIBLE       0
+#define __ISO_C_VISIBLE     1990
+#elif defined(_C99_SOURCE)  /* Localism to specify strict C99 env. */
+#define __POSIX_VISIBLE     0
+#define __XSI_VISIBLE       0
+#define __BSD_VISIBLE       0
+#define __ISO_C_VISIBLE     1999
+#else               /* Default environment: show everything. */
+#define __POSIX_VISIBLE     200809
+#define __XSI_VISIBLE       700
+#define __BSD_VISIBLE       1
+#define __ISO_C_VISIBLE     1999
+#endif
+
+/*
+ * Default values.
+ */
+#ifndef __XPG_VISIBLE
+# define __XPG_VISIBLE          700
+#endif
+#ifndef __POSIX_VISIBLE
+# define __POSIX_VISIBLE        200809
+#endif
+#ifndef __ISO_C_VISIBLE
+# define __ISO_C_VISIBLE        1999
+#endif
+#ifndef __BSD_VISIBLE
+# define __BSD_VISIBLE          1
+#endif
+
+
+/*
+ * Some of the FreeBSD sources used in Bionic need this.
+ * Originally, this is used to embed the rcs versions of each source file
+ * in the generated binary. We certainly don't want this in Bionic.
+ */
+#define __FBSDID(s) /* nothing */
+
+#define __pure2 __attribute__((__const__)) /* Android-added: used by FreeBSD libm */
+
+/*
+ * Macro to test if we're using a GNU C compiler of a specific vintage
+ * or later, for e.g. features that appeared in a particular version
+ * of GNU C.  Usage:
+ *
+ *  #if __GNUC_PREREQ__(major, minor)
+ *  ...cool feature...
+ *  #else
+ *  ...delete feature...
+ *  #endif
+ */
+#ifdef __GNUC__
+#define __GNUC_PREREQ__(x, y)                       \
+    ((__GNUC__ == (x) && __GNUC_MINOR__ >= (y)) ||          \
+     (__GNUC__ > (x)))
+#else
+#define __GNUC_PREREQ__(x, y)   0
+#endif
+
+#if defined(__cplusplus)
+#define __BEGIN_DECLS       extern "C" {
+#define __END_DECLS     }
+#define __static_cast(x,y)  static_cast<x>(y)
+#else
+#define __BEGIN_DECLS
+#define __END_DECLS
+#define __static_cast(x,y)  (x)y
+#endif
+
+#if defined(__cplusplus)
+#define __inline    inline      /* convert to C++ keyword */
+#else
+#if !defined(__GNUC__) && !defined(__lint__)
+#define __inline            /* delete GCC keyword */
+#endif /* !__GNUC__  && !__lint__ */
+#endif /* !__cplusplus */
+
+#if __GNUC_PREREQ__(3, 1)
+#define __always_inline __attribute__((__always_inline__))
+#else
+#define __always_inline
+#endif
+
+
+
+#endif /* !_SYS_CDEFS_H_ */
diff --git a/src/bsp/lk/lib/libm/k_cos.c b/src/bsp/lk/lib/libm/k_cos.c
new file mode 100644
index 0000000..b6412bf
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_cos.c
@@ -0,0 +1,79 @@
+
+/* @(#)k_cos.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * __kernel_cos( x,  y )
+ * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ *
+ * Algorithm
+ *  1. Since cos(-x) = cos(x), we need only to consider positive x.
+ *  2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
+ *  3. cos(x) is approximated by a polynomial of degree 14 on
+ *     [0,pi/4]
+ *                           4            14
+ *      cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
+ *     where the remez error is
+ *
+ *  |              2     4     6     8     10    12     14 |     -58
+ *  |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
+ *  |                                      |
+ *
+ *                 4     6     8     10    12     14
+ *  4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
+ *         cos(x) ~ 1 - x*x/2 + r
+ *     since cos(x+y) ~ cos(x) - sin(x)*y
+ *            ~ cos(x) - x*y,
+ *     a correction term is necessary in cos(x) and hence
+ *      cos(x+y) = 1 - (x*x/2 - (r - x*y))
+ *     For better accuracy, rearrange to
+ *      cos(x+y) ~ w + (tmp + (r-x*y))
+ *     where w = 1 - x*x/2 and tmp is a tiny correction term
+ *     (1 - x*x/2 == w + tmp exactly in infinite precision).
+ *     The exactness of w + tmp in infinite precision depends on w
+ *     and tmp having the same precision as x.  If they have extra
+ *     precision due to compiler bugs, then the extra precision is
+ *     only good provided it is retained in all terms of the final
+ *     expression for cos().  Retention happens in all cases tested
+ *     under FreeBSD, so don't pessimize things by forcibly clipping
+ *     any extra precision in w.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
+C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
+C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
+C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
+C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
+C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+double
+__kernel_cos(double x, double y)
+{
+    double hz,z,r,w;
+
+    z  = x*x;
+    w  = z*z;
+    r  = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
+    hz = 0.5*z;
+    w  = one-hz;
+    return w + (((one-w)-hz) + (z*r-x*y));
+}
diff --git a/src/bsp/lk/lib/libm/k_cosf.c b/src/bsp/lk/lib/libm/k_cosf.c
new file mode 100644
index 0000000..3126b47
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_cosf.c
@@ -0,0 +1,46 @@
+/* k_cosf.c -- float version of k_cos.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Debugged and optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef INLINE_KERNEL_COSDF
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
+static const double
+one =  1.0,
+C0  = -0x1ffffffd0c5e81.0p-54,  /* -0.499999997251031003120 */
+C1  =  0x155553e1053a42.0p-57,  /*  0.0416666233237390631894 */
+C2  = -0x16c087e80f1e27.0p-62,  /* -0.00138867637746099294692 */
+C3  =  0x199342e0ee5069.0p-68;  /*  0.0000243904487962774090654 */
+
+#ifdef INLINE_KERNEL_COSDF
+static __inline
+#endif
+float
+__kernel_cosdf(double x)
+{
+    double r, w, z;
+
+    /* Try to optimize for parallel evaluation as in k_tanf.c. */
+    z = x*x;
+    w = z*z;
+    r = C2+z*C3;
+    return ((one+z*C0) + w*C1) + (w*z)*r;
+}
diff --git a/src/bsp/lk/lib/libm/k_rem_pio2.c b/src/bsp/lk/lib/libm/k_rem_pio2.c
new file mode 100644
index 0000000..03559b4
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_rem_pio2.c
@@ -0,0 +1,457 @@
+
+/* @(#)k_rem_pio2.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * __kernel_rem_pio2(x,y,e0,nx,prec)
+ * double x[],y[]; int e0,nx,prec;
+ *
+ * __kernel_rem_pio2 return the last three digits of N with
+ *      y = x - N*pi/2
+ * so that |y| < pi/2.
+ *
+ * The method is to compute the integer (mod 8) and fraction parts of
+ * (2/pi)*x without doing the full multiplication. In general we
+ * skip the part of the product that are known to be a huge integer (
+ * more accurately, = 0 mod 8 ). Thus the number of operations are
+ * independent of the exponent of the input.
+ *
+ * (2/pi) is represented by an array of 24-bit integers in ipio2[].
+ *
+ * Input parameters:
+ *  x[] The input value (must be positive) is broken into nx
+ *      pieces of 24-bit integers in double precision format.
+ *      x[i] will be the i-th 24 bit of x. The scaled exponent
+ *      of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
+ *      match x's up to 24 bits.
+ *
+ *      Example of breaking a double positive z into x[0]+x[1]+x[2]:
+ *          e0 = ilogb(z)-23
+ *          z  = scalbn(z,-e0)
+ *      for i = 0,1,2
+ *          x[i] = floor(z)
+ *          z    = (z-x[i])*2**24
+ *
+ *
+ *  y[] output result in an array of double precision numbers.
+ *      The dimension of y[] is:
+ *          24-bit  precision   1
+ *          53-bit  precision   2
+ *          64-bit  precision   2
+ *          113-bit precision   3
+ *      The actual value is the sum of them. Thus for 113-bit
+ *      precison, one may have to do something like:
+ *
+ *      long double t,w,r_head, r_tail;
+ *      t = (long double)y[2] + (long double)y[1];
+ *      w = (long double)y[0];
+ *      r_head = t+w;
+ *      r_tail = w - (r_head - t);
+ *
+ *  e0  The exponent of x[0]. Must be <= 16360 or you need to
+ *              expand the ipio2 table.
+ *
+ *  nx  dimension of x[]
+ *
+ *      prec    an integer indicating the precision:
+ *          0   24  bits (single)
+ *          1   53  bits (double)
+ *          2   64  bits (extended)
+ *          3   113 bits (quad)
+ *
+ * External function:
+ *  double scalbn(), floor();
+ *
+ *
+ * Here is the description of some local variables:
+ *
+ *  jk  jk+1 is the initial number of terms of ipio2[] needed
+ *      in the computation. The minimum and recommended value
+ *      for jk is 3,4,4,6 for single, double, extended, and quad.
+ *      jk+1 must be 2 larger than you might expect so that our
+ *      recomputation test works. (Up to 24 bits in the integer
+ *      part (the 24 bits of it that we compute) and 23 bits in
+ *      the fraction part may be lost to cancelation before we
+ *      recompute.)
+ *
+ *  jz  local integer variable indicating the number of
+ *      terms of ipio2[] used.
+ *
+ *  jx  nx - 1
+ *
+ *  jv  index for pointing to the suitable ipio2[] for the
+ *      computation. In general, we want
+ *          ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
+ *      is an integer. Thus
+ *          e0-3-24*jv >= 0 or (e0-3)/24 >= jv
+ *      Hence jv = max(0,(e0-3)/24).
+ *
+ *  jp  jp+1 is the number of terms in PIo2[] needed, jp = jk.
+ *
+ *  q[] double array with integral value, representing the
+ *      24-bits chunk of the product of x and 2/pi.
+ *
+ *  q0  the corresponding exponent of q[0]. Note that the
+ *      exponent for q[i] would be q0-24*i.
+ *
+ *  PIo2[]  double precision array, obtained by cutting pi/2
+ *      into 24 bits chunks.
+ *
+ *  f[] ipio2[] in floating point
+ *
+ *  iq[]    integer array by breaking up q[] in 24-bits chunk.
+ *
+ *  fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
+ *
+ *  ih  integer. If >0 it indicates q[] is >= 0.5, hence
+ *      it also indicates the *sign* of the result.
+ *
+ */
+
+
+/*
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const int init_jk[] = {3,4,4,6}; /* initial value for jk */
+
+/*
+ * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+ *
+ *      integer array, contains the (24*i)-th to (24*i+23)-th
+ *      bit of 2/pi after binary point. The corresponding
+ *      floating value is
+ *
+ *          ipio2[i] * 2^(-24(i+1)).
+ *
+ * NB: This table must have at least (e0-3)/24 + jk terms.
+ *     For quad precision (e0 <= 16360, jk = 6), this is 686.
+ */
+static const int32_t ipio2[] = {
+    0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
+    0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
+    0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
+    0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
+    0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
+    0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
+    0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
+    0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
+    0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
+    0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
+    0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
+
+#if LDBL_MAX_EXP > 1024
+#if LDBL_MAX_EXP > 16384
+#error "ipio2 table needs to be expanded"
+#endif
+    0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
+    0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
+    0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
+    0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
+    0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
+    0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
+    0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
+    0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
+    0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
+    0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
+    0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
+    0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
+    0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
+    0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
+    0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
+    0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
+    0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
+    0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
+    0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
+    0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
+    0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
+    0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
+    0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
+    0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
+    0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
+    0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
+    0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
+    0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
+    0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
+    0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
+    0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
+    0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
+    0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
+    0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
+    0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
+    0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
+    0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
+    0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
+    0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
+    0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
+    0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
+    0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
+    0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
+    0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
+    0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
+    0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
+    0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
+    0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
+    0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
+    0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
+    0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
+    0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
+    0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
+    0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
+    0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
+    0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
+    0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
+    0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
+    0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
+    0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
+    0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
+    0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
+    0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
+    0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
+    0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
+    0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
+    0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
+    0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
+    0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
+    0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
+    0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
+    0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
+    0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
+    0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
+    0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
+    0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
+    0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
+    0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
+    0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
+    0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
+    0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
+    0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
+    0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
+    0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
+    0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
+    0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
+    0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
+    0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
+    0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
+    0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
+    0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
+    0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
+    0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
+    0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
+    0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
+    0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
+    0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
+    0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
+    0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
+    0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
+    0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
+    0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
+    0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
+    0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
+#endif
+
+};
+
+static const double PIo2[] = {
+    1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
+    7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
+    5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
+    3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
+    1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
+    1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
+    2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
+    2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
+};
+
+static const double
+zero   = 0.0,
+one    = 1.0,
+two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
+
+int
+__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
+{
+    int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
+    double z,fw,f[20],fq[20],q[20];
+
+    /* initialize jk*/
+    jk = init_jk[prec];
+    jp = jk;
+
+    /* determine jx,jv,q0, note that 3>q0 */
+    jx =  nx-1;
+    jv = (e0-3)/24;
+    if (jv<0) jv=0;
+    q0 =  e0-24*(jv+1);
+
+    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+    j = jv-jx;
+    m = jx+jk;
+    for (i=0; i<=m; i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
+
+    /* compute q[0],q[1],...q[jk] */
+    for (i=0; i<=jk; i++) {
+        for (j=0,fw=0.0; j<=jx; j++) fw += x[j]*f[jx+i-j];
+        q[i] = fw;
+    }
+
+    jz = jk;
+recompute:
+    /* distill q[] into iq[] reversingly */
+    for (i=0,j=jz,z=q[jz]; j>0; i++,j--) {
+        fw    =  (double)((int32_t)(twon24* z));
+        iq[i] =  (int32_t)(z-two24*fw);
+        z     =  q[j-1]+fw;
+    }
+
+    /* compute n */
+    z  = scalbn(z,q0);      /* actual value of z */
+    z -= 8.0*floor(z*0.125);        /* trim off integer >= 8 */
+    n  = (int32_t) z;
+    z -= (double)n;
+    ih = 0;
+    if (q0>0) { /* need iq[jz-1] to determine n */
+        i  = (iq[jz-1]>>(24-q0));
+        n += i;
+        iq[jz-1] -= i<<(24-q0);
+        ih = iq[jz-1]>>(23-q0);
+    } else if (q0==0) ih = iq[jz-1]>>23;
+    else if (z>=0.5) ih=2;
+
+    if (ih>0) { /* q > 0.5 */
+        n += 1;
+        carry = 0;
+        for (i=0; i<jz ; i++) { /* compute 1-q */
+            j = iq[i];
+            if (carry==0) {
+                if (j!=0) {
+                    carry = 1;
+                    iq[i] = 0x1000000- j;
+                }
+            } else  iq[i] = 0xffffff - j;
+        }
+        if (q0>0) {     /* rare case: chance is 1 in 12 */
+            switch (q0) {
+                case 1:
+                    iq[jz-1] &= 0x7fffff;
+                    break;
+                case 2:
+                    iq[jz-1] &= 0x3fffff;
+                    break;
+            }
+        }
+        if (ih==2) {
+            z = one - z;
+            if (carry!=0) z -= scalbn(one,q0);
+        }
+    }
+
+    /* check if recomputation is needed */
+    if (z==zero) {
+        j = 0;
+        for (i=jz-1; i>=jk; i--) j |= iq[i];
+        if (j==0) { /* need recomputation */
+            for (k=1; iq[jk-k]==0; k++); /* k = no. of terms needed */
+
+            for (i=jz+1; i<=jz+k; i++) { /* add q[jz+1] to q[jz+k] */
+                f[jx+i] = (double) ipio2[jv+i];
+                for (j=0,fw=0.0; j<=jx; j++) fw += x[j]*f[jx+i-j];
+                q[i] = fw;
+            }
+            jz += k;
+            goto recompute;
+        }
+    }
+
+    /* chop off zero terms */
+    if (z==0.0) {
+        jz -= 1;
+        q0 -= 24;
+        while (iq[jz]==0) { jz--; q0-=24;}
+    } else { /* break z into 24-bit if necessary */
+        z = scalbn(z,-q0);
+        if (z>=two24) {
+            fw = (double)((int32_t)(twon24*z));
+            iq[jz] = (int32_t)(z-two24*fw);
+            jz += 1;
+            q0 += 24;
+            iq[jz] = (int32_t) fw;
+        } else iq[jz] = (int32_t) z ;
+    }
+
+    /* convert integer "bit" chunk to floating-point value */
+    fw = scalbn(one,q0);
+    for (i=jz; i>=0; i--) {
+        q[i] = fw*(double)iq[i];
+        fw*=twon24;
+    }
+
+    /* compute PIo2[0,...,jp]*q[jz,...,0] */
+    for (i=jz; i>=0; i--) {
+        for (fw=0.0,k=0; k<=jp&&k<=jz-i; k++) fw += PIo2[k]*q[i+k];
+        fq[jz-i] = fw;
+    }
+
+    /* compress fq[] into y[] */
+    switch (prec) {
+        case 0:
+            fw = 0.0;
+            for (i=jz; i>=0; i--) fw += fq[i];
+            y[0] = (ih==0)? fw: -fw;
+            break;
+        case 1:
+        case 2:
+            fw = 0.0;
+            for (i=jz; i>=0; i--) fw += fq[i];
+            STRICT_ASSIGN(double,fw,fw);
+            y[0] = (ih==0)? fw: -fw;
+            fw = fq[0]-fw;
+            for (i=1; i<=jz; i++) fw += fq[i];
+            y[1] = (ih==0)? fw: -fw;
+            break;
+        case 3: /* painful */
+            for (i=jz; i>0; i--) {
+                fw      = fq[i-1]+fq[i];
+                fq[i]  += fq[i-1]-fw;
+                fq[i-1] = fw;
+            }
+            for (i=jz; i>1; i--) {
+                fw      = fq[i-1]+fq[i];
+                fq[i]  += fq[i-1]-fw;
+                fq[i-1] = fw;
+            }
+            for (fw=0.0,i=jz; i>=2; i--) fw += fq[i];
+            if (ih==0) {
+                y[0] =  fq[0];
+                y[1] =  fq[1];
+                y[2] =  fw;
+            } else {
+                y[0] = -fq[0];
+                y[1] = -fq[1];
+                y[2] = -fw;
+            }
+    }
+    return n&7;
+}
diff --git a/src/bsp/lk/lib/libm/k_sin.c b/src/bsp/lk/lib/libm/k_sin.c
new file mode 100644
index 0000000..8922a04
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_sin.c
@@ -0,0 +1,70 @@
+
+/* @(#)k_sin.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __kernel_sin( x, y, iy)
+ * kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
+ *
+ * Algorithm
+ *  1. Since sin(-x) = -sin(x), we need only to consider positive x.
+ *  2. Callers must return sin(-0) = -0 without calling here since our
+ *     odd polynomial is not evaluated in a way that preserves -0.
+ *     Callers may do the optimization sin(x) ~ x for tiny x.
+ *  3. sin(x) is approximated by a polynomial of degree 13 on
+ *     [0,pi/4]
+ *                   3            13
+ *      sin(x) ~ x + S1*x + ... + S6*x
+ *     where
+ *
+ *  |sin(x)         2     4     6     8     10     12  |     -58
+ *  |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
+ *  |  x                               |
+ *
+ *  4. sin(x+y) = sin(x) + sin'(x')*y
+ *          ~ sin(x) + (1-x*x/2)*y
+ *     For better accuracy, let
+ *           3      2      2      2      2
+ *      r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
+ *     then                   3    2
+ *      sin(x) = x + (S1*x + (x *(r-y/2)+y))
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
+S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
+S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
+S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
+S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
+S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
+S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+double
+__kernel_sin(double x, double y, int iy)
+{
+    double z,r,v,w;
+
+    z   =  x*x;
+    w   =  z*z;
+    r   =  S2+z*(S3+z*S4) + z*w*(S5+z*S6);
+    v   =  z*x;
+    if (iy==0) return x+v*(S1+z*r);
+    else      return x-((z*(half*y-v*r)-y)-v*S1);
+}
diff --git a/src/bsp/lk/lib/libm/k_sinf.c b/src/bsp/lk/lib/libm/k_sinf.c
new file mode 100644
index 0000000..efc3189
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_sinf.c
@@ -0,0 +1,46 @@
+/* k_sinf.c -- float version of k_sin.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef INLINE_KERNEL_SINDF
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
+static const double
+S1 = -0x15555554cbac77.0p-55,   /* -0.166666666416265235595 */
+S2 =  0x111110896efbb2.0p-59,   /*  0.0083333293858894631756 */
+S3 = -0x1a00f9e2cae774.0p-65,   /* -0.000198393348360966317347 */
+S4 =  0x16cd878c3b46a7.0p-71;   /*  0.0000027183114939898219064 */
+
+#ifdef INLINE_KERNEL_SINDF
+static __inline
+#endif
+float
+__kernel_sindf(double x)
+{
+    double r, s, w, z;
+
+    /* Try to optimize for parallel evaluation as in k_tanf.c. */
+    z = x*x;
+    w = z*z;
+    r = S3+z*S4;
+    s = z*x;
+    return (x + s*(S1+z*S2)) + s*w*r;
+}
diff --git a/src/bsp/lk/lib/libm/k_tan.c b/src/bsp/lk/lib/libm/k_tan.c
new file mode 100644
index 0000000..edb0d4b
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_tan.c
@@ -0,0 +1,133 @@
+/* @(#)k_tan.c 1.5 04/04/22 SMI */
+
+/*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* INDENT OFF */
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* __kernel_tan( x, y, k )
+ * kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ * Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
+ *
+ * Algorithm
+ *  1. Since tan(-x) = -tan(x), we need only to consider positive x.
+ *  2. Callers must return tan(-0) = -0 without calling here since our
+ *     odd polynomial is not evaluated in a way that preserves -0.
+ *     Callers may do the optimization tan(x) ~ x for tiny x.
+ *  3. tan(x) is approximated by a odd polynomial of degree 27 on
+ *     [0,0.67434]
+ *                   3             27
+ *      tan(x) ~ x + T1*x + ... + T13*x
+ *     where
+ *
+ *          |tan(x)         2     4            26   |     -59.2
+ *          |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
+ *          |  x                    |
+ *
+ *     Note: tan(x+y) = tan(x) + tan'(x)*y
+ *                ~ tan(x) + (1+x*x)*y
+ *     Therefore, for better accuracy in computing tan(x+y), let
+ *           3      2      2       2       2
+ *      r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
+ *     then
+ *                  3    2
+ *      tan(x+y) = x + (T1*x + (x *(r+y)+y))
+ *
+ *      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
+ *      tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
+ *             = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
+ */
+
+#include "math.h"
+#include "math_private.h"
+static const double xxx[] = {
+    3.33333333333334091986e-01,    /* 3FD55555, 55555563 */
+    1.33333333333201242699e-01,    /* 3FC11111, 1110FE7A */
+    5.39682539762260521377e-02,    /* 3FABA1BA, 1BB341FE */
+    2.18694882948595424599e-02,    /* 3F9664F4, 8406D637 */
+    8.86323982359930005737e-03,    /* 3F8226E3, E96E8493 */
+    3.59207910759131235356e-03,    /* 3F6D6D22, C9560328 */
+    1.45620945432529025516e-03,    /* 3F57DBC8, FEE08315 */
+    5.88041240820264096874e-04,    /* 3F4344D8, F2F26501 */
+    2.46463134818469906812e-04,    /* 3F3026F7, 1A8D1068 */
+    7.81794442939557092300e-05,    /* 3F147E88, A03792A6 */
+    7.14072491382608190305e-05,    /* 3F12B80F, 32F0A7E9 */
+    -1.85586374855275456654e-05,    /* BEF375CB, DB605373 */
+    2.59073051863633712884e-05,    /* 3EFB2A70, 74BF7AD4 */
+    /* one */    1.00000000000000000000e+00,    /* 3FF00000, 00000000 */
+    /* pio4 */   7.85398163397448278999e-01,    /* 3FE921FB, 54442D18 */
+    /* pio4lo */     3.06161699786838301793e-17 /* 3C81A626, 33145C07 */
+};
+#define one xxx[13]
+#define pio4    xxx[14]
+#define pio4lo  xxx[15]
+#define T   xxx
+/* INDENT ON */
+
+double
+__kernel_tan(double x, double y, int iy)
+{
+    double z, r, v, w, s;
+    int32_t ix, hx;
+
+    GET_HIGH_WORD(hx,x);
+    ix = hx & 0x7fffffff;           /* high word of |x| */
+    if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */
+        if (hx < 0) {
+            x = -x;
+            y = -y;
+        }
+        z = pio4 - x;
+        w = pio4lo - y;
+        x = z + w;
+        y = 0.0;
+    }
+    z = x * x;
+    w = z * z;
+    /*
+     * Break x^5*(T[1]+x^2*T[2]+...) into
+     * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
+     * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
+     */
+    r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] +
+                                            w * T[11]))));
+    v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] +
+                                     w * T[12])))));
+    s = z * x;
+    r = y + z * (s * (r + v) + y);
+    r += T[0] * s;
+    w = x + r;
+    if (ix >= 0x3FE59428) {
+        v = (double) iy;
+        return (double) (1 - ((hx >> 30) & 2)) *
+               (v - 2.0 * (x - (w * w / (w + v) - r)));
+    }
+    if (iy == 1)
+        return w;
+    else {
+        /*
+         * if allow error up to 2 ulp, simply return
+         * -1.0 / (x+r) here
+         */
+        /* compute -1.0 / (x+r) accurately */
+        double a, t;
+        z = w;
+        SET_LOW_WORD(z,0);
+        v = r - (z - x);    /* z+v = r+x */
+        t = a = -1.0 / w;   /* a = -1.0/w */
+        SET_LOW_WORD(t,0);
+        s = 1.0 + t * z;
+        return t + a * (s + t * v);
+    }
+}
diff --git a/src/bsp/lk/lib/libm/k_tanf.c b/src/bsp/lk/lib/libm/k_tanf.c
new file mode 100644
index 0000000..887bd4c
--- /dev/null
+++ b/src/bsp/lk/lib/libm/k_tanf.c
@@ -0,0 +1,66 @@
+/* k_tanf.c -- float version of k_tan.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef INLINE_KERNEL_TANDF
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
+static const double
+T[] =  {
+    0x15554d3418c99f.0p-54,   /* 0.333331395030791399758 */
+    0x1112fd38999f72.0p-55,   /* 0.133392002712976742718 */
+    0x1b54c91d865afe.0p-57,   /* 0.0533812378445670393523 */
+    0x191df3908c33ce.0p-58,   /* 0.0245283181166547278873 */
+    0x185dadfcecf44e.0p-61,   /* 0.00297435743359967304927 */
+    0x1362b9bf971bcd.0p-59,   /* 0.00946564784943673166728 */
+};
+
+#ifdef INLINE_KERNEL_TANDF
+static __inline
+#endif
+float
+__kernel_tandf(double x, int iy)
+{
+    double z,r,w,s,t,u;
+
+    z   =  x*x;
+    /*
+     * Split up the polynomial into small independent terms to give
+     * opportunities for parallel evaluation.  The chosen splitting is
+     * micro-optimized for Athlons (XP, X64).  It costs 2 multiplications
+     * relative to Horner's method on sequential machines.
+     *
+     * We add the small terms from lowest degree up for efficiency on
+     * non-sequential machines (the lowest degree terms tend to be ready
+     * earlier).  Apart from this, we don't care about order of
+     * operations, and don't need to to care since we have precision to
+     * spare.  However, the chosen splitting is good for accuracy too,
+     * and would give results as accurate as Horner's method if the
+     * small terms were added from highest degree down.
+     */
+    r = T[4]+z*T[5];
+    t = T[2]+z*T[3];
+    w = z*z;
+    s = z*x;
+    u = T[0]+z*T[1];
+    r = (x+s*u)+(s*w)*(t+w*r);
+    if (iy==1) return r;
+    else return -1.0/r;
+}
diff --git a/src/bsp/lk/lib/libm/math_private.h b/src/bsp/lk/lib/libm/math_private.h
new file mode 100644
index 0000000..86b1f85
--- /dev/null
+++ b/src/bsp/lk/lib/libm/math_private.h
@@ -0,0 +1,757 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ * $FreeBSD$
+ */
+
+#ifndef _MATH_PRIVATE_H_
+#define _MATH_PRIVATE_H_
+
+#include <sys/types.h>
+#include <endian.h>
+
+/*
+ * The original fdlibm code used statements like:
+ *  n0 = ((*(int*)&one)>>29)^1;     * index of high word *
+ *  ix0 = *(n0+(int*)&x);           * high word of x *
+ *  ix1 = *((1-n0)+(int*)&x);       * low word of x *
+ * to dig two 32 bit words out of the 64 bit IEEE floating point
+ * value.  That is non-ANSI, and, moreover, the gcc instruction
+ * scheduler gets it wrong.  We instead use the following macros.
+ * Unlike the original code, we determine the endianness at compile
+ * time, not at run time; I don't see much benefit to selecting
+ * endianness at run time.
+ */
+
+/*
+ * A union which permits us to convert between a double and two 32 bit
+ * ints.
+ */
+
+#ifdef __arm__
+#if defined(__VFP_FP__) || defined(__ARM_EABI__)
+#define IEEE_WORD_ORDER BYTE_ORDER
+#else
+#define IEEE_WORD_ORDER BIG_ENDIAN
+#endif
+#else /* __arm__ */
+#define IEEE_WORD_ORDER BYTE_ORDER
+#endif
+
+#if IEEE_WORD_ORDER == BIG_ENDIAN
+
+typedef union {
+    double value;
+    struct {
+        u_int32_t msw;
+        u_int32_t lsw;
+    } parts;
+    struct {
+        u_int64_t w;
+    } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+#if IEEE_WORD_ORDER == LITTLE_ENDIAN
+
+typedef union {
+    double value;
+    struct {
+        u_int32_t lsw;
+        u_int32_t msw;
+    } parts;
+    struct {
+        u_int64_t w;
+    } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0,ix1,d)                \
+do {                                \
+  ieee_double_shape_type ew_u;                  \
+  ew_u.value = (d);                     \
+  (ix0) = ew_u.parts.msw;                   \
+  (ix1) = ew_u.parts.lsw;                   \
+} while (0)
+
+/* Get a 64-bit int from a double. */
+#define EXTRACT_WORD64(ix,d)                    \
+do {                                \
+  ieee_double_shape_type ew_u;                  \
+  ew_u.value = (d);                     \
+  (ix) = ew_u.xparts.w;                     \
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i,d)                  \
+do {                                \
+  ieee_double_shape_type gh_u;                  \
+  gh_u.value = (d);                     \
+  (i) = gh_u.parts.msw;                     \
+} while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+
+#define GET_LOW_WORD(i,d)                   \
+do {                                \
+  ieee_double_shape_type gl_u;                  \
+  gl_u.value = (d);                     \
+  (i) = gl_u.parts.lsw;                     \
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+
+#define INSERT_WORDS(d,ix0,ix1)                 \
+do {                                \
+  ieee_double_shape_type iw_u;                  \
+  iw_u.parts.msw = (ix0);                   \
+  iw_u.parts.lsw = (ix1);                   \
+  (d) = iw_u.value;                     \
+} while (0)
+
+/* Set a double from a 64-bit int. */
+#define INSERT_WORD64(d,ix)                 \
+do {                                \
+  ieee_double_shape_type iw_u;                  \
+  iw_u.xparts.w = (ix);                     \
+  (d) = iw_u.value;                     \
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d,v)                  \
+do {                                \
+  ieee_double_shape_type sh_u;                  \
+  sh_u.value = (d);                     \
+  sh_u.parts.msw = (v);                     \
+  (d) = sh_u.value;                     \
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+
+#define SET_LOW_WORD(d,v)                   \
+do {                                \
+  ieee_double_shape_type sl_u;                  \
+  sl_u.value = (d);                     \
+  sl_u.parts.lsw = (v);                     \
+  (d) = sl_u.value;                     \
+} while (0)
+
+/*
+ * A union which permits us to convert between a float and a 32 bit
+ * int.
+ */
+
+typedef union {
+    float value;
+    /* FIXME: Assumes 32 bit int.  */
+    unsigned int word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float.  */
+
+#define GET_FLOAT_WORD(i,d)                 \
+do {                                \
+  ieee_float_shape_type gf_u;                   \
+  gf_u.value = (d);                     \
+  (i) = gf_u.word;                      \
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+
+#define SET_FLOAT_WORD(d,i)                 \
+do {                                \
+  ieee_float_shape_type sf_u;                   \
+  sf_u.word = (i);                      \
+  (d) = sf_u.value;                     \
+} while (0)
+
+/*
+ * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
+ * double.
+ */
+
+#define EXTRACT_LDBL80_WORDS(ix0,ix1,d)             \
+do {                                \
+  union IEEEl2bits ew_u;                    \
+  ew_u.e = (d);                         \
+  (ix0) = ew_u.xbits.expsign;                   \
+  (ix1) = ew_u.xbits.man;                   \
+} while (0)
+
+/*
+ * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
+ * long double.
+ */
+
+#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d)            \
+do {                                \
+  union IEEEl2bits ew_u;                    \
+  ew_u.e = (d);                         \
+  (ix0) = ew_u.xbits.expsign;                   \
+  (ix1) = ew_u.xbits.manh;                  \
+  (ix2) = ew_u.xbits.manl;                  \
+} while (0)
+
+/* Get expsign as a 16 bit int from a long double.  */
+
+#define GET_LDBL_EXPSIGN(i,d)                   \
+do {                                \
+  union IEEEl2bits ge_u;                    \
+  ge_u.e = (d);                         \
+  (i) = ge_u.xbits.expsign;                 \
+} while (0)
+
+/*
+ * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
+ * mantissa.
+ */
+
+#define INSERT_LDBL80_WORDS(d,ix0,ix1)              \
+do {                                \
+  union IEEEl2bits iw_u;                    \
+  iw_u.xbits.expsign = (ix0);                   \
+  iw_u.xbits.man = (ix1);                   \
+  (d) = iw_u.e;                         \
+} while (0)
+
+/*
+ * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
+ * comprising the mantissa.
+ */
+
+#define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2)         \
+do {                                \
+  union IEEEl2bits iw_u;                    \
+  iw_u.xbits.expsign = (ix0);                   \
+  iw_u.xbits.manh = (ix1);                  \
+  iw_u.xbits.manl = (ix2);                  \
+  (d) = iw_u.e;                         \
+} while (0)
+
+/* Set expsign of a long double from a 16 bit int.  */
+
+#define SET_LDBL_EXPSIGN(d,v)                   \
+do {                                \
+  union IEEEl2bits se_u;                    \
+  se_u.e = (d);                         \
+  se_u.xbits.expsign = (v);                 \
+  (d) = se_u.e;                         \
+} while (0)
+
+#ifdef __i386__
+/* Long double constants are broken on i386. */
+#define LD80C(m, ex, v) {                       \
+    .xbits.man = __CONCAT(m, ULL),                  \
+    .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0),  \
+}
+#else
+/* The above works on non-i386 too, but we use this to check v. */
+#define LD80C(m, ex, v) { .e = (v), }
+#endif
+
+#ifdef FLT_EVAL_METHOD
+/*
+ * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
+ */
+#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
+#else
+#define STRICT_ASSIGN(type, lval, rval) do {    \
+    volatile type __lval;           \
+                        \
+    if (sizeof(type) >= sizeof(long double))    \
+        (lval) = (rval);        \
+    else {                  \
+        __lval = (rval);        \
+        (lval) = __lval;        \
+    }                   \
+} while (0)
+#endif
+#endif /* FLT_EVAL_METHOD */
+
+/* Support switching the mode to FP_PE if necessary. */
+#if defined(__i386__) && !defined(NO_FPSETPREC)
+#define ENTERI()                \
+    long double __retval;           \
+    fp_prec_t __oprec;          \
+                        \
+    if ((__oprec = fpgetprec()) != FP_PE)   \
+        fpsetprec(FP_PE)
+#define RETURNI(x) do {             \
+    __retval = (x);             \
+    if (__oprec != FP_PE)           \
+        fpsetprec(__oprec);     \
+    RETURNF(__retval);          \
+} while (0)
+#else
+#define ENTERI(x)
+#define RETURNI(x)  RETURNF(x)
+#endif
+
+/* Default return statement if hack*_t() is not used. */
+#define      RETURNF(v)      return (v)
+
+/*
+ * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
+ * a == 0, but is slower.
+ */
+#define _2sum(a, b) do {    \
+    __typeof(a) __s, __w;   \
+                \
+    __w = (a) + (b);    \
+    __s = __w - (a);    \
+    (b) = ((a) - (__w - __s)) + ((b) - __s); \
+    (a) = __w;      \
+} while (0)
+
+/*
+ * 2sumF algorithm.
+ *
+ * "Normalize" the terms in the infinite-precision expression a + b for
+ * the sum of 2 floating point values so that b is as small as possible
+ * relative to 'a'.  (The resulting 'a' is the value of the expression in
+ * the same precision as 'a' and the resulting b is the rounding error.)
+ * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
+ * exponent overflow or underflow must not occur.  This uses a Theorem of
+ * Dekker (1971).  See Knuth (1981) 4.2.2 Theorem C.  The name "TwoSum"
+ * is apparently due to Skewchuk (1997).
+ *
+ * For this to always work, assignment of a + b to 'a' must not retain any
+ * extra precision in a + b.  This is required by C standards but broken
+ * in many compilers.  The brokenness cannot be worked around using
+ * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
+ * algorithm would be destroyed by non-null strict assignments.  (The
+ * compilers are correct to be broken -- the efficiency of all floating
+ * point code calculations would be destroyed similarly if they forced the
+ * conversions.)
+ *
+ * Fortunately, a case that works well can usually be arranged by building
+ * any extra precision into the type of 'a' -- 'a' should have type float_t,
+ * double_t or long double.  b's type should be no larger than 'a's type.
+ * Callers should use these types with scopes as large as possible, to
+ * reduce their own extra-precision and efficiciency problems.  In
+ * particular, they shouldn't convert back and forth just to call here.
+ */
+#ifdef DEBUG
+#define _2sumF(a, b) do {               \
+    __typeof(a) __w;                \
+    volatile __typeof(a) __ia, __ib, __r, __vw; \
+                            \
+    __ia = (a);                 \
+    __ib = (b);                 \
+    assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib));    \
+                            \
+    __w = (a) + (b);                \
+    (b) = ((a) - __w) + (b);            \
+    (a) = __w;                  \
+                            \
+    /* The next 2 assertions are weak if (a) is already long double. */ \
+    assert((long double)__ia + __ib == (long double)(a) + (b)); \
+    __vw = __ia + __ib;             \
+    __r = __ia - __vw;              \
+    __r += __ib;                    \
+    assert(__vw == (a) && __r == (b));      \
+} while (0)
+#else /* !DEBUG */
+#define _2sumF(a, b) do {   \
+    __typeof(a) __w;    \
+                \
+    __w = (a) + (b);    \
+    (b) = ((a) - __w) + (b); \
+    (a) = __w;      \
+} while (0)
+#endif /* DEBUG */
+
+/*
+ * Set x += c, where x is represented in extra precision as a + b.
+ * x must be sufficiently normalized and sufficiently larger than c,
+ * and the result is then sufficiently normalized.
+ *
+ * The details of ordering are that |a| must be >= |c| (so that (a, c)
+ * can be normalized without extra work to swap 'a' with c).  The details of
+ * the normalization are that b must be small relative to the normalized 'a'.
+ * Normalization of (a, c) makes the normalized c tiny relative to the
+ * normalized a, so b remains small relative to 'a' in the result.  However,
+ * b need not ever be tiny relative to 'a'.  For example, b might be about
+ * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
+ * That is usually enough, and adding c (which by normalization is about
+ * 2**53 times smaller than a) cannot change b significantly.  However,
+ * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
+ * significantly relative to b.  The caller must ensure that significant
+ * cancellation doesn't occur, either by having c of the same sign as 'a',
+ * or by having |c| a few percent smaller than |a|.  Pre-normalization of
+ * (a, b) may help.
+ *
+ * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
+ * exercise 19).  We gain considerable efficiency by requiring the terms to
+ * be sufficiently normalized and sufficiently increasing.
+ */
+#define _3sumF(a, b, c) do {    \
+    __typeof(a) __tmp;  \
+                \
+    __tmp = (c);        \
+    _2sumF(__tmp, (a)); \
+    (b) += (a);     \
+    (a) = __tmp;        \
+} while (0)
+
+/*
+ * Common routine to process the arguments to nan(), nanf(), and nanl().
+ */
+void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
+
+#ifdef _COMPLEX_H
+
+/*
+ * C99 specifies that complex numbers have the same representation as
+ * an array of two elements, where the first element is the real part
+ * and the second element is the imaginary part.
+ */
+typedef union {
+    float complex f;
+    float a[2];
+} float_complex;
+typedef union {
+    double complex f;
+    double a[2];
+} double_complex;
+typedef union {
+    long double complex f;
+    long double a[2];
+} long_double_complex;
+#define REALPART(z) ((z).a[0])
+#define IMAGPART(z) ((z).a[1])
+
+/*
+ * Inline functions that can be used to construct complex values.
+ *
+ * The C99 standard intends x+I*y to be used for this, but x+I*y is
+ * currently unusable in general since gcc introduces many overflow,
+ * underflow, sign and efficiency bugs by rewriting I*y as
+ * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
+ * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
+ * to -0.0+I*0.0.
+ */
+static __inline float complex
+cpackf(float x, float y)
+{
+    float_complex z;
+
+    REALPART(z) = x;
+    IMAGPART(z) = y;
+    return (z.f);
+}
+
+static __inline double complex
+cpack(double x, double y)
+{
+    double_complex z;
+
+    REALPART(z) = x;
+    IMAGPART(z) = y;
+    return (z.f);
+}
+
+static __inline long double complex
+cpackl(long double x, long double y)
+{
+    long_double_complex z;
+
+    REALPART(z) = x;
+    IMAGPART(z) = y;
+    return (z.f);
+}
+#endif /* _COMPLEX_H */
+
+#ifdef __GNUCLIKE_ASM
+
+/* Asm versions of some functions. */
+
+#ifdef __amd64__
+static __inline int
+irint(double x)
+{
+    int n;
+
+    asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
+    return (n);
+}
+#define HAVE_EFFICIENT_IRINT
+#endif
+
+#ifdef __i386__
+static __inline int
+irint(double x)
+{
+    int n;
+
+    asm("fistl %0" : "=m" (n) : "t" (x));
+    return (n);
+}
+#define HAVE_EFFICIENT_IRINT
+#endif
+
+#if defined(__amd64__) || defined(__i386__)
+static __inline int
+irintl(long double x)
+{
+    int n;
+
+    asm("fistl %0" : "=m" (n) : "t" (x));
+    return (n);
+}
+#define HAVE_EFFICIENT_IRINTL
+#endif
+
+#endif /* __GNUCLIKE_ASM */
+
+#ifdef DEBUG
+#if defined(__amd64__) || defined(__i386__)
+#define breakpoint()    asm("int $3")
+#else
+#include <signal.h>
+
+#define breakpoint()    raise(SIGTRAP)
+#endif
+#endif
+
+/* Write a pari script to test things externally. */
+#ifdef DOPRINT
+#include <stdio.h>
+
+#ifndef DOPRINT_SWIZZLE
+#define DOPRINT_SWIZZLE     0
+#endif
+
+#ifdef DOPRINT_LD80
+
+#define DOPRINT_START(xp) do {                      \
+    uint64_t __lx;                          \
+    uint16_t __hx;                          \
+                                    \
+    /* Hack to give more-problematic args. */           \
+    EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);              \
+    __lx ^= DOPRINT_SWIZZLE;                    \
+    INSERT_LDBL80_WORDS(*xp, __hx, __lx);               \
+    printf("x = %.21Lg; ", (long double)*xp);           \
+} while (0)
+#define DOPRINT_END1(v)                         \
+    printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo)                        \
+    printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",      \
+        (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_D64)
+
+#define DOPRINT_START(xp) do {                      \
+    uint32_t __hx, __lx;                        \
+                                    \
+    EXTRACT_WORDS(__hx, __lx, *xp);                 \
+    __lx ^= DOPRINT_SWIZZLE;                    \
+    INSERT_WORDS(*xp, __hx, __lx);                  \
+    printf("x = %.21Lg; ", (long double)*xp);           \
+} while (0)
+#define DOPRINT_END1(v)                         \
+    printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo)                        \
+    printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",      \
+        (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_F32)
+
+#define DOPRINT_START(xp) do {                      \
+    uint32_t __hx;                          \
+                                    \
+    GET_FLOAT_WORD(__hx, *xp);                  \
+    __hx ^= DOPRINT_SWIZZLE;                    \
+    SET_FLOAT_WORD(*xp, __hx);                  \
+    printf("x = %.21Lg; ", (long double)*xp);           \
+} while (0)
+#define DOPRINT_END1(v)                         \
+    printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo)                        \
+    printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",      \
+        (long double)(hi), (long double)(lo))
+
+#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
+
+#ifndef DOPRINT_SWIZZLE_HIGH
+#define DOPRINT_SWIZZLE_HIGH    0
+#endif
+
+#define DOPRINT_START(xp) do {                      \
+    uint64_t __lx, __llx;                       \
+    uint16_t __hx;                          \
+                                    \
+    EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);          \
+    __llx ^= DOPRINT_SWIZZLE;                   \
+    __lx ^= DOPRINT_SWIZZLE_HIGH;                   \
+    INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);           \
+    printf("x = %.36Lg; ", (long double)*xp);                   \
+} while (0)
+#define DOPRINT_END1(v)                         \
+    printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo)                        \
+    printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",      \
+        (long double)(hi), (long double)(lo))
+
+#endif /* DOPRINT_LD80 */
+
+#else /* !DOPRINT */
+#define DOPRINT_START(xp)
+#define DOPRINT_END1(v)
+#define DOPRINT_END2(hi, lo)
+#endif /* DOPRINT */
+
+#define RETURNP(x) do {         \
+    DOPRINT_END1(x);        \
+    RETURNF(x);         \
+} while (0)
+#define RETURNPI(x) do {        \
+    DOPRINT_END1(x);        \
+    RETURNI(x);         \
+} while (0)
+#define RETURN2P(x, y) do {     \
+    DOPRINT_END2((x), (y));     \
+    RETURNF((x) + (y));     \
+} while (0)
+#define RETURN2PI(x, y) do {        \
+    DOPRINT_END2((x), (y));     \
+    RETURNI((x) + (y));     \
+} while (0)
+#ifdef STRUCT_RETURN
+#define RETURNSP(rp) do {       \
+    if (!(rp)->lo_set)      \
+        RETURNP((rp)->hi);  \
+    RETURN2P((rp)->hi, (rp)->lo);   \
+} while (0)
+#define RETURNSPI(rp) do {      \
+    if (!(rp)->lo_set)      \
+        RETURNPI((rp)->hi); \
+    RETURN2PI((rp)->hi, (rp)->lo);  \
+} while (0)
+#endif
+#define SUM2P(x, y) ({          \
+    const __typeof (x) __x = (x);   \
+    const __typeof (y) __y = (y);   \
+                    \
+    DOPRINT_END2(__x, __y);     \
+    __x + __y;          \
+})
+
+/*
+ * ieee style elementary functions
+ *
+ * We rename functions here to improve other sources' diffability
+ * against fdlibm.
+ */
+#define __ieee754_sqrt  sqrt
+#define __ieee754_acos  acos
+#define __ieee754_acosh acosh
+#define __ieee754_log   log
+#define __ieee754_log2  log2
+#define __ieee754_atanh atanh
+#define __ieee754_asin  asin
+#define __ieee754_atan2 atan2
+#define __ieee754_exp   exp
+#define __ieee754_cosh  cosh
+#define __ieee754_fmod  fmod
+#define __ieee754_pow   pow
+#define __ieee754_lgamma lgamma
+#define __ieee754_gamma gamma
+#define __ieee754_lgamma_r lgamma_r
+#define __ieee754_gamma_r gamma_r
+#define __ieee754_log10 log10
+#define __ieee754_sinh  sinh
+#define __ieee754_hypot hypot
+#define __ieee754_j0    j0
+#define __ieee754_j1    j1
+#define __ieee754_y0    y0
+#define __ieee754_y1    y1
+#define __ieee754_jn    jn
+#define __ieee754_yn    yn
+#define __ieee754_remainder remainder
+#define __ieee754_scalb scalb
+#define __ieee754_sqrtf sqrtf
+#define __ieee754_acosf acosf
+#define __ieee754_acoshf acoshf
+#define __ieee754_logf  logf
+#define __ieee754_atanhf atanhf
+#define __ieee754_asinf asinf
+#define __ieee754_atan2f atan2f
+#define __ieee754_expf  expf
+#define __ieee754_coshf coshf
+#define __ieee754_fmodf fmodf
+#define __ieee754_powf  powf
+#define __ieee754_lgammaf lgammaf
+#define __ieee754_gammaf gammaf
+#define __ieee754_lgammaf_r lgammaf_r
+#define __ieee754_gammaf_r gammaf_r
+#define __ieee754_log10f log10f
+#define __ieee754_log2f log2f
+#define __ieee754_sinhf sinhf
+#define __ieee754_hypotf hypotf
+#define __ieee754_j0f   j0f
+#define __ieee754_j1f   j1f
+#define __ieee754_y0f   y0f
+#define __ieee754_y1f   y1f
+#define __ieee754_jnf   jnf
+#define __ieee754_ynf   ynf
+#define __ieee754_remainderf remainderf
+#define __ieee754_scalbf scalbf
+
+/* fdlibm kernel function */
+int __kernel_rem_pio2(double*,double*,int,int,int);
+
+/* double precision kernel functions */
+#ifndef INLINE_REM_PIO2
+int __ieee754_rem_pio2(double,double*);
+#endif
+double  __kernel_sin(double,double,int);
+double  __kernel_cos(double,double);
+double  __kernel_tan(double,double,int);
+double  __ldexp_exp(double,int);
+#ifdef _COMPLEX_H
+double complex __ldexp_cexp(double complex,int);
+#endif
+
+/* float precision kernel functions */
+#ifndef INLINE_REM_PIO2F
+int __ieee754_rem_pio2f(float,double*);
+#endif
+#ifndef INLINE_KERNEL_SINDF
+float   __kernel_sindf(double);
+#endif
+#ifndef INLINE_KERNEL_COSDF
+float   __kernel_cosdf(double);
+#endif
+#ifndef INLINE_KERNEL_TANDF
+float   __kernel_tandf(double,int);
+#endif
+float   __ldexp_expf(float,int);
+#ifdef _COMPLEX_H
+float complex __ldexp_cexpf(float complex,int);
+#endif
+
+/* long double precision kernel functions */
+long double __kernel_sinl(long double, long double, int);
+long double __kernel_cosl(long double, long double);
+long double __kernel_tanl(long double, long double, int);
+
+#endif /* !_MATH_PRIVATE_H_ */
diff --git a/src/bsp/lk/lib/libm/rules.mk b/src/bsp/lk/lib/libm/rules.mk
new file mode 100644
index 0000000..0239453
--- /dev/null
+++ b/src/bsp/lk/lib/libm/rules.mk
@@ -0,0 +1,42 @@
+LOCAL_DIR := $(GET_LOCAL_DIR)
+
+MODULE := $(LOCAL_DIR)
+
+MODULE_CFLAGS += -Wno-unused-variable -Wno-sign-compare -Wno-parentheses
+
+MODULE_SRCS += \
+	$(LOCAL_DIR)/k_sin.c \
+	$(LOCAL_DIR)/s_sin.c \
+	$(LOCAL_DIR)/s_sinf.c \
+	$(LOCAL_DIR)/k_cos.c \
+	$(LOCAL_DIR)/s_cos.c \
+	$(LOCAL_DIR)/s_cosf.c \
+	$(LOCAL_DIR)/k_tan.c \
+	$(LOCAL_DIR)/s_tan.c \
+	$(LOCAL_DIR)/s_tanf.c \
+	$(LOCAL_DIR)/e_sqrt.c \
+	$(LOCAL_DIR)/e_sqrtf.c \
+	$(LOCAL_DIR)/k_rem_pio2.c \
+	$(LOCAL_DIR)/s_floor.c \
+	$(LOCAL_DIR)/s_scalbn.c \
+	$(LOCAL_DIR)/s_scalbnf.c \
+	$(LOCAL_DIR)/s_copysign.c \
+	$(LOCAL_DIR)/s_copysignf.c \
+	$(LOCAL_DIR)/e_acos.c \
+	$(LOCAL_DIR)/e_acosf.c \
+	$(LOCAL_DIR)/e_asin.c \
+	$(LOCAL_DIR)/e_asinf.c \
+	$(LOCAL_DIR)/e_pow.c \
+	$(LOCAL_DIR)/e_powf.c \
+	$(LOCAL_DIR)/s_fabs.c \
+	$(LOCAL_DIR)/s_fabsf.c \
+	$(LOCAL_DIR)/e_fmod.c \
+	$(LOCAL_DIR)/e_log.c \
+	$(LOCAL_DIR)/e_exp.c \
+	$(LOCAL_DIR)/s_round.c \
+	$(LOCAL_DIR)/s_ceil.c \
+	$(LOCAL_DIR)/s_trunc.c \
+	$(LOCAL_DIR)/s_atan.c \
+	$(LOCAL_DIR)/e_atan2.c \
+
+include make/module.mk
diff --git a/src/bsp/lk/lib/libm/s_atan.c b/src/bsp/lk/lib/libm/s_atan.c
new file mode 100644
index 0000000..ed6b5c3
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_atan.c
@@ -0,0 +1,126 @@
+/* @(#)s_atan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* atan(x)
+ * Method
+ *   1. Reduce x to positive by atan(x) = -atan(-x).
+ *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
+ *      is further reduced to one of the following intervals and the
+ *      arctangent of t is evaluated by the corresponding formula:
+ *
+ *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
+ *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
+ *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
+ *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
+ *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double atanhi[] = {
+  4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
+  7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
+  9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
+  1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
+};
+
+static const double atanlo[] = {
+  2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
+  3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
+  1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
+  6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
+};
+
+static const double aT[] = {
+  3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
+ -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
+  1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
+ -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
+  9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
+ -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
+  6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
+ -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
+  4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
+ -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
+  1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
+};
+
+	static const double
+one   = 1.0,
+huge   = 1.0e300;
+
+double
+atan(double x)
+{
+	double w,s1,s2,z;
+	int32_t ix,hx,id;
+
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x44100000) {	/* if |x| >= 2^66 */
+	    u_int32_t low;
+	    GET_LOW_WORD(low,x);
+	    if(ix>0x7ff00000||
+		(ix==0x7ff00000&&(low!=0)))
+		return x+x;		/* NaN */
+	    if(hx>0) return  atanhi[3]+*(volatile double *)&atanlo[3];
+	    else     return -atanhi[3]-*(volatile double *)&atanlo[3];
+	} if (ix < 0x3fdc0000) {	/* |x| < 0.4375 */
+	    if (ix < 0x3e400000) {	/* |x| < 2^-27 */
+		if(huge+x>one) return x;	/* raise inexact */
+	    }
+	    id = -1;
+	} else {
+	x = fabs(x);
+	if (ix < 0x3ff30000) {		/* |x| < 1.1875 */
+	    if (ix < 0x3fe60000) {	/* 7/16 <=|x|<11/16 */
+		id = 0; x = (2.0*x-one)/(2.0+x);
+	    } else {			/* 11/16<=|x|< 19/16 */
+		id = 1; x  = (x-one)/(x+one);
+	    }
+	} else {
+	    if (ix < 0x40038000) {	/* |x| < 2.4375 */
+		id = 2; x  = (x-1.5)/(one+1.5*x);
+	    } else {			/* 2.4375 <= |x| < 2^66 */
+		id = 3; x  = -1.0/x;
+	    }
+	}}
+    /* end of argument reduction */
+	z = x*x;
+	w = z*z;
+    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+	s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
+	s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
+	if (id<0) return x - x*(s1+s2);
+	else {
+	    z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
+	    return (hx<0)? -z:z;
+	}
+}
+
+#if SUPPORT_LONG_DOUBLE
+#if LDBL_MANT_DIG == 53
+__weak_reference(atan, atanl);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/s_ceil.c b/src/bsp/lk/lib/libm/s_ceil.c
new file mode 100644
index 0000000..c8933c6
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_ceil.c
@@ -0,0 +1,80 @@
+/* @(#)s_ceil.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * ceil(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to ceil(x).
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double huge = 1.0e300;
+
+double
+ceil(double x)
+{
+	int32_t i0,i1,j0;
+	u_int32_t i,j;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;
+	if(j0<20) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0<0) {i0=0x80000000;i1=0;}
+		    else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
+		}
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(i0>0) i0 += (0x00100000)>>j0;
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (j0>51) {
+	    if(j0==0x400) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0) { 		/* raise inexact flag */
+		if(i0>0) {
+		    if(j0==20) i0+=1;
+		    else {
+			j = i1 + (1<<(52-j0));
+			if(j<i1) i0+=1;	/* got a carry */
+			i1 = j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	INSERT_WORDS(x,i0,i1);
+	return x;
+}
+
+
+#if SUPPORT_LONG_DOUBLE
+#if LDBL_MANT_DIG == 53
+__weak_reference(ceil, ceill);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/s_copysign.c b/src/bsp/lk/lib/libm/s_copysign.c
new file mode 100644
index 0000000..ac28b5c
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_copysign.c
@@ -0,0 +1,33 @@
+/* @(#)s_copysign.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * copysign(double x, double y)
+ * copysign(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+double
+copysign(double x, double y)
+{
+    u_int32_t hx,hy;
+    GET_HIGH_WORD(hx,x);
+    GET_HIGH_WORD(hy,y);
+    SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
+    return x;
+}
diff --git a/src/bsp/lk/lib/libm/s_copysignf.c b/src/bsp/lk/lib/libm/s_copysignf.c
new file mode 100644
index 0000000..346257b
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_copysignf.c
@@ -0,0 +1,35 @@
+/* s_copysignf.c -- float version of s_copysign.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+
+/*
+ * copysignf(float x, float y)
+ * copysignf(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+float
+copysignf(float x, float y)
+{
+        u_int32_t ix,iy;
+        GET_FLOAT_WORD(ix,x);
+        GET_FLOAT_WORD(iy,y);
+        SET_FLOAT_WORD(x,(ix&0x7fffffff)|(iy&0x80000000));
+        return x;
+}
diff --git a/src/bsp/lk/lib/libm/s_cos.c b/src/bsp/lk/lib/libm/s_cos.c
new file mode 100644
index 0000000..b71438b
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_cos.c
@@ -0,0 +1,93 @@
+/* @(#)s_cos.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* cos(x)
+ * Return cosine function of x.
+ *
+ * kernel function:
+ *  __kernel_sin        ... sine function on [-pi/4,pi/4]
+ *  __kernel_cos        ... cosine function on [-pi/4,pi/4]
+ *  __ieee754_rem_pio2  ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *  [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *  in [-pi/4 , +pi/4], and let n = k mod 4.
+ *  We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *      0          S       C         T
+ *      1          C      -S        -1/T
+ *      2         -S      -C         T
+ *      3         -C       S        -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *  TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_REM_PIO2
+#include "math_private.h"
+#include "e_rem_pio2.c"
+
+double
+cos(double x)
+{
+    double y[2],z=0.0;
+    int32_t n, ix;
+
+    /* High word of x. */
+    GET_HIGH_WORD(ix,x);
+
+    /* |x| ~< pi/4 */
+    ix &= 0x7fffffff;
+    if (ix <= 0x3fe921fb) {
+        if (ix<0x3e46a09e)          /* if x < 2**-27 * sqrt(2) */
+            if (((int)x)==0) return 1.0; /* generate inexact */
+        return __kernel_cos(x,z);
+    }
+
+    /* cos(Inf or NaN) is NaN */
+    else if (ix>=0x7ff00000) return x-x;
+
+    /* argument reduction needed */
+    else {
+        n = __ieee754_rem_pio2(x,y);
+        switch (n&3) {
+            case 0:
+                return  __kernel_cos(y[0],y[1]);
+            case 1:
+                return -__kernel_sin(y[0],y[1],1);
+            case 2:
+                return -__kernel_cos(y[0],y[1]);
+            default:
+                return  __kernel_sin(y[0],y[1],1);
+        }
+    }
+}
+#if SUPPORT_LONG_DOUBLE
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cos, cosl);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/s_cosf.c b/src/bsp/lk/lib/libm/s_cosf.c
new file mode 100644
index 0000000..a74ebdc
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_cosf.c
@@ -0,0 +1,90 @@
+/* s_cosf.c -- float version of s_cos.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_KERNEL_COSDF
+#define INLINE_KERNEL_SINDF
+#define INLINE_REM_PIO2F
+#include "math_private.h"
+#include "e_rem_pio2f.c"
+#include "k_cosf.c"
+#include "k_sinf.c"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+c1pio2 = 1*M_PI_2,          /* 0x3FF921FB, 0x54442D18 */
+c2pio2 = 2*M_PI_2,          /* 0x400921FB, 0x54442D18 */
+c3pio2 = 3*M_PI_2,          /* 0x4012D97C, 0x7F3321D2 */
+c4pio2 = 4*M_PI_2;          /* 0x401921FB, 0x54442D18 */
+
+float
+cosf(float x)
+{
+    double y;
+    int32_t n, hx, ix;
+
+    GET_FLOAT_WORD(hx,x);
+    ix = hx & 0x7fffffff;
+
+    if (ix <= 0x3f490fda) {     /* |x| ~<= pi/4 */
+        if (ix<0x39800000)      /* |x| < 2**-12 */
+            if (((int)x)==0) return 1.0; /* 1 with inexact if x != 0 */
+        return __kernel_cosdf(x);
+    }
+    if (ix<=0x407b53d1) {       /* |x| ~<= 5*pi/4 */
+        if (ix>0x4016cbe3)      /* |x|  ~> 3*pi/4 */
+            return -__kernel_cosdf(x + (hx > 0 ? -c2pio2 : c2pio2));
+        else {
+            if (hx>0)
+                return __kernel_sindf(c1pio2 - x);
+            else
+                return __kernel_sindf(x + c1pio2);
+        }
+    }
+    if (ix<=0x40e231d5) {       /* |x| ~<= 9*pi/4 */
+        if (ix>0x40afeddf)      /* |x|  ~> 7*pi/4 */
+            return __kernel_cosdf(x + (hx > 0 ? -c4pio2 : c4pio2));
+        else {
+            if (hx>0)
+                return __kernel_sindf(x - c3pio2);
+            else
+                return __kernel_sindf(-c3pio2 - x);
+        }
+    }
+
+    /* cos(Inf or NaN) is NaN */
+    else if (ix>=0x7f800000) return x-x;
+
+    /* general argument reduction needed */
+    else {
+        n = __ieee754_rem_pio2f(x,&y);
+        switch (n&3) {
+            case 0:
+                return  __kernel_cosdf(y);
+            case 1:
+                return  __kernel_sindf(-y);
+            case 2:
+                return -__kernel_cosdf(y);
+            default:
+                return  __kernel_sindf(y);
+        }
+    }
+}
diff --git a/src/bsp/lk/lib/libm/s_fabs.c b/src/bsp/lk/lib/libm/s_fabs.c
new file mode 100644
index 0000000..0dfa940
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_fabs.c
@@ -0,0 +1,31 @@
+/* @(#)s_fabs.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/s_fabs.c,v 1.7 2002/05/28 18:15:04 alfred Exp $";
+#endif
+
+/*
+ * fabs(x) returns the absolute value of x.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+double
+fabs(double x)
+{
+	u_int32_t high;
+	GET_HIGH_WORD(high,x);
+	SET_HIGH_WORD(x,high&0x7fffffff);
+        return x;
+}
diff --git a/src/bsp/lk/lib/libm/s_fabsf.c b/src/bsp/lk/lib/libm/s_fabsf.c
new file mode 100644
index 0000000..2200705
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_fabsf.c
@@ -0,0 +1,34 @@
+/* s_fabsf.c -- float version of s_fabs.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD: src/lib/msun/src/s_fabsf.c,v 1.7 2002/05/28 18:15:04 alfred Exp $";
+#endif
+
+/*
+ * fabsf(x) returns the absolute value of x.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+float
+fabsf(float x)
+{
+	u_int32_t ix;
+	GET_FLOAT_WORD(ix,x);
+	SET_FLOAT_WORD(x,ix&0x7fffffff);
+        return x;
+}
diff --git a/src/bsp/lk/lib/libm/s_floor.c b/src/bsp/lk/lib/libm/s_floor.c
new file mode 100644
index 0000000..e14a5cb
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_floor.c
@@ -0,0 +1,83 @@
+/* @(#)s_floor.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * floor(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *  Bit twiddling.
+ * Exception:
+ *  Inexact flag raised if x not equal to floor(x).
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double huge = 1.0e300;
+
+double
+floor(double x)
+{
+    int32_t i0,i1,j0;
+    u_int32_t i,j;
+    EXTRACT_WORDS(i0,i1,x);
+    j0 = ((i0>>20)&0x7ff)-0x3ff;
+    if (j0<20) {
+        if (j0<0) { /* raise inexact if x != 0 */
+            if (huge+x>0.0) { /* return 0*sign(x) if |x|<1 */
+                if (i0>=0) {i0=i1=0;}
+                else if (((i0&0x7fffffff)|i1)!=0)
+                    { i0=0xbff00000; i1=0;}
+            }
+        } else {
+            i = (0x000fffff)>>j0;
+            if (((i0&i)|i1)==0) return x; /* x is integral */
+            if (huge+x>0.0) {   /* raise inexact flag */
+                if (i0<0) i0 += (0x00100000)>>j0;
+                i0 &= (~i);
+                i1=0;
+            }
+        }
+    } else if (j0>51) {
+        if (j0==0x400) return x+x;  /* inf or NaN */
+        else return x;      /* x is integral */
+    } else {
+        i = ((u_int32_t)(0xffffffff))>>(j0-20);
+        if ((i1&i)==0) return x; /* x is integral */
+        if (huge+x>0.0) {       /* raise inexact flag */
+            if (i0<0) {
+                if (j0==20) i0+=1;
+                else {
+                    j = i1+(1<<(52-j0));
+                    if (j<(u_int32_t)i1) i0 +=1 ;   /* got a carry */
+                    i1=j;
+                }
+            }
+            i1 &= (~i);
+        }
+    }
+    INSERT_WORDS(x,i0,i1);
+    return x;
+}
+
+#if SUPPORT_LONG_DOUBLE
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(floor, floorl);
+#endif
+
+#endif
diff --git a/src/bsp/lk/lib/libm/s_round.c b/src/bsp/lk/lib/libm/s_round.c
new file mode 100644
index 0000000..dd293aa
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_round.c
@@ -0,0 +1,62 @@
+/*-
+ * Copyright (c) 2003, Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+double
+round(double x)
+{
+	double t;
+	uint32_t hx;
+
+	GET_HIGH_WORD(hx, x);
+	if ((hx & 0x7fffffff) == 0x7ff00000)
+		return (x + x);
+
+	if (!(hx & 0x80000000)) {
+		t = floor(x);
+		if (t - x <= -0.5)
+			t += 1;
+		return (t);
+	} else {
+		t = floor(-x);
+		if (t + x <= -0.5)
+			t += 1;
+		return (-t);
+	}
+}
+
+#if SUPPORT_LONG_DOUBLE
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(round, roundl);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/s_scalbn.c b/src/bsp/lk/lib/libm/s_scalbn.c
new file mode 100644
index 0000000..ef5a685
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_scalbn.c
@@ -0,0 +1,70 @@
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * scalbn (double x, int n)
+ * scalbn(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+#include <sys/cdefs.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
+
+double
+scalbn (double x, int n)
+{
+    int32_t k,hx,lx;
+    EXTRACT_WORDS(hx,lx,x);
+    k = (hx&0x7ff00000)>>20;        /* extract exponent */
+    if (k==0) {             /* 0 or subnormal x */
+        if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
+        x *= two54;
+        GET_HIGH_WORD(hx,x);
+        k = ((hx&0x7ff00000)>>20) - 54;
+        if (n< -50000) return tiny*x;   /*underflow*/
+    }
+    if (k==0x7ff) return x+x;       /* NaN or Inf */
+    k = k+n;
+    if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
+    if (k > 0)              /* normal result */
+        {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
+
+    if (k <= -54) {
+        if (n > 50000) {  /* in case integer overflow in n+k */
+            return huge*copysign(huge,x);   /*overflow*/
+        } else {
+            return tiny*copysign(tiny,x);   /*underflow*/
+        }
+    }
+    k += 54;                /* subnormal result */
+    SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
+    return x*twom54;
+}
+
+#if SUPPORT_LONG_DOUBLE
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(scalbn, ldexpl);
+__weak_reference(scalbn, scalbnl);
+#endif
+
+#endif
diff --git a/src/bsp/lk/lib/libm/s_scalbnf.c b/src/bsp/lk/lib/libm/s_scalbnf.c
new file mode 100644
index 0000000..dbbc2dd
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_scalbnf.c
@@ -0,0 +1,55 @@
+/* s_scalbnf.c -- float version of s_scalbn.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const float
+two25   =  3.355443200e+07,     /* 0x4c000000 */
+twom25  =  2.9802322388e-08,    /* 0x33000000 */
+huge   = 1.0e+30,
+tiny   = 1.0e-30;
+
+float
+scalbnf (float x, int n)
+{
+        int32_t k,ix;
+        GET_FLOAT_WORD(ix,x);
+        k = (ix&0x7f800000)>>23;                /* extract exponent */
+        if (k==0) {                             /* 0 or subnormal x */
+            if ((ix&0x7fffffff)==0) return x; /* +-0 */
+            x *= two25;
+            GET_FLOAT_WORD(ix,x);
+            k = ((ix&0x7f800000)>>23) - 25;
+            if (n< -50000) return tiny*x;       /*underflow*/
+            }
+        if (k==0xff) return x+x;                /* NaN or Inf */
+        k = k+n;
+        if (k >  0xfe) return huge*copysignf(huge,x); /* overflow  */
+        if (k > 0)                              /* normal result */
+            {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
+        if (k <= -25) {
+            if (n > 50000)      /* in case integer overflow in n+k */
+                return huge*copysignf(huge,x);  /*overflow*/
+            else return tiny*copysignf(tiny,x); /*underflow*/
+        }
+        k += 25;                                /* subnormal result */
+        SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
+        return x*twom25;
+}
+
+//__strong_reference(scalbnf, ldexpf);
diff --git a/src/bsp/lk/lib/libm/s_sin.c b/src/bsp/lk/lib/libm/s_sin.c
new file mode 100644
index 0000000..8d505c8
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_sin.c
@@ -0,0 +1,96 @@
+/* @(#)s_sin.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* sin(x)
+ * Return sine function of x.
+ *
+ * kernel function:
+ *  __kernel_sin        ... sine function on [-pi/4,pi/4]
+ *  __kernel_cos        ... cose function on [-pi/4,pi/4]
+ *  __ieee754_rem_pio2  ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *  [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *  in [-pi/4 , +pi/4], and let n = k mod 4.
+ *  We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *      0          S       C         T
+ *      1          C      -S        -1/T
+ *      2         -S      -C         T
+ *      3         -C       S        -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *  TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_REM_PIO2
+#include "math_private.h"
+#include "e_rem_pio2.c"
+
+double
+sin(double x)
+{
+    double y[2],z=0.0;
+    int32_t n, ix;
+
+    /* High word of x. */
+    GET_HIGH_WORD(ix,x);
+
+    /* |x| ~< pi/4 */
+    ix &= 0x7fffffff;
+    if (ix <= 0x3fe921fb) {
+        if (ix<0x3e500000)          /* |x| < 2**-26 */
+            {if ((int)x==0) return x;}   /* generate inexact */
+        return __kernel_sin(x,z,0);
+    }
+
+    /* sin(Inf or NaN) is NaN */
+    else if (ix>=0x7ff00000) return x-x;
+
+    /* argument reduction needed */
+    else {
+        n = __ieee754_rem_pio2(x,y);
+        switch (n&3) {
+            case 0:
+                return  __kernel_sin(y[0],y[1],1);
+            case 1:
+                return  __kernel_cos(y[0],y[1]);
+            case 2:
+                return -__kernel_sin(y[0],y[1],1);
+            default:
+                return -__kernel_cos(y[0],y[1]);
+        }
+    }
+}
+
+#if SUPPORT_LONG_DOUBLE
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sin, sinl);
+#endif
+
+#endif
diff --git a/src/bsp/lk/lib/libm/s_sinf.c b/src/bsp/lk/lib/libm/s_sinf.c
new file mode 100644
index 0000000..30c055c
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_sinf.c
@@ -0,0 +1,88 @@
+/* s_sinf.c -- float version of s_sin.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_KERNEL_COSDF
+#define INLINE_KERNEL_SINDF
+#define INLINE_REM_PIO2F
+#include "math_private.h"
+#include "e_rem_pio2f.c"
+#include "k_cosf.c"
+#include "k_sinf.c"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+s1pio2 = 1*M_PI_2,          /* 0x3FF921FB, 0x54442D18 */
+s2pio2 = 2*M_PI_2,          /* 0x400921FB, 0x54442D18 */
+s3pio2 = 3*M_PI_2,          /* 0x4012D97C, 0x7F3321D2 */
+s4pio2 = 4*M_PI_2;          /* 0x401921FB, 0x54442D18 */
+
+float
+sinf(float x)
+{
+    double y;
+    int32_t n, hx, ix;
+
+    GET_FLOAT_WORD(hx,x);
+    ix = hx & 0x7fffffff;
+
+    if (ix <= 0x3f490fda) {     /* |x| ~<= pi/4 */
+        if (ix<0x39800000)      /* |x| < 2**-12 */
+            if (((int)x)==0) return x;  /* x with inexact if x != 0 */
+        return __kernel_sindf(x);
+    }
+    if (ix<=0x407b53d1) {       /* |x| ~<= 5*pi/4 */
+        if (ix<=0x4016cbe3) {   /* |x| ~<= 3pi/4 */
+            if (hx>0)
+                return __kernel_cosdf(x - s1pio2);
+            else
+                return -__kernel_cosdf(x + s1pio2);
+        } else
+            return __kernel_sindf((hx > 0 ? s2pio2 : -s2pio2) - x);
+    }
+    if (ix<=0x40e231d5) {       /* |x| ~<= 9*pi/4 */
+        if (ix<=0x40afeddf) {   /* |x| ~<= 7*pi/4 */
+            if (hx>0)
+                return -__kernel_cosdf(x - s3pio2);
+            else
+                return __kernel_cosdf(x + s3pio2);
+        } else
+            return __kernel_sindf(x + (hx > 0 ? -s4pio2 : s4pio2));
+    }
+
+    /* sin(Inf or NaN) is NaN */
+    else if (ix>=0x7f800000) return x-x;
+
+    /* general argument reduction needed */
+    else {
+        n = __ieee754_rem_pio2f(x,&y);
+        switch (n&3) {
+            case 0:
+                return  __kernel_sindf(y);
+            case 1:
+                return  __kernel_cosdf(y);
+            case 2:
+                return  __kernel_sindf(-y);
+            default:
+                return -__kernel_cosdf(y);
+        }
+    }
+}
diff --git a/src/bsp/lk/lib/libm/s_tan.c b/src/bsp/lk/lib/libm/s_tan.c
new file mode 100644
index 0000000..36b1138
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_tan.c
@@ -0,0 +1,84 @@
+/* @(#)s_tan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* tan(x)
+ * Return tangent function of x.
+ *
+ * kernel function:
+ *  __kernel_tan        ... tangent function on [-pi/4,pi/4]
+ *  __ieee754_rem_pio2  ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *  [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *  in [-pi/4 , +pi/4], and let n = k mod 4.
+ *  We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *      0          S       C         T
+ *      1          C      -S        -1/T
+ *      2         -S      -C         T
+ *      3         -C       S        -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *  TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_REM_PIO2
+#include "math_private.h"
+#include "e_rem_pio2.c"
+
+double
+tan(double x)
+{
+    double y[2],z=0.0;
+    int32_t n, ix;
+
+    /* High word of x. */
+    GET_HIGH_WORD(ix,x);
+
+    /* |x| ~< pi/4 */
+    ix &= 0x7fffffff;
+    if (ix <= 0x3fe921fb) {
+        if (ix<0x3e400000)          /* x < 2**-27 */
+            if ((int)x==0) return x;    /* generate inexact */
+        return __kernel_tan(x,z,1);
+    }
+
+    /* tan(Inf or NaN) is NaN */
+    else if (ix>=0x7ff00000) return x-x;        /* NaN */
+
+    /* argument reduction needed */
+    else {
+        n = __ieee754_rem_pio2(x,y);
+        return __kernel_tan(y[0],y[1],1-((n&1)<<1)); /*   1 -- n even
+                            -1 -- n odd */
+    }
+}
+#if SUPPORT_LONG_DOUBLE
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(tan, tanl);
+#endif
+#endif
diff --git a/src/bsp/lk/lib/libm/s_tanf.c b/src/bsp/lk/lib/libm/s_tanf.c
new file mode 100644
index 0000000..abd60a2
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_tanf.c
@@ -0,0 +1,72 @@
+/* s_tanf.c -- float version of s_tan.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+
+#include "math.h"
+#define INLINE_KERNEL_TANDF
+#define INLINE_REM_PIO2F
+#include "math_private.h"
+#include "e_rem_pio2f.c"
+#include "k_tanf.c"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+t1pio2 = 1*M_PI_2,          /* 0x3FF921FB, 0x54442D18 */
+t2pio2 = 2*M_PI_2,          /* 0x400921FB, 0x54442D18 */
+t3pio2 = 3*M_PI_2,          /* 0x4012D97C, 0x7F3321D2 */
+t4pio2 = 4*M_PI_2;          /* 0x401921FB, 0x54442D18 */
+
+float
+tanf(float x)
+{
+    double y;
+    int32_t n, hx, ix;
+
+    GET_FLOAT_WORD(hx,x);
+    ix = hx & 0x7fffffff;
+
+    if (ix <= 0x3f490fda) {     /* |x| ~<= pi/4 */
+        if (ix<0x39800000)      /* |x| < 2**-12 */
+            if (((int)x)==0) return x;  /* x with inexact if x != 0 */
+        return __kernel_tandf(x,1);
+    }
+    if (ix<=0x407b53d1) {       /* |x| ~<= 5*pi/4 */
+        if (ix<=0x4016cbe3)     /* |x| ~<= 3pi/4 */
+            return __kernel_tandf(x + (hx>0 ? -t1pio2 : t1pio2), -1);
+        else
+            return __kernel_tandf(x + (hx>0 ? -t2pio2 : t2pio2), 1);
+    }
+    if (ix<=0x40e231d5) {       /* |x| ~<= 9*pi/4 */
+        if (ix<=0x40afeddf)     /* |x| ~<= 7*pi/4 */
+            return __kernel_tandf(x + (hx>0 ? -t3pio2 : t3pio2), -1);
+        else
+            return __kernel_tandf(x + (hx>0 ? -t4pio2 : t4pio2), 1);
+    }
+
+    /* tan(Inf or NaN) is NaN */
+    else if (ix>=0x7f800000) return x-x;
+
+    /* general argument reduction needed */
+    else {
+        n = __ieee754_rem_pio2f(x,&y);
+        /* integer parameter: 1 -- n even; -1 -- n odd */
+        return __kernel_tandf(y,1-((n&1)<<1));
+    }
+}
diff --git a/src/bsp/lk/lib/libm/s_trunc.c b/src/bsp/lk/lib/libm/s_trunc.c
new file mode 100644
index 0000000..bc996d1
--- /dev/null
+++ b/src/bsp/lk/lib/libm/s_trunc.c
@@ -0,0 +1,69 @@
+/* @(#)s_floor.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * trunc(x)
+ * Return x rounded toward 0 to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to trunc(x).
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+static const double huge = 1.0e300;
+
+double
+trunc(double x)
+{
+	int32_t i0,i1,j0;
+	u_int32_t i;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;
+	if(j0<20) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* |x|<1, so return 0*sign(x) */
+		    i0 &= 0x80000000U;
+		    i1 = 0;
+		}
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (j0>51) {
+	    if(j0==0x400) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0)		/* raise inexact flag */
+		i1 &= (~i);
+	}
+	INSERT_WORDS(x,i0,i1);
+	return x;
+}
+
+#ifdef SUPPORT_LONG_DOUBLE
+#if LDBL_MANT_DIG == 53
+__weak_reference(trunc, truncl);
+#endif
+#endif