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