/* Copyright (c) 2018, Oracle and/or its affiliates. All rights reserved.

 * Copyright (c) 2018, Cavium. All rights reserved. (By BELLSOFT)

 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.

 *

 * This code is free software; you can redistribute it and/or modify it

 * under the terms of the GNU General Public License version 2 only, as

 * published by the Free Software Foundation.

 *

 * This code is distributed in the hope that it will be useful, but WITHOUT

 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

 * version 2 for more details (a copy is included in the LICENSE file that

 * accompanied this code).

 *

 * You should have received a copy of the GNU General Public License version

 * 2 along with this work; if not, write to the Free Software Foundation,

 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.

 *

 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA

 * or visit www.oracle.com if you need additional information or have any

 * questions.

 *

 */

#include "precompiled.hpp"

#include "asm/assembler.hpp"

#include "asm/assembler.inline.hpp"

#include "runtime/stubRoutines.hpp"

#include "macroAssembler_aarch64.hpp"

// The following code is a optimized version of fdlibm sin/cos implementation

// (C code is in share/runtime/sharedRuntimeTrig.cpp) adapted for AARCH64.

// Please refer to sin/cos approximation via polynomial and

// trigonometric argument reduction techniques to the following literature:

//

// [1] Muller, Jean-Michel, Nicolas Brisebarre, Florent De Dinechin,

// Claude-Pierre Jeannerod, Vincent Lefevre, Guillaume Melquiond,

// Nathalie Revol, Damien Stehlé, and Serge Torres:

// Handbook of floating-point arithmetic.

// Springer Science & Business Media, 2009.

// [2] K. C. Ng

// Argument Reduction for Huge Arguments: Good to the Last Bit

// July 13, 1992, SunPro

//

// HOW TO READ THIS CODE:

// This code consists of several functions. Each function has following header:

// 1) Description

// 2) C-pseudo code with differences from fdlibm marked by comments starting

//        with "NOTE". Check unmodified fdlibm code in

//        share/runtime/SharedRuntimeTrig.cpp

// 3) Brief textual description of changes between fdlibm and current

//        implementation along with optimization notes (if applicable)

// 4) Assumptions, input and output

// 5) (Optional) additional notes about intrinsic implementation

// Each function is separated in blocks which follow the pseudo-code structure

//

// HIGH-LEVEL ALGORITHM DESCRIPTION:

//    - entry point: generate_dsin_dcos(...);

//    - check corner cases: NaN, INF, tiny argument.

//    - check if |x| < Pi/4. Then approximate sin/cos via polynomial (kernel_sin/kernel_cos)

//    -- else proceed to argument reduction routine (__ieee754_rem_pio2) and

//           use reduced argument to get result via kernel_sin/kernel_cos

//

// HIGH-LEVEL CHANGES BETWEEN INTRINSICS AND FDLIBM:

// 1) two_over_pi table fdlibm representation is int[], while intrinsic version

// has these int values converted to double representation to load converted

// double values directly (see stubRoutines_aarch4::_two_over_pi)

// 2) Several loops are unrolled and vectorized: see comments in code after

// labels: SKIP_F_LOAD, RECOMP_FOR1_CHECK, RECOMP_FOR2

// 3) fdlibm npio2_hw table now has "prefix" with constants used in

// calculation. These constants are loaded from npio2_hw table instead of

// constructing it in code (see stubRoutines_aarch64.cpp)

// 4) Polynomial coefficients for sin and cos are moved to table sin_coef

// and cos_coef to use the same optimization as in 3). It allows to load most of

// required constants via single instruction

//

//

//

///* __ieee754_rem_pio2(x,y)

// *

// * returns the remainder of x rem pi/2 in y[0]+y[1] (i.e. like x div pi/2)

// * x is input argument, y[] is hi and low parts of reduced argument (x)

// * uses __kernel_rem_pio2()

// */

// // use tables(see stubRoutines_aarch64.cpp): two_over_pi and modified npio2_hw

