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