Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
freebsd
GitHub Repository: freebsd/freebsd-src
 Path: blob/main/contrib/arm-optimized-routines/math/aarch64/experimental/erfinv_24u5.c
48375 views
1
/*

2
 * Double-precision inverse error function.

3
 *

4
 * Copyright (c) 2023-2024, Arm Limited.

5
 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception

6
 */

7
#include "math_config.h"

8
#include "poly_scalar_f64.h"

9
#include "test_sig.h"

10
#include "test_defs.h"

11


12
const static struct

13
{

14
  /*  We use P_N and Q_N to refer to arrays of coefficients, where P_N is the

15
      coeffs of the numerator in table N of Blair et al, and Q_N is the coeffs

16
      of the denominator.  */

17
  double P_17[7], Q_17[7], P_37[8], Q_37[8], P_57[9], Q_57[10];

18
} data = {

19
  .P_17 = { 0x1.007ce8f01b2e8p+4, -0x1.6b23cc5c6c6d7p+6, 0x1.74e5f6ceb3548p+7,

20
	    -0x1.5200bb15cc6bbp+7, 0x1.05d193233a849p+6, -0x1.148c5474ee5e1p+3,

21
	    0x1.689181bbafd0cp-3 },

22
  .Q_17 = { 0x1.d8fb0f913bd7bp+3, -0x1.6d7f25a3f1c24p+6, 0x1.a450d8e7f4cbbp+7,

23
	    -0x1.bc3480485857p+7, 0x1.ae6b0c504ee02p+6, -0x1.499dfec1a7f5fp+4,

24
	    0x1p+0 },

25
  .P_37 = { -0x1.f3596123109edp-7, 0x1.60b8fe375999ep-2, -0x1.779bb9bef7c0fp+1,

26
	    0x1.786ea384470a2p+3, -0x1.6a7c1453c85d3p+4, 0x1.31f0fc5613142p+4,

27
	    -0x1.5ea6c007d4dbbp+2, 0x1.e66f265ce9e5p-3 },

28
  .Q_37 = { -0x1.636b2dcf4edbep-7, 0x1.0b5411e2acf29p-2, -0x1.3413109467a0bp+1,

29
	    0x1.563e8136c554ap+3, -0x1.7b77aab1dcafbp+4, 0x1.8a3e174e05ddcp+4,

30
	    -0x1.4075c56404eecp+3, 0x1p+0 },

31
  .P_57 = { 0x1.b874f9516f7f1p-14, 0x1.5921f2916c1c4p-7, 0x1.145ae7d5b8fa4p-2,

32
	    0x1.29d6dcc3b2fb7p+1, 0x1.cabe2209a7985p+2, 0x1.11859f0745c4p+3,

33
	    0x1.b7ec7bc6a2ce5p+2, 0x1.d0419e0bb42aep+1, 0x1.c5aa03eef7258p-1 },

34
  .Q_57 = { 0x1.b8747e12691f1p-14, 0x1.59240d8ed1e0ap-7, 0x1.14aef2b181e2p-2,

35
	    0x1.2cd181bcea52p+1, 0x1.e6e63e0b7aa4cp+2, 0x1.65cf8da94aa3ap+3,

36
	    0x1.7e5c787b10a36p+3, 0x1.0626d68b6cea3p+3, 0x1.065c5f193abf6p+2,

37
	    0x1p+0 }

38
};

39


40
/* Inverse error function approximation, based on rational approximation as

41
   described in

42
   J. M. Blair, C. A. Edwards, and J. H. Johnson,

43
   "Rational Chebyshev approximations for the inverse of the error function",

44
   Math. Comp. 30, pp. 827--830 (1976).

45
   https://doi.org/10.1090/S0025-5718-1976-0421040-7

46
   Largest observed error is 24.46 ULP, in the extreme tail:

47
   erfinv(0x1.fd9504351b757p-1) got 0x1.ff72c1092917p+0

48
			       want 0x1.ff72c10929158p+0.  */

49
double

50
erfinv (double x)

51
{

52
  double a = fabs (x);

53


54
  if (a <= 0.75)

55
    {

56
      /* Largest observed error in this region is 6.06 ULP:

57
	 erfinv(0x1.1884650fd2d41p-2) got 0x1.fb65998cbd3fep-3

58
				     want 0x1.fb65998cbd404p-3.  */

59
      double t = x * x - 0.5625;

60
      return x * horner_6_f64 (t, data.P_17) / horner_6_f64 (t, data.Q_17);

61
    }

62


63
  if (a <= 0.9375)

64
    {

65
      /* Largest observed error in this region is 6.95 ULP:

66
	 erfinv(0x1.a8d65b94d8c6p-1) got 0x1.f08325591b54p-1

67
				    want 0x1.f08325591b547p-1.  */

68
      double t = x * x - 0.87890625;

69
      return x * horner_7_f64 (t, data.P_37) / horner_7_f64 (t, data.Q_37);

70
    }

71


72
  double t = 1.0 / (sqrt (-log (1 - a)));

73
  return horner_8_f64 (t, data.P_57)

74
	 / (copysign (t, x) * horner_9_f64 (t, data.Q_57));

75
}

76


77
#if USE_MPFR

78
# warning Not generating tests for erfinv, as MPFR has no suitable reference

79
#else

80
TEST_DISABLE_FENV (erfinv)

81
TEST_SIG (S, D, 1, erfinv, -0.99, 0.99)

82
TEST_ULP (erfinv, 24.0)

83
TEST_INTERVAL (erfinv, 0, 1, 40000)

84
TEST_INTERVAL (erfinv, -0x1p-1022, -1, 40000)

85
#endif

86


87