sysdeps/powerpc/powerpc64/power8/fpu/e_expf.S

   1 /* Optimized expf().  PowerPC64/POWER8 version.
   2    Copyright (C) 2016-2019 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <http://www.gnu.org/licenses/>.  */
  18
  19 #include <sysdep.h>
  20
  21 /* Short algorithm description:
  22  *
  23  *  Let K = 64 (table size).
  24  *       e^x  = 2^(x/log(2)) = 2^n * T[j] * (1 + P(y))
  25  *  where:
  26  *       x = m*log(2)/K + y,    y in [0.0..log(2)/K]
  27  *       m = n*K + j,           m,n,j - signed integer, j in [0..K-1]
  28  *       values of 2^(j/K) are tabulated as T[j].
  29  *
  30  *       P(y) is a minimax polynomial approximation of expf(y)-1
  31  *       on small interval [0.0..log(2)/K].
  32  *
  33  *       P(y) = P3*y*y*y*y + P2*y*y*y + P1*y*y + P0*y, calculated as
  34  *       z = y*y;    P(y) = (P3*z + P1)*z + (P2*z + P0)*y
  35  *
  36  * Special cases:
  37  *  expf(NaN) = NaN
  38  *  expf(+INF) = +INF
  39  *  expf(-INF) = 0
  40  *  expf(x) = 1 for subnormals
  41  *  for finite argument, only expf(0)=1 is exact
  42  *  expf(x) overflows if x>88.7228317260742190
  43  *  expf(x) underflows if x<-103.972076416015620
  44  */
  45
  46 #define C1 0x42ad496b           /* Single precision 125*log(2).  */
  47 #define C2 0x31800000           /* Single precision 2^(-28).  */
  48 #define SP_INF 0x7f800000       /* Single precision Inf.  */
  49 #define SP_EXP_BIAS 0x1fc0      /* Single precision exponent bias.  */
  50
  51 #define DATA_OFFSET r9
  52
  53 /* Implements the function
  54
  55    float [fp1] expf (float [fp1] x)  */
  56
  57         .machine power8
  58 ENTRY (__ieee754_expf, 4)
  59         addis   DATA_OFFSET,r2,.Lanchor@toc@ha
  60         addi    DATA_OFFSET,DATA_OFFSET,.Lanchor@toc@l
  61
  62         xscvdpspn v0,v1
  63         mfvsrd  r8,v0           /* r8 = x  */
  64         lfd     fp2,(.KLN2-.Lanchor)(DATA_OFFSET)
  65         lfd     fp3,(.P2-.Lanchor)(DATA_OFFSET)
  66         rldicl  r3,r8,32,33     /* r3 = |x|  */
  67         lis     r4,C1@ha        /* r4 = 125*log(2)  */
  68         ori     r4,r4,C1@l
  69         cmpw    r3,r4
  70         lfd     fp5,(.P3-.Lanchor)(DATA_OFFSET)
  71         lfd     fp4,(.RS-.Lanchor)(DATA_OFFSET)
  72         fmadd   fp2,fp1,fp2,fp4 /* fp2 = x * K/log(2) + (2^23 + 2^22)  */
  73         bge     L(special_paths)        /* |x| >= 125*log(2) ?  */
  74
  75         lis     r4,C2@ha
  76         ori     r4,r4,C2@l
  77         cmpw    r3,r4
  78         blt     L(small_args)   /* |x| < 2^(-28) ?  */
  79
  80         /* Main path: here if 2^(-28) <= |x| < 125*log(2) */
  81         frsp    fp6,fp2
  82         xscvdpsp v2,v2
  83         mfvsrd  r8,v2
  84         mr      r3,r8                   /* r3 = m  */
  85         rldicl  r8,r8,32,58             /* r8 = j  */
  86         lfs     fp4,(.SP_RS-.Lanchor)(DATA_OFFSET)
  87         fsubs   fp2,fp6,fp4             /* fp2 = m = x * K/log(2)  */
  88         srdi    r3,r3,32
  89         clrrwi  r3,r3,6                 /* r3 = n  */
  90         lfd     fp6,(.NLN2K-.Lanchor)(DATA_OFFSET)
  91         fmadd   fp0,fp2,fp6,fp1         /* fp0 = y = x - m*log(2)/K  */
  92         fmul    fp2,fp0,fp0             /* fp2 = z = y^2  */
  93         lfd     fp4,(.P1-.Lanchor)(DATA_OFFSET)
  94         lfd     fp6,(.P0-.Lanchor)(DATA_OFFSET)
  95         lis     r4,SP_EXP_BIAS@ha
  96         ori     r4,r4,SP_EXP_BIAS@l
  97         add     r3,r3,r4
  98         rldic   r3,r3,49,1              /* r3 = 2^n  */
  99         fmadd   fp4,fp5,fp2,fp4         /* fp4 = P3 * z + P1  */
 100         fmadd   fp6,fp3,fp2,fp6         /* fp6 = P2 * z + P0  */
 101         mtvsrd  v1,r3
 102         xscvspdp v1,v1
 103         fmul    fp4,fp4,fp2             /* fp4 = (P3 * z + P1)*z  */
 104         fmadd   fp0,fp0,fp6,fp4         /* fp0 = P(y)  */
 105         sldi    r8,r8,3                 /* Access doublewords from T[j].  */
 106         addi    r6,DATA_OFFSET,(.Ttable-.Lanchor)
 107         lfdx    fp3,r6,r8
 108         fmadd   fp0,fp0,fp3,fp3         /* fp0 = T[j] * (1 + P(y))  */
 109         fmul    fp1,fp1,fp0             /* fp1 = 2^n * T[j] * (1 + P(y))  */
 110         frsp    fp1,fp1
 111         blr
 112
 113         .align  4
 114 /* x is either underflow, overflow, infinite or NaN.  */
 115 L(special_paths):
 116         srdi    r8,r8,32
 117         rlwinm  r8,r8,3,29,29           /* r8 = 0, if x positive.
 118                                            r8 = 4, otherwise.  */
 119         addi    r6,DATA_OFFSET,(.SPRANGE-.Lanchor)
 120         lwzx    r4,r6,r8                /* r4 = .SPRANGE[signbit(x)]  */
 121         cmpw    r3,r4
 122         /* |x| <= .SPRANGE[signbit(x)]  */
 123         ble     L(near_under_or_overflow)
 124
 125         lis     r4,SP_INF@ha
 126         ori     r4,r4,SP_INF@l
 127         cmpw    r3,r4
 128         bge     L(arg_inf_or_nan)       /* |x| > Infinite ?  */
 129
 130         addi    r6,DATA_OFFSET,(.SPLARGE_SMALL-.Lanchor)
 131         lfsx    fp1,r6,r8
 132         fmuls   fp1,fp1,fp1
 133         blr
 134
 135
 136         .align  4
 137 L(small_args):
 138         /* expf(x) = 1.0, where |x| < |2^(-28)|  */
 139         lfs     fp2,(.SPone-.Lanchor)(DATA_OFFSET)
 140         fadds   fp1,fp1,fp2
 141         blr
 142
 143
 144         .align  4
 145 L(arg_inf_or_nan:)
 146         bne     L(arg_nan)
 147
 148         /* expf(+INF) = +INF
 149            expf(-INF) = 0  */
 150         addi    r6,DATA_OFFSET,(.INF_ZERO-.Lanchor)
 151         lfsx    fp1,r6,r8
 152         blr
 153
 154
 155         .align  4
 156 L(arg_nan):
 157         /* expf(NaN) = NaN  */
 158         fadd    fp1,fp1,fp1
 159         frsp    fp1,fp1
 160         blr
 161
 162         .align  4
 163 L(near_under_or_overflow):
 164         frsp    fp6,fp2
 165         xscvdpsp v2,v2
 166         mfvsrd  r8,v2
 167         mr      r3,r8                   /* r3 = m  */
 168         rldicl  r8,r8,32,58             /* r8 = j  */
 169         lfs     fp4,(.SP_RS-.Lanchor)(DATA_OFFSET)
 170         fsubs   fp2,fp6,fp4             /* fp2 = m = x * K/log(2)  */
 171         srdi    r3,r3,32
 172         clrrwi  r3,r3,6                 /* r3 = n  */
 173         lfd     fp6,(.NLN2K-.Lanchor)(DATA_OFFSET)
 174         fmadd   fp0,fp2,fp6,fp1         /* fp0 = y = x - m*log(2)/K  */
 175         fmul    fp2,fp0,fp0             /* fp2 = z = y^2  */
 176         lfd     fp4,(.P1-.Lanchor)(DATA_OFFSET)
 177         lfd     fp6,(.P0-.Lanchor)(DATA_OFFSET)
 178         ld      r4,(.DP_EXP_BIAS-.Lanchor)(DATA_OFFSET)
 179         add     r3,r3,r4
 180         rldic   r3,r3,46,1              /* r3 = 2  */
 181         fmadd   fp4,fp5,fp2,fp4         /* fp4 = P3 * z + P1  */
 182         fmadd   fp6,fp3,fp2,fp6         /* fp6 = P2 * z + P0  */
 183         mtvsrd  v1,r3
 184         fmul    fp4,fp4,fp2             /* fp4 = (P3*z + P1)*z  */
 185         fmadd   fp0,fp0,fp6,fp4         /* fp0 = P(y)  */
 186         sldi    r8,r8,3                 /* Access doublewords from T[j].  */
 187         addi    r6,DATA_OFFSET,(.Ttable-.Lanchor)
 188         lfdx    fp3,r6,r8
 189         fmadd   fp0,fp0,fp3,fp3         /* fp0 = T[j] * (1 + T[j])  */
 190         fmul    fp1,fp1,fp0             /* fp1 = 2^n * T[j] * (1 + T[j])  */
 191         frsp    fp1,fp1
 192         blr
 193 END(__ieee754_expf)
 194
 195         .section .rodata, "a",@progbits
 196 .Lanchor:
 197         .balign 8
 198 /* Table T[j] = 2^(j/K).  Double precision.  */
 199 .Ttable:
 200         .8byte  0x3ff0000000000000
 201         .8byte  0x3ff02c9a3e778061
 202         .8byte  0x3ff059b0d3158574
 203         .8byte  0x3ff0874518759bc8
 204         .8byte  0x3ff0b5586cf9890f
 205         .8byte  0x3ff0e3ec32d3d1a2
 206         .8byte  0x3ff11301d0125b51
 207         .8byte  0x3ff1429aaea92de0
 208         .8byte  0x3ff172b83c7d517b
 209         .8byte  0x3ff1a35beb6fcb75
 210         .8byte  0x3ff1d4873168b9aa
 211         .8byte  0x3ff2063b88628cd6
 212         .8byte  0x3ff2387a6e756238
 213         .8byte  0x3ff26b4565e27cdd
 214         .8byte  0x3ff29e9df51fdee1
 215         .8byte  0x3ff2d285a6e4030b
 216         .8byte  0x3ff306fe0a31b715
 217         .8byte  0x3ff33c08b26416ff
 218         .8byte  0x3ff371a7373aa9cb
 219         .8byte  0x3ff3a7db34e59ff7
 220         .8byte  0x3ff3dea64c123422
 221         .8byte  0x3ff4160a21f72e2a
 222         .8byte  0x3ff44e086061892d
 223         .8byte  0x3ff486a2b5c13cd0
 224         .8byte  0x3ff4bfdad5362a27
 225         .8byte  0x3ff4f9b2769d2ca7
 226         .8byte  0x3ff5342b569d4f82
 227         .8byte  0x3ff56f4736b527da
 228         .8byte  0x3ff5ab07dd485429
 229         .8byte  0x3ff5e76f15ad2148
 230         .8byte  0x3ff6247eb03a5585
 231         .8byte  0x3ff6623882552225
 232         .8byte  0x3ff6a09e667f3bcd
 233         .8byte  0x3ff6dfb23c651a2f
 234         .8byte  0x3ff71f75e8ec5f74
 235         .8byte  0x3ff75feb564267c9
 236         .8byte  0x3ff7a11473eb0187
 237         .8byte  0x3ff7e2f336cf4e62
 238         .8byte  0x3ff82589994cce13
 239         .8byte  0x3ff868d99b4492ed
 240         .8byte  0x3ff8ace5422aa0db
 241         .8byte  0x3ff8f1ae99157736
 242         .8byte  0x3ff93737b0cdc5e5
 243         .8byte  0x3ff97d829fde4e50
 244         .8byte  0x3ff9c49182a3f090
 245         .8byte  0x3ffa0c667b5de565
 246         .8byte  0x3ffa5503b23e255d
 247         .8byte  0x3ffa9e6b5579fdbf
 248         .8byte  0x3ffae89f995ad3ad
 249         .8byte  0x3ffb33a2b84f15fb
 250         .8byte  0x3ffb7f76f2fb5e47
 251         .8byte  0x3ffbcc1e904bc1d2
 252         .8byte  0x3ffc199bdd85529c
 253         .8byte  0x3ffc67f12e57d14b
 254         .8byte  0x3ffcb720dcef9069
 255         .8byte  0x3ffd072d4a07897c
 256         .8byte  0x3ffd5818dcfba487
 257         .8byte  0x3ffda9e603db3285
 258         .8byte  0x3ffdfc97337b9b5f
 259         .8byte  0x3ffe502ee78b3ff6
 260         .8byte  0x3ffea4afa2a490da
 261         .8byte  0x3ffefa1bee615a27
 262         .8byte  0x3fff50765b6e4540
 263         .8byte  0x3fffa7c1819e90d8
 264
 265 .KLN2:
 266         .8byte  0x40571547652b82fe      /* Double precision K/log(2).  */
 267
 268 /* Double precision polynomial coefficients.  */
 269 .P0:
 270         .8byte  0x3fefffffffffe7c6
 271 .P1:
 272         .8byte  0x3fe00000008d6118
 273 .P2:
 274         .8byte  0x3fc55550da752d4f
 275 .P3:
 276         .8byte  0x3fa56420eb78fa85
 277
 278 .RS:
 279         .8byte  0x4168000000000000      /* Double precision 2^23 + 2^22.  */
 280 .NLN2K:
 281         .8byte  0xbf862e42fefa39ef      /* Double precision -log(2)/K.  */
 282 .DP_EXP_BIAS:
 283         .8byte  0x000000000000ffc0      /* Double precision exponent bias.  */
 284
 285         .balign 4
 286 .SPone:
 287         .4byte  0x3f800000      /* Single precision 1.0.  */
 288 .SP_RS:
 289         .4byte  0x4b400000      /* Single precision 2^23 + 2^22.  */
 290
 291 .SPRANGE: /* Single precision overflow/underflow bounds.  */
 292         .4byte  0x42b17217      /* if x>this bound, then result overflows.  */
 293         .4byte  0x42cff1b4      /* if x<this bound, then result underflows.  */
 294
 295 .SPLARGE_SMALL:
 296         .4byte  0x71800000      /* 2^100.  */
 297         .4byte  0x0d800000      /* 2^-100.  */
 298
 299 .INF_ZERO:
 300         .4byte  0x7f800000      /* Single precision Inf.  */
 301         .4byte  0               /* Single precision zero.  */
 302
 303 strong_alias (__ieee754_expf, __expf_finite)