sysdeps/x86_64/fpu/multiarch/svml_d_atan22_core_sse4.S

   1 /* Function atan2 vectorized with SSE4.
   2    Copyright (C) 2021-2023 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    https://www.gnu.org/licenses/.  */
  18
  19 /*
  20  * ALGORITHM DESCRIPTION:
  21  *      For    0.0    <= x <=  7.0/16.0: atan(x) = atan(0.0) + atan(s), where s=(x-0.0)/(1.0+0.0*x)
  22  *      For  7.0/16.0 <= x <= 11.0/16.0: atan(x) = atan(0.5) + atan(s), where s=(x-0.5)/(1.0+0.5*x)
  23  *      For 11.0/16.0 <= x <= 19.0/16.0: atan(x) = atan(1.0) + atan(s), where s=(x-1.0)/(1.0+1.0*x)
  24  *      For 19.0/16.0 <= x <= 39.0/16.0: atan(x) = atan(1.5) + atan(s), where s=(x-1.5)/(1.0+1.5*x)
  25  *      For 39.0/16.0 <= x <=    inf   : atan(x) = atan(inf) + atan(s), where s=-1.0/x
  26  *      Where atan(s) ~= s+s^3*Poly11(s^2) on interval |s|<7.0/0.16.
  27  *
  28  *
  29  */
  30
  31 /* Offsets for data table __svml_datan2_data_internal
  32  */
  33 #define dPI                             0
  34 #define dPIO2                           16
  35 #define dA19                            32
  36 #define dA18                            48
  37 #define dA17                            64
  38 #define dA16                            80
  39 #define dA15                            96
  40 #define dA14                            112
  41 #define dA13                            128
  42 #define dA12                            144
  43 #define dA11                            160
  44 #define dA10                            176
  45 #define dA09                            192
  46 #define dA08                            208
  47 #define dA07                            224
  48 #define dA06                            240
  49 #define dA05                            256
  50 #define dA04                            272
  51 #define dA03                            288
  52 #define dA02                            304
  53 #define dA01                            320
  54 #define dA00                            336
  55 #define dSIGN_MASK                      352
  56 #define iCHK_WORK_SUB                   368
  57 #define iCHK_WORK_CMP                   384
  58 #define dABS_MASK                       400
  59 #define dZERO                           416
  60
  61 #include <sysdep.h>
  62
  63         .section .text.sse4, "ax", @progbits
  64 ENTRY(_ZGVbN2vv_atan2_sse4)
  65         subq    $88, %rsp
  66         cfi_def_cfa_offset(96)
  67         movaps  %xmm1, %xmm11
  68
  69         /*
  70          * #define NO_VECTOR_ZERO_ATAN2_ARGS
  71          *  Declarations
  72          * Variables
  73          * Constants
  74          *  The end of declarations
  75          *  Implementation
  76          * Get r0~=1/B
  77          * Cannot be replaced by VQRCP(D, dR0, dB);
  78          * Argument Absolute values
  79          */
  80         movups  dABS_MASK+__svml_datan2_data_internal(%rip), %xmm1
  81         movaps  %xmm0, %xmm10
  82         movaps  %xmm1, %xmm9
  83         andps   %xmm10, %xmm1
  84         andps   %xmm11, %xmm9
  85         movaps  %xmm1, %xmm4
  86         cmpnltpd %xmm9, %xmm4
  87
  88         /* Argument signs */
  89         movups  dSIGN_MASK+__svml_datan2_data_internal(%rip), %xmm5
  90         movaps  %xmm4, %xmm0
  91         movaps  %xmm5, %xmm8
  92         movaps  %xmm5, %xmm7
  93
  94         /*
  95          * 1) If y<x then a= y, b=x, PIO2=0
  96          * 2) If y>x then a=-x, b=y, PIO2=Pi/2
  97          */
  98         orps    %xmm9, %xmm5
  99         andnps  %xmm1, %xmm0
 100         andps   %xmm4, %xmm5
 101         andps   %xmm11, %xmm8
 102         movups  dPIO2+__svml_datan2_data_internal(%rip), %xmm6
 103         orps    %xmm5, %xmm0
 104         movaps  %xmm4, %xmm5
 105         andps   %xmm4, %xmm6
 106         andnps  %xmm9, %xmm5
 107         andps   %xmm1, %xmm4
 108         orps    %xmm4, %xmm5
 109         andps   %xmm10, %xmm7
 110         divpd   %xmm5, %xmm0
 111         movq    iCHK_WORK_SUB+__svml_datan2_data_internal(%rip), %xmm2
 112         xorl    %edx, %edx
 113
 114         /* Check if y and x are on main path. */
 115         pshufd  $221, %xmm9, %xmm3
 116         xorl    %eax, %eax
 117         pshufd  $221, %xmm1, %xmm13
 118         psubd   %xmm2, %xmm3
 119         psubd   %xmm2, %xmm13
 120         movdqa  %xmm3, %xmm4
 121         movq    iCHK_WORK_CMP+__svml_datan2_data_internal(%rip), %xmm12
 122         movdqa  %xmm13, %xmm14
 123         pcmpgtd %xmm12, %xmm4
 124         pcmpeqd %xmm12, %xmm3
 125         pcmpgtd %xmm12, %xmm14
 126         pcmpeqd %xmm12, %xmm13
 127
 128         /* Polynomial. */
 129         movaps  %xmm0, %xmm12
 130         por     %xmm3, %xmm4
 131         mulpd   %xmm0, %xmm12
 132
 133         /* P = A19*R2 + A18 */
 134         movups  dA19+__svml_datan2_data_internal(%rip), %xmm15
 135         movaps  %xmm11, %xmm2
 136         mulpd   %xmm12, %xmm15
 137         addpd   dA18+__svml_datan2_data_internal(%rip), %xmm15
 138
 139         /* P = P*R2 + A17 */
 140         mulpd   %xmm12, %xmm15
 141         addpd   dA17+__svml_datan2_data_internal(%rip), %xmm15
 142
 143         /* P = P*R2 + A16 */
 144         mulpd   %xmm12, %xmm15
 145         addpd   dA16+__svml_datan2_data_internal(%rip), %xmm15
 146
 147         /* P = P*R2 + A15 */
 148         mulpd   %xmm12, %xmm15
 149         addpd   dA15+__svml_datan2_data_internal(%rip), %xmm15
 150
 151         /* P = P*R2 + A14 */
 152         mulpd   %xmm12, %xmm15
 153         addpd   dA14+__svml_datan2_data_internal(%rip), %xmm15
 154
 155         /* P = P*R2 + A13 */
 156         mulpd   %xmm12, %xmm15
 157         addpd   dA13+__svml_datan2_data_internal(%rip), %xmm15
 158
 159         /* P = P*R2 + A12 */
 160         mulpd   %xmm12, %xmm15
 161         addpd   dA12+__svml_datan2_data_internal(%rip), %xmm15
 162
 163         /* P = P*R2 + A11 */
 164         mulpd   %xmm12, %xmm15
 165         addpd   dA11+__svml_datan2_data_internal(%rip), %xmm15
 166
 167         /* P = P*R2 + A10 */
 168         mulpd   %xmm12, %xmm15
 169         addpd   dA10+__svml_datan2_data_internal(%rip), %xmm15
 170
 171         /* P = P*R2 + A09 */
 172         mulpd   %xmm12, %xmm15
 173         addpd   dA09+__svml_datan2_data_internal(%rip), %xmm15
 174
 175         /* P = P*R2 + A08 */
 176         mulpd   %xmm12, %xmm15
 177         addpd   dA08+__svml_datan2_data_internal(%rip), %xmm15
 178
 179         /* P = P*R2 + A07 */
 180         mulpd   %xmm12, %xmm15
 181         addpd   dA07+__svml_datan2_data_internal(%rip), %xmm15
 182
 183         /* P = P*R2 + A06 */
 184         mulpd   %xmm12, %xmm15
 185         addpd   dA06+__svml_datan2_data_internal(%rip), %xmm15
 186
 187         /* P = P*R2 + A05 */
 188         mulpd   %xmm12, %xmm15
 189         addpd   dA05+__svml_datan2_data_internal(%rip), %xmm15
 190
 191         /* P = P*R2 + A04 */
 192         mulpd   %xmm12, %xmm15
 193         addpd   dA04+__svml_datan2_data_internal(%rip), %xmm15
 194
 195         /* P = P*R2 + A03 */
 196         mulpd   %xmm12, %xmm15
 197         addpd   dA03+__svml_datan2_data_internal(%rip), %xmm15
 198
 199         /* P = P*R2 + A02 */
 200         mulpd   %xmm12, %xmm15
 201         addpd   dA02+__svml_datan2_data_internal(%rip), %xmm15
 202
 203         /* P = P*R2 + A01 */
 204         mulpd   %xmm12, %xmm15
 205         addpd   dA01+__svml_datan2_data_internal(%rip), %xmm15
 206
 207         /* P = P*R2 */
 208         mulpd   %xmm15, %xmm12
 209
 210         /*
 211          * Reconstruction.
 212          * dP=(R+R*dP) + dPIO2
 213          */
 214         mulpd   %xmm0, %xmm12
 215         addpd   %xmm12, %xmm0
 216
 217         /* if x<0, dPI = Pi, else dPI =0 */
 218         movups  dZERO+__svml_datan2_data_internal(%rip), %xmm3
 219         por     %xmm13, %xmm14
 220         cmplepd %xmm3, %xmm2
 221         addpd   %xmm6, %xmm0
 222         andps   __svml_datan2_data_internal(%rip), %xmm2
 223         orps    %xmm8, %xmm0
 224         addpd   %xmm2, %xmm0
 225         por     %xmm14, %xmm4
 226         orps    %xmm7, %xmm0
 227         movmskps %xmm4, %ecx
 228
 229         /*  Special branch for fast (vector) processing of zero arguments  */
 230         testb   $3, %cl
 231
 232         /* Go to auxilary branch */
 233         jne     L(AUX_BRANCH)
 234         # LOE rbx rbp r12 r13 r14 r15 eax edx xmm0 xmm1 xmm2 xmm3 xmm4 xmm5 xmm6 xmm7 xmm8 xmm9 xmm10 xmm11
 235
 236         /* Return from auxilary branch
 237          * for out of main path inputs
 238          */
 239
 240 L(AUX_BRANCH_RETURN):
 241         /*
 242          *  Special branch for fast (vector) processing of zero arguments
 243          *  The end of implementation
 244          */
 245         testl   %edx, %edx
 246
 247         /* Go to special inputs processing branch */
 248         jne     L(SPECIAL_VALUES_BRANCH)
 249         # LOE rbx rbp r12 r13 r14 r15 eax edx xmm0 xmm10 xmm11
 250
 251         /* Restore registers
 252          * and exit the function
 253          */
 254
 255 L(EXIT):
 256         addq    $88, %rsp
 257         cfi_def_cfa_offset(8)
 258         ret
 259         cfi_def_cfa_offset(96)
 260
 261         /* Branch to process
 262          * special inputs
 263          */
 264
 265 L(SPECIAL_VALUES_BRANCH):
 266         movups  %xmm10, 32(%rsp)
 267         movups  %xmm11, 48(%rsp)
 268         movups  %xmm0, 64(%rsp)
 269         # LOE rbx rbp r12 r13 r14 r15 eax edx xmm0
 270
 271         movq    %r12, 16(%rsp)
 272         cfi_offset(12, -80)
 273         movl    %eax, %r12d
 274         movq    %r13, 8(%rsp)
 275         cfi_offset(13, -88)
 276         movl    %edx, %r13d
 277         movq    %r14, (%rsp)
 278         cfi_offset(14, -96)
 279         # LOE rbx rbp r15 r12d r13d
 280
 281         /* Range mask
 282          * bits check
 283          */
 284
 285 L(RANGEMASK_CHECK):
 286         btl     %r12d, %r13d
 287
 288         /* Call scalar math function */
 289         jc      L(SCALAR_MATH_CALL)
 290         # LOE rbx rbp r15 r12d r13d
 291
 292         /* Special inputs
 293          * processing loop
 294          */
 295
 296 L(SPECIAL_VALUES_LOOP):
 297         incl    %r12d
 298         cmpl    $2, %r12d
 299
 300         /* Check bits in range mask */
 301         jl      L(RANGEMASK_CHECK)
 302         # LOE rbx rbp r15 r12d r13d
 303
 304         movq    16(%rsp), %r12
 305         cfi_restore(12)
 306         movq    8(%rsp), %r13
 307         cfi_restore(13)
 308         movq    (%rsp), %r14
 309         cfi_restore(14)
 310         movups  64(%rsp), %xmm0
 311
 312         /* Go to exit */
 313         jmp     L(EXIT)
 314         cfi_offset(12, -80)
 315         cfi_offset(13, -88)
 316         cfi_offset(14, -96)
 317         # LOE rbx rbp r12 r13 r14 r15 xmm0
 318
 319         /* Scalar math fucntion call
 320          * to process special input
 321          */
 322
 323 L(SCALAR_MATH_CALL):
 324         movl    %r12d, %r14d
 325         movsd   32(%rsp, %r14, 8), %xmm0
 326         movsd   48(%rsp, %r14, 8), %xmm1
 327         call    atan2@PLT
 328         # LOE rbx rbp r14 r15 r12d r13d xmm0
 329
 330         movsd   %xmm0, 64(%rsp, %r14, 8)
 331
 332         /* Process special inputs in loop */
 333         jmp     L(SPECIAL_VALUES_LOOP)
 334         cfi_restore(12)
 335         cfi_restore(13)
 336         cfi_restore(14)
 337         # LOE rbx rbp r15 r12d r13d
 338
 339         /* Auxilary branch
 340          * for out of main path inputs
 341          */
 342
 343 L(AUX_BRANCH):
 344         /* Check if both X & Y are not NaNs:  iXYnotNAN */
 345         movaps  %xmm11, %xmm13
 346         movaps  %xmm10, %xmm12
 347         cmpordpd %xmm11, %xmm13
 348         cmpordpd %xmm10, %xmm12
 349
 350         /* Check if at least on of Y or Y is zero: iAXAYZERO */
 351         cmpeqpd %xmm3, %xmm9
 352         cmpeqpd %xmm3, %xmm1
 353
 354         /*
 355          *  Path for zero arguments (at least one of both)
 356          * Check if both args are zeros (den. is zero)
 357          */
 358         cmpeqpd %xmm3, %xmm5
 359         andps   %xmm12, %xmm13
 360         orps    %xmm1, %xmm9
 361         pshufd  $221, %xmm9, %xmm1
 362         pshufd  $221, %xmm13, %xmm9
 363
 364         /* Check if at least on of Y or Y is zero and not NaN: iAXAYZEROnotNAN */
 365         pand    %xmm9, %xmm1
 366
 367         /* Exclude from previous callout mask zero (and not NaN) arguments */
 368         movdqa  %xmm1, %xmm14
 369         pandn   %xmm4, %xmm14
 370
 371         /* Set sPIO2 to zero if den. is zero */
 372         movaps  %xmm5, %xmm4
 373         andnps  %xmm6, %xmm4
 374         andps   %xmm3, %xmm5
 375
 376         /* Res = sign(Y)*(X<0)?(PIO2+PI):PIO2 */
 377         pshufd  $221, %xmm3, %xmm3
 378         orps    %xmm5, %xmm4
 379         pshufd  $221, %xmm11, %xmm5
 380         orps    %xmm8, %xmm4
 381         pcmpgtd %xmm5, %xmm3
 382         pshufd  $80, %xmm3, %xmm6
 383         andps   %xmm2, %xmm6
 384         addpd   %xmm6, %xmm4
 385
 386         /* Go to callout */
 387         movmskps %xmm14, %edx
 388
 389         /* Merge results from main and spec path */
 390         pshufd  $80, %xmm1, %xmm2
 391         orps    %xmm7, %xmm4
 392         movdqa  %xmm2, %xmm7
 393         andps   %xmm2, %xmm4
 394         andnps  %xmm0, %xmm7
 395         andl    $3, %edx
 396         movaps  %xmm7, %xmm0
 397         orps    %xmm4, %xmm0
 398
 399         /* Return to main vector processing path */
 400         jmp     L(AUX_BRANCH_RETURN)
 401         # LOE rbx rbp r12 r13 r14 r15 eax edx xmm0 xmm10 xmm11
 402 END(_ZGVbN2vv_atan2_sse4)
 403
 404         .section .rodata, "a"
 405         .align  16
 406
 407 #ifdef __svml_datan2_data_internal_typedef
 408 typedef unsigned int VUINT32;
 409 typedef struct {
 410         __declspec(align(16)) VUINT32 dPI[2][2];
 411         __declspec(align(16)) VUINT32 dPIO2[2][2];
 412         __declspec(align(16)) VUINT32 dA19[2][2];
 413         __declspec(align(16)) VUINT32 dA18[2][2];
 414         __declspec(align(16)) VUINT32 dA17[2][2];
 415         __declspec(align(16)) VUINT32 dA16[2][2];
 416         __declspec(align(16)) VUINT32 dA15[2][2];
 417         __declspec(align(16)) VUINT32 dA14[2][2];
 418         __declspec(align(16)) VUINT32 dA13[2][2];
 419         __declspec(align(16)) VUINT32 dA12[2][2];
 420         __declspec(align(16)) VUINT32 dA11[2][2];
 421         __declspec(align(16)) VUINT32 dA10[2][2];
 422         __declspec(align(16)) VUINT32 dA09[2][2];
 423         __declspec(align(16)) VUINT32 dA08[2][2];
 424         __declspec(align(16)) VUINT32 dA07[2][2];
 425         __declspec(align(16)) VUINT32 dA06[2][2];
 426         __declspec(align(16)) VUINT32 dA05[2][2];
 427         __declspec(align(16)) VUINT32 dA04[2][2];
 428         __declspec(align(16)) VUINT32 dA03[2][2];
 429         __declspec(align(16)) VUINT32 dA02[2][2];
 430         __declspec(align(16)) VUINT32 dA01[2][2];
 431         __declspec(align(16)) VUINT32 dA00[2][2];
 432         __declspec(align(16)) VUINT32 dSIGN_MASK[2][2];
 433         __declspec(align(16)) VUINT32 iCHK_WORK_SUB[4][1];
 434         __declspec(align(16)) VUINT32 iCHK_WORK_CMP[4][1];
 435         __declspec(align(16)) VUINT32 dABS_MASK[2][2];
 436         __declspec(align(16)) VUINT32 dZERO[2][2];
 437 } __svml_datan2_data_internal;
 438 #endif
 439 __svml_datan2_data_internal:
 440         .quad   0x400921FB54442D18, 0x400921FB54442D18 // dPI
 441         .align  16
 442         .quad   0x3FF921FB54442D18, 0x3FF921FB54442D18 // dPIO2
 443         .align  16
 444         .quad   0xBEF4FDB537ABC7A3, 0xBEF4FDB537ABC7A3 // dA19
 445         .align  16
 446         .quad   0x3F2CED0A36665209, 0x3F2CED0A36665209 // dA18
 447         .align  16
 448         .quad   0xBF52E67C93954C23, 0xBF52E67C93954C23 // dA17
 449         .align  16
 450         .quad   0x3F6F5A1DAE82AFB3, 0x3F6F5A1DAE82AFB3 // dA16
 451         .align  16
 452         .quad   0xBF82B2EC618E4BAD, 0xBF82B2EC618E4BAD // dA15
 453         .align  16
 454         .quad   0x3F914F4C661116A5, 0x3F914F4C661116A5 // dA14
 455         .align  16
 456         .quad   0xBF9A5E83B081F69C, 0xBF9A5E83B081F69C // dA13
 457         .align  16
 458         .quad   0x3FA169980CB6AD4F, 0x3FA169980CB6AD4F // dA12
 459         .align  16
 460         .quad   0xBFA4EFA2E563C1BC, 0xBFA4EFA2E563C1BC // dA11
 461         .align  16
 462         .quad   0x3FA7EC0FBC50683B, 0x3FA7EC0FBC50683B // dA10
 463         .align  16
 464         .quad   0xBFAAD261EAA09954, 0xBFAAD261EAA09954 // dA09
 465         .align  16
 466         .quad   0x3FAE1749BD612DCF, 0x3FAE1749BD612DCF // dA08
 467         .align  16
 468         .quad   0xBFB11084009435E0, 0xBFB11084009435E0 // dA07
 469         .align  16
 470         .quad   0x3FB3B12A49295651, 0x3FB3B12A49295651 // dA06
 471         .align  16
 472         .quad   0xBFB745D009BADA94, 0xBFB745D009BADA94 // dA05
 473         .align  16
 474         .quad   0x3FBC71C707F7D5B5, 0x3FBC71C707F7D5B5 // dA04
 475         .align  16
 476         .quad   0xBFC2492491EE55C7, 0xBFC2492491EE55C7 // dA03
 477         .align  16
 478         .quad   0x3FC999999997EE34, 0x3FC999999997EE34 // dA02
 479         .align  16
 480         .quad   0xBFD55555555553C5, 0xBFD55555555553C5 // dA01
 481         .align  16
 482         .quad   0x3FF0000000000000, 0x3FF0000000000000 // dA00
 483         .align  16
 484         .quad   0x8000000000000000, 0x8000000000000000 // dSIGN_MASK
 485         .align  16
 486         .long   0x80300000, 0x80300000, 0x80300000, 0x80300000 // iCHK_WORK_SUB
 487         .align  16
 488         .long   0xfdd00000, 0xfdd00000, 0xfdd00000, 0xfdd00000 // iCHK_WORK_CMP
 489         .align  16
 490         .quad   0x7fffffffffffffff, 0x7fffffffffffffff // dABS_MASK
 491         .align  16
 492         .quad   0x0000000000000000, 0x0000000000000000 // dZERO
 493         .align  16
 494         .type   __svml_datan2_data_internal, @object
 495         .size   __svml_datan2_data_internal, .-__svml_datan2_data_internal