| blob | 3d48aab189d60246a9012fd6c3a24f910060fabb |
1 .file "scalb.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 // History
41 //==============================================================
42 // 02/02/00 Initial version
43 // 01/26/01 Scalb completely reworked and now standalone version
44 // 05/20/02 Cleaned up namespace and sf0 syntax
45 // 02/10/03 Reordered header: .section, .global, .proc, .align
46 // 08/06/03 Improved performance
47 //
48 // API
49 //==============================================================
50 // double = scalb (double x, double n)
51 // input floating point f8 and floating point f9
52 // output floating point f8
53 //
54 // int_type = 0 if int is 32 bits
55 // int_type = 1 if int is 64 bits
56 //
57 // Returns x* 2**n using an fma and detects overflow
58 // and underflow.
59 //
60 //
61 // Strategy:
62 // Compute biased exponent of result exp_Result = N + exp_X
63 // Break into ranges:
64 // exp_Result > 0x103fe -> Certain overflow
65 // exp_Result = 0x103fe -> Possible overflow
66 // 0x0fc01 <= exp_Result < 0x103fe -> No over/underflow (main path)
67 // 0x0fc01 - 52 <= exp_Result < 0x0fc01 -> Possible underflow
68 // exp_Result < 0x0fc01 - 52 -> Certain underflow
70 FR_Big = f6
71 FR_NBig = f7
72 FR_Floating_X = f8
73 FR_Result = f8
74 FR_Floating_N = f9
75 FR_Result2 = f9
76 FR_Result3 = f10
77 FR_Norm_X = f11
78 FR_Two_N = f12
79 FR_N_float_int = f13
80 FR_Norm_N = f14
82 GR_neg_ov_limit= r14
83 GR_big_exp = r14
84 GR_N_Biased = r15
85 GR_Big = r16
86 GR_exp_Result = r18
87 GR_pos_ov_limit= r19
88 GR_exp_sure_ou = r19
89 GR_Bias = r20
90 GR_N_as_int = r21
91 GR_signexp_X = r22
92 GR_exp_X = r23
93 GR_exp_mask = r24
94 GR_max_exp = r25
95 GR_min_exp = r26
96 GR_min_den_exp = r27
97 GR_Scratch = r28
98 GR_signexp_N = r29
99 GR_exp_N = r30
101 GR_SAVE_B0 = r32
102 GR_SAVE_GP = r33
103 GR_SAVE_PFS = r34
104 GR_Parameter_X = r35
105 GR_Parameter_Y = r36
106 GR_Parameter_RESULT = r37
107 GR_Tag = r38
109 .section .text
110 GLOBAL_IEEE754_ENTRY(scalb)
112 //
113 // Is x NAN, INF, ZERO, +-?
114 // Build the exponent Bias
115 //
116 { .mfi
117 getf.exp GR_signexp_N = FR_Floating_N // Get signexp of n
118 fclass.m p6,p0 = FR_Floating_X, 0xe7 // @snan | @qnan | @inf | @zero
119 mov GR_Bias = 0x0ffff
120 }
121 { .mfi
122 mov GR_Big = 35000 // If N this big then certain overflow
123 fcvt.fx.trunc.s1 FR_N_float_int = FR_Floating_N // Get N in significand
124 nop.i 0
125 }
126 ;;
128 { .mfi
129 getf.exp GR_signexp_X = FR_Floating_X // Get signexp of x
130 fclass.m p7,p0 = FR_Floating_N, 0x0b // Test for n=unorm
131 nop.i 0
132 }
133 //
134 // Normalize n
135 //
136 { .mfi
137 mov GR_exp_mask = 0x1ffff // Exponent mask
138 fnorm.s1 FR_Norm_N = FR_Floating_N
139 nop.i 0
140 }
141 ;;
143 //
144 // Is n NAN, INF, ZERO, +-?
145 //
146 { .mfi
147 mov GR_big_exp = 0x1003e // Exponent at which n is integer
148 fclass.m p9,p0 = FR_Floating_N, 0xe7 // @snan | @qnan | @inf | @zero
149 mov GR_max_exp = 0x103fe // Exponent of maximum double
150 }
151 //
152 // Normalize x
153 //
154 { .mfb
155 nop.m 0
156 fnorm.s1 FR_Norm_X = FR_Floating_X
157 (p7) br.cond.spnt SCALB_N_UNORM // Branch if n=unorm
158 }
159 ;;
161 SCALB_COMMON1:
162 // Main path continues. Also return here from u=unorm path.
163 // Handle special cases if x = Nan, Inf, Zero
164 { .mfb
165 nop.m 0
166 fcmp.lt.s1 p7,p0 = FR_Floating_N, f0 // Test N negative
167 (p6) br.cond.spnt SCALB_NAN_INF_ZERO
168 }
169 ;;
171 // Handle special cases if n = Nan, Inf, Zero
172 { .mfi
173 getf.sig GR_N_as_int = FR_N_float_int // Get n from significand
174 fclass.m p8,p0 = FR_Floating_X, 0x0b // Test for x=unorm
175 mov GR_exp_sure_ou = 0x1000e // Exp_N where x*2^N sure over/under
176 }
177 { .mfb
178 mov GR_min_exp = 0x0fc01 // Exponent of minimum double
179 fcvt.xf FR_N_float_int = FR_N_float_int // Convert N to FP integer
180 (p9) br.cond.spnt SCALB_NAN_INF_ZERO
181 }
182 ;;
184 { .mmi
185 and GR_exp_N = GR_exp_mask, GR_signexp_N // Get exponent of N
186 (p7) sub GR_Big = r0, GR_Big // Limit for N
187 nop.i 0
188 }
189 ;;
191 { .mib
192 cmp.lt p9,p0 = GR_exp_N, GR_big_exp // N possible non-integer?
193 cmp.ge p6,p0 = GR_exp_N, GR_exp_sure_ou // N certain over/under?
194 (p8) br.cond.spnt SCALB_X_UNORM // Branch if x=unorm
195 }
196 ;;
198 SCALB_COMMON2:
199 // Main path continues. Also return here from x=unorm path.
200 // Create biased exponent for 2**N
201 { .mmi
202 (p6) mov GR_N_as_int = GR_Big // Limit N
203 ;;
204 add GR_N_Biased = GR_Bias,GR_N_as_int
205 nop.i 0
206 }
207 ;;
209 { .mfi
210 setf.exp FR_Two_N = GR_N_Biased // Form 2**N
211 (p9) fcmp.neq.unc.s1 p9,p0 = FR_Norm_N, FR_N_float_int // Test if N an integer
212 and GR_exp_X = GR_exp_mask, GR_signexp_X // Get exponent of X
213 }
214 ;;
216 //
217 // Compute biased result exponent
218 // Branch if N is not an integer
219 //
220 { .mib
221 add GR_exp_Result = GR_exp_X, GR_N_as_int
222 mov GR_min_den_exp = 0x0fc01 - 52 // Exponent of min denorm dble
223 (p9) br.cond.spnt SCALB_N_NOT_INT
224 }
225 ;;
227 //
228 // Raise Denormal operand flag with compare
229 // Do final operation
230 //
231 { .mfi
232 cmp.lt p7,p6 = GR_exp_Result, GR_max_exp // Test no overflow
233 fcmp.ge.s0 p0,p11 = FR_Floating_X,FR_Floating_N // Dummy to set denorm
234 cmp.lt p9,p0 = GR_exp_Result, GR_min_den_exp // Test sure underflow
235 }
236 { .mfb
237 nop.m 0
238 fma.d.s0 FR_Result = FR_Two_N,FR_Norm_X,f0
239 (p9) br.cond.spnt SCALB_UNDERFLOW // Branch if certain underflow
240 }
241 ;;
243 { .mib
244 (p6) cmp.gt.unc p6,p8 = GR_exp_Result, GR_max_exp // Test sure overflow
245 (p7) cmp.ge.unc p7,p9 = GR_exp_Result, GR_min_exp // Test no over/underflow
246 (p7) br.ret.sptk b0 // Return from main path
247 }
248 ;;
250 { .bbb
251 (p6) br.cond.spnt SCALB_OVERFLOW // Branch if certain overflow
252 (p8) br.cond.spnt SCALB_POSSIBLE_OVERFLOW // Branch if possible overflow
253 (p9) br.cond.spnt SCALB_POSSIBLE_UNDERFLOW // Branch if possible underflow
254 }
255 ;;
257 // Here if possible underflow.
258 // Resulting exponent: 0x0fc01-52 <= exp_Result < 0x0fc01
259 SCALB_POSSIBLE_UNDERFLOW:
260 //
261 // Here if possible overflow.
262 // Resulting exponent: 0x103fe = exp_Result
263 SCALB_POSSIBLE_OVERFLOW:
265 // Set up necessary status fields
266 //
267 // S0 user supplied status
268 // S2 user supplied status + WRE + TD (Overflows)
269 // S3 user supplied status + FZ + TD (Underflows)
270 //
271 { .mfi
272 mov GR_pos_ov_limit = 0x103ff // Exponent for positive overflow
273 fsetc.s3 0x7F,0x41
274 nop.i 0
275 }
276 { .mfi
277 mov GR_neg_ov_limit = 0x303ff // Exponent for negative overflow
278 fsetc.s2 0x7F,0x42
279 nop.i 0
280 }
281 ;;
283 //
284 // Do final operation with s2 and s3
285 //
286 { .mfi
287 setf.exp FR_NBig = GR_neg_ov_limit
288 fma.d.s3 FR_Result3 = FR_Two_N,FR_Norm_X,f0
289 nop.i 0
290 }
291 { .mfi
292 setf.exp FR_Big = GR_pos_ov_limit
293 fma.d.s2 FR_Result2 = FR_Two_N,FR_Norm_X,f0
294 nop.i 0
295 }
296 ;;
298 // Check for overflow or underflow.
299 // Restore s3
300 // Restore s2
301 //
302 { .mfi
303 nop.m 0
304 fsetc.s3 0x7F,0x40
305 nop.i 0
306 }
307 { .mfi
308 nop.m 0
309 fsetc.s2 0x7F,0x40
310 nop.i 0
311 }
312 ;;
314 //
315 // Is the result zero?
316 //
317 { .mfi
318 nop.m 0
319 fclass.m p6, p0 = FR_Result3, 0x007
320 nop.i 0
321 }
322 { .mfi
323 nop.m 0
324 fcmp.ge.s1 p7, p8 = FR_Result2 , FR_Big
325 nop.i 0
326 }
327 ;;
329 //
330 // Detect masked underflow - Tiny + Inexact Only
331 //
332 { .mfi
333 nop.m 0
334 (p6) fcmp.neq.unc.s1 p6, p0 = FR_Result , FR_Result2
335 nop.i 0
336 }
337 ;;
339 //
340 // Is result bigger the allowed range?
341 // Branch out for underflow
342 //
343 { .mfb
344 nop.m 0
345 (p8) fcmp.le.unc.s1 p9, p10 = FR_Result2 , FR_NBig
346 (p6) br.cond.spnt SCALB_UNDERFLOW
347 }
348 ;;
350 //
351 // Branch out for overflow
352 //
353 { .bbb
354 (p7) br.cond.spnt SCALB_OVERFLOW
355 (p9) br.cond.spnt SCALB_OVERFLOW
356 br.ret.sptk b0 // Return from main path.
357 }
358 ;;
360 // Here if result overflows
361 SCALB_OVERFLOW:
362 { .mib
363 alloc r32=ar.pfs,3,0,4,0
364 addl GR_Tag = 53, r0 // Set error tag for overflow
365 br.cond.sptk __libm_error_region // Call error support for overflow
366 }
367 ;;
369 // Here if result underflows
370 SCALB_UNDERFLOW:
371 { .mib
372 alloc r32=ar.pfs,3,0,4,0
373 addl GR_Tag = 54, r0 // Set error tag for underflow
374 br.cond.sptk __libm_error_region // Call error support for underflow
375 }
376 ;;
378 SCALB_NAN_INF_ZERO:
380 //
381 // Before entry, N has been converted to a fp integer in significand of
382 // FR_N_float_int
383 //
384 // Convert N_float_int to floating point value
385 //
386 { .mfi
387 getf.sig GR_N_as_int = FR_N_float_int
388 fclass.m p6,p0 = FR_Floating_N, 0xc3 //@snan | @qnan
389 nop.i 0
390 }
391 { .mfi
392 addl GR_Scratch = 1,r0
393 fcvt.xf FR_N_float_int = FR_N_float_int
394 nop.i 0
395 }
396 ;;
398 { .mfi
399 nop.m 0
400 fclass.m p7,p0 = FR_Floating_X, 0xc3 //@snan | @qnan
401 shl GR_Scratch = GR_Scratch,63
402 }
403 ;;
405 { .mfi
406 nop.m 0
407 fclass.m p8,p0 = FR_Floating_N, 0x21 // @inf
408 nop.i 0
409 }
410 { .mfi
411 nop.m 0
412 fclass.m p9,p0 = FR_Floating_N, 0x22 // @-inf
413 nop.i 0
414 }
415 ;;
417 //
418 // Either X or N is a Nan, return result and possible raise invalid.
419 //
420 { .mfb
421 nop.m 0
422 (p6) fma.d.s0 FR_Result = FR_Floating_N,FR_Floating_X,f0
423 (p6) br.ret.spnt b0
424 }
425 ;;
427 { .mfb
428 nop.m 0
429 (p7) fma.d.s0 FR_Result = FR_Floating_N,FR_Floating_X,f0
430 (p7) br.ret.spnt b0
431 }
432 ;;
434 //
435 // If N + Inf do something special
436 // For N = -Inf, create Int
437 //
438 { .mfb
439 nop.m 0
440 (p8) fma.d.s0 FR_Result = FR_Floating_X, FR_Floating_N,f0
441 (p8) br.ret.spnt b0
442 }
443 { .mfi
444 nop.m 0
445 (p9) fnma.d.s0 FR_Floating_N = FR_Floating_N, f1, f0
446 nop.i 0
447 }
448 ;;
450 //
451 // If N==-Inf,return x/(-N)
452 //
453 { .mfb
454 cmp.ne p7,p0 = GR_N_as_int,GR_Scratch
455 (p9) frcpa.s0 FR_Result,p0 = FR_Floating_X,FR_Floating_N
456 (p9) br.ret.spnt b0
457 }
458 ;;
460 //
461 // Is N an integer.
462 //
463 { .mfi
464 nop.m 0
465 (p7) fcmp.neq.unc.s1 p7,p0 = FR_Norm_N, FR_N_float_int
466 nop.i 0
467 }
468 ;;
470 //
471 // If N not an int, return NaN and raise invalid.
472 //
473 { .mfb
474 nop.m 0
475 (p7) frcpa.s0 FR_Result,p0 = f0,f0
476 (p7) br.ret.spnt b0
477 }
478 ;;
480 //
481 // Always return x in other path.
482 //
483 { .mfb
484 nop.m 0
485 fma.d.s0 FR_Result = FR_Floating_X,f1,f0
486 br.ret.sptk b0
487 }
488 ;;
490 // Here if n not int
491 // Return NaN and raise invalid.
492 SCALB_N_NOT_INT:
493 { .mfb
494 nop.m 0
495 frcpa.s0 FR_Result,p0 = f0,f0
496 br.ret.sptk b0
497 }
498 ;;
500 // Here if n=unorm
501 SCALB_N_UNORM:
502 { .mfb
503 getf.exp GR_signexp_N = FR_Norm_N // Get signexp of normalized n
504 fcvt.fx.trunc.s1 FR_N_float_int = FR_Norm_N // Get N in significand
505 br.cond.sptk SCALB_COMMON1 // Return to main path
506 }
507 ;;
509 // Here if x=unorm
510 SCALB_X_UNORM:
511 { .mib
512 getf.exp GR_signexp_X = FR_Norm_X // Get signexp of normalized x
513 nop.i 0
514 br.cond.sptk SCALB_COMMON2 // Return to main path
515 }
516 ;;
518 GLOBAL_IEEE754_END(scalb)
519 LOCAL_LIBM_ENTRY(__libm_error_region)
521 //
522 // Get stack address of N
523 //
524 .prologue
525 { .mfi
526 add GR_Parameter_Y=-32,sp
527 nop.f 0
528 .save ar.pfs,GR_SAVE_PFS
529 mov GR_SAVE_PFS=ar.pfs
530 }
531 //
532 // Adjust sp
533 //
534 { .mfi
535 .fframe 64
536 add sp=-64,sp
537 nop.f 0
538 mov GR_SAVE_GP=gp
539 };;
541 //
542 // Store N on stack in correct position
543 // Locate the address of x on stack
544 //
545 { .mmi
546 stfd [GR_Parameter_Y] = FR_Norm_N,16
547 add GR_Parameter_X = 16,sp
548 .save b0, GR_SAVE_B0
549 mov GR_SAVE_B0=b0
550 };;
552 //
553 // Store x on the stack.
554 // Get address for result on stack.
555 //
556 .body
557 { .mib
558 stfd [GR_Parameter_X] = FR_Norm_X
559 add GR_Parameter_RESULT = 0,GR_Parameter_Y
560 nop.b 0
561 }
562 { .mib
563 stfd [GR_Parameter_Y] = FR_Result
564 add GR_Parameter_Y = -16,GR_Parameter_Y
565 br.call.sptk b0=__libm_error_support#
566 };;
568 //
569 // Get location of result on stack
570 //
571 { .mmi
572 add GR_Parameter_RESULT = 48,sp
573 nop.m 0
574 nop.i 0
575 };;
577 //
578 // Get the new result
579 //
580 { .mmi
581 ldfd FR_Result = [GR_Parameter_RESULT]
582 .restore sp
583 add sp = 64,sp
584 mov b0 = GR_SAVE_B0
585 };;
587 //
588 // Restore gp, ar.pfs and return
589 //
590 { .mib
591 mov gp = GR_SAVE_GP
592 mov ar.pfs = GR_SAVE_PFS
593 br.ret.sptk b0
594 };;
596 LOCAL_LIBM_END(__libm_error_region)
598 .type __libm_error_support#,@function
599 .global __libm_error_support#