asin-td.c

/*
 * Correctly rounded arcsine
 *
 * Author : Christoph Lauter (ENS Lyon)
 *
 * This file is part of the crlibm library developed by the Arenaire
 * project at Ecole Normale Superieure de Lyon
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or 
 * (at your option) any later version.
 * 
 * This program 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 for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
*/

#include <stdio.h>
#include <stdlib.h>
#include "crlibm.h"
#include "crlibm_private.h"
#include "triple-double.h"
#include "asin-td.h"

#define AVOID_FMA 1

void asin_accurate_lower(double *asinh, double *asinm, double *asinl, double x, double xSqh, double xSql, double sign) {
  double highPoly;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l, t7h, t7l;
  double tt1h, tt1l;
  double t8h, t8m, t8l, t9h, t9m, t9l, t10h, t10m, t10l, t11h, t11m, t11l, t12h, t12m, t12l;
  double tt8h, tt8m, tt8l, tt9h, tt9m, tt9l, tt10h, tt10m, tt10l, tt11h, tt11m, tt11l, tt12h, tt12m, tt12l;
  double xCubeh, xCubem, xCubel, tt13h, tt13m, tt13l, t13h, t13m, t13l, polyh, polym, polyl;
  double tt11hover, tt11mover, tt11lover;

#if EVAL_PERF
  crlibm_second_step_taken++;
#endif

  /* Evaluate the polynomial of degree 37
     Its coefficients start at tbl[0]

     p(x) = x + x * x^2 * (c3 + x^2 * (c5 + ...

     We receive x^2 as xSqh + xSql = x * x (exactly)
     in argument

     |x| <= 0.185 = 2^(-2.43)

     Compute monomials 27 to 37 in double precision
     monomials 13 to 25 in double-double and
     1 to 11 in triple-double precision in a 
     modified Horner form

  */

  /* Double computations */
  
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(tbl[33],xSqh,tbl[32]),xSqh,tbl[31]),xSqh,tbl[30]),xSqh,tbl[29]),xSqh,tbl[28]);		 
#else
  highPoly = tbl[28] + xSqh * (tbl[29] + xSqh * (tbl[30] + xSqh * (tbl[31] + xSqh * (tbl[32] + xSqh * tbl[33]))));
#endif

  /* Double-double computations */

  Mul12(&tt1h,&tt1l,xSqh,highPoly);
  Add22(&t1h,&t1l,tbl[27],0,tt1h,tt1l);

  MulAdd22(&t2h,&t2l,tbl[25],tbl[26],xSqh,xSql,t1h,t1l);
  MulAdd22(&t3h,&t3l,tbl[23],tbl[24],xSqh,xSql,t2h,t2l);
  MulAdd22(&t4h,&t4l,tbl[21],tbl[22],xSqh,xSql,t3h,t3l);
  MulAdd22(&t5h,&t5l,tbl[19],tbl[20],xSqh,xSql,t4h,t4l);
  MulAdd22(&t6h,&t6l,tbl[17],tbl[18],xSqh,xSql,t5h,t5l);
  MulAdd22(&t7h,&t7l,tbl[15],tbl[16],xSqh,xSql,t6h,t6l);

  /* Triple-double computations */

  Mul23(&tt8h,&tt8m,&tt8l,xSqh,xSql,t7h,t7l);                             /* 149 - 48/53 */
  Add33(&t8h,&t8m,&t8l,tbl[12],tbl[13],tbl[14],tt8h,tt8m,tt8l);           /* 145 - 43/53 */
  Mul233(&tt9h,&tt9m,&tt9l,xSqh,xSql,t8h,t8m,t8l);                        /* 139 - 39/53 */
  Add33(&t9h,&t9m,&t9l,tbl[9],tbl[10],tbl[11],tt9h,tt9m,tt9l);            /* 136 - 34/53 */
  Mul233(&tt10h,&tt10m,&tt10l,xSqh,xSql,t9h,t9m,t9l);                     /* 130 - 30/53 */
  Add33(&t10h,&t10m,&t10l,tbl[6],tbl[7],tbl[8],tt10h,tt10m,tt10l);        /* 127 - 25/53 */
  Mul233(&tt11hover,&tt11mover,&tt11lover,xSqh,xSql,t10h,t10m,t10l);      /* 121 - 21/53 */

  Renormalize3(&tt11h,&tt11m,&tt11l,tt11hover,tt11mover,tt11lover);       /* infty - 52/53 */

  Add33(&t11h,&t11m,&t11l,tbl[3],tbl[4],tbl[5],tt11h,tt11m,tt11l);        /* 149 - 47/53 */
  Mul233(&tt12h,&tt12m,&tt12l,xSqh,xSql,t11h,t11m,t11l);                  /* 143 - 43/53 */
  Add33(&t12h,&t12m,&t12l,tbl[0],tbl[1],tbl[2],tt12h,tt12m,tt12l);        /* 140 - 38/53 */

  Mul123(&xCubeh,&xCubem,&xCubel,x,xSqh,xSql);                            /* 154 - 47/53 */
 
  Mul33(&tt13h,&tt13m,&tt13l,xCubeh,xCubem,xCubel,t12h,t12m,t12l);        /* 136 - 34/53 */
  Add133(&t13h,&t13m,&t13l,x,tt13h,tt13m,tt13l);                          /* 138 - 32/53 */

  Renormalize3(&polyh,&polym,&polyl,t13h,t13m,t13l);                      /* infty - 52/53 */
  *asinh = sign * polyh;
  *asinm = sign * polym;
  *asinl = sign * polyl;
}


void  asin_accurate_middle(double *asinh, double *asinm, double *asinl, double z, int i, double sign) {
  double highPoly;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l, t7h, t7l, t8h, t8l, t9h, t9l;
  double t10h, t10m, t10l, t11h, t11m, t11l, t12h, t12m, t12l, t13h, t13m, t13l, t14h, t14m, t14l;
  double t15h, t15m, t15l, t16h, t16m, t16l;
  double tt1h, tt1l;
  double tt10h, tt10m, tt10l, tt11h, tt11m, tt11l, tt12h, tt12m, tt12l;
  double tt13h, tt13m, tt13l, tt14h, tt14m, tt14l, tt15h, tt15m, tt15l, tt16h, tt16m, tt16l;
  double polyh, polym, polyl, tt13hover, tt13mover, tt13lover;

#if EVAL_PERF
  crlibm_second_step_taken++;
#endif

  /* Evaluate the polynomial of degree 35
     Its coefficients start at tbl[i+1] 
     Evaluate degrees 35 to 20 in double precision,
     degrees 20 to 7 in double-double precision and
     finally degrees 6 to 1 in triple-double. 
     The constant coefficient is a double-double, the 
     computations are nevertheless in triple-double 
  */

  /* Double computations */
  
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(
             tbl[i+58] ,z,tbl[i+57]),z,tbl[i+56]),z,tbl[i+55]),z,tbl[i+54]),z,
             tbl[i+53]),z,tbl[i+52]),z,tbl[i+51]),z,tbl[i+50]),z,tbl[i+49]),z,
             tbl[i+48]),z,tbl[i+47]),z,tbl[i+46]),z,tbl[i+45]),z,tbl[i+44]),z,
	     tbl[i+43]),z,tbl[i+42]),z,tbl[i+41]),z,tbl[i+40]),z,tbl[i+39]);
								 
#else
  highPoly = tbl[i+39] + z * (tbl[i+40] + z * (tbl[i+41] + z * (tbl[i+42] + z * (
             tbl[i+43] + z * (tbl[i+44] + z * (tbl[i+45] + z * (tbl[i+46] + z * (
             tbl[i+47] + z * (tbl[i+48] + z * (tbl[i+49] + z * (tbl[i+50] + z * (
             tbl[i+51] + z * (tbl[i+52] + z * (tbl[i+53] + z * (tbl[i+54] + z * (
             tbl[i+55] + z * (tbl[i+56] + z * (tbl[i+57] + z * tbl[i+58]))))))))))))))))));
