Changeset 10966


Ignore:
Timestamp:
Oct 3, 2008, 11:53:46 PM (11 years ago)
Author:
gb
Message:

Provide win64 versions of acosh[f],asinh[f],atanh[f], using public-domain
code from mingw32.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/lisp-kernel/windows-calls.c

    r10946 r10966  
    456456  return gettimeofday(tp, tzp);
    457457}
     458
     459#ifdef WIN_64
     460
     461/* Make sure that the lisp calls these functions, when they do something */
     462/* This code is taken from the 32-bit mingw library and is in the
     463   public domain */
     464double
     465acosh(double x)
     466{
     467  if (isnan (x))
     468    return x;
     469
     470  if (x < 1.0)
     471    {
     472      errno = EDOM;
     473      return nan("");
     474    }
     475
     476  if (x > 0x1p32)
     477    /*  Avoid overflow (and unnecessary calculation when
     478        sqrt (x * x - 1) == x). GCC optimizes by replacing
     479        the long double M_LN2 const with a fldln2 insn.  */
     480    return log (x) + 6.9314718055994530941723E-1L;
     481
     482  /* Since  x >= 1, the arg to log will always be greater than
     483     the fyl2xp1 limit (approx 0.29) so just use logl. */
     484  return log (x + sqrt((x + 1.0) * (x - 1.0)));
     485}
     486
     487float
     488acoshf(float x)
     489{
     490  if (isnan (x))
     491    return x;
     492  if (x < 1.0f)
     493    {
     494      errno = EDOM;
     495      return nan("");
     496    }
     497
     498 if (x > 0x1p32f)
     499    /*  Avoid overflow (and unnecessary calculation when
     500        sqrt (x * x - 1) == x). GCC optimizes by replacing
     501        the long double M_LN2 const with a fldln2 insn.  */
     502    return log (x) + 6.9314718055994530941723E-1L;
     503
     504  /* Since  x >= 1, the arg to log will always be greater than
     505     the fyl2xp1 limit (approx 0.29) so just use logl. */
     506  return log (x + sqrt((x + 1.0) * (x - 1.0)));
     507}
     508
     509double
     510asinh(double x)
     511{
     512  double z;
     513  if (!isfinite (x))
     514    return x;
     515  z = fabs (x);
     516
     517  /* Avoid setting FPU underflow exception flag in x * x. */
     518#if 0
     519  if ( z < 0x1p-32)
     520    return x;
     521#endif
     522
     523  /* Use log1p to avoid cancellation with small x. Put
     524     x * x in denom, so overflow is harmless.
     525     asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
     526              = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0))  */
     527
     528  z = log1p (z + z * z / (sqrt (z * z + 1.0) + 1.0));
     529
     530  return ( x > 0.0 ? z : -z);
     531}
     532
     533float
     534asinhf(float x)
     535{
     536  float z;
     537  if (!isfinite (x))
     538    return x;
     539  z = fabsf (x);
     540
     541  /* Avoid setting FPU underflow exception flag in x * x. */
     542#if 0
     543  if ( z < 0x1p-32)
     544    return x;
     545#endif
     546
     547
     548  /* Use log1p to avoid cancellation with small x. Put
     549     x * x in denom, so overflow is harmless.
     550     asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
     551              = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0))  */
     552
     553  z = log1p (z + z * z / (sqrt (z * z + 1.0) + 1.0));
     554
     555  return ( x > 0.0 ? z : -z);
     556}
     557
     558double
     559atanh(double x)
     560{
     561  double z;
     562  if (isnan (x))
     563    return x;
     564  z = fabs (x);
     565  if (z == 1.0)
     566    {
     567      errno  = ERANGE;
     568      return (x > 0 ? INFINITY : -INFINITY);
     569    }
     570  if (z > 1.0)
     571    {
     572      errno = EDOM;
     573      return nan("");
     574    }
     575  /* Rearrange formula to avoid precision loss for small x.
     576
     577  atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x))
     578           = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0)
     579           = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x))
     580           = 0.5 * log1p ((2.0 * x ) / (1.0 - x))  */
     581  z = 0.5 * log1p ((z + z) / (1.0 - z));
     582  return x >= 0 ? z : -z;
     583}
     584
     585float
     586atanhf(float x)
     587{
     588  float z;
     589  if (isnan (x))
     590    return x;
     591  z = fabsf (x);
     592  if (z == 1.0)
     593    {
     594      errno  = ERANGE;
     595      return (x > 0 ? INFINITY : -INFINITY);
     596    }
     597  if ( z > 1.0)
     598    {
     599      errno = EDOM;
     600      return nanf("");
     601    }
     602  /* Rearrange formula to avoid precision loss for small x.
     603
     604  atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x))
     605           = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0)
     606           = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x))
     607           = 0.5 * log1p ((2.0 * x ) / (1.0 - x))  */
     608  z = 0.5 * log1p ((z + z) / (1.0 - z));
     609  return x >= 0 ? z : -z;
     610}
     611
     612#endif
    458613
    459614typedef struct {
Note: See TracChangeset for help on using the changeset viewer.