Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /* ef_rem_pio2.c -- float version of e_rem_pio2.c
  2.  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  3.  */
  4.  
  5. /*
  6.  * ====================================================
  7.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  8.  *
  9.  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10.  * Permission to use, copy, modify, and distribute this
  11.  * software is freely granted, provided that this notice
  12.  * is preserved.
  13.  * ====================================================
  14.  *
  15.  */
  16.  
  17. /* __ieee754_rem_pio2f(x,y)
  18.  *
  19.  * return the remainder of x rem pi/2 in y[0]+y[1]
  20.  * use __kernel_rem_pio2f()
  21.  */
  22.  
  23. #include "fdlibm.h"
  24.  
  25. /*
  26.  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
  27.  */
  28. #ifdef __STDC__
  29. static const __int32_t two_over_pi[] = {
  30. #else
  31. static __int32_t two_over_pi[] = {
  32. #endif
  33. 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
  34. 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
  35. 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
  36. 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
  37. 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
  38. 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
  39. 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
  40. 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
  41. 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
  42. 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
  43. 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
  44. 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
  45. 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
  46. 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
  47. 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
  48. 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
  49. 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
  50. 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
  51. 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
  52. 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
  53. 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
  54. 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
  55. };
  56.  
  57. /* This array is like the one in e_rem_pio2.c, but the numbers are
  58.    single precision and the last 8 bits are forced to 0.  */
  59. #ifdef __STDC__
  60. static const __int32_t npio2_hw[] = {
  61. #else
  62. static __int32_t npio2_hw[] = {
  63. #endif
  64. 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
  65. 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
  66. 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
  67. 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
  68. 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
  69. 0x4242c700, 0x42490f00
  70. };
  71.  
  72. /*
  73.  * invpio2:  24 bits of 2/pi
  74.  * pio2_1:   first  17 bit of pi/2
  75.  * pio2_1t:  pi/2 - pio2_1
  76.  * pio2_2:   second 17 bit of pi/2
  77.  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
  78.  * pio2_3:   third  17 bit of pi/2
  79.  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  80.  */
  81.  
  82. #ifdef __STDC__
  83. static const float
  84. #else
  85. static float
  86. #endif
  87. zero =  0.0000000000e+00, /* 0x00000000 */
  88. half =  5.0000000000e-01, /* 0x3f000000 */
  89. two8 =  2.5600000000e+02, /* 0x43800000 */
  90. invpio2 =  6.3661980629e-01, /* 0x3f22f984 */
  91. pio2_1  =  1.5707855225e+00, /* 0x3fc90f80 */
  92. pio2_1t =  1.0804334124e-05, /* 0x37354443 */
  93. pio2_2  =  1.0804273188e-05, /* 0x37354400 */
  94. pio2_2t =  6.0770999344e-11, /* 0x2e85a308 */
  95. pio2_3  =  6.0770943833e-11, /* 0x2e85a300 */
  96. pio2_3t =  6.1232342629e-17; /* 0x248d3132 */
  97.  
  98. #ifdef __STDC__
  99.         __int32_t __ieee754_rem_pio2f(float x, float *y)
  100. #else
  101.         __int32_t __ieee754_rem_pio2f(x,y)
  102.         float x,y[];
  103. #endif
  104. {
  105.         float z,w,t,r,fn;
  106.         float tx[3];
  107.         __int32_t i,j,n,ix,hx;
  108.         int e0,nx;
  109.  
  110.         GET_FLOAT_WORD(hx,x);
  111.         ix = hx&0x7fffffff;
  112.         if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
  113.             {y[0] = x; y[1] = 0; return 0;}
  114.         if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
  115.             if(hx>0) {
  116.                 z = x - pio2_1;
  117.                 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
  118.                     y[0] = z - pio2_1t;
  119.                     y[1] = (z-y[0])-pio2_1t;
  120.                 } else {                /* near pi/2, use 24+24+24 bit pi */
  121.                     z -= pio2_2;
  122.                     y[0] = z - pio2_2t;
  123.                     y[1] = (z-y[0])-pio2_2t;
  124.                 }
  125.                 return 1;
  126.             } else {    /* negative x */
  127.                 z = x + pio2_1;
  128.                 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
  129.                     y[0] = z + pio2_1t;
  130.                     y[1] = (z-y[0])+pio2_1t;
  131.                 } else {                /* near pi/2, use 24+24+24 bit pi */
  132.                     z += pio2_2;
  133.                     y[0] = z + pio2_2t;
  134.                     y[1] = (z-y[0])+pio2_2t;
  135.                 }
  136.                 return -1;
  137.             }
  138.         }
  139.         if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
  140.             t  = fabsf(x);
  141.             n  = (__int32_t) (t*invpio2+half);
  142.             fn = (float)n;
  143.             r  = t-fn*pio2_1;
  144.             w  = fn*pio2_1t;    /* 1st round good to 40 bit */
  145.             if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) { 
  146.                 y[0] = r-w;     /* quick check no cancellation */
  147.             } else {
  148.                 __uint32_t high;
  149.                 j  = ix>>23;
  150.                 y[0] = r-w;
  151.                 GET_FLOAT_WORD(high,y[0]);
  152.                 i = j-((high>>23)&0xff);
  153.                 if(i>8) {  /* 2nd iteration needed, good to 57 */
  154.                     t  = r;
  155.                     w  = fn*pio2_2;    
  156.                     r  = t-w;
  157.                     w  = fn*pio2_2t-((t-r)-w); 
  158.                     y[0] = r-w;
  159.                     GET_FLOAT_WORD(high,y[0]);
  160.                     i = j-((high>>23)&0xff);
  161.                     if(i>25)  { /* 3rd iteration need, 74 bits acc */
  162.                         t  = r; /* will cover all possible cases */
  163.                         w  = fn*pio2_3;
  164.                         r  = t-w;
  165.                         w  = fn*pio2_3t-((t-r)-w);     
  166.                         y[0] = r-w;
  167.                     }
  168.                 }
  169.             }
  170.             y[1] = (r-y[0])-w;
  171.             if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
  172.             else         return n;
  173.         }
  174.     /*
  175.      * all other (large) arguments
  176.      */
  177.         if(!FLT_UWORD_IS_FINITE(ix)) {
  178.             y[0]=y[1]=x-x; return 0;
  179.         }
  180.     /* set z = scalbn(|x|,ilogb(x)-7) */
  181.         e0      = (int)((ix>>23)-134);  /* e0 = ilogb(z)-7; */
  182.         SET_FLOAT_WORD(z, ix - ((__int32_t)e0<<23));
  183.         for(i=0;i<2;i++) {
  184.                 tx[i] = (float)((__int32_t)(z));
  185.                 z     = (z-tx[i])*two8;
  186.         }
  187.         tx[2] = z;
  188.         nx = 3;
  189.         while(tx[nx-1]==zero) nx--;     /* skip zero term */
  190.         n  =  __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
  191.         if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
  192.         return n;
  193. }
  194.