#endif

  
  /* Double-double computations */
  
  Mul12(&tt1h,&tt1l,z,highPoly);
  Add22(&t1h,&t1l,tbl[i+37],tbl[i+38],tt1h,tt1l);

  MulAdd212(&t2h,&t2l,tbl[i+35],tbl[i+36],z,t1h,t1l);
  MulAdd212(&t3h,&t3l,tbl[i+33],tbl[i+34],z,t2h,t2l);
  MulAdd212(&t4h,&t4l,tbl[i+31],tbl[i+32],z,t3h,t3l);
  MulAdd212(&t5h,&t5l,tbl[i+29],tbl[i+30],z,t4h,t4l);
  MulAdd212(&t6h,&t6l,tbl[i+27],tbl[i+28],z,t5h,t5l);
  MulAdd212(&t7h,&t7l,tbl[i+25],tbl[i+26],z,t6h,t6l);
  MulAdd212(&t8h,&t8l,tbl[i+23],tbl[i+24],z,t7h,t7l);
  MulAdd212(&t9h,&t9l,tbl[i+21],tbl[i+22],z,t8h,t8l);

  /* Triple-double computations */

  Mul123(&tt10h,&tt10m,&tt10l,z,t9h,t9l);                                          /* 154 - 47/53 */
  Add33(&t10h,&t10m,&t10l,tbl[i+18],tbl[i+19],tbl[i+20],tt10h,tt10m,tt10l);        /* 144 - 42/53 */
  Mul133(&tt11h,&tt11m,&tt11l,z,t10h,t10m,t10l);                                   /* 142 - 38/53 */
  Add33(&t11h,&t11m,&t11l,tbl[i+15],tbl[i+16],tbl[i+17],tt11h,tt11m,tt11l);        /* 136 - 33/53 */
  Mul133(&tt12h,&tt12m,&tt12l,z,t11h,t11m,t11l);                                   /* 133 - 28/53 */
  Add33(&t12h,&t12m,&t12l,tbl[i+12],tbl[i+13],tbl[i+14],tt12h,tt12m,tt12l);        /* 125 - 23/53 */
  Mul133(&tt13hover,&tt13mover,&tt13lover,z,t12h,t12m,t12l);                       /* 123 - 18/53 */

  Renormalize3(&tt13h,&tt13m,&tt13l,tt13hover,tt13mover,tt13lover);                /* infty - 52/53 */

  Add33(&t13h,&t13m,&t13l,tbl[i+9],tbl[i+10],tbl[i+11],tt13h,tt13m,tt13l);         /* 149 - 47/53 */
  Mul133(&tt14h,&tt14m,&tt14l,z,t13h,t13m,t13l);                                   /* 147 - 42/53 */
  Add33(&t14h,&t14m,&t14l,tbl[i+6],tbl[i+7],tbl[i+8],tt14h,tt14m,tt14l);           /* 139 - 37/53 */
  Mul133(&tt15h,&tt15m,&tt15l,z,t14h,t14m,t14l);                                   /* 137 - 32/53 */
  Add33(&t15h,&t15m,&t15l,tbl[i+3],tbl[i+4],tbl[i+5],tt15h,tt15m,tt15l);           /* 129 - 28/53 */
  Mul133(&tt16h,&tt16m,&tt16l,z,t15h,t15m,t15l);                                   /* 128 - 23/53 */
  Add233(&t16h,&t16m,&t16l,tbl[i+1],tbl[i+2],tt16h,tt16m,tt16l);                   /* 126 - 19/53 */

  Renormalize3(&polyh,&polym,&polyl,t16h,t16m,t16l);                               /* infty - 52/53 */
  *asinh = sign * polyh;
  *asinm = sign * polym;
  *asinl = sign * polyl;
}


void asin_accurate_higher(double *asinh, double *asinm, double *asinl, double z, double sign) {
  double highPoly;
  double tt1h, tt1l;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l, t7h, t7l, t8h, t8l;
  double tt10h, tt10m, tt10l, tt11h, tt11m, tt11l, tt12h, tt12m, tt12l, tt13h, tt13m, tt13l;
  double tt14h, tt14m, tt14l, tt15h, tt15m, tt15l, tt16h, tt16m, tt16l, tt17h, tt17m, tt17l;
  double t9h, t9l, t10h, t10m, t10l, t11h, t11m, t11l, t12h, t12m, t12l, t13h, t13m, t13l;
  double t14h, t14m, t14l, t15h, t15m, t15l, t16h, t16m, t16l, t17h, t17m, t17l;
  double tt18h, tt18m, tt18l, polyh, polym, polyl;
  double sqrtzh, sqrtzm, sqrtzl, twoZ, pTimesSh, pTimesSm, pTimesSl;
  double allhover, allmover, alllover, allh, allm, alll;
  double tt13hover, tt13mover, tt13lover, tt16hover, tt16mover, tt16lover;

#if EVAL_PERF
  crlibm_second_step_taken++;
#endif

  /* We evaluate asin(x) as 

     asin(x) = f(z) * sqrt(2*z) + Pi/2

     with z = 1 - x and

     f(z) = (asin(z) - Pi/2) / sqrt(2*z)

     f(z) is approximated by p(z)

     The polynomial p(z) is of degree 29
     Its coefficients start at tbl[TBLIDX10] 
     Coefficients for degrees 29 to 18 are in double precision,
     for degrees 17 to 9 in double-double precision and
     finally for degrees 8 to 1 in triple-double. 
     The constant coefficient (-1) is not stored in the table, 
     the computations are nevertheless in triple-double 
     We evaluate the monomials in the precision in which
     the correspondant coefficients are stored
     The coefficients' values decrease very quickly 
     so even with |z| < 2^-2.18 we can compute degree 18 
     already in double precision

     Compute than sqrt(2*z) as a triple-double
     multiply in triple-double and add Pi/2
     We will cancel no bit in the addition since
     f(z) < 0.5 * Pi/2

  */

  /* Double computations */

#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(
             tbl[TBLIDX10+53] ,z,tbl[TBLIDX10+52]),z,tbl[TBLIDX10+51]),z,
             tbl[TBLIDX10+50]),z,tbl[TBLIDX10+49]),z,tbl[TBLIDX10+48]),z,
             tbl[TBLIDX10+47]),z,tbl[TBLIDX10+46]),z,tbl[TBLIDX10+45]),z,
             tbl[TBLIDX10+44]),z,tbl[TBLIDX10+43]),z,tbl[TBLIDX10+42]);						 
#else
  highPoly = tbl[TBLIDX10+42] + z * (tbl[TBLIDX10+43] + z * (tbl[TBLIDX10+44] + z * (
             tbl[TBLIDX10+45] + z * (tbl[TBLIDX10+46] + z * (tbl[TBLIDX10+47] + z * (
             tbl[TBLIDX10+48] + z * (tbl[TBLIDX10+49] + z * (tbl[TBLIDX10+50] + z * (
             tbl[TBLIDX10+51] + z * (tbl[TBLIDX10+52] + z *  tbl[TBLIDX10+53]))))))))));
