Added stubs for additional newlib system calls
[akaros.git] / user / parlib / src / sparc / divsqrt.c
1 #include <math.h>
2
3 static inline double i2d(unsigned long long i)
4 {
5   double d;
6   asm volatile("std %2,[%1]; ldd [%1],%0" : "=f"(d) : "r"(&d),"r"(i));
7   return d;
8 }
9
10 static inline unsigned long long d2i(double d)
11 {
12   unsigned long long i;
13   asm volatile("std %2,[%1]; ldd [%1],%0" : "=r"(i) : "r"(&i),"f"(d));
14   return i;
15 }
16
17 double do_recip(double b)
18 {
19   unsigned long long i = d2i(b);
20   unsigned long long i2 = ((2046-((i>>52)&~0x800)) | ((i>>52)&0x800)) << 52;
21   unsigned long long i3 = (i >> 50) & 3;
22   static const double divlut[4] = {1.0,0.8,0.666,0.571};
23   double x = i2d(i2)*divlut[i3];
24
25   x = x*(2.0-b*x);
26   x = x*(2.0-b*x);
27   x = x*(2.0-b*x);
28   x = x*(2.0-b*x);
29   x = x*(2.0-b*x);
30
31   return x;
32 }
33
34 double do_rsqrt(double b)
35 {
36   unsigned long long i = d2i(b);
37   unsigned long long i2 = ((3069-((i>>52)&~0x800))>>1 | ((i>>52)&0x800)) << 52;
38   unsigned long long i3 = (i >> 50) & 7;
39   double x = i2d(i2);
40
41   static const double sqrtlut[8] = {1.4142,1.264,1.155,1.069, 1.0,0.894,0.816,0.756};
42   x *= sqrtlut[i3];
43
44   x = 0.5*x*(3.0-b*x*x);
45   x = 0.5*x*(3.0-b*x*x);
46   x = 0.5*x*(3.0-b*x*x);
47   x = 0.5*x*(3.0-b*x*x);
48   x = 0.5*x*(3.0-b*x*x);
49
50   return x;
51 }
52
53 double do_fdiv(double x, double y)
54 {
55   if((d2i(y) & 0x7FF0000000000000ULL) == 0)
56     return x/y;
57   return x*do_recip(y);
58 }
59
60 double do_fsqrt(double x)
61 {
62   if(d2i(x) & 0x8000000000000000ULL)
63     return sqrt(x);
64   return x*do_rsqrt(x);
65 }