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  4
  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         .long   0x7f800001
  74 pinfinity:
  75         .long   0x7f800000
  76 ninfinity:
  77         .long   0xff800000
  78 
  79 
  80         ENTRY(powl)
  81         pushl   %ebp
  82         movl    %esp,%ebp
  83         PIC_SETUP(1)
  84 
  85         fldt    8(%ebp)                 / 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    20(%ebp)                / 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(movzwl,__xpg6,eax)
 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         fucomp  %st(2)                  / y, x
 120         fnstsw  %ax                     / store status in %ax
 121         sahf                            / 80387 flags in %ax to 80386 flags
 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         fcoms   PIC_L(one)              / y , x
 148         fnstsw  %ax                     / store status in %ax
 149         sahf                            / 80387 flags in %ax to 80386 flags
 150         jne     1f
 151         fstp    %st(0)                  / x
 152         ret
 153 
 154 1:      / y is not 1
 155         / +-(|x| > 1) **  +inf is +inf
 156         / +-(|x| > 1) **  -inf is +0
 157         / +-(|x| < 1) **  +inf is +0
 158         / +-(|x| < 1) **  -inf is +inf
 159         / +-(|x| = 1) ** +-inf is NaN
 160         movb    %dl,%cl
 161         andb    $0x47,%cl
 162         cmpb    $0x05,%cl               / C3=0 C2=1 C1=0 C0=1 when +inf
 163         je      .yispinf
 164         cmpb    $0x07,%cl               / C3=0 C2=1 C1=1 C0=1 when -inf
 165         je      .yisninf
 166 
 167         / +0 ** +y (except 0, NaN)              is +0
 168         / -0 ** +y (except 0, NaN, odd int)     is +0
 169         / +0 ** -y (except 0, NaN)              is +inf (z flag)
 170         / -0 ** -y (except 0, NaN, odd int)     is +inf (z flag)
 171         / -0 ** y (odd int)                     is - (+0 ** x)
 172         movb    %dh,%cl
 173         andb    $0x47,%cl
 174         cmpb    $0x40,%cl               / C3=1 C2=0 C1=0 C0=0 when +0
 175         je      .xispzero
 176         cmpb    $0x42,%cl               / C3=1 C2=0 C1=1 C0=0 when -0
 177         je      .xisnzero
 178 
 179         / +inf ** +y (except 0, NaN)    is +inf
 180         / +inf ** -y (except 0, NaN)    is +0
 181         / -inf ** +-y (except 0, NaN)   is -0 ** -+y (NO z flag)
 182         movb    %dh,%cl
 183         andb    $0x47,%cl
 184         cmpb    $0x05,%cl               / C3=0 C2=1 C1=0 C0=1 when +inf
 185         je      .xispinf
 186         cmpb    $0x07,%cl               / C3=0 C2=1 C1=1 C0=1 when -inf
 187         je      .xisninf
 188 
 189         / x ** -1 is 1/x
 190         fcoms   PIC_L(negone)           / y , x
 191         fnstsw  %ax                     / store status in %ax
 192         sahf                            / 80387 flags in %ax to 80386 flags
 193         jne     1f
 194         fld     %st(1)                  / x , y , x
 195         fdivrs  PIC_L(one)              / 1/x , y , x
 196         jmp     .signok                 / check for over/underflow
 197 
 198 1:      / y is not -1
 199         / x ** 2 is x*x
 200         fcoms   PIC_L(two)              / y , x
 201         fnstsw  %ax                     / store status in %ax
 202         sahf                            / 80387 flags in %ax to 80386 flags
 203         jne     1f
 204         fld     %st(1)                  / x , y , x
 205         fld     %st(0)                  / x , x , y , x
 206         fmulp                           / x^2 , y , x
 207         jmp     .signok                 / check for over/underflow
 208 
 209 1:      / y is not 2
 210         / x ** 1/2 is sqrt(x)
 211         fcoms   PIC_L(half)             / y , x
 212         fnstsw  %ax                     / store status in %ax
 213         sahf                            / 80387 flags in %ax to 80386 flags
 214         jne     1f
 215         fld     %st(1)                  / x , y , x
 216         fsqrt                           / sqrt(x) , y , x
 217         jmp     .signok                 / check for over/underflow
 218 
 219 1:      / y is not 1/2
 220         / make copies of x & y
 221         fld     %st(1)                  / x , y , x
 222         fld     %st(1)                  / y , x , y , x
 223 
 224         / -x ** y (an integer) is (-1)**(y) * (+x)**(y)
 225         / x ** y (x negative & y not integer) is  NaN
 226         movl    $0,%ecx                 / track whether to flip sign of result
 227         fld     %st(1)                  / x , y , x , y , x
 228         ftst                            / compare %st(0) with 0
 229         fnstsw  %ax                     / store status in %ax
 230         sahf                            / 80387 flags in %ax to 80386 flags
 231         fstp    %st(0)                  / y , x , y , x
 232         ja      .merge                  / x > 0
 233         / x < 0
 234         call    .y_is_int
 235         cmpl    $0,%ecx
 236         jne     1f
 237         / x < 0 & y != int so x**y = NaN (i flag)
 238         fstp    %st(0)                  / x , y , x
 239         fstp    %st(0)                  / y , x
 240         fstp    %st(0)                  / x
 241         fstp    %st(0)                  / stack empty
 242         fldz
 243         fdiv    %st,%st(0)              / 0/0
 244         ret
 245 
 246 1:      / x < 0 & y = int
 247         fxch                            / x , y , y , x
 248         fchs                            / px = -x , y , y , x
 249         fxch                            / y , px , y , x
 250 .merge:
 251         / px > 0
 252         fxch                            / px , y , y , x
 253 
 254         / x**y   =   exp(y*ln(x))
 255         fyl2x                           / t=y*log2(px) , y , x
 256         fld     %st(0)                  / t , t , y , x
 257         frndint                         / [t] , t , y , x
 258         fxch                            / t , [t] , y , x
 259         fucom
 260         fnstsw  %ax                     / store status in %ax
 261         sahf                            / 80387 flags in %ax to 80386 flags
 262         je      1f                      / t is integral
 263         fsub    %st(1),%st              / t-[t] , [t] , y , x
 264         f2xm1                           / 2**(t-[t])-1 , [t] , y , x
 265         fadds   PIC_L(one)              / 2**(t-[t]) , [t] , y , x
 266         fscale                          / 2**t = px**y , [t] , y , x
 267         jmp     2f
 268 1:
 269         fstp    %st(0)                  / t=[t] , y , x
 270         fld1                            / 1 , t , y , x
 271         fscale                          / 1*2**t = x**y , t , y , x
 272 2:
 273         fstp    %st(1)                  / x**y , y , x
 274         cmpl    $1,%ecx
 275         jne     .signok
 276         fchs                            / change sign since x<0 & y=-int
 277 .signok:
 278         fstp    %st(2)                  / y , x**y
 279         fstp    %st(0)                  / x**y
 280         ret
 281 
 282 / ------------------------------------------------------------------------
 283 
 284 .xispinf:
 285         ftst                            / compare %st(0) with 0
 286         fnstsw  %ax                     / store status in %ax
 287         sahf                            / 80387 flags in %ax to 80386 flags
 288         ja      .retpinf                / y > 0
 289         jmp     .retpzero               / y < 0
 290 
 291 .xisninf:
 292         / -inf ** +-y is -0 ** -+y
 293         fchs                            / -y , x
 294         flds    PIC_L(negzero)          / -0 , -y , x
 295         fstp    %st(2)                  / -y , -0
 296         jmp     .xisnzero
 297 
 298 .yispinf:
 299         fld     %st(1)                  / x , y , x
 300         fabs                            / |x| , y , x
 301         fcomps  PIC_L(one)              / y , x
 302         fnstsw  %ax                     / store status in %ax
 303         sahf                            / 80387 flags in %ax to 80386 flags
 304         je      .retponeorinvalid       / x == -1       C99
 305         ja      .retpinf                / |x| > 1
 306         jmp     .retpzero               / |x| < 1
 307 
 308 .yisninf:
 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      .retpzero               / |x| > 1
 316         jmp     .retpinf                / |x| < 1
 317 
 318 .xispzero:
 319         / y cannot be 0 or NaN ; stack has      y , x
 320         ftst                            / compare %st(0) with 0
 321         fnstsw  %ax                     / store status in %ax
 322         sahf                            / 80387 flags in %ax to 80386 flags
 323         ja      .retpzero               / y > 0
 324         / x = +0 & y < 0 so x**y = +inf
 325         jmp     .retpinfzflag           / ret +inf & z flag
 326 
 327 .xisnzero:
 328         / y cannot be 0 or NaN ; stack has      y , x
 329         call    .y_is_int
 330         cmpl    $1,%ecx
 331         jne     1f                      / y is not an odd integer
 332         / y is an odd integer
 333         ftst                            / compare %st(0) with 0
 334         fnstsw  %ax                     / store status in %ax
 335         sahf                            / 80387 flags in %ax to 80386 flags
 336         ja      .retnzero               / y > 0
 337         / x = -0 & y < 0 (odd int)       return -inf (z flag)
 338         / x = -inf & y != 0 or NaN  return -inf (NO z flag)
 339         movb    %dh,%cl
 340         andb    $0x45,%cl
 341         cmpb    $0x05,%cl               / C3=0 C2=1 C1=? C0=1 when +-inf
 342         je      2f
 343         fdiv    %st,%st(1)              / y / x, x (raise z flag)
 344 2:
 345         fstp    %st(0)                  / x
 346         fstp    %st(0)                  / stack empty
 347         flds    PIC_L(ninfinity)        / -inf
 348         ret
 349 
 350 1:      / y is not an odd integer
 351         ftst                            / compare %st(0) with 0
 352         fnstsw  %ax                     / store status in %ax
 353         sahf                            / 80387 flags in %ax to 80386 flags
 354         ja      .retpzero               / y > 0
 355         / x = -0 & y < 0 (not odd int)   return +inf (z flag)
 356         / x = -inf & y not 0 or NaN         return +inf (NO z flag)
 357         movb    %dh,%cl
 358         andb    $0x45,%cl
 359         cmpb    $0x05,%cl               / C3=0 C2=1 C1=? C0=1 when +-inf
 360         jne     .retpinfzflag           / ret +inf & divide-by-0 flag
 361         jmp     .retpinf                / return +inf (NO z flag)
 362 
 363 .retpzero:
 364         fstp    %st(0)                  / x
 365         fstp    %st(0)                  / stack empty
 366         fldz                            / +0
 367         ret
 368 
 369 .retnzero:
 370         fstp    %st(0)                  / x
 371         fstp    %st(0)                  / stack empty
 372         flds    PIC_L(negzero)          / -0
 373         ret
 374 
 375 .retponeorinvalid:
 376         PIC_G_LOAD(movzwl,__xpg6,eax)
 377         andl    $_C99SUSv3_pow_treats_Inf_as_an_even_int,%eax
 378         cmpl    $0,%eax
 379         je      1f
 380         fstp    %st(0)                  / x
 381         fstp    %st(0)                  / stack empty
 382         fld1                            / 1
 383         ret
 384 
 385 1:
 386         fstp    %st(0)                  / x
 387         fstp    %st(0)                  / stack empty
 388         flds    PIC_L(Snan)             / Q NaN (i flag)
 389         fwait
 390         ret
 391 
 392 .retpinf:
 393         fstp    %st(0)                  / x
 394         fstp    %st(0)                  / stack empty
 395         flds    PIC_L(pinfinity)        / +inf
 396         ret
 397 
 398 .retpinfzflag:
 399         fstp    %st(0)                  / x
 400         fstp    %st(0)                  / stack empty
 401         fldz
 402         fdivrs  PIC_L(one)              / 1/0
 403         ret
 404 
 405 / Set %ecx to 2 if y is an even integer, 1 if y is an odd integer,
 406 / 0 otherwise.  Assume y is not zero.  Do not raise inexact or modify
 407 / %edx.
 408 .y_is_int:
 409         movl    28(%ebp),%eax
 410         andl    $0x7fff,%eax            / exponent of y
 411         cmpl    $0x403f,%eax
 412         jae     1f                      / |y| >= 2^64, an even int
 413         cmpl    $0x3fff,%eax
 414         jb      2f                      / |y| < 1, can't be an int
 415         movl    %eax,%ecx
 416         subl    $0x403e,%ecx
 417         negl    %ecx                    / 63 - unbiased exponent of y
 418         movl    20(%ebp),%eax
 419         bsfl    %eax,%eax               / index of least sig. 1 bit
 420         jne     3f                      / jump if 1 bit found
 421         movl    24(%ebp),%eax
 422         bsfl    %eax,%eax
 423         addl    $32,%eax                / 32 + index of least sig. 1 bit
 424 3:
 425         cmpl    %ecx,%eax
 426         jb      2f
 427         ja      1f
 428         movl    $1,%ecx
 429         ret
 430 1:
 431         movl    $2,%ecx
 432         ret
 433 2:
 434         xorl    %ecx,%ecx
 435         ret
 436         .align  4
 437         SET_SIZE(powl)