#endif
  
  /* Double-double computations */

  Mul12(&tt1h,&tt1l,z,highPoly);
  Add22(&t1h,&t1l,tbl[TBLIDX10+40],tbl[TBLIDX10+41],tt1h,tt1l);

  MulAdd212(&t2h,&t2l,tbl[TBLIDX10+38],tbl[TBLIDX10+39],z,t1h,t1l);
  MulAdd212(&t3h,&t3l,tbl[TBLIDX10+36],tbl[TBLIDX10+37],z,t2h,t2l);
  MulAdd212(&t4h,&t4l,tbl[TBLIDX10+34],tbl[TBLIDX10+35],z,t3h,t3l);
  MulAdd212(&t5h,&t5l,tbl[TBLIDX10+32],tbl[TBLIDX10+33],z,t4h,t4l);
  MulAdd212(&t6h,&t6l,tbl[TBLIDX10+30],tbl[TBLIDX10+31],z,t5h,t5l);
  MulAdd212(&t7h,&t7l,tbl[TBLIDX10+28],tbl[TBLIDX10+29],z,t6h,t6l);
  MulAdd212(&t8h,&t8l,tbl[TBLIDX10+26],tbl[TBLIDX10+27],z,t7h,t7l);
  MulAdd212(&t9h,&t9l,tbl[TBLIDX10+24],tbl[TBLIDX10+25],z,t8h,t8l);

  /* Triple-double computations */

  Mul123(&tt10h,&tt10m,&tt10l,z,t9h,t9l);                                                         /* 154 - 47/53 */
  Add33(&t10h,&t10m,&t10l,tbl[TBLIDX10+21],tbl[TBLIDX10+22],tbl[TBLIDX10+23],tt10h,tt10m,tt10l);  /* 144 - 42/53 */
  Mul133(&tt11h,&tt11m,&tt11l,z,t10h,t10m,t10l);                                                  /* 142 - 37/53 */
  Add33(&t11h,&t11m,&t11l,tbl[TBLIDX10+18],tbl[TBLIDX10+19],tbl[TBLIDX10+20],tt11h,tt11m,tt11l);  /* 134 - 32/53 */
  Mul133(&tt12h,&tt12m,&tt12l,z,t11h,t11m,t11l);                                                  /* 132 - 27/53 */
  Add33(&t12h,&t12m,&t12l,tbl[TBLIDX10+15],tbl[TBLIDX10+16],tbl[TBLIDX10+17],tt12h,tt12m,tt12l);  /* 124 - 22/53 */
  Mul133(&tt13hover,&tt13mover,&tt13lover,z,t12h,t12m,t12l);                                      /* 122 - 17/53 */

  Renormalize3(&tt13h,&tt13m,&tt13l,tt13hover,tt13mover,tt13lover);                               /* infty - 52/53 */

  Add33(&t13h,&t13m,&t13l,tbl[TBLIDX10+12],tbl[TBLIDX10+13],tbl[TBLIDX10+14],tt13h,tt13m,tt13l);  /* 149 - 47/53 */
  Mul133(&tt14h,&tt14m,&tt14l,z,t13h,t13m,t13l);                                                  /* 147 - 42/53 */
  Add33(&t14h,&t14m,&t14l,tbl[TBLIDX10+9],tbl[TBLIDX10+10],tbl[TBLIDX10+11],tt14h,tt14m,tt14l);   /* 139 - 37/53 */
  Mul133(&tt15h,&tt15m,&tt15l,z,t14h,t14m,t14l);                                                  /* 137 - 32/53 */
  Add33(&t15h,&t15m,&t15l,tbl[TBLIDX10+6],tbl[TBLIDX10+7],tbl[TBLIDX10+8],tt15h,tt15m,tt15l);     /* 129 - 27/53 */
  Mul133(&tt16hover,&tt16mover,&tt16lover,z,t15h,t15m,t15l);                                      /* 127 - 22/53 */

  Renormalize3(&tt16h,&tt16m,&tt16l,tt16hover,tt16mover,tt16lover);                               /* infty - 52/53 */

  Add33(&t16h,&t16m,&t16l,tbl[TBLIDX10+3],tbl[TBLIDX10+4],tbl[TBLIDX10+5],tt16h,tt16m,tt16l);     /* 149 - 47/53 */
  Mul133(&tt17h,&tt17m,&tt17l,z,t16h,t16m,t16l);                                                  /* 147 - 42/53 */
  Add33(&t17h,&t17m,&t17l,tbl[TBLIDX10+0],tbl[TBLIDX10+1],tbl[TBLIDX10+2],tt17h,tt17m,tt17l);     /* 139 - 37/53 */
  Mul133(&tt18h,&tt18m,&tt18l,z,t17h,t17m,t17l);                                                  /* 137 - 32/53 */
  Add133(&polyh,&polym,&polyl,-1,tt18h,tt18m,tt18l);                                              /* 136 - 30/53 */

  /* Compute sqrt(2*z) as a triple-double */

  twoZ = 2 * z;
  Sqrt13(&sqrtzh,&sqrtzm,&sqrtzl,twoZ);                                                           /* 146 - 52/53 */

  /* Multiply p(z) by sqrt(2*z) and add Pi/2 */

  Mul33(&pTimesSh,&pTimesSm,&pTimesSl,polyh,polym,polyl,sqrtzh,sqrtzm,sqrtzl);                    /* 128 - 26/53 */
  Add33(&allhover,&allmover,&alllover,PIHALFH,PIHALFM,PIHALFL,pTimesSh,pTimesSm,pTimesSl);        /* 126 - 21/53 */

  /* Renormalize and multiply by sign */
  Renormalize3(&allh,&allm,&alll,allhover,allmover,alllover);                                     /* infty - 52/53 */
  *asinh = sign * allh;
  *asinm = sign * allm;
  *asinl = sign * alll;
}








double asin_rn(double x) {
  db_number xdb;
  double sign, z, asinh, asinm, asinl;
  int i;
  double xSqh, xSql;
  double tt1h, tt1l;
  double tt6h, tt6l;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l;
  double t7h, t7l, t8h, t8l, polyh, polyl, twoZ, sqrtzh, sqrtzl;
  double pTimesSh, pTimesSl, allh, alll, highPoly, xCubeh, xCubel;
  double tmp1, tmp2, tmp3, tmp4, tmp5;

  /* Transform the argument into integer */
  xdb.d = x;

  /* Special case handling */
  
  /* Strip off the sign of argument x */
  if (xdb.i[HI] & 0x80000000) sign = -1; else sign = 1;
  xdb.i[HI] &= 0x7fffffff;

  /* asin is defined on -1 <= x <= 1, elsewhere it is NaN */
  if ((xdb.i[HI] > 0x3ff00000) || ((xdb.i[HI] == 0x3ff00000) && (xdb.i[LO] != 0x00000000))) {
    return (x-x)/0.0;    /* return NaN */
  }

  /* If |x| < 2^(-28) we have
     
     arcsin(x) = x * ( 1 + xi ) 

     with 0 <= xi < 2^(-55) 
          
     So we can decide the rounding without any computation 
  */
  if (xdb.i[HI] < 0x3e300000) {
    return x;
  }

  /* Recast x */
  x = xdb.d;

  /* Find correspondant interval and compute index to the table
     We start by filtering the two special cases around 0 and 1
  */

  if (xdb.i[HI] < BOUND1) {
    /* Special interval 0..BOUND1 
       The polynomial has no even monomials
       We must prove extra accuracy in the interval 0..sin(2^(-18))
    */

    /* Quick phase starts */

    /* Compute square of x for both quick and accurate phases */
    Mul12(&xSqh,&xSql,x,x);

    tmp4 = tbl[3];
    tmp5 = tbl[4];
    t4h = tmp4;
    t4l = tmp5;
    if (xdb.i[HI] > EXTRABOUND) {
      /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
      highPoly = FMA(FMA(FMA(FMA(tbl[23],xSqh,tbl[21]),xSqh,tbl[19]),xSqh,tbl[17]),xSqh,tbl[15]);
#else
      highPoly = tbl[15] + xSqh * (tbl[17] + xSqh * (tbl[19] + xSqh * (tbl[21] + xSqh * tbl[23])));
#endif

      /* Double-double precision evaluation */
      Mul12(&tt1h,&tt1l,xSqh,highPoly);
      Add22(&t1h,&t1l,tbl[12],tbl[13],tt1h,tt1l);
      
      MulAdd212(&t2h,&t2l,tbl[9],tbl[10],xSqh,t1h,t1l);
      MulAdd212(&t3h,&t3l,tbl[6],tbl[7],xSqh,t2h,t2l);
      MulAdd22(&t4h,&t4l,tmp4,tmp5,xSqh,xSql,t3h,t3l);
    }

    MulAdd22(&t5h,&t5l,tbl[0],tbl[1],xSqh,xSql,t4h,t4l);

    Mul122(&xCubeh,&xCubel,x,xSqh,xSql);
    Mul22(&tt6h,&tt6l,xCubeh,xCubel,t5h,t5l);
    
    Add12(tmp1,tmp2,x,tt6h);
    tmp3 = tmp2 + tt6l;
    Add12(polyh,polyl,tmp1,tmp3);

    /* Multiply by sign */
    asinh = sign * polyh;
    asinm = sign * polyl;

    /* Rounding test (on polyh+polyl, equivalently to asinh+asinm)
       The RN rounding constant is at tbl[34]
    */
    if(polyh == (polyh + (polyl * tbl[34]))) 
      return asinh;

    /* Launch accurate phase */

    asin_accurate_lower(&asinh,&asinm,&asinl,x,xSqh,xSql,sign);

    ReturnRoundToNearest3(asinh,asinm,asinl); 
  }

  if (xdb.i[HI] >= BOUND9) {
    /* Special interval BOUND9..1
       We use an asymptotic development of arcsin in sqrt(1 - x)
    */

    /* Argument reduction for quick and accurate phase
       z = 1 - x
       The operation is exact as per Sterbenz' lemma
    */

    z = 1 - x;

    /* Quick phase starts */

    /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
    highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	       tbl[TBLIDX10+42] ,z,tbl[TBLIDX10+40]),z,tbl[TBLIDX10+38]),z,
               tbl[TBLIDX10+36]),z,tbl[TBLIDX10+34]),z,tbl[TBLIDX10+32]),z,
               tbl[TBLIDX10+30]),z,tbl[TBLIDX10+28]),z,tbl[TBLIDX10+26]),z,
               tbl[TBLIDX10+24]);