//

// BEGIN __ieee754_rem_pio2 PSEUDO CODE

//

//static int __ieee754_rem_pio2(double x, double *y) {

//  double z,w,t,r,fn;

//  double tx[3];

//  int e0,i,j,nx,n,ix,hx,i0;

//

//  i0 = ((*(int*)&two24A)>>30)^1;        /* high word index */

//  hx = *(i0+(int*)&x);          /* high word of x */

//  ix = hx&0x7fffffff;

//  if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */

//    if(hx>0) {

//      z = x - pio2_1;

//      if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */

//        y[0] = z - pio2_1t;

//        y[1] = (z-y[0])-pio2_1t;

//      } else {                /* near pi/2, use 33+33+53 bit pi */

//        z -= pio2_2;

//        y[0] = z - pio2_2t;

//        y[1] = (z-y[0])-pio2_2t;

//      }

//      return 1;

//    } else {    /* negative x */

//      z = x + pio2_1;

//      if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */

//        y[0] = z + pio2_1t;

//        y[1] = (z-y[0])+pio2_1t;

//      } else {                /* near pi/2, use 33+33+53 bit pi */

//        z += pio2_2;

//        y[0] = z + pio2_2t;

//        y[1] = (z-y[0])+pio2_2t;

//      }

//      return -1;

//    }

//  }

//  if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */

//    t  = fabsd(x);

//    n  = (int) (t*invpio2+half);

//    fn = (double)n;

//    r  = t-fn*pio2_1;

//    w  = fn*pio2_1t;    /* 1st round good to 85 bit */

//    // NOTE: y[0] = r-w; is moved from if/else below to be before "if"

//    y[0] = r-w;

//    if(n<32&&ix!=npio2_hw[n-1]) {

//      // y[0] = r-w;       /* quick check no cancellation */ // NOTE: moved earlier

//    } else {

//      j  = ix>>20;

//      // y[0] = r-w; // NOTE: moved earlier

//      i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff);

//      if(i>16) {  /* 2nd iteration needed, good to 118 */

//        t  = r;

//        w  = fn*pio2_2;

//        r  = t-w;

//        w  = fn*pio2_2t-((t-r)-w);

//        y[0] = r-w;

//        i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff);

//        if(i>49)  {     /* 3rd iteration need, 151 bits acc */

//          t  = r;       /* will cover all possible cases */

//          w  = fn*pio2_3;

//          r  = t-w;

//          w  = fn*pio2_3t-((t-r)-w);

//          y[0] = r-w;

//        }

//      }

//    }

//    y[1] = (r-y[0])-w;

//    if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}

//    else         return n;

//  }

//  /*

//   * all other (large) arguments

//   */

//  // NOTE: this check is removed, because it was checked in dsin/dcos

//  // if(ix>=0x7ff00000) {          /* x is inf or NaN */

//  //  y[0]=y[1]=x-x; return 0;

//  // }

//  /* set z = scalbn(|x|,ilogb(x)-23) */

//  *(1-i0+(int*)&z) = *(1-i0+(int*)&x);

//  e0    = (ix>>20)-1046;        /* e0 = ilogb(z)-23; */

//  *(i0+(int*)&z) = ix - (e0<<20);

//

//  // NOTE: "for" loop below in unrolled. See comments in asm code

//  for(i=0;i<2;i++) {

//    tx[i] = (double)((int)(z));

//    z     = (z-tx[i])*two24A;

//  }

//

//  tx[2] = z;

//  nx = 3;

//

//  // NOTE: while(tx[nx-1]==zeroA) nx--;  is unrolled. See comments in asm code

//  while(tx[nx-1]==zeroA) nx--;  /* skip zero term */

//

//  n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);

//  if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}

//  return n;

//}

//

// END __ieee754_rem_pio2 PSEUDO CODE

//

// Changes between fdlibm and intrinsic for __ieee754_rem_pio2:

