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)