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