#else
    highPoly = tbl[TBLIDX10+24] + z * (tbl[TBLIDX10+26] + z * (tbl[TBLIDX10+28] + z * (
	       tbl[TBLIDX10+30] + z * (tbl[TBLIDX10+32] + z * (tbl[TBLIDX10+34] + z * (
	       tbl[TBLIDX10+36] + z * (tbl[TBLIDX10+38] + z * (tbl[TBLIDX10+40] + z * 
               tbl[TBLIDX10+42]))))))));
#endif
    
    /* Double-double precision evaluation */
    Mul12(&tt1h,&tt1l,z,highPoly);
    Add22(&t1h,&t1l,tbl[TBLIDX10+21],tbl[TBLIDX10+22],tt1h,tt1l);

    MulAdd212(&t2h,&t2l,tbl[TBLIDX10+18],tbl[TBLIDX10+19],z,t1h,t1l);
    MulAdd212(&t3h,&t3l,tbl[TBLIDX10+15],tbl[TBLIDX10+16],z,t2h,t2l);
    MulAdd212(&t4h,&t4l,tbl[TBLIDX10+12],tbl[TBLIDX10+13],z,t3h,t3l);
    MulAdd212(&t5h,&t5l,tbl[TBLIDX10+9],tbl[TBLIDX10+10],z,t4h,t4l);
    MulAdd212(&t6h,&t6l,tbl[TBLIDX10+6],tbl[TBLIDX10+7],z,t5h,t5l);
    MulAdd212(&t7h,&t7l,tbl[TBLIDX10+3],tbl[TBLIDX10+4],z,t6h,t6l);
    MulAdd212(&t8h,&t8l,tbl[TBLIDX10+0],tbl[TBLIDX10+1],z,t7h,t7l);
    MulAdd212(&polyh,&polyl,-1,0,z,t8h,t8l);

    /* Compute sqrt(2*z) as a double-double */

    twoZ = 2 * z;
    sqrt12(&sqrtzh,&sqrtzl,twoZ);                                                         

    /* Multiply p(z) by sqrt(2*z) and add Pi/2 */

    Mul22(&pTimesSh,&pTimesSl,polyh,polyl,sqrtzh,sqrtzl);                    
    Add22(&allh,&alll,PIHALFH,PIHALFM,pTimesSh,pTimesSl);       

    /* Multiply by sign */
    asinh = sign * allh;
    asinm = sign * alll;

    /* Rounding test 
       The RN rounding constant is at tbl[TBLIDX10+54]
    */

    if(allh == (allh + (alll * tbl[TBLIDX10+54]))) 
      return asinh;

    /* Launch accurate phase */

    asin_accurate_higher(&asinh,&asinm,&asinl,z,sign);

    ReturnRoundToNearest3(asinh,asinm,asinl); 
  }
  
  /* General 8 main intervals 
     We can already suppose that BOUND1 <= x <= BOUND9
  */
  
  if (xdb.i[HI] < BOUND5) {
    if (xdb.i[HI] < BOUND3) {
      if (xdb.i[HI] < BOUND2) i = TBLIDX2; else i = TBLIDX3;
    } else {
      if (xdb.i[HI] < BOUND4) i = TBLIDX4; else i = TBLIDX5;
    }
  } else {
    if (xdb.i[HI] < BOUND7) {
      if (xdb.i[HI] < BOUND6) i = TBLIDX6; else i = TBLIDX7;
    } else {
      if (xdb.i[HI] < BOUND8) i = TBLIDX8; else i = TBLIDX9;
    }
  }

  /* Argument reduction 
     i points to the interval midpoint value in the table
  */
  z = x - tbl[i];

  /* Quick phase starts */

  /* Double precision evaluation */

#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	     tbl[i+35] ,z,tbl[i+33]),z,tbl[i+31]),z,tbl[i+29]),z,
             tbl[i+27]),z,tbl[i+25]),z,tbl[i+23]),z,tbl[i+21]);
#else
  highPoly = tbl[i+21] + z * (tbl[i+23] + z * (tbl[i+25] + z * (
             tbl[i+27] + z * (tbl[i+29] + z * (tbl[i+31] + z * ( 
             tbl[i+33] + z *  tbl[i+35]))))));
#endif

  /* Double-double precision evaluation */
    
  Mul12(&tt1h,&tt1l,z,highPoly);
  Add22(&t1h,&t1l,tbl[i+18],tbl[i+19],tt1h,tt1l);

  MulAdd212(&t2h,&t2l,tbl[i+15],tbl[i+16],z,t1h,t1l);
  MulAdd212(&t3h,&t3l,tbl[i+12],tbl[i+13],z,t2h,t2l);
  MulAdd212(&t4h,&t4l,tbl[i+9],tbl[i+10],z,t3h,t3l);
  MulAdd212(&t5h,&t5l,tbl[i+6],tbl[i+7],z,t4h,t4l);
  MulAdd212(&t6h,&t6l,tbl[i+3],tbl[i+4],z,t5h,t5l);
  MulAdd212(&polyh,&polyl,tbl[i+1],tbl[i+2],z,t6h,t6l);

  /* Multiply by sign */
  asinh = sign * polyh;
  asinm = sign * polyl;

  /* Rounding test 
     The RN rounding constant is at tbl[i+59]
  */
    if(polyh == (polyh + (polyl * tbl[i+59]))) 
      return asinh;

  /* Launch accurate phase */

  asin_accurate_middle(&asinh,&asinm,&asinl,z,i,sign);

  ReturnRoundToNearest3(asinh,asinm,asinl); 
}