//     1. INF/NaN check for huge argument is removed in comparison with fdlibm

//     code, because this check is already done in dcos/dsin code

//     2. Most constants are now loaded from table instead of direct initialization

//     3. Two loops are unrolled

// Assumptions:

//     1. Assume |X| >= PI/4

//     2. Assume rscratch1 = 0x3fe921fb00000000  (~ PI/4)

//     3. Assume ix = r3

// Input and output:

//     1. Input: X = r0

//     2. Return n in r2, y[0] == y0 == v4, y[1] == y1 == v5

// NOTE: general purpose register names match local variable names in C code

// NOTE: fpu registers are actively reused. See comments in code about their usage

void MacroAssembler::generate__ieee754_rem_pio2(address npio2_hw,

    address two_over_pi, address pio2) {

  const long PIO2_1t = 0x3DD0B4611A626331UL;

  const long PIO2_2  = 0x3DD0B4611A600000UL;

  const long PIO2_2t = 0x3BA3198A2E037073UL;

  Label X_IS_NEGATIVE, X_IS_MEDIUM_OR_LARGE, X_IS_POSITIVE_LONG_PI, LARGE_ELSE,

      REDUCTION_DONE, X_IS_MEDIUM_BRANCH_DONE, X_IS_LARGE, NX_SET,

      X_IS_NEGATIVE_LONG_PI;

  Register X = r0, n = r2, ix = r3, jv = r4, tmp5 = r5, jx = r6,

      tmp3 = r7, iqBase = r10, ih = r11, i = r17;

    // initializing constants first

    // rscratch1 = 0x3fe921fb00000000 (see assumptions)

    movk(rscratch1, 0x3ff9, 48); // was 0x3fe921fb0..0 now it's 0x3ff921fb0..0

    mov(rscratch2, 0x4002d97c); // 3*PI/4 high word

    movk(rscratch1, 0x5440, 16); // now rscratch1 == PIO2_1

    fmovd(v1, rscratch1); // v1 = PIO2_1

    cmp(rscratch2, ix);

    br(LE, X_IS_MEDIUM_OR_LARGE);

    block_comment("if(ix<0x4002d97c) {...  /* |x| ~< 3pi/4 */ "); {

      cmp(X, zr);

      br(LT, X_IS_NEGATIVE);

      block_comment("if(hx>0) {"); {

        fsubd(v2, v0, v1); // v2 = z = x - pio2_1

        cmp(ix, rscratch1, LSR, 32);

        mov(n, 1);

        br(EQ, X_IS_POSITIVE_LONG_PI);

        block_comment("case: hx > 0 &&  ix!=0x3ff921fb {"); { /* 33+53 bit pi is good enough */

          mov(rscratch2, PIO2_1t);

          fmovd(v27, rscratch2);

          fsubd(v4, v2, v27); // v4 = y[0] = z - pio2_1t;

          fsubd(v5, v2, v4);

          fsubd(v5, v5, v27); // v5 = y[1] = (z-y[0])-pio2_1t

          b(REDUCTION_DONE);

        }

        block_comment("case: hx > 0 &*& ix==0x3ff921fb {"); { /* near pi/2, use 33+33+53 bit pi */

          bind(X_IS_POSITIVE_LONG_PI);

            mov(rscratch1, PIO2_2);

            mov(rscratch2, PIO2_2t);

            fmovd(v27, rscratch1);

            fmovd(v6, rscratch2);

            fsubd(v2, v2, v27); // z-= pio2_2

            fsubd(v4, v2, v6);  // y[0] = z - pio2_2t

            fsubd(v5, v2, v4);

            fsubd(v5, v5, v6);  // v5 = (z - y[0]) - pio2_2t

            b(REDUCTION_DONE);

        }

      }

      block_comment("case: hx <= 0)"); {

        bind(X_IS_NEGATIVE);

          faddd(v2, v0, v1); // v2 = z = x + pio2_1

          cmp(ix, rscratch1, LSR, 32);

          mov(n, -1);

          br(EQ, X_IS_NEGATIVE_LONG_PI);

          block_comment("case: hx <= 0 && ix!=0x3ff921fb) {"); { /* 33+53 bit pi is good enough */

            mov(rscratch2, PIO2_1t);

            fmovd(v27, rscratch2);

            faddd(v4, v2, v27); // v4 = y[0] = z + pio2_1t;

            fsubd(v5, v2, v4);

            faddd(v5, v5, v27); // v5 = y[1] = (z-y[0]) + pio2_1t

            b(REDUCTION_DONE);

          }

          block_comment("case: hx <= 0 && ix==0x3ff921fb"); { /* near pi/2, use 33+33+53 bit pi */

            bind(X_IS_NEGATIVE_LONG_PI);

              mov(rscratch1, PIO2_2);

              mov(rscratch2, PIO2_2t);

              fmovd(v27, rscratch1);

              fmovd(v6, rscratch2);

              faddd(v2, v2, v27); // z += pio2_2

              faddd(v4, v2, v6);  // y[0] = z + pio2_2t

              fsubd(v5, v2, v4);

              faddd(v5, v5, v6);  // v5 = (z - y[0]) + pio2_2t

              b(REDUCTION_DONE);

          }

      }

  }

  bind(X_IS_MEDIUM_OR_LARGE);

    mov(rscratch1, 0x413921fb);

    cmp(ix, rscratch1); // ix < = 0x413921fb ?

    br(GT, X_IS_LARGE);

    block_comment("|x| ~<= 2^19*(pi/2), medium size"); {

      lea(ih, ExternalAddress(npio2_hw));

      ld1(v4, v5, v6, v7, T1D, ih);

      fabsd(v31, v0);          // v31 = t = |x|

      add(ih, ih, 64);

      fmaddd(v2, v31, v5, v4); // v2 = t * invpio2 + half (invpio2 = 53 bits of 2/pi, half = 0.5)

      fcvtzdw(n, v2);          // n = (int) v2

      frintzd(v2, v2);

      fmsubd(v3, v2, v6, v31); // v3 = r = t - fn * pio2_1

      fmuld(v26, v2, v7);      // v26 = w = fn * pio2_1t

      fsubd(v4, v3, v26);      // y[0] = r - w. Calculated before branch

      cmp(n, (u1)32);

      br(GT, LARGE_ELSE);

      subw(tmp5, n, 1);        // tmp5 = n - 1

      ldrw(jv, Address(ih, tmp5, Address::lsl(2)));

      cmp(ix, jv);

      br(NE, X_IS_MEDIUM_BRANCH_DONE);

      block_comment("else block for if(n<32&&ix!=npio2_hw[n-1])"); {

        bind(LARGE_ELSE);

          fmovd(jx, v4);

          lsr(tmp5, ix, 20);                       // j = ix >> 20

          lsl(jx, jx, 1);

          sub(tmp3, tmp5, jx, LSR, 32 + 20 + 1);   // r7 = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff);

          block_comment("if(i>16)"); {

            cmp(tmp3, (u1)16);

            br(LE, X_IS_MEDIUM_BRANCH_DONE);

            // i > 16. 2nd iteration needed

            ldpd(v6, v7, Address(ih, -32));

            fmovd(v28, v3);                        // t = r

            fmuld(v29, v2, v6);                    // w = v29 = fn * pio2_2

            fsubd(v3, v28, v29);                   // r = t - w

            fsubd(v31, v28, v3);                   // v31 = (t - r)

            fsubd(v31, v29, v31);                  // v31 = w - (t - r) = - ((t - r) - w)

            fmaddd(v26, v2, v7, v31);              // v26 = w = fn*pio2_2t - ((t - r) - w)

            fsubd(v4, v3, v26);                    // y[0] = r - w

            fmovd(jx, v4);

            lsl(jx, jx, 1);

            sub(tmp3, tmp5, jx, LSR, 32 + 20 + 1); // r7 = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff);

            block_comment("if(i>49)"); {

              cmp(tmp3, (u1)49);

              br(LE, X_IS_MEDIUM_BRANCH_DONE);

              // 3rd iteration need, 151 bits acc

              ldpd(v6, v7, Address(ih, -16));

              fmovd(v28, v3);                      // save "r"

              fmuld(v29, v2, v6);                  // v29 = fn * pio2_3

              fsubd(v3, v28, v29);                 // r = r - w

              fsubd(v31, v28, v3);                 // v31 = (t - r)

              fsubd(v31, v29, v31);                // v31 = w - (t - r) = - ((t - r) - w)

              fmaddd(v26, v2, v7, v31);            // v26 = w = fn*pio2_3t - ((t - r) - w)

              fsubd(v4, v3, v26);                  // y[0] = r - w

            }

          }

      }

    block_comment("medium x tail"); {

      bind(X_IS_MEDIUM_BRANCH_DONE);

        fsubd(v5, v3, v4);                         // v5 = y[1] = (r - y[0])

        fsubd(v5, v5, v26);                        // v5 = y[1] = (r - y[0]) - w

        cmp(X, zr);

        br(GT, REDUCTION_DONE);

        fnegd(v4, v4);

        negw(n, n);

        fnegd(v5, v5);

        b(REDUCTION_DONE);

    }

  }

  block_comment("all other (large) arguments"); {

    bind(X_IS_LARGE);

      lsr(rscratch1, ix, 20);                      // ix >> 20

      movz(tmp5, 0x4170, 48);

      subw(rscratch1, rscratch1, 1046);            // e0

      fmovd(v10, tmp5);                            // init two24A value

      subw(jv, ix, rscratch1, LSL, 20);            // ix - (e0<<20)

      lsl(jv, jv, 32);

      subw(rscratch2, rscratch1, 3);

      bfm(jv, X, 0, 31);                           // jv = z

      movw(i, 24);

      fmovd(v26, jv);                              // v26 = z

      block_comment("unrolled for(i=0;i<2;i++) {tx[i] = (double)((int)(z));z = (z-tx[i])*two24A;}"); {

        // tx[0,1,2] = v6,v7,v26

        frintzd(v6, v26);                          // v6 = (double)((int)v26)

        sdivw(jv, rscratch2, i);                   // jv = (e0 - 3)/24

        fsubd(v26, v26, v6);

        sub(sp, sp, 560);

        fmuld(v26, v26, v10);

        frintzd(v7, v26);                          // v7 = (double)((int)v26)

        movw(jx, 2); // calculate jx as nx - 1, which is initially 2. Not a part of unrolled loop

        fsubd(v26, v26, v7);

      }

      block_comment("nx calculation with unrolled while(tx[nx-1]==zeroA) nx--;"); {

        add(iqBase, sp, 480);                      // base of iq[]

        fmuld(v3, v26, v10);

        br(NE, NX_SET);

        csetw(jx, NE);

      }

    bind(NX_SET);

      generate__kernel_rem_pio2(two_over_pi, pio2);

      // now we have y[0] = v4, y[1] = v5 and n = r2

      cmp(X, zr);

      br(GE, REDUCTION_DONE);

      fnegd(v4, v4);

      fnegd(v5, v5);

      negw(n, n);

  }

  bind(REDUCTION_DONE);

}

