1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 /* 22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 .file "powf.s" 30 31 / Note: 0^SNaN should not signal "invalid" but this implementation 32 / does because y is placed on the NPX stack. 33 34 / Special cases: 35 / 36 / x ** 0 is 1 37 / 1 ** y is 1 (C99) 38 / x ** NaN is NaN 39 / NaN ** y (except 0) is NaN 40 / x ** 1 is x 41 / +-(|x| > 1) ** +inf is +inf 42 / +-(|x| > 1) ** -inf is +0 43 / +-(|x| < 1) ** +inf is +0 44 / +-(|x| < 1) ** -inf is +inf 45 / (-1) ** +-inf is +1 (C99) 46 / +0 ** +y (except 0, NaN) is +0 47 / -0 ** +y (except 0, NaN, odd int) is +0 48 / +0 ** -y (except 0, NaN) is +inf (z flag) 49 / -0 ** -y (except 0, NaN, odd int) is +inf (z flag) 50 / -0 ** y (odd int) is - (+0 ** x) 51 / +inf ** +y (except 0, NaN) is +inf 52 / +inf ** -y (except 0, NaN) is +0 53 / -inf ** +-y (except 0, NaN) is -0 ** -+y (NO z flag) 54 / x ** -1 is 1/x 55 / x ** 2 is x*x 56 / -x ** y (an integer) is (-1)**(y) * (+x)**(y) 57 / x ** y (x negative & y not integer) is NaN (i flag) 58 59 #include "libm.h" 60 LIBM_ANSI_PRAGMA_WEAK(powf,function) 61 #include "libm_protos.h" 62 #include "xpg6.h" 63 64 .data 65 .align 4 66 negzero: 67 .float -0.0 68 half: 69 .float 0.5 70 one: 71 .float 1.0 72 negone: 73 .float -1.0 74 two: 75 .float 2.0 76 Snan: 77 .long 0x7f800001 78 pinfinity: 79 .long 0x7f800000 80 ninfinity: 81 .long 0xff800000 82 83 84 ENTRY(powf) 85 pushl %ebp 86 movl %esp,%ebp 87 PIC_SETUP(1) 88 89 flds 8(%ebp) / x 90 fxam / determine class of x 91 fnstsw %ax / store status in %ax 92 movb %ah,%dh / %dh <- condition code of x 93 94 flds 12(%ebp) / y , x 95 fxam / determine class of y 96 fnstsw %ax / store status in %ax 97 movb %ah,%dl / %dl <- condition code of y 98 99 call .pow_main /// LOCAL 100 PIC_WRAPUP 101 leave 102 ret 103 104 .pow_main: 105 / x ** 0 is 1 106 movb %dl,%cl 107 andb $0x45,%cl 108 cmpb $0x40,%cl / C3=1 C2=0 C1=? C0=0 when +-0 109 jne 1f 110 fstp %st(0) / x 111 fstp %st(0) / stack empty 112 fld1 / 1 113 ret 114 115 1: / y is not zero 116 PIC_G_LOAD(movzwl,__xpg6,eax) 117 andl $_C99SUSv3_pow_treats_Inf_as_an_even_int,%eax 118 cmpl $0,%eax 119 je 1f 120 121 / C99: 1 ** anything is 1 122 fld1 / 1, y, x 123 fucomp %st(2) / y, x 124 fnstsw %ax / store status in %ax 125 sahf / 80387 flags in %ax to 80386 flags 126 jp 1f / so that pow(NaN1,NaN2) returns NaN2 127 jne 1f 128 fstp %st(0) / x 129 ret 130 131 1: 132 / x ** NaN is NaN 133 movb %dl,%cl 134 andb $0x45,%cl 135 cmpb $0x01,%cl / C3=0 C2=0 C1=? C0=1 when +-NaN 136 jne 1f 137 fstp %st(1) / y 138 ret 139 140 1: / y is not NaN 141 / NaN ** y (except 0) is NaN 142 movb %dh,%cl 143 andb $0x45,%cl 144 cmpb $0x01,%cl / C3=0 C2=0 C1=? C0=1 when +-NaN 145 jne 1f 146 fstp %st(0) / x 147 ret 148 149 1: / x is not NaN 150 / x ** 1 is x 151 fcoms PIC_L(one) / y , x 152 fnstsw %ax / store status in %ax 153 sahf / 80387 flags in %ax to 80386 flags 154 jne 1f 155 fstp %st(0) / x 156 ret 157 158 1: / y is not 1 159 / +-(|x| > 1) ** +inf is +inf 160 / +-(|x| > 1) ** -inf is +0 161 / +-(|x| < 1) ** +inf is +0 162 / +-(|x| < 1) ** -inf is +inf 163 / +-(|x| = 1) ** +-inf is NaN 164 movb %dl,%cl 165 andb $0x47,%cl 166 cmpb $0x05,%cl / C3=0 C2=1 C1=0 C0=1 when +inf 167 je .yispinf 168 cmpb $0x07,%cl / C3=0 C2=1 C1=1 C0=1 when -inf 169 je .yisninf 170 171 / +0 ** +y (except 0, NaN) is +0 172 / -0 ** +y (except 0, NaN, odd int) is +0 173 / +0 ** -y (except 0, NaN) is +inf (z flag) 174 / -0 ** -y (except 0, NaN, odd int) is +inf (z flag) 175 / -0 ** y (odd int) is - (+0 ** x) 176 movb %dh,%cl 177 andb $0x47,%cl 178 cmpb $0x40,%cl / C3=1 C2=0 C1=0 C0=0 when +0 179 je .xispzero 180 cmpb $0x42,%cl / C3=1 C2=0 C1=1 C0=0 when -0 181 je .xisnzero 182 183 / +inf ** +y (except 0, NaN) is +inf 184 / +inf ** -y (except 0, NaN) is +0 185 / -inf ** +-y (except 0, NaN) is -0 ** -+y (NO z flag) 186 movb %dh,%cl 187 andb $0x47,%cl 188 cmpb $0x05,%cl / C3=0 C2=1 C1=0 C0=1 when +inf 189 je .xispinf 190 cmpb $0x07,%cl / C3=0 C2=1 C1=1 C0=1 when -inf 191 je .xisninf 192 193 / x ** -1 is 1/x 194 fcoms PIC_L(negone) / y , x 195 fnstsw %ax / store status in %ax 196 sahf / 80387 flags in %ax to 80386 flags 197 jne 1f 198 fld %st(1) / x , y , x 199 fdivrs PIC_L(one) / 1/x , y , x 200 jmp .signok / check for over/underflow 201 202 1: / y is not -1 203 / x ** 2 is square(x) 204 fcoms PIC_L(two) / y , x 205 fnstsw %ax / store status in %ax 206 sahf / 80387 flags in %ax to 80386 flags 207 jne 1f 208 fld %st(1) / x , y , x 209 fld %st(0) / x , x , y , x 210 fmulp / x^2 , y , x 211 jmp .signok / check for over/underflow 212 213 1: / y is not 2 214 / x ** 1/2 is sqrt(x) 215 fcoms PIC_L(half) / y , x 216 fnstsw %ax / store status in %ax 217 sahf / 80387 flags in %ax to 80386 flags 218 jne 1f 219 fld %st(1) / x , y , x 220 fsqrt / sqrt(x) , y , x 221 jmp .signok / check for over/underflow 222 223 1: / y is not 2 224 / make copies of x & y 225 fld %st(1) / x , y , x 226 fld %st(1) / y , x , y , x 227 228 / -x ** y (an integer) is (-1)**(y) * (+x)**(y) 229 / x ** y (x negative & y not integer) is NaN 230 movl $0,%ecx / track whether to flip sign of result 231 fld %st(1) / x , y , x , y , x 232 ftst / compare %st(0) with 0 233 fnstsw %ax / store status in %ax 234 sahf / 80387 flags in %ax to 80386 flags 235 fstp %st(0) / y , x , y , x 236 ja .merge / x > 0 237 / x < 0 238 call .y_is_int 239 cmpl $0,%ecx 240 jne 1f 241 / x < 0 & y != int so x**y = NaN (i flag) 242 fstp %st(0) / x , y , x 243 fstp %st(0) / y , x 244 fstp %st(0) / y , x 245 fstp %st(0) / y , x 246 fldz 247 fdiv %st,%st(0) / 0/0 248 ret 249 250 1: / x < 0 & y = int 251 fxch / x , y , y , x 252 fchs / px = -x , y , y , x 253 fxch / y , px , y , x 254 .merge: 255 / px > 0 256 fxch / px , y , y , x 257 258 / x**y = exp(y*ln(x)) 259 fyl2x / t=y*log2(px) , y , x 260 fld %st(0) / t , t , y , x 261 frndint / [t] , t , y , x 262 fxch / t , [t] , y , x 263 fucom 264 fnstsw %ax / store status in %ax 265 sahf / 80387 flags in %ax to 80386 flags 266 je 1f / px = int 267 fsub %st(1),%st / t-[t] , [t] , y , x 268 f2xm1 / 2**(t-[t])-1 , [t] , y , x 269 fadds PIC_L(one) / 2**(t-[t]) , [t] , y , x 270 fscale / 2**t = px**y , [t] , y , x 271 jmp 2f 272 1: 273 fstp %st(0) / t=[t] , y , x 274 fld1 / 1 , t , y , x 275 fscale / 1*2**t = x**y , t , y , x 276 2: 277 fstp %st(1) / x**y , y , x 278 cmpl $1,%ecx 279 jne .signok 280 fchs / change sign since x<0 & y=-int 281 .signok: 282 subl $4,%esp 283 fstps (%esp) / round to single precision 284 flds (%esp) / place result on NPX stack 285 addl $4,%esp 286 fstp %st(2) / y , x**y 287 fstp %st(0) / x**y 288 ret 289 290 / ------------------------------------------------------------------------ 291 292 .xispinf: 293 ftst / compare %st(0) with 0 294 fnstsw %ax / store status in %ax 295 sahf / 80387 flags in %ax to 80386 flags 296 ja .retpinf / y > 0 297 jmp .retpzero / y < 0 298 299 .xisninf: 300 / -inf ** +-y is -0 ** -+y 301 fchs / -y , x 302 flds PIC_L(negzero) / -0 , -y , x 303 fstp %st(2) / -y , -0 304 jmp .xisnzero 305 306 .yispinf: 307 fld %st(1) / x , y , x 308 fabs / |x| , y , x 309 fcomps PIC_L(one) / y , x 310 fnstsw %ax / store status in %ax 311 sahf / 80387 flags in %ax to 80386 flags 312 je .retponeorinvalid / x == -1 C99 313 ja .retpinf / |x| > 1 314 jmp .retpzero / |x| < 1 315 316 .yisninf: 317 fld %st(1) / x , y , x 318 fabs / |x| , y , x 319 fcomps PIC_L(one) / y , x 320 fnstsw %ax / store status in %ax 321 sahf / 80387 flags in %ax to 80386 flags 322 je .retponeorinvalid / x == -1 C99 323 ja .retpzero / |x| > 1 324 jmp .retpinf / |x| < 1 325 326 .xispzero: 327 / y cannot be 0 or NaN ; stack has y , x 328 ftst / compare %st(0) with 0 329 fnstsw %ax / store status in %ax 330 sahf / 80387 flags in %ax to 80386 flags 331 ja .retpzero / y > 0 332 / x = +0 & y < 0 so x**y = +inf 333 jmp .retpinfzflag / ret +inf & z flag 334 335 .xisnzero: 336 / y cannot be 0 or NaN ; stack has y , x 337 call .y_is_int 338 cmpl $1,%ecx 339 jne 1f / y is not an odd integer 340 / y is an odd integer 341 ftst / compare %st(0) with 0 342 fnstsw %ax / store status in %ax 343 sahf / 80387 flags in %ax to 80386 flags 344 ja .retnzero / y > 0 345 / x = -0 & y < 0 (odd int) return -inf (z flag) 346 / x = -inf & y != 0 or NaN return -inf (NO z flag) 347 movb %dh,%cl 348 andb $0x45,%cl 349 cmpb $0x05,%cl / C3=0 C2=1 C1=? C0=1 when +-inf 350 je 2f 351 fdiv %st,%st(1) / y / x, x (raise z flag) 352 2: 353 fstp %st(0) / x 354 fstp %st(0) / stack empty 355 flds PIC_L(ninfinity) / -inf 356 ret 357 358 1: / y is not an odd integer 359 ftst / compare %st(0) with 0 360 fnstsw %ax / store status in %ax 361 sahf / 80387 flags in %ax to 80386 flags 362 ja .retpzero / y > 0 363 / x = -0 & y < 0 (not odd int) return +inf (z flag) 364 / x = -inf & y not 0 or NaN return +inf (NO z flag) 365 movb %dh,%cl 366 andb $0x45,%cl 367 cmpb $0x05,%cl / C3=0 C2=1 C1=? C0=1 when +-inf 368 jne .retpinfzflag / ret +inf & divide-by-0 flag 369 jmp .retpinf / return +inf (NO z flag) 370 371 .retpzero: 372 fstp %st(0) / x 373 fstp %st(0) / stack empty 374 fldz / +0 375 ret 376 377 .retnzero: 378 fstp %st(0) / x 379 fstp %st(0) / stack empty 380 flds PIC_L(negzero) / -0 381 ret 382 383 .retponeorinvalid: 384 PIC_G_LOAD(movzwl,__xpg6,eax) 385 andl $_C99SUSv3_pow_treats_Inf_as_an_even_int,%eax 386 cmpl $0,%eax 387 je 1f 388 fstp %st(0) / x 389 fstp %st(0) / stack empty 390 fld1 / 1 391 ret 392 393 1: 394 fstp %st(0) / x 395 fstp %st(0) / stack empty 396 flds PIC_L(Snan) / Q NaN (i flag) 397 fwait 398 ret 399 400 .retpinf: 401 fstp %st(0) / x 402 fstp %st(0) / stack empty 403 flds PIC_L(pinfinity) / +inf 404 ret 405 406 .retpinfzflag: 407 fstp %st(0) / x 408 fstp %st(0) / stack empty 409 fldz 410 fdivrs PIC_L(one) / 1/0 411 ret 412 413 / Set %ecx to 2 if y is an even integer, 1 if y is an odd integer, 414 / 0 otherwise. Assume y is not zero. Do not raise inexact or modify 415 / %edx. 416 .y_is_int: 417 movl 12(%ebp),%eax 418 andl $0x7fffffff,%eax / |y| 419 cmpl $0x4b800000,%eax 420 jae 1f / |y| >= 2^24, an even int 421 cmpl $0x3f800000,%eax 422 jb 2f / |y| < 1, can't be an int 423 movl %eax,%ecx 424 sarl $23,%ecx 425 subl $150,%ecx 426 negl %ecx / 23 - unbiased exponent of y 427 bsfl %eax,%eax / index of least sig. 1 bit 428 cmpl %ecx,%eax 429 jb 2f 430 ja 1f 431 movl $1,%ecx 432 ret 433 1: 434 movl $2,%ecx 435 ret 436 2: 437 xorl %ecx,%ecx 438 ret 439 .align 4 440 SET_SIZE(powf)