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)