///*

// * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)

// * double x[],y[]; int e0,nx,prec; int ipio2[];

// *

// * __kernel_rem_pio2 return the last three digits of N with

// *              y = x - N*pi/2

// * so that |y| < pi/2.

// *

// * The method is to compute the integer (mod 8) and fraction parts of

// * (2/pi)*x without doing the full multiplication. In general we

// * skip the part of the product that are known to be a huge integer (

// * more accurately, = 0 mod 8 ). Thus the number of operations are

// * independent of the exponent of the input.

// *

// * NOTE: 2/pi int representation is converted to double

// * // (2/pi) is represented by an array of 24-bit integers in ipio2[].

// *

// * Input parameters:

// *      x[]     The input value (must be positive) is broken into nx

// *              pieces of 24-bit integers in double precision format.

// *              x[i] will be the i-th 24 bit of x. The scaled exponent

// *              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0

// *              match x's up to 24 bits.

// *

// *              Example of breaking a double positive z into x[0]+x[1]+x[2]:

// *                      e0 = ilogb(z)-23

// *                      z  = scalbn(z,-e0)

// *              for i = 0,1,2

// *                      x[i] = floor(z)

// *                      z    = (z-x[i])*2**24

// *

// *

// *      y[]     ouput result in an array of double precision numbers.

