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

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

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

     4  *

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

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

     7  * published by the Free Software Foundation.

     8  *

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

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

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

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

    13  * accompanied this code).

    14  *

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

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

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

    18  *

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

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

    21  * questions.

    22  *

    23  */

    24 

    25 #include "precompiled.hpp"

    26 #include "asm/assembler.hpp"

    27 #include "asm/assembler.inline.hpp"

    28 #include "runtime/stubRoutines.hpp"

    29 #include "macroAssembler_aarch64.hpp"

    30 

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

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

    33 

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

    35 // trigonometric argument reduction techniques to the following literature:

    36 //

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

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

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

    40 // Handbook of floating-point arithmetic.

    41 // Springer Science & Business Media, 2009.

    42 // [2] K. C. Ng

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

    44 // July 13, 1992, SunPro

    45 //

    46 // HOW TO READ THIS CODE:

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

    48 // 1) Description

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

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

    51 //        share/runtime/SharedRuntimeTrig.cpp

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

    53 //        implementation along with optimization notes (if applicable)

    54 // 4) Assumptions, input and output

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

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

    57 //

    58 // HIGH-LEVEL ALGORITHM DESCRIPTION:

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

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

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

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

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

    64 //

    65 // HIGH-LEVEL CHANGES BETWEEN INTRINSICS AND FDLIBM:

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

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

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

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

    70 // labels: SKIP_F_LOAD, RECOMP_FOR1_CHECK, RECOMP_FOR2

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

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

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

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

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

    76 // required constants via single instruction

    77 //

    78 //

    79 //

    80 ///* __ieee754_rem_pio2(x,y)

    81 // *

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

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

    84 // * uses __kernel_rem_pio2()

    85 // */

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

    87 //

    88 // BEGIN __ieee754_rem_pio2 PSEUDO CODE

    89 //

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

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

    92 //  double tx[3];

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

    94 //

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

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

    97 //  ix = hx&0x7fffffff;

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

    99 //    if(hx>0) {

   100 //      z = x - pio2_1;

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

   102 //        y[0] = z - pio2_1t;

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

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

   105 //        z -= pio2_2;

   106 //        y[0] = z - pio2_2t;

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

   108 //      }

   109 //      return 1;

   110 //    } else {    /* negative x */

   111 //      z = x + pio2_1;

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

   113 //        y[0] = z + pio2_1t;

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

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

   116 //        z += pio2_2;

   117 //        y[0] = z + pio2_2t;

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

   119 //      }

   120 //      return -1;

   121 //    }

   122 //  }

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

   124 //    t  = fabsd(x);

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

   126 //    fn = (double)n;

   127 //    r  = t-fn*pio2_1;

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

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

   130 //    y[0] = r-w;

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

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

   133 //    } else {

   134 //      j  = ix>>20;

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

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

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

   138 //        t  = r;

   139 //        w  = fn*pio2_2;

   140 //        r  = t-w;

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

   142 //        y[0] = r-w;

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

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

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

   146 //          w  = fn*pio2_3;

   147 //          r  = t-w;

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

   149 //          y[0] = r-w;

   150 //        }

   151 //      }

   152 //    }

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

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

   155 //    else         return n;

   156 //  }

   157 //  /*

   158 //   * all other (large) arguments

   159 //   */

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

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

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

   163 //  // }

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

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

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

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

   168 //

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

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

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

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

   173 //  }

   174 //

   175 //  tx[2] = z;

   176 //  nx = 3;

   177 //

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

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

   180 //

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

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

   183 //  return n;

   184 //}

   185 //

   186 // END __ieee754_rem_pio2 PSEUDO CODE

   187 //

   188 // Changes between fdlibm and intrinsic for __ieee754_rem_pio2:

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

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

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

   192 //     3. Two loops are unrolled

   193 // Assumptions:

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

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

   196 //     3. Assume ix = r3

   197 // Input and output:

   198 //     1. Input: X = r0

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

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

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

   202 void MacroAssembler::generate__ieee754_rem_pio2(address npio2_hw,

   203     address two_over_pi, address pio2) {

   204   const long PIO2_1t = 0x3DD0B4611A626331UL;

   205   const long PIO2_2  = 0x3DD0B4611A600000UL;

   206   const long PIO2_2t = 0x3BA3198A2E037073UL;

   207   Label X_IS_NEGATIVE, X_IS_MEDIUM_OR_LARGE, X_IS_POSITIVE_LONG_PI, LARGE_ELSE,

   208       REDUCTION_DONE, X_IS_MEDIUM_BRANCH_DONE, X_IS_LARGE, NX_SET,

   209       X_IS_NEGATIVE_LONG_PI;

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

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

   212     // initializing constants first

   213     // rscratch1 = 0x3fe921fb00000000 (see assumptions)

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

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

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

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

   218     cmp(rscratch2, ix);

   219     br(LE, X_IS_MEDIUM_OR_LARGE);

   220 

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

   222       cmp(X, zr);

   223       br(LT, X_IS_NEGATIVE);

   224 

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

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

   227         cmp(ix, rscratch1, LSR, 32);

   228         mov(n, 1);

   229         br(EQ, X_IS_POSITIVE_LONG_PI);

   230 

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

   232           mov(rscratch2, PIO2_1t);

   233           fmovd(v27, rscratch2);

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

   235           fsubd(v5, v2, v4);

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

   237           b(REDUCTION_DONE);

   238         }

   239 

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

   241           bind(X_IS_POSITIVE_LONG_PI);

   242             mov(rscratch1, PIO2_2);

   243             mov(rscratch2, PIO2_2t);

   244             fmovd(v27, rscratch1);

   245             fmovd(v6, rscratch2);

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

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

   248             fsubd(v5, v2, v4);

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

   250             b(REDUCTION_DONE);

   251         }

   252       }

   253 

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

   255         bind(X_IS_NEGATIVE);

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

   257           cmp(ix, rscratch1, LSR, 32);

   258           mov(n, -1);

   259           br(EQ, X_IS_NEGATIVE_LONG_PI);

   260 

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

   262             mov(rscratch2, PIO2_1t);

   263             fmovd(v27, rscratch2);

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

   265             fsubd(v5, v2, v4);

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

   267             b(REDUCTION_DONE);

   268           }

   269 

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

   271             bind(X_IS_NEGATIVE_LONG_PI);

   272               mov(rscratch1, PIO2_2);

   273               mov(rscratch2, PIO2_2t);

   274               fmovd(v27, rscratch1);

   275               fmovd(v6, rscratch2);

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

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

   278               fsubd(v5, v2, v4);

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

   280               b(REDUCTION_DONE);

   281           }

   282       }

   283   }

   284   bind(X_IS_MEDIUM_OR_LARGE);

   285     mov(rscratch1, 0x413921fb);

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

   287     br(GT, X_IS_LARGE);

   288 

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

   290       lea(ih, ExternalAddress(npio2_hw));

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

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

   293       add(ih, ih, 64);

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

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

   296       frintzd(v2, v2);

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

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

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

   300       cmp(n, 32);

   301       br(GT, LARGE_ELSE);

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

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

   304       cmp(ix, jv);

   305       br(NE, X_IS_MEDIUM_BRANCH_DONE);

   306 

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

   308         bind(LARGE_ELSE);

   309           fmovd(jx, v4);

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

   311           lsl(jx, jx, 1);

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

   313 

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

   315             cmp(tmp3, 16);

   316             br(LE, X_IS_MEDIUM_BRANCH_DONE);

   317             // i > 16. 2nd iteration needed

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

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

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

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

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

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

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

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

   326             fmovd(jx, v4);

   327             lsl(jx, jx, 1);

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

   329 

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

   331               cmp(tmp3, 49);

   332               br(LE, X_IS_MEDIUM_BRANCH_DONE);

   333               // 3rd iteration need, 151 bits acc

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

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

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

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

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

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

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

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

   342             }

   343           }

   344       }

   345     block_comment("medium x tail"); {

   346       bind(X_IS_MEDIUM_BRANCH_DONE);

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

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

   349         cmp(X, zr);

   350         br(GT, REDUCTION_DONE);

   351         fnegd(v4, v4);

   352         negw(n, n);

   353         fnegd(v5, v5);

   354         b(REDUCTION_DONE);

   355     }

   356   }

   357 

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

   359     bind(X_IS_LARGE);

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

   361       movz(tmp5, 0x4170, 48);

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

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

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

   365       lsl(jv, jv, 32);

   366       subw(rscratch2, rscratch1, 3);

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

   368       movw(i, 24);

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

   370 

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

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

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

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

   375         fsubd(v26, v26, v6);

   376         sub(sp, sp, 560);

   377         fmuld(v26, v26, v10);

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

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

   380         fsubd(v26, v26, v7);

   381       }

   382 

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

   384         fcmpd(v26, 0.0d);                          // if NE then jx == 2. else it's 1 or 0

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

   386         fmuld(v3, v26, v10);

   387         br(NE, NX_SET);

   388         fcmpd(v7, 0.0d);                           // v7 == 0 => jx = 0. Else jx = 1

   389         csetw(jx, NE);

   390       }

   391     bind(NX_SET);

   392       generate__kernel_rem_pio2(two_over_pi, pio2);

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

   394       cmp(X, zr);

   395       br(GE, REDUCTION_DONE);

   396       fnegd(v4, v4);

   397       fnegd(v5, v5);

   398       negw(n, n);

   399   }

   400   bind(REDUCTION_DONE);

   401 }

   402 

   403 ///*

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

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

   406 // *

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

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

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

   410 // *

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

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

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

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

   415 // * independent of the exponent of the input.

   416 // *

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

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

   419 // *

   420 // * Input parameters:

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

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

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

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

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

   426 // *

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

   428 // *                      e0 = ilogb(z)-23

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

   430 // *              for i = 0,1,2

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

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

   433 // *

   434 // *

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

   436 // *              The dimension of y[] is:

   437 // *                      24-bit  precision       1

   438 // *                      53-bit  precision       2

   439 // *                      64-bit  precision       2

   440 // *                      113-bit precision       3

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

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

   443 // *

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

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

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

   447 // *              r_head = t+w;

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

   449 // *

   450 // *      e0      The exponent of x[0]

   451 // *

   452 // *      nx      dimension of x[]

   453 // *

   454 // *      prec    an interger indicating the precision:

   455 // *                      0       24  bits (single)

   456 // *                      1       53  bits (double)

   457 // *                      2       64  bits (extended)

   458 // *                      3       113 bits (quad)

   459 // *

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

   461 // *      //ipio2[]

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

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

   464 // *      //        floating value is

   465 // *

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

   467 // *

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

   469 // *

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

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

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

   473 // *

   474 // *      jz      local integer variable indicating the number of

   475 // *              terms of ipio2[] used.

   476 // *

   477 // *      jx      nx - 1

   478 // *

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

   480 // *              computation. In general, we want

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

   482 // *              is an integer. Thus

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

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

   485 // *

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

   487 // *

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

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

   490 // *

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

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

   493 // *

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

   495 // *              into 24 bits chunks.

   496 // *

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

   498 // *

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

   500 // *

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

   502 // *

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

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

   505 // *

   506 // */

   507 //

   508 // Use PIo2 table(see stubRoutines_aarch64.cpp)

   509 //

   510 // BEGIN __kernel_rem_pio2 PSEUDO CODE

   511 //

   512 //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) {

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

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

   515 //

   516 //  /* initialize jk*/

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

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

   519 //

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

   521 //  jx =  nx-1;

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

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

   524 //

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

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

   527 //

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

   529 //  //       allows the use of wider loads/stores

   530 //  for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : /* NOTE: converted to double */ ipio2[j]; //(double) ipio2[j];

   531 //

   532 //  // NOTE: unrolled and vectorized "for". See comments in asm code

   533 //  /* compute q[0],q[1],...q[jk] */

   534 //  for (i=0;i<=jk;i++) {

   535 //    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;

   536 //  }

   537 //

   538 //  jz = jk;

   539 //recompute:

   540 //  /* distill q[] into iq[] reversingly */

   541 //  for(i=0,j=jz,z=q[jz];j>0;i++,j--) {

   542 //    fw    =  (double)((int)(twon24* z));

   543 //    iq[i] =  (int)(z-two24B*fw);

   544 //    z     =  q[j-1]+fw;

   545 //  }

   546 //

   547 //  /* compute n */

   548 //  z  = scalbnA(z,q0);           /* actual value of z */

   549 //  z -= 8.0*floor(z*0.125);              /* trim off integer >= 8 */

   550 //  n  = (int) z;

   551 //  z -= (double)n;

   552 //  ih = 0;

   553 //  if(q0>0) {    /* need iq[jz-1] to determine n */

   554 //    i  = (iq[jz-1]>>(24-q0)); n += i;

   555 //    iq[jz-1] -= i<<(24-q0);

   556 //    ih = iq[jz-1]>>(23-q0);

   557 //  }

   558 //  else if(q0==0) ih = iq[jz-1]>>23;

   559 //  else if(z>=0.5) ih=2;

   560 //

   561 //  if(ih>0) {    /* q > 0.5 */

   562 //    n += 1; carry = 0;

   563 //    for(i=0;i<jz ;i++) {        /* compute 1-q */

   564 //      j = iq[i];

   565 //      if(carry==0) {

   566 //        if(j!=0) {

   567 //          carry = 1; iq[i] = 0x1000000- j;

   568 //        }

   569 //      } else  iq[i] = 0xffffff - j;

   570 //    }

   571 //    if(q0>0) {          /* rare case: chance is 1 in 12 */

   572 //      switch(q0) {

   573 //      case 1:

   574 //        iq[jz-1] &= 0x7fffff; break;

   575 //      case 2:

   576 //        iq[jz-1] &= 0x3fffff; break;

   577 //      }

   578 //    }

   579 //    if(ih==2) {

   580 //      z = one - z;

   581 //      if(carry!=0) z -= scalbnA(one,q0);

   582 //    }

   583 //  }

   584 //

   585 //  /* check if recomputation is needed */

   586 //  if(z==zeroB) {

   587 //    j = 0;

   588 //    for (i=jz-1;i>=jk;i--) j |= iq[i];

   589 //    if(j==0) { /* need recomputation */

   590 //      for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */

   591 //

   592 //      for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */

   593 //        f[jx+i] = /* NOTE: converted to double */ ipio2[jv+i]; //(double) ipio2[jv+i];

   594 //        for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];

   595 //        q[i] = fw;

   596 //      }

   597 //      jz += k;

   598 //      goto recompute;

   599 //    }

   600 //  }

   601 //

   602 //  /* chop off zero terms */

   603 //  if(z==0.0) {

   604 //    jz -= 1; q0 -= 24;

   605 //    while(iq[jz]==0) { jz--; q0-=24;}

   606 //  } else { /* break z into 24-bit if necessary */

   607 //    z = scalbnA(z,-q0);

   608 //    if(z>=two24B) {

   609 //      fw = (double)((int)(twon24*z));

   610 //      iq[jz] = (int)(z-two24B*fw);

   611 //      jz += 1; q0 += 24;

   612 //      iq[jz] = (int) fw;

   613 //    } else iq[jz] = (int) z ;

   614 //  }

   615 //

   616 //  /* convert integer "bit" chunk to floating-point value */

   617 //  fw = scalbnA(one,q0);

   618 //  for(i=jz;i>=0;i--) {

   619 //    q[i] = fw*(double)iq[i]; fw*=twon24;

   620 //  }

   621 //

   622 //  /* compute PIo2[0,...,jp]*q[jz,...,0] */

   623 //  for(i=jz;i>=0;i--) {

   624 //    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];

   625 //    fq[jz-i] = fw;

   626 //  }

   627 //

   628 //  // NOTE: switch below is eliminated, because prec is always 2 for doubles

   629 //  /* compress fq[] into y[] */

   630 //  //switch(prec) {

   631 //  //case 0:

   632 //  //  fw = 0.0;

   633 //  //  for (i=jz;i>=0;i--) fw += fq[i];

   634 //  //  y[0] = (ih==0)? fw: -fw;

   635 //  //  break;

   636 //  //case 1:

   637 //  //case 2:

   638 //    fw = 0.0;

   639 //    for (i=jz;i>=0;i--) fw += fq[i];

   640 //    y[0] = (ih==0)? fw: -fw;

   641 //    fw = fq[0]-fw;

   642 //    for (i=1;i<=jz;i++) fw += fq[i];

   643 //    y[1] = (ih==0)? fw: -fw;

   644 //  //  break;

   645 //  //case 3:       /* painful */

   646 //  //  for (i=jz;i>0;i--) {

   647 //  //    fw      = fq[i-1]+fq[i];

   648 //  // fq[i]  += fq[i-1]-fw;

   649 //  //    fq[i-1] = fw;

   650 //  //  }

   651 //  //  for (i=jz;i>1;i--) {

   652 //  //    fw      = fq[i-1]+fq[i];

   653 //  //    fq[i]  += fq[i-1]-fw;

   654 //  //    fq[i-1] = fw;

   655 //  //  }

   656 //  //  for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];

   657 //  //  if(ih==0) {

   658 //  //    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;

   659 //  //  } else {

   660 //  //    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;

   661 //  //  }

   662 //  //}

   663 //  return n&7;

   664 //}

   665 //

   666 // END __kernel_rem_pio2 PSEUDO CODE

   667 //

   668 // Changes between fdlibm and intrinsic:

   669 //     1. One loop is unrolled and vectorized (see comments in code)

   670 //     2. One loop is split into 2 loops (see comments in code)

   671 //     3. Non-double code is removed(last switch). Sevaral variables became

   672 //         constants because of that (see comments in code)

   673 //     4. Use of jx, which is nx-1 instead of nx

   674 // Assumptions:

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

   676 // Input and output:

   677 //     1. Input: X = r0, jx == nx - 1 == r6, e0 == rscratch1

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

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

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

   681 void MacroAssembler::generate__kernel_rem_pio2(address two_over_pi, address pio2) {

   682   Label Q_DONE, JX_IS_0, JX_IS_2, COMP_INNER_LOOP, RECOMP_FOR2, Q0_ZERO_CMP_LT,

   683       RECOMP_CHECK_DONE_NOT_ZERO, Q0_ZERO_CMP_DONE, COMP_FOR, Q0_ZERO_CMP_EQ,

   684       INIT_F_ZERO, RECOMPUTE, IH_FOR_INCREMENT, IH_FOR_STORE, RECOMP_CHECK_DONE,

   685       Z_IS_LESS_THAN_TWO24B, Z_IS_ZERO, FW_Y1_NO_NEGATION,

   686       RECOMP_FW_UPDATED, Z_ZERO_CHECK_DONE, FW_FOR1, IH_AFTER_SWITCH, IH_HANDLED,

   687       CONVERTION_FOR, FW_Y0_NO_NEGATION, FW_FOR1_DONE, FW_FOR2, FW_FOR2_DONE,

   688       IH_FOR, SKIP_F_LOAD, RECOMP_FOR1, RECOMP_FIRST_FOR, INIT_F_COPY,

   689       RECOMP_FOR1_CHECK;

   690   Register tmp2 = r1, n = r2, jv = r4, tmp5 = r5, jx = r6,

   691       tmp3 = r7, iqBase = r10, ih = r11, tmp4 = r12, tmp1 = r13,

   692       jz = r14, j = r15, twoOverPiBase = r16, i = r17, qBase = r18;

   693     // jp = jk == init_jk[prec] = init_jk[2] == {2,3,4,6}[2] == 4

   694     // jx = nx - 1

   695     lea(twoOverPiBase, ExternalAddress(two_over_pi));

   696     cmpw(jv, zr);

   697     addw(tmp4, jx, 4); // tmp4 = m = jx + jk = jx + 4. jx is in {0,1,2} so m is in [4,5,6]

   698     cselw(jv, jv, zr, GE);

   699     fmovd(v26, 0.0d);

   700     addw(tmp5, jv, 1);                    // jv+1

   701     subsw(j, jv, jx);

   702     add(qBase, sp, 320);                  // base of q[]

   703     msubw(rscratch1, i, tmp5, rscratch1); // q0 =  e0-24*(jv+1)

   704     // use double f[20], fq[20], q[20], iq[20] on stack, which is

   705     // (20 + 20 + 20) x 8 + 20 x 4 = 560 bytes. From lower to upper addresses it

   706     // will contain f[20], fq[20], q[20], iq[20]

   707     // now initialize f[20] indexes 0..m (inclusive)

   708     // for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : /* NOTE: converted to double */ ipio2[j]; // (double) ipio2[j];

   709     mov(tmp5, sp);

   710 

   711     block_comment("for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : /* NOTE: converted to double */ ipio2[j]; // (double) ipio2[j];"); {

   712         eorw(i, i, i);

   713         br(GE, INIT_F_COPY);

   714       bind(INIT_F_ZERO);

   715         stpq(v26, v26, Address(post(tmp5, 32)));

   716         addw(i, i, 4);

   717         addsw(j, j, 4);

   718         br(LT, INIT_F_ZERO);

   719         subw(i, i, j);

   720         movw(j, zr);

   721       bind(INIT_F_COPY);

   722         add(tmp1, twoOverPiBase, j, LSL, 3); // ipio2[j] start address

   723         ld1(v18, v19, v20, v21, T16B, tmp1);

   724         add(tmp5, sp, i, ext::uxtx, 3);

   725         st1(v18, v19, v20, v21, T16B, tmp5);

   726     }

   727     // v18..v21 can actually contain f[0..7]

   728     cbz(i, SKIP_F_LOAD); // i == 0 => f[i] == f[0] => already loaded

   729     ld1(v18, v19, v20, v21, T2D, Address(sp)); // load f[0..7]

   730   bind(SKIP_F_LOAD);

   731     // calculate 2^q0 and 2^-q0, which we'll need further.

   732     // q0 is exponent. So, calculate biased exponent(q0+1023)

   733     negw(tmp4, rscratch1);

   734     addw(tmp5, rscratch1, 1023);

   735     addw(tmp4, tmp4, 1023);

   736     // Unroll following for(s) depending on jx in [0,1,2]

   737     // for (i=0;i<=jk;i++) {

   738     //   for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;

   739     // }

   740     // Unrolling for jx == 0 case:

   741     //   q[0] = x[0] * f[0]

   742     //   q[1] = x[0] * f[1]

   743     //   q[2] = x[0] * f[2]

   744     //   q[3] = x[0] * f[3]

   745     //   q[4] = x[0] * f[4]

   746     //

   747     // Vectorization for unrolled jx == 0 case:

   748     //   {q[0], q[1]} = {f[0], f[1]} * x[0]

   749     //   {q[2], q[3]} = {f[2], f[3]} * x[0]

   750     //   q[4] = f[4] * x[0]

   751     //

   752     // Unrolling for jx == 1 case:

   753     //   q[0] = x[0] * f[1] + x[1] * f[0]

   754     //   q[1] = x[0] * f[2] + x[1] * f[1]

   755     //   q[2] = x[0] * f[3] + x[1] * f[2]

   756     //   q[3] = x[0] * f[4] + x[1] * f[3]

   757     //   q[4] = x[0] * f[5] + x[1] * f[4]

   758     //

   759     // Vectorization for unrolled jx == 1 case:

   760     //   {q[0], q[1]} = {f[0], f[1]} * x[1]

   761     //   {q[2], q[3]} = {f[2], f[3]} * x[1]

   762     //   q[4] = f[4] * x[1]

   763     //   {q[0], q[1]} += {f[1], f[2]} * x[0]

   764     //   {q[2], q[3]} += {f[3], f[4]} * x[0]

   765     //   q[4] += f[5] * x[0]

   766     //

   767     // Unrolling for jx == 2 case:

   768     //   q[0] = x[0] * f[2] + x[1] * f[1] + x[2] * f[0]

   769     //   q[1] = x[0] * f[3] + x[1] * f[2] + x[2] * f[1]

   770     //   q[2] = x[0] * f[4] + x[1] * f[3] + x[2] * f[2]

   771     //   q[3] = x[0] * f[5] + x[1] * f[4] + x[2] * f[3]

   772     //   q[4] = x[0] * f[6] + x[1] * f[5] + x[2] * f[4]

   773     //

   774     // Vectorization for unrolled jx == 2 case:

   775     //   {q[0], q[1]} = {f[0], f[1]} * x[2]

   776     //   {q[2], q[3]} = {f[2], f[3]} * x[2]

   777     //   q[4] = f[4] * x[2]

   778     //   {q[0], q[1]} += {f[1], f[2]} * x[1]

   779     //   {q[2], q[3]} += {f[3], f[4]} * x[1]

   780     //   q[4] += f[5] * x[1]

   781     //   {q[0], q[1]} += {f[2], f[3]} * x[0]

   782     //   {q[2], q[3]} += {f[4], f[5]} * x[0]

   783     //   q[4] += f[6] * x[0]

   784   block_comment("unrolled and vectorized computation of q[0]..q[jk]"); {

   785       cmpw(jx, 1);

   786       lsl(tmp5, tmp5, 52);                     // now it's 2^q0 double value

   787       lsl(tmp4, tmp4, 52);                     // now it's 2^-q0 double value

   788       br(LT, JX_IS_0);

   789       add(i, sp, 8);

   790       ldpq(v26, v27, i);                       // load f[1..4]

   791       br(GT, JX_IS_2);

   792       // jx == 1

   793       fmulxvs(v28, T2D, v18, v7);              // f[0,1] * x[1]

   794       fmulxvs(v29, T2D, v19, v7);              // f[2,3] * x[1]

   795       fmuld(v30, v20, v7);                     // f[4] * x[1]

   796       fmlavs(v28, T2D, v26, v6, 0);

   797       fmlavs(v29, T2D, v27, v6, 0);

   798       fmlavs(v30, T2D, v6, v20, 1);            // v30 += f[5] * x[0]

   799       b(Q_DONE);

   800     bind(JX_IS_2);

   801       fmulxvs(v28, T2D, v18, v3);              // f[0,1] * x[2]

   802       fmulxvs(v29, T2D, v19, v3);              // f[2,3] * x[2]

   803       fmuld(v30, v20, v3);                     // f[4] * x[2]

   804       fmlavs(v28, T2D, v26, v7, 0);

   805       fmlavs(v29, T2D, v27, v7, 0);

   806       fmlavs(v30, T2D, v7, v20, 1);            // v30 += f[5] * x[1]

   807       fmlavs(v28, T2D, v19, v6, 0);

   808       fmlavs(v29, T2D, v20, v6, 0);

   809       fmlavs(v30, T2D, v6, v21, 0);            // v30 += f[6] * x[0]

   810       b(Q_DONE);

   811     bind(JX_IS_0);

   812       fmulxvs(v28, T2D, v18, v6);              // f[0,1] * x[0]

   813       fmulxvs(v29, T2D, v19, v6);              // f[2,3] * x[0]

   814       fmuld(v30, v20, v6);                     // f[4] * x[0]

   815     bind(Q_DONE);

   816       st1(v28, v29, v30, T2D, Address(qBase)); // save calculated q[0]...q[jk]

   817   }

   818   movz(i, 0x3E70, 48);

   819   movw(jz, 4);

   820   fmovd(v17, i);                               // v17 = twon24

   821   fmovd(v30, tmp5);                            // 2^q0

   822   fmovd(v21, 0.125d);

   823   fmovd(v20, 8.0d);

   824   fmovd(v22, tmp4);                            // 2^-q0

   825 

   826   block_comment("recompute loop"); {

   827     bind(RECOMPUTE);

   828       //  for(i=0,j=jz,z=q[jz];j>0;i++,j--) {

   829       //    fw    =  (double)((int)(twon24* z));

   830       //    iq[i] =  (int)(z-two24A*fw);

   831       //    z     =  q[j-1]+fw;

   832       //  }

   833       block_comment("distill q[] into iq[] reversingly"); {

   834           eorw(i, i, i);

   835           movw(j, jz);

   836           add(tmp2, qBase, jz, LSL, 3);                    // q[jz] address

   837           ldrd(v18, post(tmp2, -8));                       // z = q[j] and moving address to q[j-1]

   838         bind(RECOMP_FIRST_FOR);

   839           ldrd(v27, post(tmp2, -8));

   840           fmuld(v29, v17, v18);                            // twon24*z

   841           frintzd(v29, v29);                               // (double)(int)

   842           fmsubd(v28, v10, v29, v18);                      // v28 = z-two24A*fw

   843           fcvtzdw(tmp1, v28);                              // (int)(z-two24A*fw)

   844           strw(tmp1, Address(iqBase, i, Address::lsl(2)));

   845           faddd(v18, v27, v29);

   846           add(i, i, 1);

   847           subs(j, j, 1);

   848           br(GT, RECOMP_FIRST_FOR);

   849       }

   850       // compute n

   851       fmuld(v18, v18, v30);

   852       fmuld(v2, v18, v21);

   853       frintmd(v2, v2);                                     // v2 = floor(v2) == rounding towards -inf

   854       fmsubd(v18, v2, v20, v18);                           // z -= 8.0*floor(z*0.125);

   855       movw(ih, 2);

   856       frintzd(v2, v18);                                    // v2 = (double)((int)z)

   857       fcvtzdw(n, v18);                                     // n  = (int) z;

   858       fsubd(v18, v18, v2);                                 // z -= (double)n;

   859 

   860       block_comment("q0-dependent initialization"); {

   861           cmpw(rscratch1, 0);                              // if (q0 > 0)

   862           br(LT, Q0_ZERO_CMP_LT);

   863           subw(j, jz, 1);                                  // j = jz - 1

   864           ldrw(tmp2, Address(iqBase, j, Address::lsl(2))); // tmp2 = iq[jz-1]

   865           br(EQ, Q0_ZERO_CMP_EQ);

   866           movw(tmp4, 24);

   867           subw(tmp4, tmp4, rscratch1);                     // == 24 - q0

   868           lsrvw(i, tmp2, tmp4);                            // i = iq[jz-1] >> (24-q0)

   869           lslvw(tmp5, i, tmp4);

   870           subw(tmp2, tmp2, tmp5);                          // iq[jz-1] -= i<<(24-q0);

   871           strw(tmp2, Address(iqBase, j, Address::lsl(2))); // store iq[jz-1]

   872           subw(rscratch2, tmp4, 1);                        // == 23 - q0

   873           addw(n, n, i);                                   // n+=i

   874           lsrvw(ih, tmp2, rscratch2);                      // ih = iq[jz-1] >> (23-q0)

   875           b(Q0_ZERO_CMP_DONE);

   876         bind(Q0_ZERO_CMP_EQ);

   877           lsr(ih, tmp2, 23);                               // ih = iq[z-1] >> 23

   878           b(Q0_ZERO_CMP_DONE);

   879         bind(Q0_ZERO_CMP_LT);

   880           fmovd(v4, 0.5d);

   881           fcmpd(v18, v4);

   882           cselw(ih, zr, ih, LT);                           // if (z<0.5) ih = 0

   883       }

   884     bind(Q0_ZERO_CMP_DONE);

   885       cmpw(ih, zr);

   886       br(LE, IH_HANDLED);

   887 

   888     block_comment("if(ih>) {"); {

   889       // use rscratch2 as carry

   890 

   891       block_comment("for(i=0;i<jz ;i++) {...}"); {

   892           addw(n, n, 1);

   893           eorw(i, i, i);

   894           eorw(rscratch2, rscratch2, rscratch2);

   895         bind(IH_FOR);

   896           ldrw(j, Address(iqBase, i, Address::lsl(2)));    // j = iq[i]

   897           movw(tmp3, 0x1000000);

   898           subw(tmp3, tmp3, rscratch2);

   899           cbnzw(rscratch2, IH_FOR_STORE);

   900           cbzw(j, IH_FOR_INCREMENT);

   901           movw(rscratch2, 1);

   902         bind(IH_FOR_STORE);

   903           subw(tmp3, tmp3, j);

   904           strw(tmp3, Address(iqBase, i, Address::lsl(2))); // iq[i] = 0xffffff - j

   905         bind(IH_FOR_INCREMENT);

   906           addw(i, i, 1);

   907           cmpw(i, jz);

   908           br(LT, IH_FOR);

   909       }

   910 

   911       block_comment("if(q0>0) {"); {

   912         cmpw(rscratch1, zr);

   913         br(LE, IH_AFTER_SWITCH);

   914         // tmp3 still has iq[jz-1] value. no need to reload

   915         // now, zero high tmp3 bits (rscratch1 number of bits)

   916         movw(j, -1);

   917         subw(i, jz, 1);                                    // set i to jz-1

   918         lsrv(j, j, rscratch1);

   919         andw(tmp3, tmp3, j, LSR, 8);                       // we have 24-bit-based constants

   920         strw(tmp3, Address(iqBase, i, Address::lsl(2)));   // save iq[jz-1]

   921       }

   922       bind(IH_AFTER_SWITCH);

   923         cmpw(ih, 2);

   924         br(NE, IH_HANDLED);

   925 

   926         block_comment("if(ih==2) {"); {

   927           fmovd(v25, 1.0d);

   928           fsubd(v18, v25, v18);                            // z = one - z;

   929           cbzw(rscratch2, IH_HANDLED);

   930           fsubd(v18, v18, v30);                            // z -= scalbnA(one,q0);

   931         }

   932     }

   933     bind(IH_HANDLED);

   934       // check if recomputation is needed

   935       fcmpd(v18, 0.0d);

   936       br(NE, RECOMP_CHECK_DONE_NOT_ZERO);

   937 

   938       block_comment("if(z==zeroB) {"); {

   939 

   940         block_comment("for (i=jz-1;i>=jk;i--) j |= iq[i];"); {

   941             subw(i, jz, 1);

   942             eorw(j, j, j);

   943             b(RECOMP_FOR1_CHECK);

   944           bind(RECOMP_FOR1);

   945             ldrw(tmp1, Address(iqBase, i, Address::lsl(2)));

   946             orrw(j, j, tmp1);

   947             subw(i, i, 1);

   948           bind(RECOMP_FOR1_CHECK);

   949             cmpw(i, 4);

   950             br(GE, RECOMP_FOR1);

   951         }

   952         cbnzw(j, RECOMP_CHECK_DONE);

   953 

   954         block_comment("if(j==0) {"); {

   955             // for(k=1;iq[jk-k]==0;k++); // let's unroll it. jk == 4. So, read

   956             // iq[3], iq[2], iq[1], iq[0] until non-zero value

   957             ldp(tmp1, tmp3, iqBase);               // iq[0..3]

   958             movw(j, 2);

   959             cmp(tmp3, zr);

   960             csel(tmp1, tmp1, tmp3, EQ);            // set register for further consideration

   961             cselw(j, j, zr, EQ);                   // set initial k. Use j as k

   962             cmp(zr, tmp1, LSR, 32);

   963             addw(i, jz, 1);

   964             csincw(j, j, j, NE);

   965 

   966           block_comment("for(i=jz+1;i<=jz+k;i++) {...}"); {

   967               addw(jz, i, j); // i = jz+1, j = k-1. j+i = jz+k (which is a new jz)

   968             bind(RECOMP_FOR2);

   969               addw(tmp1, jv, i);

   970               ldrd(v29, Address(twoOverPiBase, tmp1, Address::lsl(3)));

   971               addw(tmp2, jx, i);

   972               strd(v29, Address(sp, tmp2, Address::lsl(3)));

   973               // f[jx+i] = /* NOTE: converted to double */ ipio2[jv+i]; //(double) ipio2[jv+i];

   974               // since jx = 0, 1 or 2 we can unroll it:

   975               // for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];

   976               // f[jx+i-j] == (for first iteration) f[jx+i], which is already v29

   977               add(tmp2, sp, tmp2, ext::uxtx, 3); // address of f[jx+i]

   978               ldpd(v4, v5, Address(tmp2, -16)); // load f[jx+i-2] and f[jx+i-1]

   979               fmuld(v26, v6, v29); // initial fw

   980               cbzw(jx, RECOMP_FW_UPDATED);

   981               fmaddd(v26, v7, v5, v26);

   982               cmpw(jx, 1);

   983               br(EQ, RECOMP_FW_UPDATED);

   984               fmaddd(v26, v3, v4, v26);

   985             bind(RECOMP_FW_UPDATED);

   986               strd(v26, Address(qBase, i, Address::lsl(3))); // q[i] = fw;

   987               addw(i, i, 1);

   988               cmpw(i, jz);                                   // jz here is "old jz" + k

   989               br(LE, RECOMP_FOR2);

   990           }

   991             b(RECOMPUTE);

   992         }

   993       }

   994     }

   995     bind(RECOMP_CHECK_DONE);

   996       // chop off zero terms

   997       fcmpd(v18, 0.0d);

   998       br(EQ, Z_IS_ZERO);

   999 

  1000       block_comment("else block of if(z==0.0) {"); {

  1001         bind(RECOMP_CHECK_DONE_NOT_ZERO);

  1002           fmuld(v18, v18, v22);

  1003           fcmpd(v18, v10);                                   // v10 is stil two24A

  1004           br(LT, Z_IS_LESS_THAN_TWO24B);

  1005           fmuld(v1, v18, v17);                               // twon24*z

  1006           frintzd(v1, v1);                                   // v1 = (double)(int)(v1)

  1007           fmaddd(v2, v10, v1, v18);

  1008           fcvtzdw(tmp3, v1);                                 // (int)fw

  1009           fcvtzdw(tmp2, v2);                                 // double to int

  1010           strw(tmp2, Address(iqBase, jz, Address::lsl(2)));

  1011           addw(rscratch1, rscratch1, 24);

  1012           addw(jz, jz, 1);

  1013           strw(tmp3, Address(iqBase, jz, Address::lsl(2)));  // iq[jz] = (int) fw

  1014           b(Z_ZERO_CHECK_DONE);

  1015         bind(Z_IS_LESS_THAN_TWO24B);

  1016           fcvtzdw(tmp3, v18);                                // (int)z

  1017           strw(tmp3, Address(iqBase, jz, Address::lsl(2)));  // iq[jz] = (int) z

  1018           b(Z_ZERO_CHECK_DONE);

  1019       }

  1020 

  1021       block_comment("if(z==0.0) {"); {

  1022         bind(Z_IS_ZERO);

  1023           subw(jz, jz, 1);

  1024           ldrw(tmp1, Address(iqBase, jz, Address::lsl(2)));

  1025           subw(rscratch1, rscratch1, 24);

  1026           cbz(tmp1, Z_IS_ZERO);

  1027       }

  1028       bind(Z_ZERO_CHECK_DONE);

  1029         // convert integer "bit" chunk to floating-point value

  1030         // v17 = twon24

  1031         // update v30, which was scalbnA(1.0, <old q0>);

  1032         addw(tmp2, rscratch1, 1023); // biased exponent

  1033         lsl(tmp2, tmp2, 52); // put at correct position

  1034         mov(i, jz);

  1035         fmovd(v30, tmp2);

  1036 

  1037         block_comment("for(i=jz;i>=0;i--) {q[i] = fw*(double)iq[i]; fw*=twon24;}"); {

  1038           bind(CONVERTION_FOR);

  1039             ldrw(tmp1, Address(iqBase, i, Address::lsl(2)));

  1040             scvtfwd(v31, tmp1);

  1041             fmuld(v31, v31, v30);

  1042             strd(v31, Address(qBase, i, Address::lsl(3)));

  1043             fmuld(v30, v30, v17);

  1044             subsw(i, i, 1);

  1045             br(GE, CONVERTION_FOR);

  1046         }

  1047         add(rscratch2, sp, 160); // base for fq

  1048         // reusing twoOverPiBase

  1049         lea(twoOverPiBase, ExternalAddress(pio2));

  1050 

  1051       block_comment("compute PIo2[0,...,jp]*q[jz,...,0]. for(i=jz;i>=0;i--) {...}"); {

  1052           movw(i, jz);

  1053           movw(tmp2, zr); // tmp2 will keep jz - i == 0 at start

  1054         bind(COMP_FOR);

  1055           // for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];

  1056           fmovd(v30, 0.0d);

  1057           add(tmp5, qBase, i, LSL, 3); // address of q[i+k] for k==0

  1058           movw(tmp3, 4);

  1059           movw(tmp4, zr);              // used as k

  1060           cmpw(tmp2, 4);

  1061           add(tmp1, qBase, i, LSL, 3); // used as q[i] address

  1062           cselw(tmp3, tmp2, tmp3, LE); // min(jz - i, jp)

  1063 

  1064           block_comment("for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];"); {

  1065             bind(COMP_INNER_LOOP);

  1066               ldrd(v18, Address(tmp1, tmp4, Address::lsl(3)));          // q[i+k]

  1067               ldrd(v19, Address(twoOverPiBase, tmp4, Address::lsl(3))); // PIo2[k]

  1068               fmaddd(v30, v18, v19, v30);                               // fw += PIo2[k]*q[i+k];

  1069               addw(tmp4, tmp4, 1);                                      // k++

  1070               cmpw(tmp4, tmp3);

  1071               br(LE, COMP_INNER_LOOP);

  1072           }

  1073           strd(v30, Address(rscratch2, tmp2, Address::lsl(3)));         // fq[jz-i]

  1074           add(tmp2, tmp2, 1);

  1075           subsw(i, i, 1);

  1076           br(GE, COMP_FOR);

  1077       }

  1078 

  1079       block_comment("switch(prec) {...}. case 2:"); {

  1080         // compress fq into y[]

  1081         // remember prec == 2

  1082 

  1083         block_comment("for (i=jz;i>=0;i--) fw += fq[i];"); {

  1084             fmovd(v4, 0.0d);

  1085             mov(i, jz);

  1086           bind(FW_FOR1);

  1087             ldrd(v1, Address(rscratch2, i, Address::lsl(3)));

  1088             subsw(i, i, 1);

  1089             faddd(v4, v4, v1);

  1090             br(GE, FW_FOR1);

  1091         }

  1092         bind(FW_FOR1_DONE);

  1093           // v1 contains fq[0]. so, keep it so far

  1094           fsubd(v5, v1, v4); // fw = fq[0] - fw

  1095           cbzw(ih, FW_Y0_NO_NEGATION);

  1096           fnegd(v4, v4);

  1097         bind(FW_Y0_NO_NEGATION);

  1098 

  1099         block_comment("for (i=1;i<=jz;i++) fw += fq[i];"); {

  1100             movw(i, 1);

  1101               cmpw(jz, 1);

  1102             br(LT, FW_FOR2_DONE);

  1103           bind(FW_FOR2);

  1104             ldrd(v1, Address(rscratch2, i, Address::lsl(3)));

  1105             addw(i, i, 1);

  1106             cmp(i, jz);

  1107             faddd(v5, v5, v1);

  1108             br(LE, FW_FOR2);

  1109         }

  1110         bind(FW_FOR2_DONE);

  1111           cbz(ih, FW_Y1_NO_NEGATION);

  1112           fnegd(v5, v5);

  1113         bind(FW_Y1_NO_NEGATION);

  1114           add(sp, sp, 560);

  1115       }

  1116 }

  1117 

  1118 ///* __kernel_sin( x, y, iy)

  1119 // * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854

  1120 // * Input x is assumed to be bounded by ~pi/4 in magnitude.

  1121 // * Input y is the tail of x.

  1122 // * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).

  1123 // *

  1124 // * Algorithm

  1125 // *      1. Since sin(-x) = -sin(x), we need only to consider positive x.

  1126 // *      2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.

  1127 // *      3. sin(x) is approximated by a polynomial of degree 13 on

  1128 // *         [0,pi/4]

  1129 // *                               3            13

  1130 // *              sin(x) ~ x + S1*x + ... + S6*x

  1131 // *         where

  1132 // *

  1133 // *      |sin(x)         2     4     6     8     10     12  |     -58

  1134 // *      |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2

  1135 // *      |  x                                               |