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

   1 /* Optimized sinf().  PowerPC64/POWER8 version.
   2    Copyright (C) 2016 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 #define _ERRNO_H        1
  21 #include <bits/errno.h>
  22
  23 #define FRAMESIZE (FRAME_MIN_SIZE+16)
  24
  25 #define FLOAT_EXPONENT_SHIFT    23
  26 #define FLOAT_EXPONENT_BIAS     127
  27 #define INTEGER_BITS            3
  28
  29 #define PI_4            0x3f490fdb      /* PI/4 */
  30 #define NINEPI_4        0x40e231d6      /* 9 * PI/4 */
  31 #define TWO_PN5         0x3d000000      /* 2^-5 */
  32 #define TWO_PN27        0x32000000      /* 2^-27 */
  33 #define INFINITY        0x7f800000
  34 #define TWO_P23         0x4b000000      /* 2^27 */
  35 #define FX_FRACTION_1_28 0x9249250      /* 0x100000000 / 28 + 1 */
  36
  37         /* Implements the function
  38
  39            float [fp1] sinf (float [fp1] x)  */
  40
  41         .machine power8
  42 EALIGN(__sinf, 4, 0)
  43         addis   r9,r2,L(anchor)@toc@ha
  44         addi    r9,r9,L(anchor)@toc@l
  45
  46         lis     r4,PI_4@h
  47         ori     r4,r4,PI_4@l
  48
  49         xscvdpspn v0,v1
  50         mfvsrd  r8,v0
  51         rldicl  r3,r8,32,33             /* Remove sign bit.  */
  52
  53         cmpw    r3,r4
  54         bge     L(greater_or_equal_pio4)
  55
  56         lis     r4,TWO_PN5@h
  57         ori     r4,r4,TWO_PN5@l
  58
  59         cmpw    r3,r4
  60         blt     L(less_2pn5)
  61
  62         /* Chebyshev polynomial of the form:
  63          * x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
  64
  65         lfd     fp9,(L(S0)-L(anchor))(r9)
  66         lfd     fp10,(L(S1)-L(anchor))(r9)
  67         lfd     fp11,(L(S2)-L(anchor))(r9)
  68         lfd     fp12,(L(S3)-L(anchor))(r9)
  69         lfd     fp13,(L(S4)-L(anchor))(r9)
  70
  71         fmul    fp2,fp1,fp1             /* x^2 */
  72         fmul    fp3,fp2,fp1             /* x^3 */
  73
  74         fmadd   fp4,fp2,fp13,fp12       /* S3+x^2*S4 */
  75         fmadd   fp4,fp2,fp4,fp11        /* S2+x^2*(S3+x^2*S4) */
  76         fmadd   fp4,fp2,fp4,fp10        /* S1+x^2*(S2+x^2*(S3+x^2*S4)) */
  77         fmadd   fp4,fp2,fp4,fp9         /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))) */
  78         fmadd   fp1,fp3,fp4,fp1         /* x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))) */
  79         frsp    fp1,fp1                 /* Round to single precision.  */
  80
  81         blr
  82
  83         .balign 16
  84 L(greater_or_equal_pio4):
  85         lis     r4,NINEPI_4@h
  86         ori     r4,r4,NINEPI_4@l
  87         cmpw    r3,r4
  88         bge     L(greater_or_equal_9pio4)
  89
  90         /* Calculate quotient of |x|/(PI/4).  */
  91         lfd     fp2,(L(invpio4)-L(anchor))(r9)
  92         fabs    fp1,fp1                 /* |x| */
  93         fmul    fp2,fp1,fp2             /* |x|/(PI/4) */
  94         fctiduz fp2,fp2
  95         mfvsrd  r3,v2                   /* n = |x| mod PI/4 */
  96
  97         /* Now use that quotient to find |x| mod (PI/2).  */
  98         addi    r7,r3,1
  99         rldicr  r5,r7,2,60              /* ((n+1) >> 1) << 3 */
 100         addi    r6,r9,(L(pio2_table)-L(anchor))
 101         lfdx    fp4,r5,r6
 102         fsub    fp1,fp1,fp4
 103
 104         .balign 16
 105 L(reduced):
 106         /* Now we are in the range -PI/4 to PI/4.  */
 107
 108         /* Work out if we are in a positive or negative primary interval.  */
 109         rldicl  r4,r7,62,63             /* ((n+1) >> 2) & 1 */
 110
 111         /* We are operating on |x|, so we need to add back the original
 112            sign.  */
 113         rldicl  r8,r8,33,63             /* (x >> 31) & 1, ie the sign bit.  */
 114         xor     r4,r4,r8                /* 0 if result should be positive,
 115                                            1 if negative.  */
 116
 117         /* Load a 1.0 or -1.0.  */
 118         addi    r5,r9,(L(ones)-L(anchor))
 119         sldi    r4,r4,3
 120         lfdx    fp0,r4,r5
 121
 122         /* Are we in the primary interval of sin or cos?  */
 123         andi.   r4,r7,0x2
 124         bne     L(cos)
 125
 126         /* Chebyshev polynomial of the form:
 127            x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
 128
 129         lfd     fp9,(L(S0)-L(anchor))(r9)
 130         lfd     fp10,(L(S1)-L(anchor))(r9)
 131         lfd     fp11,(L(S2)-L(anchor))(r9)
 132         lfd     fp12,(L(S3)-L(anchor))(r9)
 133         lfd     fp13,(L(S4)-L(anchor))(r9)
 134
 135         fmul    fp2,fp1,fp1             /* x^2 */
 136         fmul    fp3,fp2,fp1             /* x^3 */
 137
 138         fmadd   fp4,fp2,fp13,fp12       /* S3+x^2*S4 */
 139         fmadd   fp4,fp2,fp4,fp11        /* S2+x^2*(S3+x^2*S4) */
 140         fmadd   fp4,fp2,fp4,fp10        /* S1+x^2*(S2+x^2*(S3+x^2*S4)) */
 141         fmadd   fp4,fp2,fp4,fp9         /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))) */
 142         fmadd   fp4,fp3,fp4,fp1         /* x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))) */
 143         fmul    fp4,fp4,fp0             /* Add in the sign.  */
 144         frsp    fp1,fp4                 /* Round to single precision.  */
 145
 146         blr
 147
 148         .balign 16
 149 L(cos):
 150         /* Chebyshev polynomial of the form:
 151            1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).  */
 152
 153         lfd     fp9,(L(C0)-L(anchor))(r9)
 154         lfd     fp10,(L(C1)-L(anchor))(r9)
 155         lfd     fp11,(L(C2)-L(anchor))(r9)
 156         lfd     fp12,(L(C3)-L(anchor))(r9)
 157         lfd     fp13,(L(C4)-L(anchor))(r9)
 158
 159         fmul    fp2,fp1,fp1             /* x^2 */
 160         lfd     fp3,(L(DPone)-L(anchor))(r9)
 161
 162         fmadd   fp4,fp2,fp13,fp12       /* C3+x^2*C4 */
 163         fmadd   fp4,fp2,fp4,fp11        /* C2+x^2*(C3+x^2*C4) */
 164         fmadd   fp4,fp2,fp4,fp10        /* C1+x^2*(C2+x^2*(C3+x^2*C4)) */
 165         fmadd   fp4,fp2,fp4,fp9         /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))) */
 166         fmadd   fp4,fp2,fp4,fp3         /* 1.0 + x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))) */
 167         fmul    fp4,fp4,fp0             /* Add in the sign.  */
 168         frsp    fp1,fp4                 /* Round to single precision.  */
 169
 170         blr
 171
 172         .balign 16
 173 L(greater_or_equal_9pio4):
 174         lis     r4,INFINITY@h
 175         ori     r4,r4,INFINITY@l
 176         cmpw    r3,r4
 177         bge     L(inf_or_nan)
 178
 179         lis     r4,TWO_P23@h
 180         ori     r4,r4,TWO_P23@l
 181         cmpw    r3,r4
 182         bge     L(greater_or_equal_2p23)
 183
 184         fabs    fp1,fp1                 /* |x| */
 185
 186         /* Calculate quotient of |x|/(PI/4).  */
 187         lfd     fp2,(L(invpio4)-L(anchor))(r9)
 188
 189         lfd     fp3,(L(DPone)-L(anchor))(r9)
 190         lfd     fp4,(L(DPhalf)-L(anchor))(r9)
 191         fmul    fp2,fp1,fp2             /* |x|/(PI/4) */
 192         friz    fp2,fp2                 /* n = floor(|x|/(PI/4)) */
 193
 194         /* Calculate (n + 1) / 2.  */
 195         fadd    fp2,fp2,fp3             /* n + 1 */
 196         fmul    fp3,fp2,fp4             /* (n + 1) / 2 */
 197         friz    fp3,fp3
 198
 199         lfd     fp4,(L(pio2hi)-L(anchor))(r9)
 200         lfd     fp5,(L(pio2lo)-L(anchor))(r9)
 201
 202         fmul    fp6,fp4,fp3
 203         fadd    fp6,fp6,fp1
 204         fmadd   fp1,fp5,fp3,fp6
 205
 206         fctiduz fp2,fp2
 207         mfvsrd  r7,v2                   /* n + 1 */
 208
 209         b       L(reduced)
 210
 211         .balign 16
 212 L(inf_or_nan):
 213         bne     L(skip_errno_setting)   /* Is a NAN?  */
 214
 215         /* We delayed the creation of the stack frame, as well as the saving of
 216            the link register, because only at this point, we are sure that
 217            doing so is actually needed.  */
 218
 219         stfd    fp1,-8(r1)
 220
 221         /* Save the link register.  */
 222         mflr    r0
 223         std     r0,16(r1)
 224         cfi_offset(lr, 16)
 225
 226         /* Create the stack frame.  */
 227         stdu    r1,-FRAMESIZE(r1)
 228         cfi_adjust_cfa_offset(FRAMESIZE)
 229
 230         bl      JUMPTARGET(__errno_location)
 231         nop
 232
 233         /* Restore the stack frame.  */
 234         addi    r1,r1,FRAMESIZE
 235         cfi_adjust_cfa_offset(-FRAMESIZE)
 236         /* Restore the link register.  */
 237         ld      r0,16(r1)
 238         mtlr    r0
 239
 240         lfd     fp1,-8(r1)
 241
 242         /* errno = EDOM */
 243         li      r4,EDOM
 244         stw     r4,0(r3)
 245
 246 L(skip_errno_setting):
 247         fsub    fp1,fp1,fp1             /* x - x */
 248         blr
 249
 250         .balign 16
 251 L(greater_or_equal_2p23):
 252         fabs    fp1,fp1
 253
 254         srwi    r4,r3,FLOAT_EXPONENT_SHIFT
 255         subi    r4,r4,FLOAT_EXPONENT_BIAS
 256
 257         /* We reduce the input modulo pi/4, so we need 3 bits of integer
 258            to determine where in 2*pi we are. Index into our array
 259            accordingly.  */
 260         addi r4,r4,INTEGER_BITS
 261
 262         /* To avoid an expensive divide, for the range we care about (0 - 127)
 263            we can transform x/28 into:
 264
 265            x/28 = (x * ((0x100000000 / 28) + 1)) >> 32
 266
 267            mulhwu returns the top 32 bits of the 64 bit result, doing the
 268            shift for us in the same instruction. The top 32 bits are undefined,
 269            so we have to mask them.  */
 270
 271         lis     r6,FX_FRACTION_1_28@h
 272         ori     r6,r6,FX_FRACTION_1_28@l
 273         mulhwu  r5,r4,r6
 274         clrldi  r5,r5,32
 275
 276         /* Get our pointer into the invpio4_table array.  */
 277         sldi    r4,r5,3
 278         addi    r6,r9,(L(invpio4_table)-L(anchor))
 279         add     r4,r4,r6
 280
 281         lfd     fp2,0(r4)
 282         lfd     fp3,8(r4)
 283         lfd     fp4,16(r4)
 284         lfd     fp5,24(r4)
 285
 286         fmul    fp6,fp2,fp1
 287         fmul    fp7,fp3,fp1
 288         fmul    fp8,fp4,fp1
 289         fmul    fp9,fp5,fp1
 290
 291         /* Mask off larger integer bits in highest double word that we don't
 292            care about to avoid losing precision when combining with smaller
 293            values.  */
 294         fctiduz fp10,fp6
 295         mfvsrd  r7,v10
 296         rldicr  r7,r7,0,(63-INTEGER_BITS)
 297         mtvsrd  v10,r7
 298         fcfidu  fp10,fp10               /* Integer bits.  */
 299
 300         fsub    fp6,fp6,fp10            /* highest -= integer bits */
 301
 302         /* Work out the integer component, rounded down. Use the top two
 303            limbs for this.  */
 304         fadd    fp10,fp6,fp7            /* highest + higher */
 305
 306         fctiduz fp10,fp10
 307         mfvsrd  r7,v10
 308         andi.   r0,r7,1
 309         fcfidu  fp10,fp10
 310
 311         /* Subtract integer component from highest limb.  */
 312         fsub    fp12,fp6,fp10
 313
 314         beq     L(even_integer)
 315
 316         /* Our integer component is odd, so we are in the -PI/4 to 0 primary
 317            region. We need to shift our result down by PI/4, and to do this
 318            in the mod (4/PI) space we simply subtract 1.  */
 319         lfd     fp11,(L(DPone)-L(anchor))(r9)
 320         fsub    fp12,fp12,fp11
 321
 322         /* Now add up all the limbs in order.  */
 323         fadd    fp12,fp12,fp7
 324         fadd    fp12,fp12,fp8
 325         fadd    fp12,fp12,fp9
 326
 327         /* And finally multiply by pi/4.  */
 328         lfd     fp13,(L(pio4)-L(anchor))(r9)
 329         fmul    fp1,fp12,fp13
 330
 331         addi    r7,r7,1
 332         b       L(reduced)
 333
 334 L(even_integer):
 335         lfd     fp11,(L(DPone)-L(anchor))(r9)
 336
 337         /* Now add up all the limbs in order.  */
 338         fadd    fp12,fp12,fp7
 339         fadd    fp12,r12,fp8
 340         fadd    fp12,r12,fp9
 341
 342         /* We need to check if the addition of all the limbs resulted in us
 343            overflowing 1.0.  */
 344         fcmpu   0,fp12,fp11
 345         bgt     L(greater_than_one)
 346
 347         /* And finally multiply by pi/4.  */
 348         lfd     fp13,(L(pio4)-L(anchor))(r9)
 349         fmul    fp1,fp12,fp13
 350
 351         addi    r7,r7,1
 352         b       L(reduced)
 353
 354 L(greater_than_one):
 355         /* We did overflow 1.0 when adding up all the limbs. Add 1.0 to our
 356            integer, and subtract 1.0 from our result. Since that makes the
 357            integer component odd, we need to subtract another 1.0 as
 358            explained above.  */
 359         addi    r7,r7,1
 360
 361         lfd     fp11,(L(DPtwo)-L(anchor))(r9)
 362         fsub    fp12,fp12,fp11
 363
 364         /* And finally multiply by pi/4.  */
 365         lfd     fp13,(L(pio4)-L(anchor))(r9)
 366         fmul    fp1,fp12,fp13
 367
 368         addi    r7,r7,1
 369         b       L(reduced)
 370
 371         .balign 16
 372 L(less_2pn5):
 373         lis     r4,TWO_PN27@h
 374         ori     r4,r4,TWO_PN27@l
 375
 376         cmpw    r3,r4
 377         blt     L(less_2pn27)
 378
 379         /* A simpler Chebyshev approximation is close enough for this range:
 380            x+x^3*(SS0+x^2*SS1).  */
 381
 382         lfd     fp10,(L(SS0)-L(anchor))(r9)
 383         lfd     fp11,(L(SS1)-L(anchor))(r9)
 384
 385         fmul    fp2,fp1,fp1             /* x^2 */
 386         fmul    fp3,fp2,fp1             /* x^3 */
 387
 388         fmadd   fp4,fp2,fp11,fp10       /* SS0+x^2*SS1 */
 389         fmadd   fp1,fp3,fp4,fp1         /* x+x^3*(SS0+x^2*SS1) */
 390
 391         frsp    fp1,fp1                 /* Round to single precision.  */
 392
 393         blr
 394
 395         .balign 16
 396 L(less_2pn27):
 397         cmpwi   r3,0
 398         beq     L(zero)
 399
 400         /* Handle some special cases:
 401
 402            sinf(subnormal) raises inexact/underflow
 403            sinf(min_normalized) raises inexact/underflow
 404            sinf(normalized) raises inexact.  */
 405
 406         lfd     fp2,(L(small)-L(anchor))(r9)
 407
 408         fmul    fp2,fp1,fp2             /* x * small */
 409         fsub    fp1,fp1,fp2             /* x - x * small */
 410
 411         frsp    fp1,fp1
 412
 413         blr
 414
 415         .balign 16
 416 L(zero):
 417         blr
 418
 419 END (__sinf)
 420
 421         .section .rodata, "a"
 422
 423         .balign 8
 424
 425 L(anchor):
 426
 427         /* Chebyshev constants for sin, range -PI/4 - PI/4.  */
 428 L(S0):  .8byte  0xbfc5555555551cd9
 429 L(S1):  .8byte  0x3f81111110c2688b
 430 L(S2):  .8byte  0xbf2a019f8b4bd1f9
 431 L(S3):  .8byte  0x3ec71d7264e6b5b4
 432 L(S4):  .8byte  0xbe5a947e1674b58a
 433
 434         /* Chebyshev constants for sin, range 2^-27 - 2^-5.  */
 435 L(SS0): .8byte  0xbfc555555543d49d
 436 L(SS1): .8byte  0x3f8110f475cec8c5
 437
 438         /* Chebyshev constants for cos, range -PI/4 - PI/4.  */
 439 L(C0):  .8byte  0xbfdffffffffe98ae
 440 L(C1):  .8byte  0x3fa55555545c50c7
 441 L(C2):  .8byte  0xbf56c16b348b6874
 442 L(C3):  .8byte  0x3efa00eb9ac43cc0
 443 L(C4):  .8byte  0xbe923c97dd8844d7
 444
 445 L(invpio2):
 446         .8byte  0x3fe45f306dc9c883      /* 2/PI */
 447
 448 L(invpio4):
 449         .8byte  0x3ff45f306dc9c883      /* 4/PI */
 450
 451 L(invpio4_table):
 452         .8byte  0x0000000000000000
 453         .8byte  0x3ff45f306c000000
 454         .8byte  0x3e3c9c882a000000
 455         .8byte  0x3c54fe13a8000000
 456         .8byte  0x3aaf47d4d0000000
 457         .8byte  0x38fbb81b6c000000
 458         .8byte  0x3714acc9e0000000
 459         .8byte  0x3560e4107c000000
 460         .8byte  0x33bca2c756000000
 461         .8byte  0x31fbd778ac000000
 462         .8byte  0x300b7246e0000000
 463         .8byte  0x2e5d2126e8000000
 464         .8byte  0x2c97003248000000
 465         .8byte  0x2ad77504e8000000
 466         .8byte  0x290921cfe0000000
 467         .8byte  0x274deb1cb0000000
 468         .8byte  0x25829a73e0000000
 469         .8byte  0x23fd1046be000000
 470         .8byte  0x2224baed10000000
 471         .8byte  0x20709d338e000000
 472         .8byte  0x1e535a2f80000000
 473         .8byte  0x1cef904e64000000
 474         .8byte  0x1b0d639830000000
 475         .8byte  0x1964ce7d24000000
 476         .8byte  0x17b908bf16000000
 477
 478 L(pio4):
 479         .8byte  0x3fe921fb54442d18      /* PI/4 */
 480
 481 /* PI/2 as a sum of two doubles. We only use 32 bits of the upper limb
 482    to avoid losing significant bits when multiplying with up to
 483    (2^22)/(pi/2).  */
 484 L(pio2hi):
 485         .8byte  0xbff921fb54400000
 486
 487 L(pio2lo):
 488         .8byte  0xbdd0b4611a626332
 489
 490 L(pio2_table):
 491         .8byte  0
 492         .8byte  0x3ff921fb54442d18      /* 1 * PI/2 */
 493         .8byte  0x400921fb54442d18      /* 2 * PI/2 */
 494         .8byte  0x4012d97c7f3321d2      /* 3 * PI/2 */
 495         .8byte  0x401921fb54442d18      /* 4 * PI/2 */
 496         .8byte  0x401f6a7a2955385e      /* 5 * PI/2 */
 497         .8byte  0x4022d97c7f3321d2      /* 6 * PI/2 */
 498         .8byte  0x4025fdbbe9bba775      /* 7 * PI/2 */
 499         .8byte  0x402921fb54442d18      /* 8 * PI/2 */
 500         .8byte  0x402c463abeccb2bb      /* 9 * PI/2 */
 501         .8byte  0x402f6a7a2955385e      /* 10 * PI/2 */
 502
 503 L(small):
 504         .8byte  0x3cd0000000000000      /* 2^-50 */
 505
 506 L(ones):
 507         .8byte  0x3ff0000000000000      /* +1.0 */
 508         .8byte  0xbff0000000000000      /* -1.0 */
 509
 510 L(DPhalf):
 511         .8byte  0x3fe0000000000000      /* 0.5 */
 512
 513 L(DPone):
 514         .8byte  0x3ff0000000000000      /* 1.0 */
 515
 516 L(DPtwo):
 517         .8byte  0x4000000000000000      /* 2.0 */
 518
 519 weak_alias(__sinf, sinf)