| blob | 60be3b3ffcf2f2680ced03671168a85f6297038e |
1 .file "scalb.s"
3 // Copyright (c) 2000, 2001, Intel Corporation
4 // All rights reserved.
5 //
6 // Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
7 // and Ping Tak Peter Tang of the Computational Software Lab, Intel Corporation.
8 //
9 // WARRANTY DISCLAIMER
10 //
11 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
12 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
13 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
14 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
15 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
16 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
17 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
18 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
19 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
20 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
21 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
22 //
23 // Intel Corporation is the author of this code, and requests that all
24 // problem reports or change requests be submitted to it directly at
25 // http://developer.intel.com/opensource.
26 //
27 // History
28 //==============================================================
29 // 2/02/00 Initial version
30 // 1/26/01 Scalb completely reworked and now standalone version
31 //
32 // API
33 //==============================================================
34 // double = scalb (double x, double n)
35 // input floating point f8 and floating point f9
36 // output floating point f8
37 //
38 // Returns x* 2**n using an fma and detects overflow
39 // and underflow.
40 //
41 //
43 #include "libm_support.h"
45 FR_Floating_X = f8
46 FR_Result = f8
47 FR_Floating_N = f9
48 FR_Result2 = f9
49 FR_Norm_N = f10
50 FR_Result3 = f11
51 FR_Norm_X = f12
52 FR_N_float_int = f13
53 FR_Two_N = f14
54 FR_Two_to_Big = f15
55 FR_Big = f6
56 FR_NBig = f7
58 GR_N_Biased = r15
59 GR_Big = r16
60 GR_NBig = r17
61 GR_Scratch = r18
62 GR_Scratch1 = r19
63 GR_Bias = r20
64 GR_N_as_int = r21
66 GR_SAVE_B0 = r32
67 GR_SAVE_GP = r33
68 GR_SAVE_PFS = r34
69 GR_Parameter_X = r35
70 GR_Parameter_Y = r36
71 GR_Parameter_RESULT = r37
72 GR_Tag = r38
74 .align 32
75 .global scalb
77 .section .text
78 .proc scalb
79 .align 32
81 scalb:
82 #ifdef _LIBC
83 .global __ieee754_scalb
84 .type __ieee754_scalb,@function
85 __ieee754_scalb:
86 #endif
88 //
89 // Is x NAN, INF, ZERO, +-?
90 //
91 { .mfi
92 alloc r32=ar.pfs,0,3,4,0
93 fclass.m.unc p7,p0 = FR_Floating_X, 0xe7 //@snan | @qnan | @inf | @zero
94 addl GR_Scratch = 0x019C3F,r0
95 }
96 //
97 // Is y a NAN, INF, ZERO, +-?
98 //
99 { .mfi
100 nop.m 999
101 fclass.m.unc p6,p0 = FR_Floating_N, 0xe7 //@snan | @qnan | @inf | @zero
102 addl GR_Scratch1 = 0x063BF,r0
103 }
104 ;;
106 //
107 // Convert N to a fp integer
108 // Normalize x
109 //
110 { .mfi
111 nop.m 0
112 fnorm.s1 FR_Norm_N = FR_Floating_N
113 nop.i 999
114 }
115 { .mfi
116 nop.m 999
117 fnorm.s1 FR_Norm_X = FR_Floating_X
118 nop.i 999
119 };;
121 //
122 // Create 2*big
123 // Create 2**-big
124 // Normalize x
125 // Branch on special values.
126 //
127 { .mib
128 setf.exp FR_Big = GR_Scratch
129 nop.i 0
130 (p6) br.cond.spnt L(SCALB_NAN_INF_ZERO)
131 }
132 { .mib
133 setf.exp FR_NBig = GR_Scratch1
134 nop.i 0
135 (p7) br.cond.spnt L(SCALB_NAN_INF_ZERO)
136 };;
138 //
139 // Convert N to a fp integer
140 // Create -35000
141 //
142 { .mfi
143 addl GR_Scratch = 1,r0
144 fcvt.fx.trunc.s1 FR_N_float_int = FR_Norm_N
145 addl GR_NBig = -35000,r0
146 }
147 ;;
149 //
150 // Put N if a GP register
151 // Convert N_float_int to floating point value
152 // Create 35000
153 // Build the exponent Bias
154 //
155 { .mii
156 getf.sig GR_N_as_int = FR_N_float_int
157 shl GR_Scratch = GR_Scratch,63
158 addl GR_Big = 35000,r0
159 }
160 { .mfi
161 addl GR_Bias = 0x0FFFF,r0
162 fcvt.xf FR_N_float_int = FR_N_float_int
163 nop.i 0
164 };;
166 //
167 // Catch those fp values that are beyond 2**64-1
168 // Is N > 35000
169 // Is N < -35000
170 //
171 { .mfi
172 cmp.ne.unc p9,p10 = GR_N_as_int,GR_Scratch
173 nop.f 0
174 nop.i 0
175 }
176 { .mmi
177 cmp.ge.unc p6, p0 = GR_N_as_int, GR_Big
178 cmp.le.unc p8, p0 = GR_N_as_int, GR_NBig
179 nop.i 0
180 };;
182 //
183 // Is N really an int, only for those non-int indefinites?
184 // Create exp bias.
185 //
186 { .mfi
187 add GR_N_Biased = GR_Bias,GR_N_as_int
188 (p9) fcmp.neq.unc.s1 p7,p0 = FR_Norm_N, FR_N_float_int
189 nop.i 0
190 };;
192 //
193 // Branch and return if N is not an int.
194 // Main path, create 2**N
195 //
196 { .mfi
197 setf.exp FR_Two_N = GR_N_Biased
198 nop.i 999
199 }
200 { .mfb
201 nop.m 0
202 (p7) frcpa f8,p11 = f0,f0
203 (p7) br.ret.spnt b0
204 };;
206 //
207 // Set denormal on denormal input x and denormal input N
208 //
209 { .mfi
210 nop.m 999
211 (p10)fcmp.ge.s1 p6,p8 = FR_Norm_N,f0
212 nop.i 0
213 };;
214 { .mfi
215 nop.m 999
216 fcmp.ge.s0 p0,p11 = FR_Floating_X,f0
217 nop.i 999
218 }
219 { .mfi
220 nop.m 999
221 fcmp.ge.s0 p12,p13 = FR_Floating_N,f0
222 nop.i 0
223 };;
225 //
226 // Adjust 2**N if N was very small or very large
227 //
229 { .mfi
230 nop.m 0
231 (p6) fma.s1 FR_Two_N = FR_Big,f1,f0
232 nop.i 0
233 }
234 { .mlx
235 nop.m 999
236 (p0) movl GR_Scratch = 0x00000000000303FF
237 };;
238 { .mfi
239 nop.m 0
240 (p8) fma.s1 FR_Two_N = FR_NBig,f1,f0
241 nop.i 0
242 }
243 { .mlx
244 nop.m 999
245 (p0) movl GR_Scratch1= 0x00000000000103FF
246 };;
248 // Set up necessary status fields
249 //
250 // S0 user supplied status
251 // S2 user supplied status + WRE + TD (Overflows)
252 // S3 user supplied status + FZ + TD (Underflows)
253 //
254 { .mfi
255 nop.m 999
256 (p0) fsetc.s3 0x7F,0x41
257 nop.i 999
258 }
259 { .mfi
260 nop.m 999
261 (p0) fsetc.s2 0x7F,0x42
262 nop.i 999
263 };;
265 //
266 // Do final operation
267 //
268 { .mfi
269 setf.exp FR_NBig = GR_Scratch
270 fma.d.s0 FR_Result = FR_Two_N,FR_Norm_X,f0
271 nop.i 999
272 }
273 { .mfi
274 nop.m 999
275 fma.d.s3 FR_Result3 = FR_Two_N,FR_Norm_X,f0
276 nop.i 999
277 };;
278 { .mfi
279 setf.exp FR_Big = GR_Scratch1
280 fma.d.s2 FR_Result2 = FR_Two_N,FR_Norm_X,f0
281 nop.i 999
282 };;
284 // Check for overflow or underflow.
285 //
286 // S0 user supplied status
287 // S2 user supplied status + WRE + TD (Overflow)
288 // S3 user supplied status + FZ + TD (Underflow)
289 //
290 //
291 // Restore s3
292 // Restore s2
293 //
294 { .mfi
295 nop.m 0
296 fsetc.s3 0x7F,0x40
297 nop.i 999
298 }
299 { .mfi
300 nop.m 0
301 fsetc.s2 0x7F,0x40
302 nop.i 999
303 };;
305 //
306 // Is the result zero?
307 //
308 { .mfi
309 nop.m 999
310 fclass.m.unc p6, p0 = FR_Result3, 0x007
311 nop.i 999
312 }
313 { .mfi
314 addl GR_Tag = 53, r0
315 fcmp.ge.unc.s1 p7, p8 = FR_Result2 , FR_Big
316 nop.i 0
317 };;
319 //
320 // Detect masked underflow - Tiny + Inexact Only
321 //
322 { .mfi
323 nop.m 999
324 (p6) fcmp.neq.unc.s1 p6, p0 = FR_Result , FR_Result2
325 nop.i 999
326 };;
328 //
329 // Is result bigger the allowed range?
330 // Branch out for underflow
331 //
332 { .mfb
333 (p6) addl GR_Tag = 54, r0
334 (p8) fcmp.le.unc.s1 p9, p10 = FR_Result2 , FR_NBig
335 (p6) br.cond.spnt L(SCALB_UNDERFLOW)
336 };;
338 //
339 // Branch out for overflow
340 //
341 { .mbb
342 nop.m 0
343 (p7) br.cond.spnt L(SCALB_OVERFLOW)
344 (p9) br.cond.spnt L(SCALB_OVERFLOW)
345 };;
347 //
348 // Return from main path.
349 //
350 { .mfb
351 nop.m 999
352 nop.f 0
353 br.ret.sptk b0;;
354 }
356 L(SCALB_NAN_INF_ZERO):
358 //
359 // Convert N to a fp integer
360 //
361 { .mfi
362 addl GR_Scratch = 1,r0
363 fcvt.fx.trunc.s1 FR_N_float_int = FR_Norm_N
364 nop.i 999
365 }
366 { .mfi
367 nop.m 0
368 fclass.m.unc p6,p0 = FR_Floating_N, 0xc3 //@snan | @qnan
369 nop.i 0
370 };;
371 { .mfi
372 nop.m 0
373 fclass.m.unc p7,p0 = FR_Floating_X, 0xc3 //@snan | @qnan
374 shl GR_Scratch = GR_Scratch,63
375 };;
376 { .mfi
377 nop.m 0
378 fclass.m.unc p8,p0 = FR_Floating_N, 0x21 // @inf
379 nop.i 0
380 }
381 { .mfi
382 nop.m 0
383 fclass.m.unc p9,p0 = FR_Floating_N, 0x22 // @-inf
384 nop.i 0
385 };;
387 //
388 // Either X or N is a Nan, return result and possible raise invalid.
389 //
390 { .mfb
391 nop.m 0
392 (p6) fma.d.s0 FR_Result = FR_Floating_N,FR_Floating_X,f0
393 (p6) br.ret.spnt b0
394 };;
395 { .mfb
396 getf.sig GR_N_as_int = FR_N_float_int
397 (p7) fma.d.s0 FR_Result = FR_Floating_N,FR_Floating_X,f0
398 (p7) br.ret.spnt b0
399 };;
401 //
402 // If N + Inf do something special
403 // For N = -Inf, create Int
404 //
405 { .mfb
406 nop.m 0
407 (p8) fma.d.s0 FR_Result = FR_Floating_X, FR_Floating_N,f0
408 (p8) br.ret.spnt b0
409 }
410 { .mfi
411 nop.m 0
412 (p9) fnma.d.s0 FR_Floating_N = FR_Floating_N, f1, f0
413 nop.i 0
414 };;
416 //
417 // If N==-Inf,return x/(-N)
418 //
419 { .mfb
420 nop.m 0
421 (p9) frcpa.s0 FR_Result,p6 = FR_Floating_X,FR_Floating_N
422 (p9) br.ret.spnt b0
423 };;
425 //
426 // Convert N_float_int to floating point value
427 //
428 { .mfi
429 cmp.ne.unc p9,p0 = GR_N_as_int,GR_Scratch
430 fcvt.xf FR_N_float_int = FR_N_float_int
431 nop.i 0
432 };;
434 //
435 // Is N an integer.
436 //
437 { .mfi
438 nop.m 0
439 (p9) fcmp.neq.unc.s1 p7,p0 = FR_Norm_N, FR_N_float_int
440 nop.i 0
441 };;
443 //
444 // If N not an int, return NaN and raise invalid.
445 //
446 { .mfb
447 nop.m 0
448 (p7) frcpa.s0 FR_Result,p6 = f0,f0
449 (p7) br.ret.spnt b0
450 };;
452 //
453 // Always return x in other path.
454 //
455 { .mfb
456 nop.m 0
457 fma.d.s0 FR_Result = FR_Floating_X,f1,f0
458 br.ret.sptk b0
459 };;
461 .endp scalb
462 ASM_SIZE_DIRECTIVE(scalb)
463 #ifdef _LIBC
464 ASM_SIZE_DIRECTIVE(__ieee754_scalb)
465 #endif
466 .proc __libm_error_region
467 __libm_error_region:
469 L(SCALB_OVERFLOW):
470 L(SCALB_UNDERFLOW):
472 //
473 // Get stack address of N
474 //
475 .prologue
476 { .mfi
477 add GR_Parameter_Y=-32,sp
478 nop.f 0
479 .save ar.pfs,GR_SAVE_PFS
480 mov GR_SAVE_PFS=ar.pfs
481 }
482 //
483 // Adjust sp
484 //
485 { .mfi
486 .fframe 64
487 add sp=-64,sp
488 nop.f 0
489 mov GR_SAVE_GP=gp
490 };;
492 //
493 // Store N on stack in correct position
494 // Locate the address of x on stack
495 //
496 { .mmi
497 stfd [GR_Parameter_Y] = FR_Norm_N,16
498 add GR_Parameter_X = 16,sp
499 .save b0, GR_SAVE_B0
500 mov GR_SAVE_B0=b0
501 };;
503 //
504 // Store x on the stack.
505 // Get address for result on stack.
506 //
507 .body
508 { .mib
509 stfd [GR_Parameter_X] = FR_Norm_X
510 add GR_Parameter_RESULT = 0,GR_Parameter_Y
511 nop.b 0
512 }
513 { .mib
514 stfd [GR_Parameter_Y] = FR_Result
515 add GR_Parameter_Y = -16,GR_Parameter_Y
516 br.call.sptk b0=__libm_error_support#
517 };;
519 //
520 // Get location of result on stack
521 //
522 { .mmi
523 nop.m 0
524 nop.m 0
525 add GR_Parameter_RESULT = 48,sp
526 };;
528 //
529 // Get the new result
530 //
531 { .mmi
532 ldfd FR_Result = [GR_Parameter_RESULT]
533 .restore sp
534 add sp = 64,sp
535 mov b0 = GR_SAVE_B0
536 };;
538 //
539 // Restore gp, ar.pfs and return
540 //
541 { .mib
542 mov gp = GR_SAVE_GP
543 mov ar.pfs = GR_SAVE_PFS
544 br.ret.sptk b0
545 };;
547 .endp __libm_error_region
548 ASM_SIZE_DIRECTIVE(__libm_error_region)
550 .type __libm_error_support#,@function
551 .global __libm_error_support#