| blob | 36cfd1e6677b8d4387b5edbeee7068a5eff83db8 |
1 .file "hypot.s"
4 // Copyright (c) 2000 - 2003, Intel Corporation
5 // All rights reserved.
6 //
7 // Contributed 2000 by the Intel Numerics Group, Intel Corporation
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are
11 // met:
12 //
13 // * Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
15 //
16 // * Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
19 //
20 // * The name of Intel Corporation may not be used to endorse or promote
21 // products derived from this software without specific prior written
22 // permission.
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35 //
36 // Intel Corporation is the author of this code, and requests that all
37 // problem reports or change requests be submitted to it directly at
38 // http://www.intel.com/software/products/opensource/libraries/num.htm.
39 //
40 //*********************************************************************
41 //
42 // History:
43 // 02/02/00 hand-optimized
44 // 04/04/00 Unwind support added
45 // 06/20/00 new version
46 // 08/15/00 Bundle added after call to __libm_error_support to properly
47 // set [the previously overwritten] GR_Parameter_RESULT.
48 // 05/20/02 Cleaned up namespace and sf0 syntax
49 // 02/10/03 Reordered header: .section, .global, .proc, .align
50 // 04/17/03 Added missing mutex directive
51 //
52 //*********************************************************************
53 // ___________
54 // Function: hypot(x,y) = |(x^2 + y^2) = for double precision values
55 // x and y
56 // Also provides cabs functionality.
57 //
58 //*********************************************************************
59 //
60 // Resources Used:
61 //
62 // Floating-Point Registers: f8 (Input and Return Value)
63 // f9 (Input)
64 // f6 -f15, f32-f34
65 //
66 // General Purpose Registers:
67 // r2,r3,r29 (Scratch)
68 // r32-r36 (Locals)
69 // r37-r40 (Used to pass arguments to error handling routine)
70 //
71 // Predicate Registers: p6 - p10
72 //
73 //*********************************************************************
74 //
75 // IEEE Special Conditions:
76 //
77 // All faults and exceptions should be raised correctly.
78 // Overflow can occur.
79 // hypot(Infinity and anything) = +Infinity
80 // hypot(QNaN and anything) = QNaN
81 // hypot(SNaN and anything ) = QNaN
82 //
83 //*********************************************************************
84 //
85 // Implementation:
86 // x2 = x * x in double-extended
87 // y2 = y * y in double-extended
88 // temp = x2 + y2 in double-extended
89 // sqrt(temp) rounded to double
90 //
91 //*********************************************************************
93 GR_SAVE_PFS = r33
94 GR_SAVE_B0 = r34
95 GR_SAVE_GP = r35
96 GR_Parameter_X = r36
97 GR_Parameter_Y = r37
98 GR_Parameter_RESULT = r38
99 GR_Parameter_TAG = r39
101 FR_X = f32
102 FR_Y = f33
103 FR_RESULT = f8
105 .section .text
107 LOCAL_LIBM_ENTRY(cabs)
108 LOCAL_LIBM_END(cabs)
110 GLOBAL_IEEE754_ENTRY(hypot)
112 {.mfi
113 alloc r32= ar.pfs,0,4,4,0
114 // Compute x*x
115 fma.s1 f10=f8,f8,f0
116 // r2=bias-1
117 mov r2=0xfffe
118 }
119 {.mfi
120 // 63/8
121 mov r3=0x40fc //0000
122 // y*y
123 fma.s1 f11=f9,f9,f0
124 // r29=429/16
125 mov r29=0x41d68;; //000
126 }
128 { .mfi
129 nop.m 0
130 // Check if x is an Inf - if so return Inf even
131 // if y is a NaN (C9X)
132 fclass.m.unc p7, p6 = f8, 0x023
133 shl r3=r3,16
134 }
135 {.mfi
136 nop.m 0
137 // if possible overflow, copy f8 to f32
138 // set Denormal, if necessary
139 // (p8)
140 fma.d.s0 f32=f8,f1,f0
141 nop.i 0;;
142 }
143 { .mfi
144 nop.m 0
145 // Check if y is an Inf - if so return Inf even
146 // if x is a NaN (C9X)
147 fclass.m.unc p8, p9 = f9, 0x023
148 shl r29=r29,12
149 }
150 { .mfb
151 // f7=0.5
152 setf.exp f7=r2
153 // For x=inf, multiply y by 1 to raise invalid on y an SNaN
154 // (p7) fma.s0 f9=f9,f1,f0
155 // copy f9 to f33; set Denormal, if necessary
156 fma.d.s0 f33=f9,f1,f0
157 nop.b 0;;
158 }
159 {.mfb
160 // f13=63/8
161 setf.s f13=r3
162 // is y Zero ?
163 (p6) fclass.m p6,p0=f9,0x7
164 nop.b 0
165 }
166 {.mlx
167 nop.m 0
168 movl r2=0x408c0000;;
169 }
171 {.mfi
172 // f34=429/16
173 setf.s f34=r29
174 // is x Zero ?
175 (p9) fclass.m p9,p0=f8,0x7
176 // 231/16
177 mov r3=0x4167;; //0000
178 }
179 {.mfi
180 nop.m 0
181 // a=x2+y2
182 fma.s1 f12=f10,f1,f11
183 nop.i 0;;
184 }
185 {.mfi
186 nop.m 0
187 // y not NaN ?
188 (p9) fclass.m p8,p0=f9,0x3f
189 shl r3=r3,16
190 }
191 {.mfi
192 nop.m 0
193 // f6=2
194 fma.s1 f6=f1,f1,f1
195 nop.i 0;;
196 }
199 {.mfi
200 nop.m 0
201 // x not NaN ?
202 (p6) fclass.m p7,p0=f8,0x3f
203 nop.i 0;;
204 }
205 {.mfi
206 // f9=35/8
207 setf.s f9=r2
208 nop.f 0
209 // 2*emax-2
210 mov r2=0x107fb;;
211 }
213 .pred.rel "mutex",p7,p8
214 {.mfb
215 nop.m 0
216 // if f8=Infinity or f9=Zero, return |f8|
217 (p7) fmerge.s f8=f0,f32
218 (p7) br.ret.spnt b0
219 }
220 {.mfb
221 nop.m 0
222 // if f9=Infinity or f8=Zero, return |f9|
223 (p8) fmerge.s f8=f0,f33
224 (p8) br.ret.spnt b0;;
225 }
228 {.mfi
229 // f10 =231/16
230 setf.s f10=r3
231 // z0=frsqrta(a)
232 frsqrta.s1 f8,p6=f12
233 nop.i 0;;
234 }
236 { .mfi
237 nop.m 0
238 // Identify Natvals, Infs, NaNs, and Zeros
239 // and return result
240 fclass.m.unc p7, p0 = f12, 0x1E7
241 nop.i 0;;
242 }
243 {.mfb
244 // get exponent of x^2+y^2
245 getf.exp r3=f12
246 // if special case, set f8
247 (p7) mov f8=f12
248 (p7) br.ret.spnt b0;;
249 }
252 {.mfi
253 nop.m 0
254 // S0=a*z0
255 (p6) fma.s1 f14=f12,f8,f0
256 nop.i 0
257 }
258 {.mfi
259 nop.m 0
260 // H0=0.5*z0
261 (p6) fma.s1 f15=f8,f7,f0
262 nop.i 0;;
263 }
266 {.mfi
267 nop.m 0
268 // f6=5/2
269 fma.s1 f6=f7,f1,f6
270 nop.i 0
271 }
272 {.mfi
273 nop.m 0
274 // f11=3/2
275 fma.s1 f11=f7,f1,f1
276 nop.i 0;;
277 }
279 {.mfi
280 nop.m 0
281 // d=0.5-S0*H0
282 (p6) fnma.s1 f7=f14,f15,f7
283 nop.i 0;;
284 }
286 {.mfi
287 nop.m 0
288 // P67=231/16+429/16*d
289 (p6) fma.s1 f10=f34,f7,f10
290 nop.i 0
291 }
292 {.mfi
293 nop.m 0
294 // P45=63/8*d+35/8
295 (p6) fma.s1 f9=f13,f7,f9
296 nop.i 0;;
297 }
298 {.mfi
299 nop.m 0
300 // P23=5/2*d+3/2
301 (p6) fma.s1 f11=f6,f7,f11
302 nop.i 0
303 }
304 {.mfi
305 nop.m 0
306 // d2=d*d
307 (p6) fma.s1 f13=f7,f7,f0
308 nop.i 0;;
309 }
311 {.mfi
312 nop.m 0
313 // P47=d2*P67+P45
314 (p6) fma.s1 f10=f10,f13,f9
315 nop.i 0
316 }
317 {.mfi
318 nop.m 0
319 // P13=d*P23+1
320 (p6) fma.s1 f11=f11,f7,f1
321 nop.i 0;;
322 }
323 {.mfi
324 nop.m 0
325 // d3=d2*d
326 (p6) fma.s1 f13=f13,f7,f0
327 nop.i 0;;
328 }
330 {.mfi
331 nop.m 0
332 // T0=d*S0
333 (p6) fma.s1 f15=f7,f14,f0
334 nop.i 0
335 }
336 {.mfi
337 // Is x^2 + y^2 well less than the overflow
338 // threshold?
339 (p6) cmp.lt.unc p7, p8 = r3,r2
340 // P=P13+d3*P47
341 (p6) fma.s1 f10=f13,f10,f11
342 nop.i 0;;
343 }
345 {.mfb
346 nop.m 0
347 // S=P*T0+S0
348 fma.d.s0 f8=f10,f15,f14
349 // No overflow in this case
350 (p7) br.ret.sptk b0;;
351 }
353 { .mfi
354 nop.m 0
355 (p8) fsetc.s2 0x7F,0x42
356 // Possible overflow path, must detect by
357 // Setting widest range exponent with prevailing
358 // rounding mode.
359 nop.i 0 ;;
360 }
363 { .mfi
364 // bias+0x400 (bias+EMAX+1)
365 (p8) mov r2=0x103ff
366 // S=P*T0+S0
367 (p8) fma.d.s2 f12=f10,f15,f14
368 nop.i 0 ;;
369 }
370 { .mfi
371 (p8) setf.exp f11 = r2
372 (p8) fsetc.s2 0x7F,0x40
373 // Restore Original Mode in S2
374 nop.i 0 ;;
375 }
376 { .mfi
377 nop.m 0
378 (p8) fcmp.lt.unc.s1 p9, p10 = f12, f11
379 nop.i 0 ;;
380 }
381 { .mib
382 nop.m 0
383 mov GR_Parameter_TAG = 46
384 // No overflow
385 (p9) br.ret.sptk b0;;
386 }
387 GLOBAL_IEEE754_END(hypot)
389 LOCAL_LIBM_ENTRY(__libm_error_region)
390 .prologue
391 { .mfi
392 add GR_Parameter_Y=-32,sp // Parameter 2 value
393 nop.f 0
394 .save ar.pfs,GR_SAVE_PFS
395 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
396 }
397 { .mfi
398 .fframe 64
399 add sp=-64,sp // Create new stack
400 nop.f 0
401 mov GR_SAVE_GP=gp // Save gp
402 };;
403 { .mmi
404 stfd [GR_Parameter_Y] = FR_Y,16 // Save Parameter 2 on stack
405 add GR_Parameter_X = 16,sp // Parameter 1 address
406 .save b0, GR_SAVE_B0
407 mov GR_SAVE_B0=b0 // Save b0
408 };;
409 .body
410 { .mib
411 stfd [GR_Parameter_X] = FR_X // Store Parameter 1 on stack
412 add GR_Parameter_RESULT = 0,GR_Parameter_Y
413 nop.b 0 // Parameter 3 address
414 }
415 { .mib
416 stfd [GR_Parameter_Y] = FR_RESULT // Store Parameter 3 on stack
417 add GR_Parameter_Y = -16,GR_Parameter_Y
418 br.call.sptk b0=__libm_error_support# // Call error handling function
419 };;
420 { .mmi
421 nop.m 0
422 nop.m 0
423 add GR_Parameter_RESULT = 48,sp
424 };;
425 { .mmi
426 ldfd f8 = [GR_Parameter_RESULT] // Get return result off stack
427 .restore sp
428 add sp = 64,sp // Restore stack pointer
429 mov b0 = GR_SAVE_B0 // Restore return address
430 };;
431 { .mib
432 mov gp = GR_SAVE_GP // Restore gp
433 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
434 br.ret.sptk b0 // Return
435 };;
436 LOCAL_LIBM_END(__libm_error_region#)
437 .type __libm_error_support#,@function
438 .global __libm_error_support#