double asin_ru(double x) {
  db_number xdb;
  double sign, z, asinh, asinm, asinl;
  int i;
  double xSqh, xSql;
  double tt1h, tt1l;
  double tt6h, tt6l;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l;
  double t7h, t7l, t8h, t8l, polyh, polyl, twoZ, sqrtzh, sqrtzl;
  double pTimesSh, pTimesSl, allh, alll, highPoly, xCubeh, xCubel;
  double tmp1, tmp2, tmp3, tmp4, tmp5;

  /* Transform the argument into integer */
  xdb.d = x;

  /* Special case handling */
  
  /* Strip off the sign of argument x */
  if (xdb.i[HI] & 0x80000000) sign = -1; else sign = 1;
  xdb.i[HI] &= 0x7fffffff;

  /* asin is defined on -1 <= x <= 1, elsewhere it is NaN */
  if ((xdb.i[HI] > 0x3ff00000) || ((xdb.i[HI] == 0x3ff00000) && (xdb.i[LO] != 0x00000000))) {
    return (x-x)/0.0;    /* return NaN */
  }

  /* If |x| < 2^(-28) we have
     
     arcsin(x) = x * ( 1 + xi ) 

     with 0 <= xi < 2^(-55) 
          
     So we can decide the rounding without any computation 
  */
  if (xdb.i[HI] < 0x3e300000) {
    /* If x == 0 then we got the algebraic result arcsin(0) = 0
       If x < 0 then the truncation rest is negative but less than 
       1 ulp; we round upwards by returning x
    */
    if (x <= 0) return x;
    /* Otherwise the rest is positive, less than 1 ulp and the
       image is not algebraic 
       We return x + 1ulp
    */
    xdb.l++;
    return xdb.d;
  }

  /* Recast x */
  x = xdb.d;

  /* Find correspondant interval and compute index to the table
     We start by filtering the two special cases around 0 and 1
  */

  if (xdb.i[HI] < BOUND1) {
    /* Special interval 0..BOUND1 
       The polynomial has no even monomials
       We must prove extra accuracy in the interval 0..sin(2^(-18))
    */

    /* Quick phase starts */

    /* Compute square of x for both quick and accurate phases */
    Mul12(&xSqh,&xSql,x,x);

    tmp4 = tbl[3];
    tmp5 = tbl[4];
    t4h = tmp4;
    t4l = tmp5;
    if (xdb.i[HI] > EXTRABOUND) {
      /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
      highPoly = FMA(FMA(FMA(FMA(tbl[23],xSqh,tbl[21]),xSqh,tbl[19]),xSqh,tbl[17]),xSqh,tbl[15]);
#else
      highPoly = tbl[15] + xSqh * (tbl[17] + xSqh * (tbl[19] + xSqh * (tbl[21] + xSqh * tbl[23])));
#endif

      /* Double-double precision evaluation */
      Mul12(&tt1h,&tt1l,xSqh,highPoly);
      Add22(&t1h,&t1l,tbl[12],tbl[13],tt1h,tt1l);
      
      MulAdd212(&t2h,&t2l,tbl[9],tbl[10],xSqh,t1h,t1l);
      MulAdd212(&t3h,&t3l,tbl[6],tbl[7],xSqh,t2h,t2l);
      MulAdd22(&t4h,&t4l,tmp4,tmp5,xSqh,xSql,t3h,t3l);
    }

    MulAdd22(&t5h,&t5l,tbl[0],tbl[1],xSqh,xSql,t4h,t4l);

    Mul122(&xCubeh,&xCubel,x,xSqh,xSql);
    Mul22(&tt6h,&tt6l,xCubeh,xCubel,t5h,t5l);
    
    Add12(tmp1,tmp2,x,tt6h);
    tmp3 = tmp2 + tt6l;
    Add12(polyh,polyl,tmp1,tmp3);

    /* Multiply by sign */
    asinh = sign * polyh;
    asinm = sign * polyl;

    /* Rounding test 
       The RU rounding constant is at tbl[35]
    */
    TEST_AND_RETURN_RU(asinh, asinm, tbl[35]);

    /* Launch accurate phase */

    asin_accurate_lower(&asinh,&asinm,&asinl,x,xSqh,xSql,sign);

    ReturnRoundUpwards3(asinh,asinm,asinl); 
  }

  if (xdb.i[HI] > BOUND9) {
    /* Special interval BOUND9..1
       We use an asymptotic development of arcsin in sqrt(1 - x)
    */

    /* Argument reduction for quick and accurate phase
       z = 1 - x
       The operation is exact as per Sterbenz' lemma
    */

    z = 1 - x;

    /* Quick phase starts */

    /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
    highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	       tbl[TBLIDX10+42] ,z,tbl[TBLIDX10+40]),z,tbl[TBLIDX10+38]),z,
               tbl[TBLIDX10+36]),z,tbl[TBLIDX10+34]),z,tbl[TBLIDX10+32]),z,
               tbl[TBLIDX10+30]),z,tbl[TBLIDX10+28]),z,tbl[TBLIDX10+26]),z,
               tbl[TBLIDX10+24]);
#else
    highPoly = tbl[TBLIDX10+24] + z * (tbl[TBLIDX10+26] + z * (tbl[TBLIDX10+28] + z * (
	       tbl[TBLIDX10+30] + z * (tbl[TBLIDX10+32] + z * (tbl[TBLIDX10+34] + z * (
	       tbl[TBLIDX10+36] + z * (tbl[TBLIDX10+38] + z * (tbl[TBLIDX10+40] + z * 
               tbl[TBLIDX10+42]))))))));
#endif
    
    /* Double-double precision evaluation */
    Mul12(&tt1h,&tt1l,z,highPoly);
    Add22(&t1h,&t1l,tbl[TBLIDX10+21],tbl[TBLIDX10+22],tt1h,tt1l);

    MulAdd212(&t2h,&t2l,tbl[TBLIDX10+18],tbl[TBLIDX10+19],z,t1h,t1l);
    MulAdd212(&t3h,&t3l,tbl[TBLIDX10+15],tbl[TBLIDX10+16],z,t2h,t2l);
    MulAdd212(&t4h,&t4l,tbl[TBLIDX10+12],tbl[TBLIDX10+13],z,t3h,t3l);
    MulAdd212(&t5h,&t5l,tbl[TBLIDX10+9],tbl[TBLIDX10+10],z,t4h,t4l);
    MulAdd212(&t6h,&t6l,tbl[TBLIDX10+6],tbl[TBLIDX10+7],z,t5h,t5l);
    MulAdd212(&t7h,&t7l,tbl[TBLIDX10+3],tbl[TBLIDX10+4],z,t6h,t6l);
    MulAdd212(&t8h,&t8l,tbl[TBLIDX10+0],tbl[TBLIDX10+1],z,t7h,t7l);
    MulAdd212(&polyh,&polyl,-1,0,z,t8h,t8l);

    /* Compute sqrt(2*z) as a double-double */

    twoZ = 2 * z;
    sqrt12(&sqrtzh,&sqrtzl,twoZ);                                                         

    /* Multiply p(z) by sqrt(2*z) and add Pi/2 */

    Mul22(&pTimesSh,&pTimesSl,polyh,polyl,sqrtzh,sqrtzl);                    
    Add22(&allh,&alll,PIHALFH,PIHALFM,pTimesSh,pTimesSl);       

    /* Multiply by sign */
    asinh = sign * allh;
    asinm = sign * alll;

    /* Rounding test 
       The RU rounding constant is at tbl[TBLIDX10+55]
    */
    TEST_AND_RETURN_RU(asinh, asinm, tbl[TBLIDX10+55]);

    /* Launch accurate phase */

    asin_accurate_higher(&asinh,&asinm,&asinl,z,sign);

    ReturnRoundUpwards3(asinh,asinm,asinl); 
  }
  
  /* General 8 main intervals 
     We can already suppose that BOUND1 <= x <= BOUND9
  */
  
  if (xdb.i[HI] < BOUND5) {
    if (xdb.i[HI] < BOUND3) {
      if (xdb.i[HI] < BOUND2) i = TBLIDX2; else i = TBLIDX3;
    } else {
      if (xdb.i[HI] < BOUND4) i = TBLIDX4; else i = TBLIDX5;
    }
  } else {
    if (xdb.i[HI] < BOUND7) {
      if (xdb.i[HI] < BOUND6) i = TBLIDX6; else i = TBLIDX7;
    } else {
      if (xdb.i[HI] < BOUND8) i = TBLIDX8; else i = TBLIDX9;
    }
  }

  /* Argument reduction 
     i points to the interval midpoint value in the table
  */
  z = x - tbl[i];

  /* Quick phase starts */

  /* Double precision evaluation */