// *              The dimension of y[] is:

// *                      24-bit  precision       1

// *                      53-bit  precision       2

// *                      64-bit  precision       2

// *                      113-bit precision       3

// *              The actual value is the sum of them. Thus for 113-bit

// *              precsion, one may have to do something like:

// *

// *              long double t,w,r_head, r_tail;

// *              t = (long double)y[2] + (long double)y[1];

// *              w = (long double)y[0];

// *              r_head = t+w;

// *              r_tail = w - (r_head - t);

// *

// *      e0      The exponent of x[0]

// *

// *      nx      dimension of x[]

// *

// *      prec    an interger indicating the precision:

// *                      0       24  bits (single)

// *                      1       53  bits (double)

// *                      2       64  bits (extended)

// *                      3       113 bits (quad)

// *

// *      NOTE: ipio2[] array below is converted to double representation

// *      //ipio2[]

// *      //        integer array, contains the (24*i)-th to (24*i+23)-th

// *      //        bit of 2/pi after binary point. The corresponding

// *      //        floating value is

// *

// *                      ipio2[i] * 2^(-24(i+1)).

// *

// * Here is the description of some local variables:

// *

// *      jk      jk+1 is the initial number of terms of ipio2[] needed

// *              in the computation. The recommended value is 2,3,4,

