开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 通讯/手机编程 > J2ME > 查看源码
MicroDouble.java
资源名称：GMapsTest.rar　[点击查看]
上传用户：hygd004
上传日期：2022-07-01
资源大小：246k
文件大小：135k
源码类别：

J2ME

开发平台：

Java

MicroDouble.java：源码内容
  1.           lo = negate(LN2_LO);
  2.           k = -1;
  3.         }
  4.       } else {
  5.         long tmp = mul(INV_LN2, d);
  6.         if (xsb == 0) {
  7.           tmp = add(tmp, ONE_HALF);
  8.         } else {
  9.           tmp = sub(tmp, ONE_HALF);
  10.         }
  11.         k = intValue(add(mul(INV_LN2, d), 
  12.                 (xsb == 0) ? ONE_HALF : negate(ONE_HALF)));
  13.         long t = intToDouble(k);
  14.         hi = sub(d, mul(t, LN2_HI)); //  t*ln2_hi is exact here 
  15.         lo = mul(t, LN2_LO);
  16.       }
  17.       d = sub(hi, lo);
  18.       c = sub(sub(hi, d), lo);
  19.     } 
  20.     else if(hx < 0x3c900000) { // when |x|<2**-54, return x 
  21.       return d;
  22.     }
  23.     else {
  24.       k = 0;
  25.       c = 0;
  26.     }
  28.     // x is now in primary range 
  29.     long hfx = scalbn(d, -1);
  30.     long hxs = mul(d, hfx);
  31.     long r1 = add(ONE, mul(hxs, add(Q1, mul(hxs, add(Q2, mul(hxs, 
  32.             add(Q3, mul(hxs, add(Q4, mul(hxs, Q5))))))))));
  33.     long t = sub(THREE, mul(r1, hfx));
  34.     long e = mul(hxs, div(sub(r1, t), sub(SIX, mul(d, t))));
  35.     if(k==0) return sub(d, sub(mul(d, e), hxs)); // c is 0
  36.     else {
  37.       e = sub(sub(mul(d, sub(e, c)), c), hxs);
  38.       if (k == -1) return sub(scalbn(sub(d, e), -1), ONE_HALF);
  39.       if(k==1) 
  40.         if (lt(d, negate(ONE_FOURTH))) return negate(scalbn(sub(e, add(d, 
  41.                 ONE_HALF)), 1));
  42.         else return add(ONE, scalbn(sub(d, e), 1));
  43.         if (k <= -2 || k>56) { // suffice to return exp(x)-1 
  44.           y = sub(ONE, sub(e, d));
  45.           y = setHI(y, getHI(y) + (k << 20)); // add k to y's exponent 
  46.           return sub(y, ONE);
  47.         }
  48.         t = ONE;
  49.         if(k<20) {
  50.           t = setHI(t, 0x3ff00000 - (0x200000>>k)); // t=1-2^-k 
  51.           y = sub(t, sub(e, d));
  52.           y = setHI(y, getHI(y) + (k << 20)); // add k to y's exponent 
  53.        } else {
  54.          t = setHI(t, ((0x3ff-k)<<20)); // 2^-k 
  55.          y = add(sub(d, add(e, t)), ONE);
  56.          y = setHI(y, getHI(y) + (k << 20)); //add k to y's exponent 
  57.        }
  58.     }
  59.     return y;
  60.   }
  62.   private static final long BP[]       = { ONE, THREE_HALVES };
  63.   private static final long DP_HI[]    = { ZERO, 0x3fe2b80340000000L}; // 5.84962487220764160156e-01 
  64.   private static final long DP_LO[]    = { ZERO, 0x3e4cfdeb43cfd006L}; // 1.35003920212974897128e-08 
  65.   // poly coefs for (3/2)*(log(x)-2s-2/3*s**3 
  66.   private static final long L1         = 0x3fe3333333333303L; // 5.99999999999994648725e-01 
  67.   private static final long L2         = 0x3fdb6db6db6fabffL; //  4.28571428578550184252e-01 
  68.   private static final long L3         = 0x3fd55555518f264dL; //  3.33333329818377432918e-01 
  69.   private static final long L4         = 0x3fd17460a91d4101L; //  2.72728123808534006489e-01 
  70.   private static final long L5         = 0x3fcd864a93c9db65L; //  2.30660745775561754067e-01 
  71.   private static final long L6         = 0x3fca7e284a454eefL; //  2.06975017800338417784e-01 
  72.   private static final long LN2_HI_B   = 0x3fe62e4300000000L; //  6.93147182464599609375e-01 
  73.   private static final long LN2_LO_B   = 0xbe205c610ca86c39L; // -1.90465429995776804525e-09) 0xbe205c610ca86c39 
  74.   private static final long OVT        = 0x3c971547652b82feL; // 8.0085662595372944372e-17) // -(1024-log2(ovfl+.5ulp)) 
  75.   private static final long CP         = 0x3feec709dc3a03fdL; //  9.61796693925975554329e-01 = 2/(3ln2)
  76.   private static final long CP_HI      = 0x3feec709e0000000L; //  9.61796700954437255859e-01 = (float)cp
  77.   private static final long CP_LO      = 0xbe3e2fe0145b01f5L; // -7.02846165095275826516e-09 = tail of cp_h
  78.   private static final long INV_LN2_HI = 0x3ff7154760000000L; //  1.44269502162933349609e+00 = 24b 1/ln2
  79.   private static final long INV_LN2_LO = 0x3e54ae0bf85ddf44L; //  1.92596299112661746887e-08 = 1/ln2 tail
  81.   /**
  82.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#pow(double, double)">Math.pow(double, double)</a>.
  83.    */
  84.   public static long pow(long d1, long d2) {
  85.     if (isZero(d2)) {
  86.       return ONE;
  87.     } else if (isNaN(d1) || isNaN(d2)) {
  88.       return NaN;
  89.     }
  90.       
  91.     int i0 = (int) ((ONE  >> 29) ^ 1);
  92.     int i1 = 1 - i0;    
  93.     int hx = getHI(d1);
  94.     int lx = getLO(d1);
  95.     int hy = getHI(d2);
  96.     int ly = getLO(d2);
  97.     int ix = hx & 0x7fffffff;
  98.     int iy = hy & 0x7fffffff;
  100.     // determine if y is an odd int when x < 0
  101.     // yisint = 0 ... y is not an integer
  102.     // yisint = 1 ... y is an odd int
  103.     // yisint = 2 ... y is an even int
  104.     int yisint  = 0;
  105.     if (hx < 0) {
  106.       if (iy >= 0x43400000)
  107.         yisint = 2; // even integer y
  108.       else if(iy >= 0x3ff00000) {
  109.         int k = (iy>>20)-0x3ff; // exponent
  110.         if (k > 20) {
  111.           int j = ly>>>(52-k);
  112.           if ((j<<(52-k)) == ly)
  113.             yisint = 2-(j&1);
  114.         } else if (ly == 0) {
  115.           int j = iy >> (20-k);
  116.           if ((j<<(20-k)) == iy)
  117.             yisint = 2-(j&1);
  118.         }
  119.       }
  120.     }
  122.     // special value of y
  123.     if (ly==0) {
  124.       if (iy == 0x7ff00000) { // y is +-inf
  125.         if (((ix-0x3ff00000)|lx) == 0)
  126.                 return  NaN; // inf**+-1 is NaN 
  127.         else if (ix >= 0x3ff00000) // (|x|>1)**+-inf = inf,0 
  128.                 return (hy>=0)? POSITIVE_INFINITY : ZERO;
  129.         else // (|x|<1)**-,+inf = inf,0 
  130.                 return (hy<0) ? POSITIVE_INFINITY : ZERO;
  131.       } 
  132.       if (iy == 0x3ff00000) { // y is  +-1 
  133.         if (hy < 0)
  134.           return div(ONE, d1);
  135.         else
  136.           return d1;
  137.       }
  138.       if (hy == 0x40000000)
  139.         return mul(d1, d1); // y is  2
  140.       if (hy == 0x3fe00000) { // y is  0.5
  141.         if (hx>=0) // x >= +0 
  142.           return sqrt(d1);
  143.       }
  144.     }
  146.     long ax = abs(d1);
  147.     // special value of x
  148.     if (lx == 0) {
  149.       if (ix==0x7ff00000 || ix==0 || ix==0x3ff00000) {
  150.         long z = ax; // x is +-0,+-inf,+-1
  151.         if (hy < 0 )
  152.           z = div(ONE, z); // z = (1/|x|)
  153.         if (hx < 0) {
  154.           if (((ix-0x3ff00000)|yisint) == 0) {
  155.             z = NaN; // (-1)**non-int is NaN
  156.           } else if (yisint==1) 
  157.             z = negate(z); // (x<0)**odd = -(|x|**odd)
  158.         }
  159.         return z;
  160.       }
  161.     }
  163.   int n = (hx >> 31) + 1;
  165.     /* (x<0)**(non-int) is NaN */
  166.   if ((n | yisint) == 0) {
  167.       return NaN;
  168.     }
  170.     boolean negative = ((n | (yisint-1)) == 0);
  172.     // |y| is huge
  173.     long t1, t2;
  174.     if (iy > 0x41e00000) { // if |y| > 2**31
  175.       if (iy > 0x43f00000){ // if |y| > 2**64, must o/uflow 
  176.         if (ix <= 0x3fefffff) {
  177.           return ((hy < 0) ? POSITIVE_INFINITY : ZERO);
  178.         }
  179.         return ((hy > 0) ? POSITIVE_INFINITY : ZERO);
  180.       }
  181.       // over/underflow if x is not close to one
  182.       if (ix < 0x3fefffff) {
  183.         if (negative) {
  184.           return ((hy < 0) ? NEGATIVE_INFINITY : NEGATIVE_ZERO);
  185.         } else {
  186.           return ((hy < 0) ? POSITIVE_INFINITY : ZERO);
  187.         }
  188.       }
  189.       if (ix > 0x3ff00000) {
  190.         if (negative) {
  191.           return ((hy > 0) ? NEGATIVE_INFINITY : NEGATIVE_ZERO);
  192.         } else {
  193.           return ((hy > 0) ? POSITIVE_INFINITY : ZERO);
  194.         }
  195.       }
  196.       // now |1-x| is tiny <= 2**-20, suffice to compute
  197.       // log(x) by x-x^2/2+x^3/3-x^4/4
  198.       long t = sub(ax, ONE); // t has 20 trailing zeros
  199.       long w = mul(mul(t, t), sub(ONE_HALF, mul(t, 
  200.               sub(ONE_THIRD, mul(t, ONE_FOURTH)))));
  201.       long u = mul(INV_LN2_HI, t); // INV_LN2_HI has 21 sig. bits
  202.       long v = sub(mul(t, INV_LN2_LO), mul(w, INV_LN2));
  203.       t1 = setLO(add(u, v), 0);
  204.       t2 = sub(v, sub(t1, u));
  205.     } else {
  206.       n = 0;
  207.       // take care subnormal number
  208.       if (ix<0x00100000) {
  209.         ax = scalbn(ax, 53);
  210.         n -= 53;
  211.         ix = getHI(ax);
  212.       }
  213.       n  += ((ix)>>20) - 0x3ff;
  214.       int j  = ix & 0x000fffff;
  215.       // determine interval 
  216.       ix = j|0x3ff00000; // normalize ix
  217.       int k;
  218.       if (j <= 0x3988E)
  219.         k=0; // |x|<sqrt(3/2)
  220.       else if (j < 0xBB67A)
  221.         k=1; // |x|<sqrt(3)
  222.       else {
  223.         k = 0;
  224.         n+=1;
  225.         ix -= 0x00100000;
  226.       }
  227.       ax = setHI(ax, ix);
  229.       // compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5)
  230.       long u = sub(ax, BP[k]); // bp[0]=1.0, bp[1]=1.5
  231.       long v = div(ONE, add(ax, BP[k]));
  232.       long ss = mul(u, v);
  233.       long s_h = setLO(ss, 0);
  234.       // t_h=ax+bp[k] High 
  235.       long t_h = setHI(ZERO, ((ix>>1)|0x20000000)+0x00080000+(k<<18)); 
  236.       long t_l = sub(ax, sub(t_h, BP[k]));
  237.       long s_l = mul(v, sub(sub(u, mul(s_h, t_h)), mul(
  238.               s_h, t_l)));
  239.       // compute log(ax)
  240.       long s2 = mul(ss, ss);
  241.       long r = mul(mul(s2, s2), add(L1, mul(s2, add(L2, mul(
  242.               s2, add(L3, mul(s2, add(L4, mul(s2, add(L5, mul(
  243.               s2, L6)))))))))));
  244.       r = add(r, mul(s_l, add(s_h, ss)));
  245.       s2 = mul(s_h, s_h);
  246.       t_h = setLO(add(add(THREE, s2), r), 0);
  247.       t_l = sub(r, sub(sub(t_h, THREE), s2));
  248.       // u+v = ss*(1+...) 
  249.       u = mul(s_h, t_h);
  250.       v = add(mul(s_l, t_h), mul(t_l, ss));
  251.       // 2/(3log2)*(ss+...) 
  252.       long p_h = setLO(add(u, v), 0);
  253.       long p_l = sub(v, sub(p_h, u));
  254.       long z_h = mul(CP_HI, p_h); // cp_h+cp_l = 2/(3*log2) 
  255.       long z_l = add(add(mul(CP_LO, p_h), mul(p_l, CP)), DP_LO[k]);
  256.       // log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l 
  257.       long t = intToDouble(n);
  258.       t1 = setLO(add(add(add(z_h, z_l), DP_HI[k]), t), 0);
  259.       t2 = sub(z_l, sub(sub(sub(t1, t), DP_HI[k]), z_h));
  260.     }
  262.     // split up y into y1+y2 and compute (y1+y2)*(t1+t2) 
  263.     long y1 = setLO(d2, 0);
  264.     long p_l = add(mul(sub(d2, y1), t1), mul(d2, t2));
  265.     long p_h = mul(y1, t1);
  266.     long z = add(p_l, p_h);
  267.     int j = getHI(z);
  268.     int i = getLO(z);
  269.     if (j >= 0x40900000) { // z >= 1024
  270.       if ((((j-0x40900000)|i) != 0) // if z > 1024
  271.           || gt(add(p_l, OVT), sub(z, p_h))) {
  272.         return negative ? NEGATIVE_INFINITY : POSITIVE_INFINITY;
  273.       }
  274.     } else if((j&0x7fffffff) >= 0x4090cc00) { // z <= -1075
  275.       if ((((j-0xc090cc00)|i) != 0) // z < -1075
  276.           || le(p_l, sub(z, p_h))) {
  277.         return negative ? NEGATIVE_ZERO : ZERO;
  278.       }
  279.     }
  280.     // compute 2**(p_h+p_l)
  281.     i = j&0x7fffffff;
  282.     int k = (i>>20)-0x3ff;
  283.     n = 0;
  284.     if (i > 0x3fe00000) { // if |z| > 0.5, set n = [z+0.5]
  285.       n = j + (0x00100000>>(k+1));
  286.       k = ((n&0x7fffffff)>>20) - 0x3ff; // new k for n
  287.       long t = ZERO;
  288.       t = setHI(t, (n&~(0x000fffff>>k)));
  289.       n = ((n&0x000fffff)|0x00100000)>>(20-k);
  290.       if (j < 0) n = -n;
  291.       p_h = sub(p_h, t);
  292.     }
  293.     long t = setLO(add(p_l, p_h), 0);
  294.     long u = mul(t, LN2_HI_B);
  295.     long v = add(mul(sub(p_l, sub(t, p_h)), LN2), mul(t, 
  296.             LN2_LO_B));
  297.     z = add(u, v);
  298.     long w = sub(v, sub(z, u));
  299.     t = mul(z, z);
  300.     t1 = sub(z, mul(t, add(P1, mul(t, add(P2, mul(t, add(P3,
  301.             mul(t, add(P4, mul(t, P5))))))))));
  302.     long r = sub(div(mul(z, t1), sub(t1, TWO)), add(
  303.             w, mul(z, w)));
  304.     z = sub(ONE, sub(r, z));
  305.     j = getHI(z);
  306.     j += (n<<20);
  307.     if ((j>>20) <= 0) {
  308.       z = scalbn(z,n); //subnormal output
  309.     } else {
  310.       i = getHI(z);
  311.       i += (n << 20);
  312.       z = setHI(z, i);
  313.     }
  314.     if (negative) {
  315.       z = negate(z);
  316.     }
  317.     return z;
  318.   }
  319.   
  320.   /**
  321.    * Returns the logarithm of a <code>double</code> value using a specified
  322.    * base.  For most arguments, the return value is computed as:
  323.    * <code>log(d) / log(base)</code>.  If <code>base</code> is <code>E</code> or 
  324.    * <code>10</code> the dedicated log function is used.  If <code>base</code>
  325.    * is zero, infinite, <code>NaN</code>, or negative, <code>NaN</code> is
  326.    * returned.
  327.    *
  328.    * @param   d   a <code>double</code> value greater than <code>0.0</code>.
  329.    * @param   base   a <code>double</code> value greater than <code>0.0</code>.
  330.    * @return  the value log<sub><code>base</code></sub>&nbsp;<code>d</code>
  331.    */
  332.   public static long log(long d, long base) {
  333.     if (base == E) {
  334.       return log(d);
  335.     } else if (base == TEN) {
  336.       return log10(d);
  337.     } else if (isZero(base) || isInfinite(base) || isNaN(base) 
  338.             || unpackSign(base)) {
  339.       return NaN;
  340.     }
  341.     return div(log(d), log(base));
  342.   }
  344.   private static final long IVLN10    = 0x3FDBCB7B1526E50EL; // 4.34294481903251816668e-01
  345.   private static final long LOG10_2HI = 0x3FD34413509F6000L; // 3.01029995663611771306e-01
  346.   private static final long LOG10_2LO = 0x3D59FEF311F12B36L; // 3.69423907715893078616e-13
  348.   /**
  349.    * Returns the base 10 logarithm of a <code>double</code>
  350.    * value.  
  351.    *
  352.    * @param   d   a <code>double</code> value greater than <code>0.0</code>.
  353.    * @return  the value log<sub>10</sub>&nbsp;<code>d</code>
  354.    */
  355.   public static long log10(long d) {
  356.     if (isZero(d)) {
  357.       return NEGATIVE_INFINITY;
  358.     } else if (isNaN(d) || unpackSign(d)) {
  359.       return NaN;
  360.     } else if (d == POSITIVE_INFINITY) {
  361.       return d;
  362.     }
  363.     int n = ilogb(d);
  364.     if (n < 0) {
  365.       n++;
  366.     }
  367.     d = scalbn(d, -n);
  368.     long dn = intToDouble(n);
  369.     return add(mul(dn, LOG10_2HI), add(mul(dn, LOG10_2LO), 
  370.             mul(IVLN10, log(d))));
  371.   }
  372.   
  373.   private static final long LG1 = 0x3fe5555555555593L;  // 6.666666666666735130e-01
  374.   private static final long LG2 = 0x3fd999999997fa04L;  // 3.999999999940941908e-01
  375.   private static final long LG3 = 0x3fd2492494229359L;  // 2.857142874366239149e-01
  376.   private static final long LG4 = 0x3fcc71c51d8e78afL;  // 2.222219843214978396e-01
  377.   private static final long LG5 = 0x3fc7466496cb03deL;  // 1.818357216161805012e-01
  378.   private static final long LG6 = 0x3fc39a09d078c69fL;  // 1.531383769920937332e-01
  379.   private static final long LG7 = 0x3fc2f112df3e5244L;  // 1.479819860511658591e-01
  381.   /**
  382.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#log(double)">Math.log(double)</a>.
  383.    */
  384.   public static long log(long d) {
  385.     if (isZero(d)) {
  386.       return NEGATIVE_INFINITY;
  387.     } else if (isNaN(d) || unpackSign(d)) {
  388.       return NaN;
  389.     } else if (d == POSITIVE_INFINITY) {
  390.       return d;
  391.     }
  393.     int hx = getHI(d); // high word of x 
  394.     int k=0;
  395.     if (hx < 0x00100000) { // x < 2**-1022
  396.       k -= 54;
  397.       d = scalbn(d, 54);
  398.       hx = getHI(d);
  399.     } 
  400.     k += (hx>>20)-1023;
  401.     hx &= 0x000fffff;
  402.     int i = (hx+0x95f64)&0x100000;
  403.     d = setHI(d, hx|(i^0x3ff00000)); // normalize x or x/2
  404.     k += (i>>20);
  405.     long f = sub(d, ONE);
  406.     if ((0x000fffff&(2+hx))<3) { // |f| < 2**-20
  407.       if (isZero(f)) {
  408.         if (k == 0) {
  409.           return ZERO;
  410.         }
  411.         long dk = intToDouble(k);
  412.         return add(mul(dk, LN2_HI), mul(dk, LN2_LO));
  413.       }
  414.       long R = mul(mul(f, f), sub(ONE_HALF, 
  415.               mul(ONE_THIRD, f)));
  416.       if (k == 0) {
  417.         return sub(f, R);
  418.       } else {
  419.         long dk = intToDouble(k);
  420.         return sub(mul(dk, LN2_HI), sub(sub(R, 
  421.                 mul(dk, LN2_LO)), f));
  422.       }
  423.     }
  424.     long dk = intToDouble(k);
  425.     long s = div(f, add(TWO, f));
  426.     long z = mul(s, s);
  427.     long w = mul(z, z);
  428.     long R = add(mul(w, add(LG2, mul(w, add(LG4, mul(w, LG6))))),
  429.             mul(z, add(LG1, mul(w, add(LG3, mul(w, add(LG5, 
  430.             mul(w, LG7))))))));
  431.     i = (hx - 0x6147a) | (0x6b851 - hx);
  432.     if (i > 0) {
  433.       long hfsq = mul(scalbn(f, -1), f);
  434.       if (k == 0) {
  435.         return sub(f, sub(hfsq, mul(s, add(hfsq, R))));
  436.       } else {
  437.         return sub(mul(dk, LN2_HI), sub(sub(hfsq, 
  438.                add(mul(s, add(hfsq, R)), mul(dk, LN2_LO))), f));
  439.       }
  440.     } else if (k==0) {
  441.       return sub(f, mul(s, sub(f, R)));
  442.     } 
  443.     return sub(mul(dk, LN2_HI), sub(sub(mul(s, 
  444.             sub(f, R)), mul(dk, LN2_LO)), f));
  445.   }
  447.   private static final long LP1 = 0x3FE5555555555593L; // 6.666666666666735130e-01
  448.   private static final long LP2 = 0x3FD999999997FA04L; // 3.999999999940941908e-01
  449.   private static final long LP3 = 0x3FD2492494229359L; // 2.857142874366239149e-01
  450.   private static final long LP4 = 0x3FCC71C51D8E78AFL; // 2.222219843214978396e-01
  451.   private static final long LP5 = 0x3FC7466496CB03DEL; // 1.818357216161805012e-01
  452.   private static final long LP6 = 0x3FC39A09D078C69FL; // 1.531383769920937332e-01
  453.   private static final long LP7 = 0x3FC2F112DF3E5244L; // 1.479819860511658591e-01
  455.   /**
  456.    * Returns the natural logarithm of <code>1 + d</code>, computed in a way 
  457.    * that is accurate even when the value of d is close to zero.
  458.    *
  459.    * @param   d   a <code>double</code> value greater than <code>-1.0</code>.
  460.    * @return  the value ln&nbsp;<code>d</code>, the natural logarithm of
  461.    *          <code>d + 1</code>.
  462.    * @see #log(long)
  463.    * @see #expm1(long)
  464.    */
  465.   public static long log1p(long d) {
  466.     int hx = getHI(d);
  467.     int ax = hx & 0x7fffffff;
  468.     int k = 1;
  469.     int hu = 0;
  470.     long f = 0;
  471.     if (hx < 0x3FDA827A) { // x < 0.41422 
  472.       if(ax>=0x3ff00000) { // x <= -1.0 
  473.         if (d == NEGATIVE_ONE) return NEGATIVE_INFINITY; // log1p(-1)=+inf 
  474.         else return NaN; // log1p(x<-1)=NaN
  475.       }
  476.       if(ax<0x3e200000) { // |x| < 2**-29 
  477.         if(ax<0x3c900000) // |x| < 2**-54 
  478.           return d;
  479.         else
  480.           return sub(d, scalbn(mul(d, d), -1));
  481.       }
  482.       if(hx>0||hx<=0xbfd2bec3) {
  483.         k=0;f=d;hu=1;} // -0.2929<x<0.41422 
  484.     } 
  485.     if (hx >= 0x7ff00000) return d;
  486.     long u;
  487.     long c = 0;
  488.     if(k!=0) {
  489.       if(hx<0x43400000) {
  490.         u = add(ONE, d);
  491.         hu = getHI(u); // high word of u 
  492.         k  = (hu>>20)-1023;
  493.         c  = (k > 0) ? sub(ONE, sub(u, d)) : sub(d, sub(u, 
  494.                 ONE)); // correction term 
  495.         c  = div(c, u);
  496.       } else {
  497.         u  = d;
  498.         hu = getHI(u); // high word of u 
  499.         k  = (hu>>20)-1023;
  500.         c  = ZERO;
  501.       }
  502.       hu &= 0x000fffff;
  503.       if(hu<0x6a09e) {
  504.         u = setHI(u, hu|0x3ff00000); // normalize u 
  505.       } else {
  506.         k += 1;
  507.         u = setHI(u, hu|0x3fe00000); // normalize u/2 
  508.         hu = (0x00100000-hu)>>2;
  509.       }
  510.       f = add(u, NEGATIVE_ONE);
  511.     }
  512.     long hfsq= scalbn(mul(f, f), -1);
  513.     long R;
  514.     long dk = intToDouble(k);
  515.     if(hu==0) { // |f| < 2**-20 
  516.       if(isZero(f)) {
  517.         if(k==0) {
  518.           return ZERO;  
  519.         } else {
  520.           c = add(c, mul(dk, LN2_LO)); 
  521.           return add(mul(dk, LN2_HI), c);
  522.         }
  523.       }
  524.       R = mul(hfsq, sub(ONE, div(scalbn(f, 2), THREE)));
  525.       if(k==0) {
  526.         return sub(f, R); 
  527.       } else {
  528.         return sub(mul(dk, LN2_HI), sub(sub(
  529.                 R, add(mul(dk, LN2_LO), c)), f));
  530.       }
  531.     }
  532.     long s = div(f, add(TWO, f));
  533.     long z = mul(s, s);
  534.     R = mul(z, add(LP1, mul(z, add(LP2, mul(z, add(LP3, mul(
  535.             z, add(LP4, mul(z, add(LP5, mul(z, add(LP6, mul(z, 
  536.             LP7)))))))))))));
  537.     if(k==0) {
  538.       return sub(f, sub(hfsq, mul(s, add(hfsq, R))));
  539.     } else {
  540.       return sub(mul(dk, LN2_HI), sub(sub(hfsq, add(
  541.               mul(s, add(hfsq, R)), add(mul(dk, LN2_LO), c))), f));
  542.     }
  543.   }
  544.   
  545.   /**
  546.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#sin(double)">Math.sin(double)</a>.
  547.    */
  548.   public static long sin(long d) {
  549.     int ix = getHI(d) & 0x7fffffff;
  550.     if (ix <= 0x3fe921fb) {
  551.       // |x| ~< pi/4 
  552.       return kernelSin(d, ZERO, 0);
  553.     } else if (ix >= 0x7ff00000) {
  554.       // sin(Inf or NaN) is NaN 
  555.       return NaN;
  556.     } else {
  557.       // argument reduction needed 
  558.       long[] y = new long[2];
  559.       int n = remPio2(d, y);
  560.       switch(n&3) {
  561.         case 0:
  562.           return  kernelSin(y[0], y[1], 1);
  563.         case 1:
  564.           return  kernelCos(y[0], y[1]);
  565.         case 2:
  566.           return negate(kernelSin(y[0], y[1], 1));
  567.         default:
  568.           return negate(kernelCos(y[0], y[1]));
  569.       }
  570.     }
  571.   }
  573.   /**
  574.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#cos(double)">Math.cos(double)</a>.
  575.    */
  576.   public static long cos(long d) {
  577.     // High word of x.
  578.     int ix = getHI(d) & 0x7fffffff;
  580.     // |x| ~< pi/4
  581.     if (ix <= 0x3fe921fb) {
  582.       return kernelCos(d, ZERO);
  583.     } else if (ix >= 0x7ff00000) {
  584.       // cos(Inf or NaN) is NaN
  585.       return NaN;
  586.     } else {
  587.       // argument reduction needed
  588.       long y[] = new long[2];
  589.       int n = remPio2(d,y);
  590.       switch(n&3) {
  591.         case 0: 
  592.           return kernelCos(y[0],y[1]);
  593.         case 1: 
  594.           return negate(kernelSin(y[0],y[1],1));
  595.         case 2: 
  596.           return negate(kernelCos(y[0],y[1]));
  597.         default:
  598.           return kernelSin(y[0],y[1],1);
  599.       }
  600.     }
  601.   }
  603.   /**
  604.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#tan(double)">Math.tan(double)</a>.
  605.    */
  606.   public static long tan(long d) {
  607.     // |x| ~< pi/4 
  608.     int ix = getHI(d) & 0x7fffffff;
  609.     if (ix <= 0x3fe921fb) {
  610.       return kernelTan(d, ZERO, 1);
  611.     } else if (ix >= 0x7ff00000) {
  612.       // tan(Inf or NaN) is NaN
  613.       return NaN;
  614.     } else {
  615.       // argument reduction needed 
  616.       long[] y = new long [2];
  617.       int n = remPio2(d, y);
  618.       // 1 -- n even -1 -- n odd
  619.       return kernelTan(y[0], y[1], 1-((n&1)<<1));
  620.     }
  621.   }
  622.         
  623.   private static final long PIO4    = 0x3fe921fb54442d18L; // 7.85398163397448278999e-01 
  624.   private static final long PIO4_LO = 0x3c81a62633145c07L; // 3.06161699786838301793e-17 
  625.   private static final long T0      = 0x3fd5555555555563L; // 3.33333333333334091986e-01 
  626.   private static final long T1      = 0x3fc111111110fe7aL; // 1.33333333333201242699e-01
  627.   private static final long T2      = 0x3faba1ba1bb341feL; // 5.39682539762260521377e-02 
  628.   private static final long T3      = 0x3f9664f48406d637L; // 2.18694882948595424599e-02 
  629.   private static final long T4      = 0x3f8226e3e96e8493L; // 8.86323982359930005737e-03 
  630.   private static final long T5      = 0x3f6d6d22c9560328L; // 3.59207910759131235356e-03 
  631.   private static final long T6      = 0x3f57dbc8fee08315L; // 1.45620945432529025516e-03 
  632.   private static final long T7      = 0x3f4344d8f2f26501L; // 5.88041240820264096874e-04 
  633.   private static final long T8      = 0x3f3026f71a8d1068L; // 2.46463134818469906812e-04 
  634.   private static final long T9      = 0x3f147e88a03792a6L; // 7.81794442939557092300e-05 
  635.   private static final long T10     = 0x3f12b80f32f0a7e9L; // 7.14072491382608190305e-05 
  636.   private static final long T11     = 0xbef375cbdb605373L; // -1.85586374855275456654e-05 
  637.   private static final long T12     = 0x3efb2a7074bf7ad4L; // 2.59073051863633712884e-05 
  638.         
  639.   private static long kernelTan(long x, long y, int iy) {
  640.     int hx = getHI(x); // high word of x
  641.     int ix = hx & 0x7fffffff;
  642.     if (ix < 0x3e300000) { // x < 2**-28 
  643.       if (intValue(x) == 0) {
  644.         if (((ix | getLO(x)) | (iy + 1)) == 0) {
  645.           return POSITIVE_INFINITY;
  646.         } else if (iy == 1) {
  647.           return x;
  648.         } else {
  649.           /* compute -1 / (x+y) carefully */
  650.           long w = add(x, y);
  651.           long z = setLO(w, 0);
  652.           long v = sub(y, sub(z, x));
  653.           long a = div(NEGATIVE_ONE, w);
  654.           long t = setLO(a, 0);
  655.           long s = add(ONE, mul(t, z));
  656.           return add(t, mul(a, add(s, mul(t, v))));
  657.         }
  658.       }
  659.     }
  660.     if (ix>=0x3FE59428) {
  661.       // |x|>=0.6744 
  662.       if (hx<0) {
  663.         x = negate(x);
  664.         y = negate(y);
  665.       }
  666.       x = add(sub(PIO4, x), sub(PIO4_LO, y));
  667.       y = ZERO;
  668.     }
  669.     long z = mul(x, x);
  670.     long w = mul(z, z);
  672.     // Break x^5*(T[1]+x^2*T[2]+...) into
  673.     // x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + 
  674.     // x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
  675.     long r = add(T1, mul(w, add(T3, mul(w, add(T5, mul(w, add(
  676.             T7, mul(w, add(T9, mul(w, T11))))))))));
  677.     long v = mul(z, add(T2, mul(w, add(T4, mul(w, add(T6, 
  678.             mul(w, add(T8, mul(w, add(T10, mul(w, T12)))))))))));
  679.     long s = mul(z, x);
  680.     r = add(add(y, mul(z, add(mul(s, add(r, v)), y))), mul(T0, s));
  681.     w = add(x, r);
  682.     if (ix >= 0x3FE59428) {
  683.       v = intToDouble(iy);
  684.       return mul(intToDouble(1-((hx>>30)&2)), sub(v, mul(
  685.               TWO, sub(x, sub(div(mul(w, w), add(
  686.               w, v)), r)))));
  687.     }
  688.     if (iy == 1) {
  689.       return w;
  690.     } else { // if allow error up to 2 ulp, simply return -1.0/(x+r) here 
  691.       // compute -1.0/(x+r) accurately 
  692.       z = setLO(w, 0);
  693.       v = sub(r, sub(z, x)); // z+v = r+x
  694.       long a = div(NEGATIVE_ONE, w); // a = -1.0/w 
  695.       long t = setLO(a, 0);
  696.       s = add(ONE, mul(t, z));
  697.       return add(t, mul(a, add(s, mul(t, v))));
  698.     }
  699.   }
  700.         
  701.   private static long S1 = 0xBFC5555555555549L; // -1.66666666666666324348e-01
  702.   private static long S2 = 0x3F8111111110F8A6L; // 8.33333333332248946124e-03 
  703.   private static long S3 = 0xBF2A01A019C161D5L; // -1.98412698298579493134e-04 
  704.   private static long S4 = 0x3EC71DE357B1FE7DL; // 2.75573137070700676789e-06  
  705.   private static long S5 = 0xBE5AE5E68A2B9CEBL; // -2.50507602534068634195e-08
  706.   private static long S6 = 0x3DE5D93A5ACFD57CL; // 1.58969099521155010221e-10
  708.   private static long kernelSin(long x, long y, int iy) {
  709.     int ix = getHI(x) & 0x7fffffff; // high word of x
  710.     if (ix < 0x3e400000) { // |x| < 2**-27 
  711.       return x;
  712.     }
  713.     long z = mul(x, x);
  714.     long v = mul(z, x);
  715.     long r = add(S2, mul(z, add(S3, mul(z, add(S4, mul(z, add(S5,
  716.             mul(z, S6))))))));
  717.     if (iy == 0) {
  718.       return add(x, mul(v, add(S1, mul(z, r))));
  719.     }
  720.     return sub(x, sub(sub(mul(z, 
  721.             sub(mul(ONE_HALF, y), mul(v, r))), y), mul(v, S1)));
  722.   }
  724.   private static final long C1  = 0x3fa555555555554cL; //  4.16666666666666019037e-02 
  725.   private static final long C2  = 0xbf56c16c16c15177L; // -1.38888888888741095749e-03 
  726.   private static final long C3  = 0x3efa01a019cb1590L; //  2.48015872894767294178e-05 
  727.   private static final long C4  = 0xbe927e4f809c52adL; // -2.75573143513906633035e-07 
  728.   private static final long C5  = 0x3e21ee9ebdb4b1c4L; //  2.08757232129817482790e-09 
  729.   private static final long C6  = 0xbda8fae9be8838d4L; // -1.13596475577881948265e-11 
  731.   private static long kernelCos(long x, long y) {
  732.     int ix = getHI(x) & 0x7fffffff; // ix = |x|'s high word
  733.     if (ix < 0x3e400000) {
  734.       // if x < 2**27
  735.       return ONE;
  736.     }
  737.     long z = mul(x, x);
  738.     long r = mul(z, add(C1, mul(z, add(C2, mul(z, add(C3, 
  739.             mul(z, add(C4, mul(z, add(C5, mul(z, C6)))))))))));
  740.     if (ix < 0x3FD33333) {
  741.       // if |x| < 0.3
  742.       return sub(ONE, sub(mul(ONE_HALF, z), sub(
  743.               mul(z, r), mul(x, y))));
  744.     }
  745.     long qx;
  746.     if (ix > 0x3fe90000) { // x > 0.78125
  747.       qx = 0x3fd2000000000000L; // 0.28125
  748.     } else {
  749.       qx = set(ix-0x00200000, 0); // x/4 
  750.     }
  751.     long hz = sub(mul(ONE_HALF, z), qx);
  752.     long a = sub(ONE, qx);
  753.     return sub(a, (sub(hz, sub(mul(z, r), mul(x, y)))));
  754.   }
  755.         
  756.   private static final long PIO2[] = {
  757.           0x3ff921fb40000000L,  // 1.57079625129699707031e+00 
  758.           0x3e74442d00000000L,  // 7.54978941586159635335e-08 
  759.           0x3cf8469880000000L,  // 5.39030252995776476554e-15 
  760.           0x3b78cc5160000000L,  // 3.28200341580791294123e-22 
  761.           0x39f01b8380000000L,  // 1.27065575308067607349e-29 
  762.           0x387a252040000000L,  // 1.22933308981111328932e-36 
  763.           0x36e3822280000000L,  // 2.73370053816464559624e-44 
  764.           0x3569f31d00000000L   // 2.16741683877804819444e-51 
  765.   };
  766.   
  767.    // Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
  768.   private static final int[] TWO_OVER_PI = {
  769.           0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, 
  770.           0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, 
  771.           0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, 
  772.           0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, 
  773.           0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, 
  774.           0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, 
  775.           0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, 
  776.           0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, 
  777.           0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, 
  778.           0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, 
  779.           0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b};
  780.           
  781.   private static final int[] NPIO2_HW = {
  782.           0x3ff921fb, 0x400921fb, 0x4012d97c, 0x401921fb, 0x401f6a7a, 0x4022d97c,
  783.           0x4025fdbb, 0x402921fb, 0x402c463a, 0x402f6a7a, 0x4031475c, 0x4032d97c,
  784.           0x40346b9c, 0x4035fdbb, 0x40378fdb, 0x403921fb, 0x403ab41b, 0x403c463a,
  785.           0x403dd85a, 0x403f6a7a, 0x40407e4c, 0x4041475c, 0x4042106c, 0x4042d97c,
  786.           0x4043a28c, 0x40446b9c, 0x404534ac, 0x4045fdbb, 0x4046c6cb, 0x40478fdb,
  787.           0x404858eb, 0x404921fb};
  788.           
  789.   private static final long TWO24    = 0x4170000000000000L; // 1.67772160000000000000e+07 
  791.   private static final long INV_PIO2 = 0x3fe45f306dc9c883L; // 6.36619772367581382433e-01 53 bits of 2/pi 
  792.   
  793.   private static final long PIO2_1   = 0x3ff921fb54400000L; // 1.57079632673412561417e+00 first  33 bit of pi/2 
  794.   private static final long PIO2_1T  = 0x3dd0b4611a626331L; // 6.07710050650619224932e-11 pi/2 - pio2_1 
  795.   private static final long PIO2_2   = 0x3dd0b4611a600000L; // 6.07710050630396597660e-11 second 33 bit of pi/2 
  796.   private static final long PIO2_2T  = 0x3ba3198a2e037073L; // 2.02226624879595063154e-21 pi/2 - (pio2_1+pio2_2) 
  797.   private static final long PIO2_3   = 0x3ba3198a2e000000L; // 2.02226624871116645580e-21 third  33 bit of pi/2 
  798.   private static final long PIO2_3T  = 0x397b839a252049c1L; // 8.47842766036889956997e-32 pi/2 - (pio2_1+pio2_2+pio2_3) 
  800.   // Return the remainder of x % pi/2 in y[0]+y[1].  This is rather complex.
  801.   // Perhaps it would make sense to do something simpler.
  802.   private static int remPio2(long x, long[] y) {
  803.     int hx = getHI(x); // high word of x
  804.     int ix = hx&0x7fffffff;
  805.     if (ix <= 0x3fe921fb) {
  806.       // |x| ~<= pi/4 , no need for reduction 
  807.       y[0] = x;
  808.       y[1] = ZERO;
  809.       return 0;
  810.     }
  812.     if (ix < 0x4002d97c) {
  813.       // |x| < 3pi/4, special case with n=+-1
  814.       long a = PIO2_1;
  815.       long b = (ix == 0x3ff921fb) ? PIO2_2T : PIO2_1T;
  816.       if (hx > 0) {
  817.         a = negate(a);
  818.         b = negate(b);
  819.       }
  820.       long z = add(x, a);
  821.       y[0] = add(z, b);
  822.       y[1] = add(sub(z, y[0]), b);
  823.       return (hx > 0) ? 1 : -1;
  824.     }
  826.     if (ix <= 0x413921fb) { // |x| ~<= 2^19*(pi/2), medium size 
  827.       long t = abs(x);
  828.       long fn = rint(mul(t, INV_PIO2));
  829.       int n = intValue(fn);
  830.       long r = sub(t, mul(fn, PIO2_1));
  831.       long w = mul(fn, PIO2_1T); // 1st round good to 85 bit
  832.       if ((n < 32) && (ix != NPIO2_HW[n-1])) {
  833.         y[0] = sub(r, w); // quick check no cancellation 
  834.       } else {
  835.         int j = ix >> 20;
  836.         y[0] = sub(r, w);
  837.         int i = j - (((getHI(y[0])) >> 20) & 0x7ff);
  838.         if (i > 16) { // 2nd iteration needed, good to 118 
  839.           t  = r;
  840.           w  = mul(fn, PIO2_2);
  841.           r  = sub(t, w);
  842.           w = sub(mul(fn, PIO2_2T), sub(sub(t, r), w));
  843.           y[0] = sub(r, w);
  844.           i = j - (((getHI(y[0])) >> 20) & 0x7ff);
  845.           if (i > 49) { // 3rd iteration need, 151 bits acc
  846.             t  = r;     // will cover all possible cases
  847.             w  = mul(fn, PIO2_3);
  848.             r  = sub(t, w);
  849.             w = sub(mul(fn, PIO2_3T), sub(sub(t, r), w));
  850.             y[0] = sub(r, w);
  851.           }
  852.         }
  853.       }
  854.       y[1] = sub(sub(r, y[0]), w);
  855.       if (hx < 0) {
  856.         y[0] = negate(y[0]);
  857.         y[1] = negate(y[1]);
  858.         return -n;
  859.       } else {
  860.         return n;
  861.       }
  862.     }
  864.     // all other (large) arguments
  865.     if (ix >= 0x7ff00000) {
  866.       // x is inf or NaN 
  867.       y[0] = y[1] = NaN;
  868.       return 0;
  869.     }
  870.     // set z = scalbn(|x|,ilogb(x)-23)
  871.     long z = getLO(x);
  872.     int e0 = (int) ((ix >> 20) - 1046);
  873.     z = setHI(z, ix - (e0 << 20));
  874.     long[] tx = new long[3];
  875.     for (int i=0; i<2; i++) {
  876.       tx[i] = intToDouble(intValue(z));
  877.       z     = scalbn(sub(z, tx[i]), 24);
  878.     }
  879.     tx[2] = z;
  880.     int nx = 3;
  881.     while (isZero(tx[nx-1])) nx--; // skip zero term 
  882.     int n = kernelRemPio2(tx, y, e0, nx);
  883.     if (hx < 0) {
  884.       y[0] = negate(y[0]);
  885.       y[1] = negate(y[1]);
  886.       return -n;
  887.     }
  888.     return n;
  889.   }
  891.   // Work even harder to compute mod pi/2.  This is extremely complex.
  892.   // Perhaps it would make sense to do something simpler.
  893.   private static int kernelRemPio2(long[] x, long[] y, int e0, int nx) {
  894.     // initialize jk
  895.     int jk = 4;
  896.     int jp = jk;
  898.     // determine jx,jv,q0, note that 3>q0 
  899.     int jx =  nx - 1;
  900.     int jv = (e0-3)/24;
  901.     if (jv < 0)  jv = 0;
  902.     int q0 =  e0-24*(jv+1);
  904.     // set up f[0] to f[jx+jk] where f[jx+jk] = two_over_pi[jv+jk]
  905.     int j = jv - jx;
  906.     int m = jx + jk;
  907.     long[] f = new long[20];
  908.     for (int i=0; i<=m; i++, j++) {
  909.       f[i] = (j<0) ? ZERO : intToDouble(TWO_OVER_PI[j]);
  910.     }
  912.     // compute q[0],q[1],...q[jk]
  913.     long[] q = new long[20];
  914.     for (int i=0; i<=jk; i++) {
  915.       long fw = ZERO;
  916.       for (j=0; j<=jx; j++) {
  917.         fw = add(fw, mul(x[j], f[jx+i-j]));
  918.       }
  919.       q[i] = fw;
  920.     }
  921.     
  922.     int jz = jk;
  924.     int n, ih;
  925.     long z;
  926.     int[] iq = new int[20];
  927.     boolean recompute;
  928.     do {
  929.       recompute = false;
  930.       // distill q[] into iq[] reversingly 
  931.       int i;
  932.       for (i=0, j=jz, z=q[jz]; j>0; i++, j--) {
  933.         long fw = intToDouble(intValue(scalbn(z, -24)));
  934.         iq[i] = intValue(sub(z, scalbn(fw, 24)));
  935.         z = add(q[j-1], fw);
  936.       }
  938.       // compute n
  939.       z = scalbn(z, q0); // actual value of z
  940.       z = sub(z, scalbn(floor(scalbn(z, -3)), 3)); // trim off integer >= 8
  941.       z = sub(z, mul(EIGHT, floor(mul(z, ONE_EIGHTH)))); // trim off integer >= 8
  942.       n = intValue(z);
  943.       z = sub(z, intToDouble(n));
  944.       ih = 0;
  945.       if (q0 > 0) { // need iq[jz-1] to determine n
  946.         i  = (iq[jz-1]>>(24-q0));
  947.         n += i;
  948.         iq[jz-1] -= i<<(24-q0);
  949.         ih = iq[jz-1]>>(23-q0);
  950.       } else if (q0==0) {
  951.         ih = iq[jz-1]>>23;
  952.       } else if(ge(z, ONE_HALF)) {
  953.         ih=2;
  954.       }
  955.       if (ih>0) { // q > 0.5
  956.         n += 1;
  957.         int carry = 0;
  958.         for(i=0; i<jz; i++) { // compute 1-q 
  959.           j = iq[i];
  960.           if (carry == 0) {
  961.             if (j != 0) {
  962.               carry = 1;
  963.               iq[i] = 0x1000000- j;
  964.             }
  965.           } else  iq[i] = 0xffffff - j;
  966.         }
  967.         if (q0 > 0) { // rare case: chance is 1 in 12
  968.           switch (q0) {
  969.           case 1:
  970.             iq[jz-1] &= 0x7fffff; break;
  971.           case 2:
  972.             iq[jz-1] &= 0x3fffff; break;
  973.           }
  974.         }
  975.         if (ih == 2) {
  976.           z = sub(ONE, z);
  977.           if (carry != 0) {
  978.             z = sub(z, scalbn(ONE,q0));
  979.           }
  980.         }
  981.       }
  982.       
  983.       // check if recomputation is needed 
  984.       if (isZero(z)) {
  985.         j = 0;
  986.         for (i = jz-1;  i >= jk;  i--)
  987.           j |= iq[i];
  988.         if (j == 0) { // need recomputation
  989.           int k;
  990.           for (k=1; iq[jk-k]==0; k++); // k = no. of terms needed
  991.           for (i=jz+1; i<=jz+k; i++) { // add q[jz+1] to q[jz+k]
  992.             f[jx+i] = intToDouble(TWO_OVER_PI[jv+i]);
  993.             long fw = ZERO;
  994.             for (j=0; j<=jx; j++)
  995.               fw = add(fw, mul(x[j], f[jx+i-j]));
  996.             q[i] = fw;
  997.           }
  998.           jz += k;
  999.           recompute = true;
  1000.         }
  1001.       }
  1002.     } while (recompute);
  1004.     // chop off zero terms 
  1005.     if (isZero(z)) {
  1006.       jz--;
  1007.       q0 -= 24;
  1008.       while (iq[jz] == 0) {
  1009.         jz--;
  1010.         q0 -= 24;
  1011.       }
  1012.     } else { // break z into 24-bit if necessary
  1013.       z = scalbn(z,-q0);
  1014.       if (ge(z, TWO24)) {
  1015.         long fw = intToDouble(intValue(scalbn(z, -24)));
  1016.         iq[jz] = intValue(sub(z, scalbn(fw, 24)));
  1017.         jz++;
  1018.         q0 += 24;
  1019.         iq[jz] = intValue(fw);
  1020.       } else iq[jz] = intValue(z);
  1021.     }
  1023.     // convert integer "bit" chunk to floating-point value
  1024.     long fw = scalbn(ONE, q0);
  1025.     for (int i=jz; i >= 0; i--) {
  1026.       q[i] = mul(fw, intToDouble(iq[i]));
  1027.       fw = scalbn(fw, -24);
  1028.     }
  1030.     // compute PIo2[0,...,jp]*q[jz,...,0] 
  1031.     long[] fq = new long[20];
  1032.     for (int i=jz; i>=0; i--) {
  1033.       fw = ZERO;
  1034.       for (int k=0; (k<=jp) && (k<=(jz-i)); k++)
  1035.         fw = add(fw, mul(PIO2[k], q[i+k]));
  1036.       fq[jz-i] = fw;
  1037.     }
  1039.     // compress fq[] into y[]
  1040.     fw = ZERO;
  1041.     for (int i=jz; i>=0; i--)
  1042.       fw = add(fw, fq[i]);
  1043.     y[0] = (ih==0)? fw : negate(fw);
  1044.     fw = sub(fq[0], fw);
  1045.     for (int i=1; i<=jz; i++)
  1046.       fw = add(fw, fq[i]);
  1047.     y[1] = (ih==0) ? fw: negate(fw); 
  1048.     return n&7;
  1049.   }
  1050.   
  1051.   private static final long PIO2_HI = 0x3FF921FB54442D18L;  // 1.57079632679489655800e+00
  1052.   private static final long PIO2_LO = 0x3C91A62633145C07L;  // 6.12323399573676603587e-17 
  1053.   private static final long PIO4_HI = 0x3FE921FB54442D18L;  // 7.85398163397448278999e-01
  1054.   // coefficient for R(x^2) 
  1055.   private static final long PS0     = 0x3fc5555555555555L;  //  1.66666666666666657415e-01 
  1056.   private static final long PS1     = 0xbfd4d61203eb6f7dL;  // -3.25565818622400915405e-01 
  1057.   private static final long PS2     = 0x3fc9c1550e884455L;  //  2.01212532134862925881e-01 
  1058.   private static final long PS3     = 0xbfa48228b5688f3bL;  // -4.00555345006794114027e-02 
  1059.   private static final long PS4     = 0x3f49efe07501b288L;  //  7.91534994289814532176e-04 
  1060.   private static final long PS5     = 0x3f023de10dfdf709L;  //  3.47933107596021167570e-05 
  1061.   private static final long QS1     = 0xc0033a271c8a2d4bL;  // -2.40339491173441421878e+00 
  1062.   private static final long QS2     = 0x40002ae59c598ac8L;  //  2.02094576023350569471e+00 
  1063.   private static final long QS3     = 0xbfe6066c1b8d0159L;  // -6.88283971605453293030e-01 
  1064.   private static final long QS4     = 0x3fb3b8c5b12e9282L;  //  7.70381505559019352791e-02 
  1066.   private static long pOverQ(long t) {
  1067.     return div(mul(t, add(PS0, mul(t, add(PS1, mul(t, add(PS2,
  1068.             mul(t, add(PS3, mul(t, add(PS4, mul(t, PS5))))))))))),
  1069.             add(ONE, mul(t, add(QS1, mul(t, add(QS2, mul(t,  
  1070.             add(QS3, mul(t, QS4)))))))));
  1071.   }
  1073.   /**
  1074.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#asin(double)">Math.asin(double)</a>.
  1075.    */
  1076.   public static final long asin(long d) {
  1077.     int hx = getHI(d);
  1078.     int ix = hx & 0x7fffffff;
  1079.     if (ix>= 0x3ff00000) { // |x|>= 1 
  1080.       if(((ix-0x3ff00000)|getLO(d))==0)
  1081.         // asin(1)=+-pi/2 with inexact
  1082.         return copySign(PIO2_HI, d);
  1083.       return NaN;  // asin(|x|>1) is NaN
  1084.     } else if (ix<0x3fe00000) { // |x|<0.5 
  1085.       if (ix<0x3e400000) { // if |x| < 2**-27
  1086.         return d;
  1087.       }
  1088.       long t = mul(d, d);
  1089.       long w = pOverQ(t);
  1090.       return add(d, mul(d, w));
  1091.     }
  1092.     // 1> |x|>= 0.5
  1093.     long w = sub(ONE, abs(d));
  1094.     long t = scalbn(w, -1);
  1095.     long s = sqrt(t);
  1096.     if(ix>=0x3FEF3333) { // if |x| > 0.975
  1097.       w = pOverQ(t);
  1098.       t = sub(PIO2_HI, sub(scalbn(add(s, mul(s, w)), 1), 
  1099.               PIO2_LO));
  1100.     } else {
  1101.       w  = setLO(s, 0);
  1102.       long c = div(sub(t, mul(w, w)), add(s, w));
  1103.       long r = pOverQ(t);
  1104.       long p = sub(mul(scalbn(s, 1), r), sub(PIO2_LO,
  1105.               scalbn(c, 1)));
  1106.       long q = sub(PIO4_HI, scalbn(w, 1));
  1107.       t = sub(PIO4_HI, sub(p, q));
  1108.     }    
  1109.     return ((hx>0) ? t : negate(t));
  1110.   }
  1112.   /**
  1113.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#acos(double)">Math.acos(double)</a>.
  1114.    */
  1115.   public static long acos(long d) {
  1116.     int hx = getHI(d);
  1117.     int ix = hx&0x7fffffff;
  1118.     if (ix >= 0x3ff00000) { // |x| >= 1 
  1119.       if (((ix-0x3ff00000)|getLO(d)) == 0) { // |x|==1
  1120.         if (hx>0) 
  1121.           return ZERO;
  1122.         else
  1123.           return PI; // acos(-1)= pi
  1124.       }
  1125.       return NaN; // acos(|x|>1) is NaN
  1126.     }
  1127.     if (ix < 0x3fe00000) { // |x| < 0.5
  1128.       if (ix <= 0x3c600000)
  1129.         return PIO2_HI; // if|x|<2**-57
  1130.       long z = mul(d, d);
  1131.       long r = pOverQ(z);
  1132.       return sub(PIO2_HI, sub(d, sub(PIO2_LO, mul(d, r))));
  1133.     } else if (hx<0) { // x < -0.5 
  1134.       long z = scalbn(add(ONE, d), -1);
  1135.       long s = sqrt(z);
  1136.       long r = pOverQ(z);
  1137.       long w = sub(mul(r, s), PIO2_LO);
  1138.       return sub(PI, scalbn(add(s, w), 1));
  1139.     } else { // x > 0.5 
  1140.       long z = scalbn(sub(ONE, d), -1);
  1141.       long s = sqrt(z);
  1142.       long df = setLO(s, 0);
  1143.       long c = div(sub(z, mul(df, df)), add(s, df));
  1144.       long r = pOverQ(z);
  1145.       long w = add(mul(r, s), c);
  1146.       return scalbn(add(df, w), 1);
  1147.     }
  1148.   }
  1149.   
  1150.   private static final long atanhi[] = {
  1151.           0x3fddac670561bb4fL,  // 4.63647609000806093515e-01 atan(0.5)hi 
  1152.           0x3fe921fb54442d18L,  // 7.85398163397448278999e-01 atan(1.0)hi
  1153.           0x3fef730bd281f69bL,  // 9.82793723247329054082e-01 atan(1.5)hi
  1154.           0x3ff921fb54442d18L   // 1.57079632679489655800e+00 atan(inf)hi
  1155.   };
  1157.   private static final long atanlo[] = {
  1158.           0x3c7a2b7f222f65e2L,  // 2.26987774529616870924e-17 atan(0.5)lo
  1159.           0x3c81a62633145c07L,  // 3.06161699786838301793e-17 atan(1.0)lo
  1160.           0x3c7007887af0cbbdL,  // 1.39033110312309984516e-17 atan(1.5)lo  
  1161.           0x3c91a62633145c07L   // 6.12323399573676603587e-17 atan(inf)lo  
  1162.   };
  1164.   private static final long AT0  = 0x3fd555555555550dL; //  3.33333333333329318027e-01 
  1165.   private static final long AT1  = 0xbfc999999998ebc4L; // -1.99999999998764832476e-01 
  1166.   private static final long AT2  = 0x3fc24924920083ffL; //  1.42857142725034663711e-01 
  1167.   private static final long AT3  = 0xbfbc71c6fe231671L; // -1.11111104054623557880e-01 
  1168.   private static final long AT4  = 0x3fb745cdc54c206eL; //  9.09088713343650656196e-02 
  1169.   private static final long AT5  = 0xbfb3b0f2af749a6dL; // -7.69187620504482999495e-02 
  1170.   private static final long AT6  = 0x3fb10d66a0d03d51L; //  6.66107313738753120669e-02 
  1171.   private static final long AT7  = 0xbfadde2d52defd9aL; // -5.83357013379057348645e-02 
  1172.   private static final long AT8  = 0x3fa97b4b24760debL; //  4.97687799461593236017e-02 
  1173.   private static final long AT9  = 0xbfa2b4442c6a6c2fL; // -3.65315727442169155270e-02 
  1174.   private static final long AT10 = 0x3f90ad3ae322da11L; //  1.62858201153657823623e-02 
  1175.         
  1176.   /**
  1177.    * Mimics <a href="http://java.sun.com/j2se/1.4.2/docs/api/java/lang/Math.html#atan(double)">Math.atan(double)</a>.
  1178.    */
  1179.   public static long atan(long d) {
  1180.     int hx = getHI(d);
  1181.     int ix = hx&0x7fffffff;
  1182.     int id;
  1183.     if(ix>=0x44100000) { // if |x| >= 2^66 
  1184.       if(ix>0x7ff00000||
  1185.         (ix==0x7ff00000&&(getLO(d)!=0)))
  1186.         return NaN;
  1187.       return (hx > 0) ? atanhi[3] : negate(atanhi[3]);
  1188.     } if (ix < 0x3fdc0000) { // |x| < 0.4375
  1189.       if (ix < 0x3e200000) { // |x| < 2^-29 
  1190.         return d;
  1191.       }
  1192.       id = -1;
  1193.     } else {
  1194.       d = abs(d);
  1195.       if (ix < 0x3ff30000) { // |x| < 1.1875 
  1196.         if (ix < 0x3fe60000) { // 7/16 <=|x|<11/16 
  1197.           id = 0;
  1198.           d = div(sub(scalbn(d, 1), ONE), add(TWO, d));
  1199.         } else { // 11/16<=|x|< 19/16 
  1200.           id = 1;
  1201.           d = div(sub(d, ONE), add(d, ONE));
  1202.         }
  1203.       } else {
  1204.         if (ix < 0x40038000) { // |x| < 2.4375 
  1205.           id = 2;
  1206.           d = div(sub(d, THREE_HALVES), add(ONE,
  1207.                   mul(THREE_HALVES, d)));
  1208.         } else { // 2.4375 <= |x| < 2^66 
  1209.           id = 3;
  1210.           d = div(NEGATIVE_ONE, d);
  1211.         }
  1212.       }
  1213.     }
  1214.     // end of argument reduction 
  1215.     long z = mul(d, d);
  1216.     long w = mul(z, z);
  1217.     // break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly 
  1218.     long s1 = mul(z, add(AT0, mul(w, add(AT2, mul(w, add(AT4, 
  1219.             mul(w, add(AT6, mul(w, add(AT8, mul(w, AT10)))))))))));
  1220.     long s2 = mul(w, add(AT1, mul(w, add(AT3, mul(w, add(AT5,
  1221.             mul(w, add(AT7, mul(w, AT9)))))))));
  1222.     if (id<0) {
  1223.       return sub(d, mul(d, add(s1, s2)));
  1224.     } else {
  1225.       z = sub(atanhi[id], sub(sub(mul(d, add(s1, s2)), 
  1226.               atanlo[id]), d));
  1227.       return (hx < 0) ? negate(z): z;
  1228.     }
  1229.   }
  1230.   
  1231.   /**
  1232.    * Returns the hyperbolic cosine of an angle.
  1233.    *
  1234.    * @param   d   an angle, in radians.
  1235.    * @return  the hyperbolic cosine of the argument.
  1236.    */
  1237.   public static long cosh(long d) {
  1238.     if (isNaN(d)) {
  1239.       return NaN;
  1240.     } else if (isInfinite(d)) {
  1241.       return POSITIVE_INFINITY;
  1242.     }
  1243.     
  1244.     int ix = getHI(d) & 0x7fffffff;
  1245.  
  1246.     // |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) 
  1247.     if(ix<0x3fd62e43) {
  1248.       long t = expm1(abs(d));
  1249.       long w = add(ONE, t);
  1250.       if (ix<0x3c800000) return w; // cosh(tiny) = 1 
  1251.       return add(ONE, div(mul(t, t), add(w, w)));
  1252.     }
  1254.     // |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2;
  1255.     if (ix < 0x40360000) {
  1256.       long t = exp(abs(d));
  1257.       return add(scalbn(t, -1), div(ONE_HALF, t));
  1258.     }
  1260.     // |x| in [22, log(maxdouble)] return half*exp(|x|) 
  1261.     if (ix < 0x40862E42)  return scalbn(exp(abs(d)), -1);
  1263.     // |x| in [log(maxdouble), overflowthresold]
  1264.     if (abs(d) <= 0x408633CE8fb9f87dL) {
  1265.       long w = exp(scalbn(abs(d), -1));
  1266.       long t = scalbn(w, -1);
  1267.       return mul(t, w);
  1268.     }
  1270.     // |x| > overflowthresold, cosh(x) overflow 
  1271.     return POSITIVE_INFINITY;
  1272.   }
  1274.   /**
  1275.    * Returns the hyperbolic sine of an angle.
  1276.    *
  1277.    * @param   d   an angle, in radians.
  1278.    * @return  the hyperbolic sine of the argument.
  1279.    */
  1280.   public static long sinh(long d) {
  1281.     if (isNaN(d)) {
  1282.       return NaN;
  1283.     } else if (isInfinite(d)) {
  1284.       return POSITIVE_INFINITY;
  1285.     }
  1286.     
  1287.     int jx = getHI(d);
  1288.     int ix = jx & 0x7fffffff;
  1289.     long h = ONE_HALF;
  1290.     if (jx < 0) h = negate(h);
  1292.     // |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) 
  1293.     if (ix < 0x40360000) { // |x|<22 
  1294.       if (ix<0x3e300000) // |x|<2**-28
  1295.         return d; // sinh(tiny) = tiny with inexact
  1296.       long t = expm1(abs(d));
  1297.       if(ix<0x3ff00000) return mul(h, sub(scalbn(t, 1), mul(t, 
  1298.               div(t, add(t, ONE)))));
  1299.       return mul(h, add(t, div(t, add(t, ONE))));
  1300.     }
  1302.     // |x| in [22, log(maxdouble)] return 0.5*exp(|x|) 
  1303.     if (ix < 0x40862E42)  return mul(h, exp(abs(d)));
  1305.     // |x| in [log(maxdouble), overflowthresold] 
  1306.     if (abs(d) <= 0x408633CE8fb9f87dL) {
  1307.       long w = exp(scalbn(abs(d), -1));
  1308.       long t = mul(h, w);
  1309.       return mul(t, w);
  1310.     }
  1312.     // |x| > overflowthresold, sinh(x) overflow 
  1313.     return copySign(POSITIVE_INFINITY, d);
  1314.   }
  1316.   /**
  1317.    * Returns the hyperbolic tangent of an angle.
  1318.    *
  1319.    * @param   d   an angle, in radians.
  1320.    * @return  the hyperbolic tangent of the argument.
  1321.    */
  1322.   public static long tanh(long d) {
  1323.     if (isNaN(d)) {
  1324.       return NaN;
  1325.     } else if (isInfinite(d)) {
  1326.       return copySign(POSITIVE_INFINITY, d);
  1327.     }
  1328.     
  1329.     int jx = getHI(d);
  1330.     int ix = jx & 0x7fffffff;
  1332.     // |x| < 22 
  1333.     long z;
  1334.     if (ix < 0x40360000) { // |x|<22 
  1335.         if (ix<0x3c800000) // |x|<2**-55 
  1336.           return d; // tanh(small) = small 
  1337.         if (ix>=0x3ff00000) { // |x|>=1  
  1338.           long t = expm1(scalbn(abs(d), 1));
  1339.           z = sub(ONE, div(TWO, add(t, TWO)));
  1340.         } else {
  1341.           long t = expm1(negate(mul(TWO, abs(d))));
  1342.           z = negate(div(t, add(t, TWO)));
  1343.         }
  1344.     // |x| > 22, return +-1 
  1345.     } else {
  1346.       z = ONE;
  1347.     }
  1348.     return (jx>=0) ? z : negate(z);
  1349.   }
  1350.   
  1351.   /**
  1352.    * Returns the arc hyperbolic cosine of an angle.
  1353.    *
  1354.    * @param   d   the value whose arc hyperbolic cosine is to be returned.
  1355.    * @return  the arc hyperbolic cosine of the argument.
  1356.    */
  1357.   public static long acosh(long d) {
  1358.     if (isNaN(d) || lt(d, ONE)) {
  1359.       return NaN;
  1360.     } else if (d == POSITIVE_INFINITY) {
  1361.       return d;
  1362.     } else if (d == ONE) {
  1363.       return ZERO;
  1364.     }
  1365.     int hx = getHI(d);
  1366.     if (hx > 0x41b00000) {
  1367.       return add(log(d), LN2);
  1368.     } else if (hx > 0x40000000) {
  1369.       long t = mul(d, d);
  1370.       return log(sub(scalbn(d, 1), div(ONE, add(d, 
  1371.               sqrt(sub(t, ONE))))));
  1372.     }
  1373.     long t = sub(d, ONE);
  1374.     return log1p(add(t, sqrt(add(scalbn(t, 1), 
  1375.             mul(t, t)))));
  1376.   }
  1378.   /**
  1379.    * Returns the arc hyperbolic sine of an angle.
  1380.    *
  1381.    * @param   d   the value whose arc hyperbolic sine is to be returned.
  1382.    * @return  the arc hyperbolic sine of the argument.
  1383.    */
  1384.   public static long asinh(long d) {
  1385.     int hx = getHI(d);
  1386.     int ix = hx&0x7fffffff;
  1387.     if(ix>=0x7ff00000) return d; // x is inf or NaN 
  1388.     if(ix< 0x3e300000) { // |x|<2**-28
  1389.       return d;
  1390.     } 
  1391.     long w;
  1392.     if(ix>0x41b00000) { /* |x| > 2**28 */
  1393.       w = add(log(abs(d)), LN2);
  1394.     } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
  1395.       long t = abs(d);
  1396.       w = log(add(scalbn(t, 1), div(ONE, add(sqrt(add(
  1397.               mul(d, d), ONE)), t))));
  1398.     } else { /* 2.0 > |x| > 2**-28 */
  1399.       long t = mul(d, d);
  1400.       w = log1p(add(abs(d), div(t, add(ONE, sqrt(add(
  1401.               ONE, t))))));
  1402.     }
  1403.     if(hx>0) return w; else return negate(w);
  1404.   }
  1405.   
  1406.   /**
  1407.    * Returns the arc hyperbolic tangent of an angle.
  1408.    *
  1409.    * @param   d   the value whose arc hyperbolic tangent is to be returned.
  1410.    * @return  the arc hyperbolic tangent of the argument.
  1411.    */
  1412.   public static long atanh(long d) {
  1413.     if (isNaN(d) || gt(d, ONE) || lt(d, NEGATIVE_ONE)) {
  1414.       return NaN;
  1415.     } 
  1416.     boolean negate = unpackSign(d);
  1417.     d = abs(d);
  1418.     if (d == ONE) {
  1419.       d = POSITIVE_INFINITY;
  1420.     } else if (lt(d, ONE_HALF)) {
  1421.       long t = add(d, d);
  1422.       d = scalbn(log1p(add(t, div(mul(t, d), sub(ONE, d)))), -1);
  1423.     } else {
  1424.       d = scalbn(log1p(div(add(d, d), sub(ONE, d))), -1);
  1425.     }
  1426.     if (negate) {
  1427.       d = negate(d);
  1428.     }
  1429.     return d;
  1430.   }
  1431.   
  1432.   /**
  1433.    * Returns the difference of d2 and d1, expressed as a percentage of d1.
  1434.    *
  1435.    * @param   d1   the "starting" value 
  1436.    * @param   d2   the "final" value
  1437.    * @return  100 * ((d2 - d1) / d1)
  1438.    */
  1439.   public static long percentChange(long d1, long d2) {
  1440.     return mul(div(sub(d2, d1), d1), ONE_HUNDRED);
  1441.   }
  1442.   
  1443.   /**
  1444.    * Returns d2 expressed as a percentage of d1
  1445.    *
  1446.    * @param   d1   the "base" value 
  1447.    * @param   d2   the other value
  1448.    * @return  100 * (d2 / d1)
  1449.    */
  1450.   public static long percentTotal(long d1, long d2) {
  1451.     return mul(div(d2, d1), ONE_HUNDRED);
  1452.   }
  1454.   /**
  1455.    * Returns the factorial of d.  If d is not a mathematical integer greater 
  1456.    * than or equal to zero, the return value is NaN.  Use the gamma function
  1457.    * for non-integer values.
  1458.    * <p>
  1459.    * This is a naive implentation.  TODO: make this better.
  1460.    *
  1461.    * @param d a <code>double</code> value.
  1462.    * @return d! 
  1463.    */
  1464.   public static long factorial(long d) {
  1465.     if (isZero(d)) {
  1466.       return ONE;
  1467.     } else if (! isPositiveInteger(d)) {
  1468.       return NaN;
  1469.     }
  1470.     return factorial(ONE, d, ONE);
  1471.   }
  1472.   
  1473.   private static boolean isPositiveInteger(long d) {
  1474.     return (! (unpackSign(d) || isNaN(d) || isInfinite(d) || isZero(d) 
  1475.             || (rint(d) != d)));
  1476.   }
  1477.   
  1478.   private static long factorial(long base, long d1, long d2) {
  1479.     while ((d1 != d2) && (base != POSITIVE_INFINITY)) {
  1480.       base = mul(base, d1);
  1481.       d1 = add(d1, NEGATIVE_ONE);
  1482.     }
  1483.     return base;
  1484.   }
  1485.   
  1486.   /**
  1487.    * Return the number of ways of obtaining an ordered subset of d2 elements 
  1488.    * from a set of d1 elements.  If d1 and d2 are not mathematical integers
  1489.    * than or equal to zero, or if d1 is less than d2, the return value is NaN.
  1490.    *
  1491.    * @param d1 a <code>double</code> value
  1492.    * @param d2 a <code>double</code> value
  1493.    * @return d1! / (d1 - d2)!
  1494.    * @see #factorial(long)
  1495.    * @see #combinations(long, long)
  1496.    */
  1497.   public static long permutations(long d1, long d2) {
  1498.     if (! (isPositiveInteger(d1) && isPositiveInteger(d2) && ge(d1, d2))) {
  1499.       return NaN;
  1500.     }
  1501.     return factorial(ONE, d1, sub(d1, d2));
  1502.   }
  1504.   /**
  1505.    * Return the number of ways of obtaining an unordered subset of d2 elements 
  1506.    * from a set of d1 elements.  Also known as the binomial coefficient.
  1507.    * If d1 and d2 are not mathematical integers greater
  1508.    * than or equal to zero, or if d1 is less than d2, the return value is NaN.
  1509.    *
  1510.    * @param d1 a <code>double</code> value
  1511.    * @param d2 a <code>double</code> value
  1512.    * @return d1! / (d2! * (d1 - d2)!)
  1513.    * @see #factorial(long)
  1514.    * @see #permutations(long, long)
  1515.    */
  1516.   public static long combinations(long d1, long d2) {
  1517.     if (! (isPositiveInteger(d1) && isPositiveInteger(d2) && ge(d1, d2))) {
  1518.       return NaN;
  1519.     }
  1520.     long d3 = sub(d1, d2);
  1521.     if (gt(d3, d2)) {
  1522.       long tmp = d3;
  1523.       d3 = d2;
  1524.       d2 = tmp;
  1525.     }
  1526.     if (isZero(d3)) {
  1527.       d3 = ONE;
  1528.     } else {
  1529.       d3 = div(ONE, factorial(ONE, d3, ONE));
  1530.     }
  1531.     return factorial(d3, d1, d2);
  1532.   }
  1534.   
  1535.   /**
  1536.    * Returns the complete gamma function of d.  
  1537.    *
  1538.    * @param d a <code>double</code> value
  1539.    * @return Gamma(d)
  1540.    * @see #lgamma(long)
  1541.    */
  1542.   public static long gamma(long d) {
  1543.     if (isNaN(d) || isZero(d) 
  1544.         || (lt(d, ZERO) && ((d == NEGATIVE_INFINITY) || (rint(d) == d)))) {
  1545.       return NaN;
  1546.     }
  1547.     return lgamma(d, true);
  1548.   }
  1549.   
  1550.   /**
  1551.    * Returns the natural logarithm of the absolute value of the gamma function
  1552.    * of d.
  1553.    *
  1554.    * @param d a <code>double</code> value
  1555.    * @return Log(|Gamma(d)|)
  1556.    * @see #gamma(long)
  1557.    */
  1558.   public static long lgamma(long d) {
  1559.     return lgamma(d, false);
  1560.   }
  1561.   
  1562.   private static final long A0  = 0x3FB3C467E37DB0C8L; // 7.72156649015328655494e-02
  1563.   private static final long A1  = 0x3FD4A34CC4A60FADL; // 3.22467033424113591611e-01
  1564.   private static final long A2  = 0x3FB13E001A5562A7L; // 6.73523010531292681824e-02
  1565.   private static final long A3  = 0x3F951322AC92547BL; // 2.05808084325167332806e-02
  1566.   private static final long A4  = 0x3F7E404FB68FEFE8L; // 7.38555086081402883957e-03
  1567.   private static final long A5  = 0x3F67ADD8CCB7926BL; // 2.89051383673415629091e-03
  1568.   private static final long A6  = 0x3F538A94116F3F5DL; // 1.19270763183362067845e-03
  1569.   private static final long A7  = 0x3F40B6C689B99C00L; // 5.10069792153511336608e-04
  1570.   private static final long A8  = 0x3F2CF2ECED10E54DL; // 2.20862790713908385557e-04
  1571.   private static final long A9  = 0x3F1C5088987DFB07L; // 1.08011567247583939954e-04
  1572.   private static final long A10 = 0x3EFA7074428CFA52L; // 2.52144565451257326939e-05
  1573.   private static final long A11 = 0x3F07858E90A45837L; // 4.48640949618915160150e-05
  1574.   private static final long TC = 0x3FF762D86356BE3FL; // 1.46163214496836224576e+00
  1575.   private static final long TF = 0xBFBF19B9BCC38A42L; // -1.21486290535849611461e-01
  1576.   // tt = -(tail of tf)
  1577.   private static final long TT = 0xBC50C7CAA48A971FL; // -3.63867699703950536541e-18
  1578.   private static final long TB0 = 0x3FDEF72BC8EE38A2L; // 4.83836122723810047042e-01
  1579.   private static final long TB1 = 0xBFC2E4278DC6C509L; // -1.47587722994593911752e-01
  1580.   private static final long TB2 = 0x3FB08B4294D5419BL; // 6.46249402391333854778e-02
  1581.   private static final long TB3 = 0xBFA0C9A8DF35B713L; // -3.27885410759859649565e-02
  1582.   private static final long TB4 = 0x3F9266E7970AF9ECL; // 1.79706750811820387126e-02
  1583.   private static final long TB5 = 0xBF851F9FBA91EC6AL; // -1.03142241298341437450e-02
  1584.   private static final long TB6 = 0x3F78FCE0E370E344L; // 6.10053870246291332635e-03
  1585.   private static final long TB7 = 0xBF6E2EFFB3E914D7L; // -3.68452016781138256760e-03
  1586.   private static final long TB8 = 0x3F6282D32E15C915L; // 2.25964780900612472250e-03
  1587.   private static final long TB9 = 0xBF56FE8EBF2D1AF1L; // -1.40346469989232843813e-03
  1588.   private static final long TB10 = 0x3F4CDF0CEF61A8E9L; // 8.81081882437654011382e-04
  1589.   private static final long TB11 = 0xBF41A6109C73E0ECL; // -5.38595305356740546715e-04
  1590.   private static final long TB12 = 0x3F34AF6D6C0EBBF7L; // 3.15632070903625950361e-04
  1591.   private static final long TB13 = 0xBF347F24ECC38C38L; // -3.12754168375120860518e-04
  1592.   private static final long TB14 = 0x3F35FD3EE8C2D3F4L; // 3.35529192635519073543e-04
  1593.   private static final long U0  = 0xBFB3C467E37DB0C8L; // -7.72156649015328655494e-02
  1594.   private static final long U1  = 0x3FE4401E8B005DFFL; // 6.32827064025093366517e-01
  1595.   private static final long U2  = 0x3FF7475CD119BD6FL; // 1.45492250137234768737e+00
  1596.   private static final long U3  = 0x3FEF497644EA8450L; // 9.77717527963372745603e-01
  1597.   private static final long U4  = 0x3FCD4EAEF6010924L; // 2.28963728064692451092e-01
  1598.   private static final long U5  = 0x3F8B678BBF2BAB09L; // 1.33810918536787660377e-02
  1599.   private static final long V1  = 0x4003A5D7C2BD619CL; // 2.45597793713041134822e+00, /* 
  1600.   private static final long V2  = 0x40010725A42B18F5L; // 2.12848976379893395361e+00, /* 
  1601.   private static final long V3  = 0x3FE89DFBE45050AFL; // 7.69285150456672783825e-01, /* 
  1602.   private static final long V4  = 0x3FBAAE55D6537C88L; // 1.04222645593369134254e-01, /* 
  1603.   private static final long V5  = 0x3F6A5ABB57D0CF61L; // 3.21709242282423911810e-03, /* 
  1604.   private static final long SB0  = 0xBFB3C467E37DB0C8L; // 7.72156649015328655494e-02, /* 
  1605.   private static final long SB1  = 0x3FCB848B36E20878L; // 2.14982415960608852501e-01, /* 
  1606.   private static final long SB2  = 0x3FD4D98F4F139F59L; // 3.25778796408930981787e-01, /* 
  1607.   private static final long SB3  = 0x3FC2BB9CBEE5F2F7L; // 1.46350472652464452805e-01, /* 
  1608.   private static final long SB4  = 0x3F9B481C7E939961L; // 2.66422703033638609560e-02, /* 
  1609.   private static final long SB5  = 0x3F5E26B67368F239L; // 1.84028451407337715652e-03, /* 
  1610.   private static final long SB6  = 0x3F00BFECDD17E945L; // 3.19475326584100867617e-05, /* 
  1611.   private static final long R1  = 0x3FF645A762C4AB74L; // 1.39200533467621045958e+00, /* 
  1612.   private static final long R2  = 0x3FE71A1893D3DCDCL; // 7.21935547567138069525e-01, /* 
  1613.   private static final long R3  = 0x3FC601EDCCFBDF27L; // 1.71933865632803078993e-01, /*
  1614.   private static final long R4  = 0x3F9317EA742ED475L; // 1.86459191715652901344e-02, /* 
  1615.   private static final long R5  = 0x3F497DDACA41A95BL; // 7.77942496381893596434e-04, /* 
  1616.   private static final long R6  = 0x3EDEBAF7A5B38140L; // 7.32668430744625636189e-06, /* 
  1617.   private static final long W0  = 0x3FDACFE390C97D69L; // 4.18938533204672725052e-01, /* 
  1618.   private static final long W1  = 0x3FB555555555553BL; // 8.33333333333329678849e-02, /* 
  1619.   private static final long W2  = 0xBF66C16C16B02E5CL; // -2.77777777728775536470e-03, /* 
  1620.   private static final long W3  = 0x3F4A019F98CF38B6L; // 7.93650558643019558500e-04, /* 
  1621.   private static final long W4  = 0xBF4380CB8C0FE741L; // -5.95187557450339963135e-04, /* 
  1622.   private static final long W5  = 0x3F4B67BA4CDAD5D1L; // 8.36339918996282139126e-04, /* 
  1623.   private static final long W6  = 0xBF5AB89D0B9E43E4L; // -1.63092934096575273989e-03; /* 
  1625.   /**
  1626.    * Return gamma(x) or ln(|gamma(x)|).  First ln(|gamma(x)|) is computed.  Then,
  1627.    * if exp is true, the exponential of that value is computed, and the sign
  1628.    * is restored.
  1629.    */
  1630.   private static long lgamma(long x, boolean exp) {
  1631.     int hx = getHI(x);
  1632.     int lx = getLO(x);
  1634.     // purge off +-inf, NaN, +-0, and negative arguments 
  1635.     boolean negative = false;
  1636.     int ix = hx&0x7fffffff;
  1637.     if(ix>=0x7ff00000) return mul(x, x); 
  1638.     if((ix|lx)==0) return POSITIVE_INFINITY; // one/zero
  1639.     if(ix<0x3b900000) { // |x|<2**-70, return -log(|x|)
  1640.       if(hx<0) {
  1641.         negative = true;
  1642.         x = negate(x);
  1643.       }
  1644.       x = negate(log(x));
  1645.       if (exp) {
  1646.         x = exp(x);
  1647.         if (negative) {
  1648.           x = negate(x);
  1649.         }
  1650.       }
  1651.       return x;
  1652.     }
  1653.     long t, nadj;
  1654.     if(hx<0) {
  1655.         if(ix>=0x43300000) // |x|>=2**52, must be -integer 
  1656.             return POSITIVE_INFINITY; // one/zero;
  1657.         t = sinPi(x);
  1658.         if(isZero(t)) return POSITIVE_INFINITY; // one/zero; // -integer 
  1659.         nadj = log(div(PI, abs(mul(t, x))));
  1660.         if (lt(t, ZERO)) {
  1661.           negative = true;
  1662.         }
  1663.         x = negate(x);
  1664.     } else {
  1665.       nadj = ZERO;
  1666.     }
  1668.     // purge off 1 and 2 
  1669.     long r, y;
  1670.     int i;
  1671.     if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = ZERO;
  1672.     // for x < 2.0 
  1673.     else if(ix<0x40000000) {
  1674.       if(ix<=0x3feccccc) { // lgamma(x) = lgamma(x+1)-log(x) 
  1675.         r = negate(log(x));
  1676.         if(ix>=0x3FE76944) {
  1677.           y = sub(ONE, x); 
  1678.           i=0;
  1679.         } else if(ix>=0x3FCDA661) {
  1680.           y= sub(x, sub(TC, ONE)); 
  1681.           i=1;
  1682.         } else {
  1683.           y = x; 
  1684.           i=2;
  1685.         }
  1686.       } else {
  1687.         r = ZERO;
  1688.         if(ix>=0x3FFBB4C3) {
  1689.           y=sub(TWO, x);
  1690.           i=0;
  1691.         } // [1.7316,2]
  1692.         else if(ix>=0x3FF3B4C4) {
  1693.           y=sub(x, TC);
  1694.           i=1;
  1695.         } // [1.23,1.73] 
  1696.         else {
  1697.           y=sub(x, ONE);
  1698.           i=2;
  1699.         }
  1700.       }
  1701.       long w, z, p1, p2, p3, p;
  1702.       switch(i) {
  1703.         case 0:
  1704.           z = mul(y, y);
  1705.           p1 = add(A0, mul(z, add(A2, mul(z, add(A4, mul(
  1706.                   z, add(A6, mul(z, add(A8, mul(z, A10))))))))));
  1707.           p2 = mul(z, add(A1, mul(z, add(A3, mul(z, add(
  1708.                   A5, mul(z, add(A7, mul(z, add(A9, mul(z, A11)))))))))));
  1709.           p = add(mul(y, p1), p2);
  1710.           r = add(r, sub(p, scalbn(y, -1))); 
  1711.           break;
  1712.         case 1:
  1713.           z = mul(y, y);
  1714.           w = mul(z, y);
  1715.           p1 = add(TB0, mul(w, add(TB3, mul(w, add(TB6, mul(w,
  1716.                   add(TB9, mul(w, TB12)))))))); // parallel comp 
  1717.           p2 = add(TB1, mul(w, add(TB4, mul(w, add(TB7, mul(w,
  1718.                   add(TB10, mul(w, TB13))))))));
  1719.           p3 = add(TB2, mul(w, add(TB5, mul(w, add(TB8, mul(w,
  1720.                   add(TB11, mul(w, TB14))))))));
  1721.           p = sub(mul(z, p1), sub(TT, mul(w, add(p2, 
  1722.                   mul(y, p3)))));
  1723.           r = add(r, add(TF, p));
  1724.           break;
  1725.         case 2:
  1726.           p1 = mul(y, add(U0, mul(y, add(U1, mul(y, add(U2,
  1727.                   mul(y, add(U3, mul(y, add(U4, mul(y, U5)))))))))));
  1728.           p2 = add(ONE, mul(y, add(V1, mul(y, add(V2, mul(y,
  1729.                   add(V3, mul(y, add(V4, mul(y, V5))))))))));
  1730.           r = add(r, add(negate(scalbn(y, -1)), div(p1, p2)));
  1731.       }
  1732.     }
  1733.     else if(ix<0x40200000) { // x < 8.0 
  1734.       i = intValue(x);
  1735.       t = ZERO;
  1736.       y = sub(x, intToDouble(i));
  1737.       long p = mul(y, add(SB0, mul(y, add(SB1, mul(y, add(SB2, 
  1738.               mul(y, add(SB3, mul(y, add(SB4, mul(y, add(SB5,
  1739.               mul(y, SB6)))))))))))));
  1740.       long q = add(ONE, mul(y, add(R1, mul(y, add(R2, mul(y, add(R3,
  1741.               mul(y, add(R4, mul(y, add(R5, mul(y, R6))))))))))));
  1742.       r = add(scalbn(y, -1), div(p, q));
  1743.       long z = ONE; // lgamma(1+s) = log(s) + lgamma(s) 
  1744.       switch(i) {
  1745.         case 7: z = mul(z, add(y, SIX)); // FALLTHRU 
  1746.         case 6: z = mul(z, add(y, FIVE)); // FALLTHRU 
  1747.         case 5: z = mul(z, add(y, FOUR)); // FALLTHRU 
  1748.         case 4: z = mul(z, add(y, THREE)); // FALLTHRU 
  1749.         case 3: z = mul(z, add(y, TWO)); // FALLTHRU 
  1750.                 r = add(r, log(z)); 
  1751.                 break;
  1752.         }
  1753.     // 8.0 <= x < 2**58 
  1754.     } else if (ix < 0x43900000) {
  1755.         t = log(x);
  1756.         long z = div(ONE, x);
  1757.         y = mul(z, z);
  1758.         long w = add(W0, mul(z, add(W1, mul(y, add(W2, mul(y, 
  1759.                 add(W3, mul(y, add(W4, mul(y, add(W5, mul(y, 
  1760.                 W6))))))))))));
  1761.         r = add(mul(sub(x, ONE_HALF), sub(t, ONE)), w);
  1762.     } else {
  1763.     // 2**58 <= x <= inf 
  1764.       r = mul(x, sub(log(x), ONE));
  1765.     }
  1766.     if(hx<0) r = sub(nadj, r);
  1767.     if (exp) {
  1768.       r = exp(r);
  1769.       if (negative) {
  1770.         r = negate(r);
  1771.       }
  1772.     }
  1773.     return r;
  1774.   }
  1775.  
  1776.   private static final long TWO52 = 0x4330000000000000L; // 4.50359962737049600000e+15
  1778.   /** used by lgamma */
  1779.   private static long sinPi(long x) {
  1780.     int ix = 0x7fffffff & getHI(x);
  1781.     if(ix<0x3fd00000) return kernelSin(mul(PI, x), ZERO, 0);
  1782.     long y = negate(x); // x is assume negative
  1784.     // argument reduction, make sure inexact flag not raised if input is 
  1785.     // an integer
  1786.     long z = floor(y);
  1787.     int n;
  1788.     if (ne(z, y)) { // inexact anyway
  1789.       y = scalbn(y, -1);
  1790.       y = scalbn(sub(y, floor(y)), 1); // y = |x| mod 2.0
  1791.       n = intValue(scalbn(y, 2));
  1792.     } else {
  1793.       if(ix>=0x43400000) {
  1794.         y = ZERO; n = 0; // y must be even
  1795.       } else {
  1796.         if(ix<0x43300000) z = add(y, TWO52); // exact
  1797.         n = getLO(z) & 1; // lower word of z 
  1798.         y = intToDouble(n);
  1799.         n <<= 2;
  1800.       }
  1801.     }
  1802.     switch (n) {
  1803.       case 0:   y =  kernelSin(mul(PI, y),ZERO,0); break;
  1804.       case 1:   
  1805.       case 2:   y =  kernelCos(mul(PI, sub(ONE_HALF, y)),ZERO); break;
  1806.       case 3:  
  1807.       case 4:   y =  kernelSin(mul(PI, sub(ONE, y)),ZERO,0); break;
  1808.       case 5:
  1809.       case 6:   y = negate(kernelCos(mul(PI, sub(y, THREE_HALVES)),ZERO)); break;
  1810.       default:  y =  kernelSin(mul(PI, sub(y, TWO)),ZERO,0); break;
  1811.     }
  1812.     return negate(y);
  1813.   }
  1814.      
  1815.   
  1816.   /////////////////////////////////////////////////////////////////////////////
  1817.   // Instance members
  1818.   /////////////////////////////////////////////////////////////////////////////
  1820.   private final long value;
  1822.   
  1823.   /////////////////////////////////////////////////////////////////////////////
  1824.   // Constructors
  1825.   /////////////////////////////////////////////////////////////////////////////
  1827.   /**
  1828.    * Constructs a newly-allocated <code>MicroDouble</code> object that represents 
  1829.    * the argument. 
  1830.    *
  1831.    * @param d the <code>double</code> value to be represented by the <code>MicroDouble</code>.
  1832.    */
  1833.   public MicroDouble(long d) {
  1834.     // canonicalize NaN values so that hashCode() and equals() can be simpler
  1835.     if (isNaN(d)) {
  1836.       d = NaN;
  1837.     }
  1838.     value = d;
  1839.   }
  1840.   
  1841.   /**
  1842.    * Constructs a newly-allocated <code>MicroDouble</code> object that represents 
  1843.    * the argument.
  1844.    *
  1845.    * @param s a <code>String</code> to be converted to a <code>MicroDouble</code>.
  1846.    * @throws NumberFormatException if the <code>String</code> does not contain a
  1847.    *         parsable number.
  1848.    * @see #parseDouble(String)
  1849.    */
  1850.   public MicroDouble(String s) {
  1851.     this(parseDouble(s));
  1852.   }
  1854.   
  1855.   /////////////////////////////////////////////////////////////////////////////
  1856.   // Instance methods
  1857.   /////////////////////////////////////////////////////////////////////////////
  1859.   /**
  1860.    * Returns the <code>double</code> value of this <code>MicroDouble</code>
  1861.    * object.
  1862.    */
  1863.   public long doubleValue() {
  1864.     return value;
  1865.   }
  1867.   /**
  1868.    * Returns a String object representing this MicroDouble's value.
  1869.    * Equivalent to <code>toString(doubleValue())</code>.
  1870.    *
  1871.    * @see #toString(long)
  1872.    */
  1873.   public String toString() {
  1874.     return toString(value);
  1875.   }
  1876.     
  1877.   /**
  1878.    * Returns a hash code for this <code>MicroDouble</code> object.
  1879.    * Equivalent to <code>(int) (doubleValue() ^ (doubleValue >>> 32))</code>
  1880.    */
  1881.   public int hashCode() {
  1882.     return ((int) value) ^ ((int) (value >>> 32));
  1883.   }
  1884.   
  1885.   /**
  1886.    * Compares this object against the specified object.
  1887.    * Equivalent to 
  1888.    * <code>((obj instanceof MicroDouble) && (compare(((MicroDouble) obj).doubleValue(), doubleValue()) == 0))</code>
  1889.    * 
  1890.    * @see #compare(long, long)
  1891.    */
  1892.   public boolean equals(Object obj) {
  1893.     return ((obj instanceof MicroDouble) 
  1894.             && (((MicroDouble) obj).value == value));
  1895.   }
  1896.   
  1897. }