#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	     tbl[i+35] ,z,tbl[i+33]),z,tbl[i+31]),z,tbl[i+29]),z,
             tbl[i+27]),z,tbl[i+25]),z,tbl[i+23]),z,tbl[i+21]);
#else
  highPoly = tbl[i+21] + z * (tbl[i+23] + z * (tbl[i+25] + z * (
             tbl[i+27] + z * (tbl[i+29] + z * (tbl[i+31] + z * ( 
             tbl[i+33] + z *  tbl[i+35]))))));
#endif

  /* Double-double precision evaluation */
    
  Mul12(&tt1h,&tt1l,z,highPoly);
  Add22(&t1h,&t1l,tbl[i+18],tbl[i+19],tt1h,tt1l);

  MulAdd212(&t2h,&t2l,tbl[i+15],tbl[i+16],z,t1h,t1l);
  MulAdd212(&t3h,&t3l,tbl[i+12],tbl[i+13],z,t2h,t2l);
  MulAdd212(&t4h,&t4l,tbl[i+9],tbl[i+10],z,t3h,t3l);
  MulAdd212(&t5h,&t5l,tbl[i+6],tbl[i+7],z,t4h,t4l);
  MulAdd212(&t6h,&t6l,tbl[i+3],tbl[i+4],z,t5h,t5l);
  MulAdd212(&polyh,&polyl,tbl[i+1],tbl[i+2],z,t6h,t6l);

  /* Multiply by sign */
  asinh = sign * polyh;
  asinm = sign * polyl;

  /* Rounding test 
     The RU rounding constant is at tbl[i+60]
  */
  TEST_AND_RETURN_RU(asinh, asinm, tbl[i+60]);

  /* Launch accurate phase */

  asin_accurate_middle(&asinh,&asinm,&asinl,z,i,sign);

  ReturnRoundUpwards3(asinh,asinm,asinl); 
}

double asin_rd(double x) {
  db_number xdb;
  double sign, z, asinh, asinm, asinl;
  int i;
  double xSqh, xSql;
  double tt1h, tt1l;
  double tt6h, tt6l;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l;
  double t7h, t7l, t8h, t8l, polyh, polyl, twoZ, sqrtzh, sqrtzl;
  double pTimesSh, pTimesSl, allh, alll, highPoly, xCubeh, xCubel;
  double tmp1, tmp2, tmp3, tmp4, tmp5;

  /* Transform the argument into integer */
  xdb.d = x;

  /* Special case handling */
  
  /* Strip off the sign of argument x */
  if (xdb.i[HI] & 0x80000000) sign = -1; else sign = 1;
  xdb.i[HI] &= 0x7fffffff;

  /* asin is defined on -1 <= x <= 1, elsewhere it is NaN */
  if ((xdb.i[HI] > 0x3ff00000) || ((xdb.i[HI] == 0x3ff00000) && (xdb.i[LO] != 0x00000000))) {
    return (x-x)/0.0;    /* return NaN */
  }

  /* If |x| < 2^(-28) we have
     
     arcsin(x) = x * ( 1 + xi ) 

     with 0 <= xi < 2^(-55) 
          
     So we can decide the rounding without any computation 
  */
  if (xdb.i[HI] < 0x3e300000) {
    /* If x == 0 then we got the algebraic result arcsin(0) = 0
       If x > 0 then the truncation rest is positive but less than 
       1 ulp; we round downwards by returning x
    */
    if (x >= 0) return x;
    /* Otherwise the rest is negative, less than 1 ulp and the
       image is not algebraic 
       We return x - 1ulp
       We stripped off the sign, so we add 1 ulp to -x (in xdb.d) and multiply by -1
    */
    xdb.l++;
    return -1 * xdb.d;
  }

  /* Recast x */
  x = xdb.d;

  /* Find correspondant interval and compute index to the table
     We start by filtering the two special cases around 0 and 1
  */

  if (xdb.i[HI] < BOUND1) {
    /* Special interval 0..BOUND1 
       The polynomial has no even monomials
       We must prove extra accuracy in the interval 0..sin(2^(-18))
    */

    /* Quick phase starts */

    /* Compute square of x for both quick and accurate phases */
    Mul12(&xSqh,&xSql,x,x);

    tmp4 = tbl[3];
    tmp5 = tbl[4];
    t4h = tmp4;
    t4l = tmp5;
    if (xdb.i[HI] > EXTRABOUND) {
      /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
      highPoly = FMA(FMA(FMA(FMA(tbl[23],xSqh,tbl[21]),xSqh,tbl[19]),xSqh,tbl[17]),xSqh,tbl[15]);
#else
      highPoly = tbl[15] + xSqh * (tbl[17] + xSqh * (tbl[19] + xSqh * (tbl[21] + xSqh * tbl[23])));
#endif

      /* Double-double precision evaluation */
      Mul12(&tt1h,&tt1l,xSqh,highPoly);
      Add22(&t1h,&t1l,tbl[12],tbl[13],tt1h,tt1l);
      
      MulAdd212(&t2h,&t2l,tbl[9],tbl[10],xSqh,t1h,t1l);
      MulAdd212(&t3h,&t3l,tbl[6],tbl[7],xSqh,t2h,t2l);
      MulAdd22(&t4h,&t4l,tmp4,tmp5,xSqh,xSql,t3h,t3l);
    }

    MulAdd22(&t5h,&t5l,tbl[0],tbl[1],xSqh,xSql,t4h,t4l);

    Mul122(&xCubeh,&xCubel,x,xSqh,xSql);
    Mul22(&tt6h,&tt6l,xCubeh,xCubel,t5h,t5l);
    
    Add12(tmp1,tmp2,x,tt6h);
    tmp3 = tmp2 + tt6l;
    Add12(polyh,polyl,tmp1,tmp3);

    /* Multiply by sign */
    asinh = sign * polyh;
    asinm = sign * polyl;

    /* Rounding test 
       The RD rounding constant is at tbl[35]
    */
    TEST_AND_RETURN_RD(asinh, asinm, tbl[35]);

    /* Launch accurate phase */

    asin_accurate_lower(&asinh,&asinm,&asinl,x,xSqh,xSql,sign);

    ReturnRoundDownwards3(asinh,asinm,asinl); 
  }

  if (xdb.i[HI] > BOUND9) {
    /* Special interval BOUND9..1
       We use an asymptotic development of arcsin in sqrt(1 - x)
    */

    /* Argument reduction for quick and accurate phase
       z = 1 - x
       The operation is exact as per Sterbenz' lemma
    */

    z = 1 - x;

    /* Quick phase starts */

    /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
    highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	       tbl[TBLIDX10+42] ,z,tbl[TBLIDX10+40]),z,tbl[TBLIDX10+38]),z,
               tbl[TBLIDX10+36]),z,tbl[TBLIDX10+34]),z,tbl[TBLIDX10+32]),z,
               tbl[TBLIDX10+30]),z,tbl[TBLIDX10+28]),z,tbl[TBLIDX10+26]),z,
               tbl[TBLIDX10+24]);
#else
    highPoly = tbl[TBLIDX10+24] + z * (tbl[TBLIDX10+26] + z * (tbl[TBLIDX10+28] + z * (
	       tbl[TBLIDX10+30] + z * (tbl[TBLIDX10+32] + z * (tbl[TBLIDX10+34] + z * (
	       tbl[TBLIDX10+36] + z * (tbl[TBLIDX10+38] + z * (tbl[TBLIDX10+40] + z * 
               tbl[TBLIDX10+42]))))))));
