src/sf-trig.c

   1 #include <gnumeric-config.h>
   2 #include <sf-trig.h>
   3 #include <mathfunc.h>
   4
   5 /* ------------------------------------------------------------------------- */
   6
   7 static gnm_float
   8 gnm_cot_helper (volatile gnm_float *x)
   9 {
  10         gnm_float s = gnm_sin (*x);
  11         gnm_float c = gnm_cos (*x);
  12
  13         if (s == 0)
  14                 return gnm_nan;
  15         else
  16                 return c / s;
  17 }
  18
  19 /**
  20  * gnm_cot:
  21  * @x: an angle in radians
  22  *
  23  * Returns: The co-tangent of the given angle.
  24  */
  25 gnm_float
  26 gnm_cot (gnm_float x)
  27 {
  28         /* See http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59089 */
  29         return gnm_cot_helper (&x);
  30 }
  31
  32 /**
  33  * gnm_acot:
  34  * @x: a number
  35  *
  36  * Returns: The inverse co-tangent of the given number.
  37  */
  38 gnm_float
  39 gnm_acot (gnm_float x)
  40 {
  41         if (gnm_finite (x)) {
  42                 if (x == 0)
  43                         return M_PIgnum / 2;
  44                 return gnm_atan (1 / x);
  45         } else {
  46                 /* +inf -> +0 */
  47                 /* -Inf -> -0 */
  48                 /* +-NaN -> +-NaN */
  49                 return 1 / x;
  50         }
  51 }
  52
  53 /**
  54  * gnm_coth:
  55  * @x: a number.
  56  *
  57  * Returns: The hyperbolic co-tangent of the given number.
  58  */
  59 gnm_float
  60 gnm_coth (gnm_float x)
  61 {
  62         return 1 / gnm_tanh (x);
  63 }
  64
  65 /**
  66  * gnm_acoth:
  67  * @x: a number
  68  *
  69  * Returns: The inverse hyperbolic co-tangent of the given number.
  70  */
  71 gnm_float
  72 gnm_acoth (gnm_float x)
  73 {
  74         return (gnm_abs (x) > 2)
  75                 ? gnm_log1p (2 / (x - 1)) / 2
  76                 : gnm_log ((x - 1) / (x + 1)) / -2;
  77 }
  78
  79 /* ------------------------------------------------------------------------- */
  80
  81
  82 static gnm_float
  83 reduce_pi_simple (gnm_float x, int *pk, int kbits)
  84 {
  85         static const gnm_float two_over_pi = GNM_const(0.63661977236758134307553505349);
  86         static const gnm_float pi_over_two[] = {
  87                 +0x1.921fb5p+0,
  88                 +0x1.110b46p-26,
  89                 +0x1.1a6263p-54,
  90                 +0x1.8a2e03p-81,
  91                 +0x1.c1cd12p-107
  92         };
  93         int i;
  94         gnm_float k = gnm_floor (x * gnm_ldexp (two_over_pi, kbits - 2) + 0.5);
  95         gnm_float xx = 0;
  96
  97         g_assert (k < (1 << 26));
  98         *pk = (int)k;
  99
 100         if (k == 0)
 101                 return x;
 102
 103         x -= k * gnm_ldexp (pi_over_two[0], 2 - kbits);
 104         for (i = 1; i < 5; i++) {
 105                 gnm_float dx = k * gnm_ldexp (pi_over_two[i], 2 - kbits);
 106                 gnm_float s = x - dx;
 107                 xx += x - s - dx;
 108                 x = s;
 109         }
 110
 111         return x + xx;
 112 }
 113
 114 /*
 115  * Add the 64-bit number p at *dst and backwards.  This would be very
 116  * simple and fast in assembler.  In C it's a bit of a mess.
 117  */
 118 static inline void
 119 add_at (guint32 *dst, guint64 p)
 120 {
 121         unsigned h = p >> 63;
 122
 123         p += *dst;
 124         *dst-- = p & 0xffffffffu;
 125         p >>= 32;
 126         if (p) {
 127                 p += *dst;
 128                 *dst-- = p & 0xffffffffu;
 129                 if (p >> 32)
 130                         while (++*dst == 0)
 131                                 dst--;
 132         } else if (h) {
 133                 /* p overflowed, pass carry on. */
 134                 dst--;
 135                 while (++*dst == 0)
 136                         dst--;
 137         }
 138 }
 139
 140 static gnm_float
 141 reduce_pi_full (gnm_float x, int *pk, int kbits)
 142 {
 143         /* FIXME?  This table isn't big enough for long double's range */
 144         static const guint32 one_over_two_pi[] = {
 145                 0x28be60dbu,
 146                 0x9391054au,
 147                 0x7f09d5f4u,
 148                 0x7d4d3770u,
 149                 0x36d8a566u,
 150                 0x4f10e410u,
 151                 0x7f9458eau,
 152                 0xf7aef158u,
 153                 0x6dc91b8eu,
 154                 0x909374b8u,
 155                 0x01924bbau,
 156                 0x82746487u,
 157                 0x3f877ac7u,
 158                 0x2c4a69cfu,
 159                 0xba208d7du,
 160                 0x4baed121u,
 161                 0x3a671c09u,
 162                 0xad17df90u,
 163                 0x4e64758eu,
 164                 0x60d4ce7du,
 165                 0x272117e2u,
 166                 0xef7e4a0eu,
 167                 0xc7fe25ffu,
 168                 0xf7816603u,
 169                 0xfbcbc462u,
 170                 0xd6829b47u,
 171                 0xdb4d9fb3u,
 172                 0xc9f2c26du,
 173                 0xd3d18fd9u,
 174                 0xa797fa8bu,
 175                 0x5d49eeb1u,
 176                 0xfaf97c5eu,
 177                 0xcf41ce7du,
 178                 0xe294a4bau,
 179                 0x9afed7ecu,
 180                 0x47e35742u
 181         };
 182         static const guint32 pi_over_four[] = {
 183                 0xc90fdaa2u,
 184                 0x2168c234u,
 185                 0xc4c6628bu,
 186                 0x80dc1cd1u
 187         };
 188
 189         gnm_float m;
 190         guint32 w2, w1, w0, wm1, wm2;
 191         int e, neg;
 192         unsigned di, i, j;
 193         guint32 r[6], r4[4];
 194         gnm_float rh, rl, l48, h48;
 195
 196         m = gnm_frexp (x, &e);
 197         if (e >= GNM_MANT_DIG) {
 198                 di = ((unsigned)e - GNM_MANT_DIG) / 32u;
 199                 e -= di * 32;
 200         } else
 201                 di = 0;
 202         m = gnm_ldexp (m, e - 64);
 203         w2  = (guint32)gnm_floor (m); m = gnm_ldexp (m - w2, 32);
 204         w1  = (guint32)gnm_floor (m); m = gnm_ldexp (m - w1, 32);
 205         w0  = (guint32)gnm_floor (m); m = gnm_ldexp (m - w0, 32);
 206         wm1 = (guint32)gnm_floor (m); m = gnm_ldexp (m - wm1, 32);
 207         wm2 = (guint32)gnm_floor (m);
 208
 209         /*
 210          * r[0] is an integer overflow area to be ignored.  We will not create
 211          * any carry into r[-1] because 5/(2pi) < 1.
 212          */
 213         r[0] = 0;
 214
 215         for (i = 0; i < 5; i++) {
 216                 g_assert (i + 2 + di < G_N_ELEMENTS (one_over_two_pi));
 217                 r[i + 1] = 0;
 218                 if (wm2 && i > 1)
 219                         add_at (&r[i + 1], (guint64)wm2 * one_over_two_pi[i - 2]);
 220                 if (wm1 && i > 0)
 221                         add_at (&r[i + 1], (guint64)wm1 * one_over_two_pi[i - 1]);
 222                 if (w0)
 223                         add_at (&r[i + 1], (guint64)w0  * one_over_two_pi[i     + di]);
 224                 if (w1)
 225                         add_at (&r[i + 1], (guint64)w1  * one_over_two_pi[i + 1 + di]);
 226                 if (w2)
 227                         add_at (&r[i + 1], (guint64)w2  * one_over_two_pi[i + 2 + di]);
 228
 229                 /*
 230                  * We're done at i==3 unless the first 31-kbits bits, not counting
 231                  * those ending up in sign and *pk, are all zeros or ones.
 232                  */
 233                 if (i == 3 && ((r[1] + 1u) & (0x7fffffffu >> kbits)) > 1)
 234                         break;
 235         }
 236
 237         *pk = kbits ? (r[1] >> (32 - kbits)) : 0;
 238         if ((neg = ((r[1] >> (31 - kbits)) & 1))) {
 239                 (*pk)++;
 240                 /* Two-complement negation */
 241                 for (j = 1; j <= i; j++) r[j] ^= 0xffffffffu;
 242                 add_at (&r[i], 1);
 243         }
 244         r[1] &= (0xffffffffu >> kbits);
 245
 246         j = 1;
 247         if (r[j] == 0) j++;
 248         r4[0] = r4[1] = r4[2] = r4[3] = 0;
 249         add_at (&r4[1], (guint64)r[j    ] * pi_over_four[0]);
 250         add_at (&r4[2], (guint64)r[j    ] * pi_over_four[1]);
 251         add_at (&r4[2], (guint64)r[j + 1] * pi_over_four[0]);
 252         add_at (&r4[3], (guint64)r[j    ] * pi_over_four[2]);
 253         add_at (&r4[3], (guint64)r[j + 1] * pi_over_four[1]);
 254         add_at (&r4[3], (guint64)r[j + 2] * pi_over_four[0]);
 255
 256         h48 = gnm_ldexp (((guint64)r4[0] << 16) | (r4[1] >> 16),
 257                          -32 * j + (kbits + 1) - 16);
 258         l48 = gnm_ldexp (((guint64)(r4[1] & 0xffff) << 32) | r4[2],
 259                          -32 * j + (kbits + 1) - 64);
 260
 261         rh = h48 + l48;
 262         rl = h48 - rh + l48;
 263
 264         if (neg) {
 265                 rh = -rh;
 266                 rl = -rl;
 267         }
 268
 269         return gnm_ldexp (rh + rl, 2 - kbits);
 270 }
 271
 272 /**
 273   * gnm_reduce_pi:
 274   * @x: number of reduce
 275   * @e: scale between -1 and 8, inclusive.
 276   * @k: (out): location to return lower @e+1 bits of reduction count
 277   *
 278   * This function performs range reduction for trigonometric functions.
 279   *
 280   * Returns: a value, xr, such that x = xr + j * Pi/2^@e for some integer
 281   * number j and |xr| <=  Pi/2^(@e+1).  The lower @e+1 bits of j will be
 282   * returned in @k.
 283   */
 284 gnm_float
 285 gnm_reduce_pi (gnm_float x, int e, int *k)
 286 {
 287         gnm_float xr;
 288         void *state;
 289
 290         g_return_val_if_fail (e >= -1 && e <= 8, x);
 291         g_return_val_if_fail (k != NULL, x);
 292
 293         if (!gnm_finite (x)) {
 294                 *k = 0;
 295                 return x * gnm_nan;
 296         }
 297
 298         /*
 299          * We aren't actually using quads, but we rely somewhat on
 300          * proper ieee double semantics.
 301          */
 302         state = gnm_quad_start ();
 303
 304         if (gnm_abs (x) < (1u << (27 - e)))
 305                 xr = reduce_pi_simple (gnm_abs (x), k, e + 1);
 306         else
 307                 xr = reduce_pi_full (gnm_abs (x), k, e + 1);
 308
 309         if (x < 0) {
 310                 xr = -xr;
 311                 *k = -*k;
 312         }
 313         *k &= ((1 << (e + 1)) - 1);
 314
 315         gnm_quad_end (state);
 316
 317         return xr;
 318 }
 319
 320
 321
 322
 323 #ifdef GNM_REDUCES_TRIG_RANGE
 324
 325 gnm_float
 326 gnm_sin (gnm_float x)
 327 {
 328         int km4;
 329         gnm_float xr = gnm_reduce_pi (x, 1, &km4);
 330
 331         switch (km4) {
 332         default:
 333         case 0: return +sin (xr);
 334         case 1: return +cos (xr);
 335         case 2: return -sin (xr);
 336         case 3: return -cos (xr);
 337         }
 338 }
 339
 340 gnm_float
 341 gnm_cos (gnm_float x)
 342 {
 343         int km4;
 344         gnm_float xr = gnm_reduce_pi (x, 1, &km4);
 345
 346         switch (km4) {
 347         default:
 348         case 0: return +cos (xr);
 349         case 1: return -sin (xr);
 350         case 2: return -cos (xr);
 351         case 3: return +sin (xr);
 352         }
 353 }
 354
 355 gnm_float
 356 gnm_tan (gnm_float x)
 357 {
 358         int km4;
 359         gnm_float xr = gnm_reduce_pi (x, 1, &km4);
 360
 361         switch (km4) {
 362         default:
 363         case 0: case 2: return +tan (xr);
 364         case 1: case 3: return -cos (xr) / sin (xr);
 365         }
 366 }
 367
 368 #endif
 369
 370 /* ------------------------------------------------------------------------- */