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