#endif
    
    /* Double-double precision evaluation */
    Mul12(&tt1h,&tt1l,z,highPoly);
    Add22(&t1h,&t1l,tbl[TBLIDX10+21],tbl[TBLIDX10+22],tt1h,tt1l);

    MulAdd212(&t2h,&t2l,tbl[TBLIDX10+18],tbl[TBLIDX10+19],z,t1h,t1l);
    MulAdd212(&t3h,&t3l,tbl[TBLIDX10+15],tbl[TBLIDX10+16],z,t2h,t2l);
    MulAdd212(&t4h,&t4l,tbl[TBLIDX10+12],tbl[TBLIDX10+13],z,t3h,t3l);
    MulAdd212(&t5h,&t5l,tbl[TBLIDX10+9],tbl[TBLIDX10+10],z,t4h,t4l);
    MulAdd212(&t6h,&t6l,tbl[TBLIDX10+6],tbl[TBLIDX10+7],z,t5h,t5l);
    MulAdd212(&t7h,&t7l,tbl[TBLIDX10+3],tbl[TBLIDX10+4],z,t6h,t6l);
    MulAdd212(&t8h,&t8l,tbl[TBLIDX10+0],tbl[TBLIDX10+1],z,t7h,t7l);
    MulAdd212(&polyh,&polyl,-1,0,z,t8h,t8l);

    /* Compute sqrt(2*z) as a double-double */

    twoZ = 2 * z;
    sqrt12(&sqrtzh,&sqrtzl,twoZ);                                                         

    /* Multiply p(z) by sqrt(2*z) and add Pi/2 */

    Mul22(&pTimesSh,&pTimesSl,polyh,polyl,sqrtzh,sqrtzl);                    
    Add22(&allh,&alll,PIHALFH,PIHALFM,pTimesSh,pTimesSl);       

    /* Multiply by sign */
    asinh = sign * allh;
    asinm = sign * alll;

    /* Rounding test 
       The RD rounding constant is at tbl[TBLIDX10+55]
    */
    TEST_AND_RETURN_RD(asinh, asinm, tbl[TBLIDX10+55]);

    /* Launch accurate phase */

    asin_accurate_higher(&asinh,&asinm,&asinl,z,sign);

    ReturnRoundDownwards3(asinh,asinm,asinl); 
  }
  
  /* General 8 main intervals 
     We can already suppose that BOUND1 <= x <= BOUND9
  */
  
  if (xdb.i[HI] < BOUND5) {
    if (xdb.i[HI] < BOUND3) {
      if (xdb.i[HI] < BOUND2) i = TBLIDX2; else i = TBLIDX3;
    } else {
      if (xdb.i[HI] < BOUND4) i = TBLIDX4; else i = TBLIDX5;
    }
  } else {
    if (xdb.i[HI] < BOUND7) {
      if (xdb.i[HI] < BOUND6) i = TBLIDX6; else i = TBLIDX7;
    } else {
      if (xdb.i[HI] < BOUND8) i = TBLIDX8; else i = TBLIDX9;
    }
  }

  /* Argument reduction 
     i points to the interval midpoint value in the table
  */
  z = x - tbl[i];

  /* Quick phase starts */

  /* Double precision evaluation */

#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	     tbl[i+35] ,z,tbl[i+33]),z,tbl[i+31]),z,tbl[i+29]),z,
             tbl[i+27]),z,tbl[i+25]),z,tbl[i+23]),z,tbl[i+21]);
#else
  highPoly = tbl[i+21] + z * (tbl[i+23] + z * (tbl[i+25] + z * (
             tbl[i+27] + z * (tbl[i+29] + z * (tbl[i+31] + z * ( 
             tbl[i+33] + z *  tbl[i+35]))))));
#endif

  /* Double-double precision evaluation */
    
  Mul12(&tt1h,&tt1l,z,highPoly);
  Add22(&t1h,&t1l,tbl[i+18],tbl[i+19],tt1h,tt1l);

  MulAdd212(&t2h,&t2l,tbl[i+15],tbl[i+16],z,t1h,t1l);
  MulAdd212(&t3h,&t3l,tbl[i+12],tbl[i+13],z,t2h,t2l);
  MulAdd212(&t4h,&t4l,tbl[i+9],tbl[i+10],z,t3h,t3l);
  MulAdd212(&t5h,&t5l,tbl[i+6],tbl[i+7],z,t4h,t4l);
  MulAdd212(&t6h,&t6l,tbl[i+3],tbl[i+4],z,t5h,t5l);
  MulAdd212(&polyh,&polyl,tbl[i+1],tbl[i+2],z,t6h,t6l);

  /* Multiply by sign */
  asinh = sign * polyh;
  asinm = sign * polyl;

  /* Rounding test 
     The RD rounding constant is at tbl[i+60]
  */
  TEST_AND_RETURN_RD(asinh, asinm, tbl[i+60]);

  /* Launch accurate phase */

  asin_accurate_middle(&asinh,&asinm,&asinl,z,i,sign);

  ReturnRoundDownwards3(asinh,asinm,asinl); 
}

double asin_rz(double x) {
  db_number xdb;
  double sign, z, asinh, asinm, asinl;
  int i;
  double xSqh, xSql;
  double tt1h, tt1l;
  double tt6h, tt6l;
  double t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l;
  double t7h, t7l, t8h, t8l, polyh, polyl, twoZ, sqrtzh, sqrtzl;
  double pTimesSh, pTimesSl, allh, alll, highPoly, xCubeh, xCubel;
  double tmp1, tmp2, tmp3, tmp4, tmp5;

  /* Transform the argument into integer */
  xdb.d = x;

  /* Special case handling */
  
  /* Strip off the sign of argument x */
  if (xdb.i[HI] & 0x80000000) sign = -1; else sign = 1;
  xdb.i[HI] &= 0x7fffffff;

  /* asin is defined on -1 <= x <= 1, elsewhere it is NaN */
  if ((xdb.i[HI] > 0x3ff00000) || ((xdb.i[HI] == 0x3ff00000) && (xdb.i[LO] != 0x00000000))) {
    return (x-x)/0.0;    /* return NaN */
  }

  /* If |x| < 2^(-28) we have
     
     arcsin(x) = x * ( 1 + xi ) 

     with 0 <= xi < 2^(-55) 
          
     So we can decide the rounding without any computation 
  */
  if (xdb.i[HI] < 0x3e300000) {
    /* If x == 0 the result is algebraic and equal to 0
       If x < 0 the truncation rest is negative and less than 1 ulp, we return x
       If x > 0 the truncation rest is positive and less than 1 ulp, we return x
    */
    return x;
  }  

  /* Recast x */
  x = xdb.d;

  /* Find correspondant interval and compute index to the table
     We start by filtering the two special cases around 0 and 1
  */

  if (xdb.i[HI] < BOUND1) {
    /* Special interval 0..BOUND1 
       The polynomial has no even monomials
       We must prove extra accuracy in the interval 0..sin(2^(-18))
    */

    /* Quick phase starts */

    /* Compute square of x for both quick and accurate phases */
    Mul12(&xSqh,&xSql,x,x);

    tmp4 = tbl[3];
    tmp5 = tbl[4];
    t4h = tmp4;
    t4l = tmp5;
    if (xdb.i[HI] > EXTRABOUND) {
      /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
      highPoly = FMA(FMA(FMA(FMA(tbl[23],xSqh,tbl[21]),xSqh,tbl[19]),xSqh,tbl[17]),xSqh,tbl[15]);
#else
      highPoly = tbl[15] + xSqh * (tbl[17] + xSqh * (tbl[19] + xSqh * (tbl[21] + xSqh * tbl[23])));
#endif

      /* Double-double precision evaluation */
      Mul12(&tt1h,&tt1l,xSqh,highPoly);
      Add22(&t1h,&t1l,tbl[12],tbl[13],tt1h,tt1l);
      
      MulAdd212(&t2h,&t2l,tbl[9],tbl[10],xSqh,t1h,t1l);
      MulAdd212(&t3h,&t3l,tbl[6],tbl[7],xSqh,t2h,t2l);
      MulAdd22(&t4h,&t4l,tmp4,tmp5,xSqh,xSql,t3h,t3l);
    }

    MulAdd22(&t5h,&t5l,tbl[0],tbl[1],xSqh,xSql,t4h,t4l);

    Mul122(&xCubeh,&xCubel,x,xSqh,xSql);
    Mul22(&tt6h,&tt6l,xCubeh,xCubel,t5h,t5l);
    
    Add12(tmp1,tmp2,x,tt6h);
    tmp3 = tmp2 + tt6l;
    Add12(polyh,polyl,tmp1,tmp3);

    /* Multiply by sign */
    asinh = sign * polyh;
    asinm = sign * polyl;

    /* Rounding test 
       The RZ rounding constant is at tbl[35]
    */
    TEST_AND_RETURN_RZ(asinh, asinm, tbl[35]);

    /* Launch accurate phase */

    asin_accurate_lower(&asinh,&asinm,&asinl,x,xSqh,xSql,sign);

    ReturnRoundTowardsZero3(asinh,asinm,asinl); 
  }

  if (xdb.i[HI] > BOUND9) {
    /* Special interval BOUND9..1
       We use an asymptotic development of arcsin in sqrt(1 - x)
    */

    /* Argument reduction for quick and accurate phase
       z = 1 - x
       The operation is exact as per Sterbenz' lemma
    */

    z = 1 - x;

    /* Quick phase starts */

    /* Double precision evaluation */
#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
    highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	       tbl[TBLIDX10+42] ,z,tbl[TBLIDX10+40]),z,tbl[TBLIDX10+38]),z,
               tbl[TBLIDX10+36]),z,tbl[TBLIDX10+34]),z,tbl[TBLIDX10+32]),z,
               tbl[TBLIDX10+30]),z,tbl[TBLIDX10+28]),z,tbl[TBLIDX10+26]),z,
               tbl[TBLIDX10+24]);
