diff --git a/libm/src/e_pow.c b/libm/src/e_pow.c
index d213132..b0a3f53 100644
--- a/libm/src/e_pow.c
+++ b/libm/src/e_pow.c
@@ -61,6 +61,14 @@
 #include "math.h"
 #include "math_private.h"
 
+#if defined(KRAIT_NEON_OPTIMIZATION) || defined(SPARROW_NEON_OPTIMIZATION)
+#if defined(KRAIT_NO_AAPCS_VFP_MODE)
+double pow_neon(double x, double y);
+#else
+double pow_neon(double x, double y, int32_t lx, int32_t hx) __attribute__((pcs("aapcs-vfp")));
+#endif
+#endif
+
 static const double
 bp[] = {1.0, 1.5,},
 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
@@ -108,12 +116,32 @@
 	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;	
+    if (ly == 0) {
+        if (hy == ly) {
+            /* y==0.0, x**0 = 1 */
+            return one;
+        }
+        else if (iy > 0x7ff00000) {
+            /* y is NaN, return x+y (NaN) */
+            return x+y;
+        }
+    }
+    else if (iy >= 0x7ff00000) {
+        /* y is NaN, return x+y (NaN) */
+        return x+y;
+    }
+
+    if (lx == 0) {
+        if (ix > 0x7ff00000) {
+            /* x is NaN, return x+y (NaN) */
+            return x+y;
+        }
+    }
+    else if (ix >= 0x7ff00000) {
+        /* x is NaN, return x+y (NaN) */
+        return x+y;
+    }
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
@@ -201,6 +229,14 @@
 	    t1 = u+v;
 	    SET_LOW_WORD(t1,0);
 	    t2 = v-(t1-u);
+#if defined(KRAIT_NEON_OPTIMIZATION) || defined(SPARROW_NEON_OPTIMIZATION)
+	} else if (ix <= 0x40100000 && iy <= 0x40100000 && hy > 0 && hx > 0) {
+#if defined(KRAIT_NO_AAPCS_VFP_MODE)
+		return pow_neon(x,y);
+#else
+		return pow_neon(x,y,lx,hx);
+#endif
+#endif
 	} else {
 	    double ss,s2,s_h,s_l,t_h,t_l;
 	    n = 0;
diff --git a/libm/src/k_cos.c b/libm/src/k_cos.c
index 00916d7..b8cdf8f 100644
--- a/libm/src/k_cos.c
+++ b/libm/src/k_cos.c
@@ -69,6 +69,17 @@
 double
 __kernel_cos(double x, double y)
 {
+#if defined(KRAIT_NEON_OPTIMIZATION)
+	double hz,z,zz,r,w,k;
+
+	z  = x*x;
+	zz = z*z;
+	k = x*y;
+	hz = (float)0.5*z;
+	r  = z*(z*(C1+z*(C2+z*((C3+z*C4)+zz*(C5+z*C6)))));
+	w  = one-hz;
+	return w + (((one-w)-hz) + (r-k));
+#else
 	double hz,z,r,w;
 
 	z  = x*x;
@@ -76,4 +87,5 @@
 	hz = (float)0.5*z;
 	w  = one-hz;
 	return w + (((one-w)-hz) + (z*r-x*y));
+#endif
 }
diff --git a/libm/src/k_sin.c b/libm/src/k_sin.c
index ae06a9d..ee641d4 100644
--- a/libm/src/k_sin.c
+++ b/libm/src/k_sin.c
@@ -60,6 +60,16 @@
 double
 __kernel_sin(double x, double y, int iy)
 {
+#if defined(KRAIT_NEON_OPTIMIZATION)
+	double z,zz,r,v;
+
+	z	=  x*x;
+	zz  =  z*z;
+	v	=  z*x;
+	r	=  S2+z*((S3+z*S4)+zz*(S5+z*S6));
+	if(iy==0) return x+v*(S1+z*r);
+	else      return x-((z*(half*y-v*r)-y)-v*S1);
+#else
 	double z,r,v;
 
 	z	=  x*x;
@@ -67,4 +77,5 @@
 	r	=  S2+z*(S3+z*(S4+z*(S5+z*S6)));
 	if(iy==0) return x+v*(S1+z*r);
 	else      return x-((z*(half*y-v*r)-y)-v*S1);
+#endif
 }
diff --git a/libm/src/math_private.h b/libm/src/math_private.h
index 5f6e088..7cda2e9 100644
--- a/libm/src/math_private.h
+++ b/libm/src/math_private.h
@@ -257,11 +257,19 @@
 #define __ieee754_ldexpf ldexpf
 
 /* fdlibm kernel function */
+#if defined(KRAIT_NEON_OPTIMIZATION)
+int	__ieee754_rem_pio2(double,double*) __attribute__((pcs("aapcs-vfp")));
+double	__kernel_sin(double,double,int) __attribute__((pcs("aapcs-vfp")));
+double	__kernel_cos(double,double) __attribute__((pcs("aapcs-vfp")));
+double	__kernel_tan(double,double,int) __attribute__((pcs("aapcs-vfp")));
+int	__kernel_rem_pio2(double*,double*,int,int,int,const int*) __attribute__((pcs("aapcs-vfp")));
+#else
 int	__ieee754_rem_pio2(double,double*);
 double	__kernel_sin(double,double,int);
 double	__kernel_cos(double,double);
 double	__kernel_tan(double,double,int);
 int	__kernel_rem_pio2(double*,double*,int,int,int,const int*);
+#endif
 
 /* float versions of fdlibm kernel functions */
 int	__ieee754_rem_pio2f(float,float*);