// *              6 for single, double, extended,and quad.

// *

// *      jz      local integer variable indicating the number of

// *              terms of ipio2[] used.

// *

// *      jx      nx - 1

// *

// *      jv      index for pointing to the suitable ipio2[] for the

// *              computation. In general, we want

// *                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8

// *              is an integer. Thus

// *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv

// *              Hence jv = max(0,(e0-3)/24).

// *

// *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.

// *

// *      q[]     double array with integral value, representing the

// *              24-bits chunk of the product of x and 2/pi.

// *

// *      q0      the corresponding exponent of q[0]. Note that the

// *              exponent for q[i] would be q0-24*i.

// *

// *      PIo2[]  double precision array, obtained by cutting pi/2

// *              into 24 bits chunks.

// *

// *      f[]     ipio2[] in floating point

// *

// *      iq[]    integer array by breaking up q[] in 24-bits chunk.

// *

// *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]

// *

// *      ih      integer. If >0 it indicates q[] is >= 0.5, hence

// *              it also indicates the *sign* of the result.

// *

// */

//

// Use PIo2 table(see stubRoutines_aarch64.cpp)

//

// BEGIN __kernel_rem_pio2 PSEUDO CODE

//

//static int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, /* NOTE: converted to double */ const double *ipio2 // const int *ipio2) {

//  int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;

//  double z,fw,f[20],fq[20],q[20];

//

//  /* initialize jk*/

//  // jk = init_jk[prec]; // NOTE: prec==2 for double. jk is always 4.

//  jp = jk; // NOTE: always 4

//

//  /* determine jx,jv,q0, note that 3>q0 */

//  jx =  nx-1;

//  jv = (e0-3)/24; if(jv<0) jv=0;

//  q0 =  e0-24*(jv+1);

//

//  /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */

//  j = jv-jx; m = jx+jk;

//

//  // NOTE: split into two for-loops: one with zeroB and one with ipio2[j]. It