overby.me
this repo has no description
1// double hypot(double x, double y)
2// double cabs(double complex z)
3//
4// Returns sqrt(x*x + y*y), computed without spurious overflow or underflow.
5//
6// -- Stephen Canon, August 2009
7
8#if !defined __i386__ && !defined __x86_64__
9#error This implementation is for i386 and x86_64 only
10#endif
11
12.text
13.globl _hypot
14.globl _cabs
15
16#if defined __i386__
17
18#define ABSHI $0x7fffffff
19#define INFHI $0x7ff00000
20
21// Entry point --------------------------------------------------------
22
23.align 4
24_cabs:
25_hypot:
26 mov 8(%esp), %eax
27 mov 16(%esp), %edx
28 mov ABSHI, %ecx
29
30// Detect special cases -----------------------------------------------
31
32 and %ecx, %eax // |x|
33 jz L_xHiIsZero //
34L_returnFromXHiIsZero:
35 and %ecx, %edx // |y|
36 jz L_yHiIsZero //
37L_returnFromYHiIsZero:
38 cmp INFHI, %eax // we need to detect if either x or y is
39 jae L_xIsSpecial // an infinity and return +infinity if so;
40 cmp INFHI, %edx // this is because hypot(inf, NaN) needs
41 jae L_returnAbsY // to return +inf, not NaN.
42
43// Mainline computation -----------------------------------------------
44
45 fldl 4(%esp) // x87 stack: { x }
46 fmul %st(0), %st(0) // { x*x }
47 fldl 12(%esp) // { y x*x }
48 fmul %st(0), %st(0) // { y*y x*x }
49 faddp // { y*y+x*x }
50 fsqrt // { sqrt(y*y + x*x) }
51 fstpl 4(%esp) // round to double
52 fldl 4(%esp)
53 ret
54
55// Special case handling ----------------------------------------------
56
57L_xHiIsZero:
58 cmpl $0, 4(%esp) // is the low word of x zero?
59 jnz L_returnFromXHiIsZero // if not, jump back to mainline
60L_returnAbsY:
61 and %ecx, 16(%esp)
62 fldl 12(%esp)
63 ret
64
65L_yHiIsZero:
66 cmpl $0, 12(%esp) // is the low word of y zero?
67 jnz L_returnFromYHiIsZero // if not, jump back to mainline
68L_returnAbsX:
69 and %ecx, 8(%esp)
70 fldl 4(%esp)
71 ret
72
73L_xIsSpecial:
74 cmp INFHI, %edx // check if y is infinity
75 jnz L_returnAbsX
76 cmpl $0, 12(%esp)
77 jz L_returnAbsY
78 jmp L_returnAbsX
79
80#else // defined __i386__
81
82#define ABSMASK $0x7fffffffffffffff
83#define INFINITY $0x7ff0000000000000
84
85// Entry point --------------------------------------------------------
86
87.align 4
88_cabs:
89_hypot:
90 movd %xmm0, %rax
91 movd %xmm1, %rdx
92 mov ABSMASK, %rcx
93
94// Detect special cases -----------------------------------------------
95
96 and %rcx, %rax // |x|
97 jz L_returnAbsY //
98 and %rcx, %rdx // |y|
99 jz L_returnAbsX //
100 mov INFINITY, %rdi //
101
102 movsd %xmm0, -8(%rsp) // store x to the red zone
103 movsd %xmm1, -16(%rsp) // store y to the red zone
104
105 cmp %rdi, %rax // we need to detect if either x or y is
106 jae L_xIsSpecial // an infinity and return +infinity if so;
107 cmp %rdi, %rdx // this is because hypot(inf, NaN) needs
108 jae L_returnAbsY // to return +inf, not NaN.
109
110// Mainline computation -----------------------------------------------
111
112 fldl -8(%rsp) // x87 stack: { x }
113 fmul %st(0), %st(0) // { x*x }
114 fldl -16(%rsp) // { y x*x }
115 fmul %st(0), %st(0) // { y*y x*x }
116 faddp // { y*y+x*x }
117 fsqrt // { sqrt(y*y + x*x) }
118 fstpl -8(%rsp) // round result to double
119 movsd -8(%rsp), %xmm0 // and load in xmm register to return
120 ret
121
122// Special case handling ----------------------------------------------
123
124L_xIsSpecial: // x is either inf or NaN
125 cmp %rdi, %rdx
126 jz L_returnAbsY // if y is infinity, return |y|
127L_returnAbsX:
128 movd %rax, %xmm0 // otherwise, return |x|
129 ret
130L_returnAbsY:
131 and %rcx, %rdx
132 movd %rdx, %xmm0 // return |y|
133 ret
134
135#endif // defined __i386__