#else
    highPoly = tbl[TBLIDX10+24] + z * (tbl[TBLIDX10+26] + z * (tbl[TBLIDX10+28] + z * (
	       tbl[TBLIDX10+30] + z * (tbl[TBLIDX10+32] + z * (tbl[TBLIDX10+34] + z * (
	       tbl[TBLIDX10+36] + z * (tbl[TBLIDX10+38] + z * (tbl[TBLIDX10+40] + z * 
               tbl[TBLIDX10+42]))))))));
#endif
    
    /* Double-double precision evaluation */
    Mul12(&tt1h,&tt1l,z,highPoly);
    Add22(&t1h,&t1l,tbl[TBLIDX10+21],tbl[TBLIDX10+22],tt1h,tt1l);

    MulAdd212(&t2h,&t2l,tbl[TBLIDX10+18],tbl[TBLIDX10+19],z,t1h,t1l);
    MulAdd212(&t3h,&t3l,tbl[TBLIDX10+15],tbl[TBLIDX10+16],z,t2h,t2l);
    MulAdd212(&t4h,&t4l,tbl[TBLIDX10+12],tbl[TBLIDX10+13],z,t3h,t3l);
    MulAdd212(&t5h,&t5l,tbl[TBLIDX10+9],tbl[TBLIDX10+10],z,t4h,t4l);
    MulAdd212(&t6h,&t6l,tbl[TBLIDX10+6],tbl[TBLIDX10+7],z,t5h,t5l);
    MulAdd212(&t7h,&t7l,tbl[TBLIDX10+3],tbl[TBLIDX10+4],z,t6h,t6l);
    MulAdd212(&t8h,&t8l,tbl[TBLIDX10+0],tbl[TBLIDX10+1],z,t7h,t7l);
    MulAdd212(&polyh,&polyl,-1,0,z,t8h,t8l);

    /* Compute sqrt(2*z) as a double-double */

    twoZ = 2 * z;
    sqrt12(&sqrtzh,&sqrtzl,twoZ);                                                         

    /* Multiply p(z) by sqrt(2*z) and add Pi/2 */

    Mul22(&pTimesSh,&pTimesSl,polyh,polyl,sqrtzh,sqrtzl);                    
    Add22(&allh,&alll,PIHALFH,PIHALFM,pTimesSh,pTimesSl);       

    /* Multiply by sign */
    asinh = sign * allh;
    asinm = sign * alll;

    /* Rounding test 
       The RZ rounding constant is at tbl[TBLIDX10+55]
    */
    TEST_AND_RETURN_RZ(asinh, asinm, tbl[TBLIDX10+55]);

    /* Launch accurate phase */

    asin_accurate_higher(&asinh,&asinm,&asinl,z,sign);

    ReturnRoundTowardsZero3(asinh,asinm,asinl); 
  }
  
  /* General 8 main intervals 
     We can already suppose that BOUND1 <= x <= BOUND9
  */
  
  if (xdb.i[HI] < BOUND5) {
    if (xdb.i[HI] < BOUND3) {
      if (xdb.i[HI] < BOUND2) i = TBLIDX2; else i = TBLIDX3;
    } else {
      if (xdb.i[HI] < BOUND4) i = TBLIDX4; else i = TBLIDX5;
    }
  } else {
    if (xdb.i[HI] < BOUND7) {
      if (xdb.i[HI] < BOUND6) i = TBLIDX6; else i = TBLIDX7;
    } else {
      if (xdb.i[HI] < BOUND8) i = TBLIDX8; else i = TBLIDX9;
    }
  }

  /* Argument reduction 
     i points to the interval midpoint value in the table
  */
  z = x - tbl[i];

  /* Quick phase starts */

  /* Double precision evaluation */

#if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
  highPoly = FMA(FMA(FMA(FMA(FMA(FMA(FMA(
	     tbl[i+35] ,z,tbl[i+33]),z,tbl[i+31]),z,tbl[i+29]),z,
             tbl[i+27]),z,tbl[i+25]),z,tbl[i+23]),z,tbl[i+21]);
#else
  highPoly = tbl[i+21] + z * (tbl[i+23] + z * (tbl[i+25] + z * (
             tbl[i+27] + z * (tbl[i+29] + z * (tbl[i+31] + z * ( 
             tbl[i+33] + z *  tbl[i+35]))))));
#endif

  /* Double-double precision evaluation */
    
  Mul12(&tt1h,&tt1l,z,highPoly);
  Add22(&t1h,&t1l,tbl[i+18],tbl[i+19],tt1h,tt1l);

  MulAdd212(&t2h,&t2l,tbl[i+15],tbl[i+16],z,t1h,t1l);
  MulAdd212(&t3h,&t3l,tbl[i+12],tbl[i+13],z,t2h,t2l);
  MulAdd212(&t4h,&t4l,tbl[i+9],tbl[i+10],z,t3h,t3l);
  MulAdd212(&t5h,&t5l,tbl[i+6],tbl[i+7],z,t4h,t4l);
  MulAdd212(&t6h,&t6l,tbl[i+3],tbl[i+4],z,t5h,t5l);
  MulAdd212(&polyh,&polyl,tbl[i+1],tbl[i+2],z,t6h,t6l);

  /* Multiply by sign */
  asinh = sign * polyh;
  asinm = sign * polyl;

  /* Rounding test 
     The RZ rounding constant is at tbl[i+60]
  */
  TEST_AND_RETURN_RZ(asinh, asinm, tbl[i+60]);

  /* Launch accurate phase */

  asin_accurate_middle(&asinh,&asinm,&asinl,z,i,sign);

  ReturnRoundTowardsZero3(asinh,asinm,asinl); 
}