Test Construction for Mathematical Functions

1
Test Construction for Mathematical Functions

• Victor Kuliamin
• Institute for System Programming, Russian Academy
  of Sciences,
• Moscow

2
Mathematical Modeling

• Astrophysics
• Geosciences
• Biosciences
• Social Sciences
• Confidence?

3
Mathematical Libraries

• Floating-point numbers arithmetics
• Basic math functions
• Elementary
• Special
• Specialized libraries
• Linear algebra
• Number theory
• Numeric calculus
• Dynamic systems
• Optimization

sqrt, pow, exp, log, sin, atan, cosh,
erf, tgamma, j0, y1,
4
How they should work?

• Intuition everybody knows sin (?)
• Standards
• IEEE 754 (Floating-point arithmetics)FP numbers,
  operations ,-,,/, sqrt
• ISO 9899 (C language and libraries)28 real
  functions
• IEEE 1003.1 (POSIX)34 real 16 complex
  functions
• ISO 10697.1-3 (Language independent
  arithmetics)Elementary real and complex functions

5
IEEE 754 Floating-point Numbers
k1
n-1
0
1
k
n, k
0
0
1
1
1
1
1
1
0
1
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
S
E
M
sign
exponent
mantissa
B 2(k1) 1

• Normal E gt 0 E lt 2k 1 X
  (1)S2(EB)(1M/2(nk1))
• Denormal E 0 X
  (1)S2(B1)(M/2(nk1))
• Special E 2k 1
• M 0 ?, ?
• M ? 0 NaN

2(1)1.1012 13/16
0, -0
n 32, k 8 float n 64, k 11 double n
79, k 15 ext. double n 128, k 15
quadruple
1/0 ?, (1)/0 ?
0/0 NaN
6
Standard Requirements IEEE 754

• , , , /, sqrt, remainder
• Correct rounding 4 rounding modes
• to ?
• to ?
• to 0
• to the nearest
• Exception flags
• INVALID Incorrect arguments
• DIVISION-BY-ZERO Infinite result
• OVERFLOW Too big result
• UNDERFLOW Too small (or denormal) result
• INEXACT Inexact result

0
7
Standard Requirements ISO C POSIX

• ISO/IEC 9899 (C language) 28 real functions
• Exact values sin(0) 0, exp(0) 1,
• DIVISION-BY-ZERO flag log(0), atanh(1),
  pow(0,x), ?(-n)
• NaN results and INVALID flag outside of domains
• IEEE 1003.1 (POSIX) 34 real 16 complex
• All IEEE 754 flags (except for INEXACT) for real
  functions
• Error settings Domain error INVALID, Range
  error OVERFLOW or UNDERFLOW
• If x is denormalf(x) x for each f(x)x in 0
  (sin, asin, sinh, expm1)

Contradiction with 3 rounding modes!
8
Standard Requirements ISO 10697

• Real and complex elementary functions (no erf,
  gamma, bessel)
• Only symmetric rounding modes (no rounding to ?
  or to ?)
• Preservation of sign
• Preservation of monotonicity
• Inaccuracy 0.5-2.0 ulp
• Evenness and oddity
• Exact values cosh(0) 1, log(1) 0,
• Asymptotics near 0 cos(x) 1, sin(x) x,
• Relations expm1 lt exp, cosh gt sinh, atan lt
  ?( p/2 ) ,

for sin, cos, tan small arguments only
9
Summary of Requirements

• Domain boundaries and poles ( flags)
• Exact values, limits and asymptotics
• Preservation of sign and monotonicity
• SymmetriesEvenness, periodicity, others
  ?(1x) x?(x)
• Relations and range boundaries
• Correct rounding (according to mode)
• Computational accuracy
• Interoperability and portabilityof libraries and
  applications

10
Contradictions with Correct Rounding

• Correct rounding
• Oddity (sym. with x, 1/x)

POSIX f(x) x for denormal x and f(x)x in 0
Range boundaries
11
Test Data Construction

• Bit structure of FP numbers
• Boundaries
• 0, -0, ?, -?, NaN
• Least and greatest positive and negative, normal
  and denormal
• Mantissa patterns
• FFFFFFFFFFFFF16 FFFFF0000000016
  555550000FFFF16
• Both arguments and values of function
• Intervals of uniform function behavior
• Points hard to compute correctly rounded function
  value

12
Intervals and their boundaries
max

• Poles and overflow points
• Zeroes and extremes
• Tangents and asymtotics horizontal and diagonal

0
13
Table Maker Dilemma (TMD)

• Rounding to the nearest
• f x.xxxxxxxxxx0111111111...1xxx...
• f x.xxxxxxxxxx1000000000...0xxx...
• Rounding to 0, ?, -?
• f x.xxxxxxxxxx0000000000...0xxx...
• f x.xxxxxxxxxx1111111111...1xxx...

!
?

• tan(1.11011111111111111111111111111111111111111111
  0001111122-22)
• 1.111000000000000000000000000000000000000000
  0101010001 0 178 01022-22
• sin(1.11100000000000000000000000000000000000000111
  0000100022-19)
• 1.110111111111111111111111111111111111110000
  0010111000 067 1110122-19

14
Number of Hard Points
Eval. ? 0, ?, -? ? N
54 0.5 1
53 1 1 2
52 2 4 4
51 4 6 6
50 8 10 12
49 16 19 21
48 32 32 37
47 64 70 67
46 128 142 106
45 256 280 239
44 512 547 518
43 1024 1073 996
42 2048 2103 1985
41 4096 4187 4040
40 8192 8325 8142

• Probabilistic evaluation
• Uniform independent bits distribution
• Total N 2(n-k-1) values
• N2-m have m consecutive equal bits
• Real data for sin on exponent -16

15
Hard Points Calculation Techniques
Feasible only for single precision numbers

• Exhaustive search
• Continued fractions (Kahan, 1983)
• Dyadic method (Tang, 1989 Kahan, 1994)
• Reduced search (Lefevre, 1997)
• Lattice reduction (Gonnet, 2002 Stehle, Lefevre,
  Zimmermann, 2003)
• Integer secants method (2007)

X Np X M2m 2(n k 1) lt M lt 2(n
k) ? p (2mM)/N
sqrt(N2m) M ½ 2(n-k-1) lt M, N lt 2(n-k)
? 2(m2)N (2M 1)2 j ? (2M 1)2
j (mod 2(m2))
33864178045159811206438920823311565991202393932998
38035242121518428537554064774221620930267583474709
60206804568602636298927181441186370849986972132271
59466226343020116976329729079225588927108306160340
38541342154669787134871905353772776431251615694251
273653 p/2 1.011010101100010110110010011000
1011001010000111111110 1857 01122849 tan(1.0110
1010110001011011001001100010110010100001111111112
2849) -1.110110011011101010011010011110010111010
1011000110101 10102260

• j 15
• sqrt(1.0010010101100101011001011100101011011100101
  1111101002)
• 1.0001001000001111100110011001111010011001001101
  110100 0 150 0002

16
Integer Secants Method

• Feasible near a tangent with integer linear
  coefficient
• The most hard points are near intersections of
  parallel secants with function graph

F(x) f(x) ax b c1x2 c2x3 c3x4
F(x) c1(G(x) )2, G(x) x d1x2 d2x3
G(x) y ? x H(y), H is reversed
series xm H(sqrt(m/c12z)) ? F(xm) m/2z
2z
17
Achievements

• Hard points
• float (single precision)
• sqrt, cbrt, exp, sin, cos
• double
• Some hard points with 48 additional bits can be
  found in crlibm tests http//lipforge.ens-lyon.fr
  /projects/crlibm
• All with 46 additional bits for sqrt
• Neighborhood of 0, with 40 additional bits for
  sin, asin, cos, acos, tan, atan, sinh, asinh,
  cosh, tanh, atanh, exp, exp2, expm1, log1p, j0
• Neighborhood of /-?, with 40 additional bits
  exp, atan, tanh
• extended double
• All with 53 additional bits for sqrt
• Neighborhood of 0, with 51 additional bits for
  sin, exp
• Test suites developed
• double sqrt, atan, exp, sin

18
Tested Libraries
ID Processor arch Library OS
x86_64 glibc x86_64 glibc 2.3.4 Linux
i686 glibc i686 glibc 2.5,2.7 Linux
ia64 ia64 glibc 2.3.6,2.4 Linux
ppc32 ppc32 glibc 2.3.5 Linux
ppc64 ppc64 glibc 2.7 Linux
s390 s390 glibc 2.4 Linux
i686 VC8 i686 Visual C 8 Windows XP
x86_64 VC6 x86_64 Visual C 6 Windows XP
sparc UltraSparc III Sun libc Solaris 10
19
Test Results sqrt
ID Errors everywhere INVALID flag
x86_64 glibc No comput errors
i686 glibc ulps 1, 0, 0, 0(12,7)
ia64 No comput errors
ppc32 No comput errors
ppc64 No comput errors
s390 No comput errors
i686 VC8 ulps 0, 1, 1, 1(40,42,42), errno
x86_64 VC6 ulps 0, 1, 1, 1(40,42,42), errno
sparc Not-NaN for negative, no comput errors
20
Test Results atan
ID Errors everywhere INVALID flag
x86_64 glibc ulps 1, 22, 18, 18(71,38,38, 6)
i686 glibc ulps 1, 1, 1, 1(27,20,20,20), up gtp /2
ia64 ulps 1, 1, 1, 1(26 for all), up gtp /2, non-POSIX
ppc32 ulps 1, 22, 18, 18(5,38,38, 6)
ppc64 ulps 1, 22, 18, 18(5,38,38, 6)
s390 ulps 1, 22, 18, 18(5,38,38, 6)
i686 VC8 ulps 1, 2, 2, 1(30,51,51,25)
x86_64 VC6 ulps 1, 2, 2, 2(26,51,51,50)
sparc ulps 1, 1, 1, 1(22,22,22,25), up gtp /2
21
Test Results exp
ID Errors everywhere INV, (ex ia64) - UND, OVER, -inf
x86_64 glibc up, down, to 0 huge, up -gt1, down,0 all, inf1
i686 glibc ulps 1, 1, 1, 1(9,12,15,16), inf1
ia64 ulps 1, 1, 1, 1(3,4,4,4), inf1
ppc32 up, down, to 0 huge, up -gt1, down,0 all, inf1
ppc64 up, down, to 0 huge, up all, down -0, inf1
s390 up, down, to 0 huge, up all, down -0, UO, inf1
i686 VC8 ulps 16, 17, 16, 16(71,99,40,40), inf2, 0
x86_64 VC6 ulps 1, 2, 1, 1(41,90,13,13), inf2, 0x
sparc ulps 1, 1, 1, 1(6,10,10,10), inf -gt float
exp( -16.9666900121946469 ) 1.34027189065658032e
300
exp( 706.945585650006592 ) -1.26136053650993277e
308
exp( 0.524284463190085037 )
1.11022302462515654e-16
22
Test Results sin
ID Errors everywhere INV, DOMAIN
x86_64 glibc ulps to nearest 1 (0,3), other - huge
i686 glibc errors increase to /-inf
ia64 ulps 1, 1, 1, 1(5,5,5,6)
ppc32 ulps to nearest 1 (0,3), other huge
ppc64 ulps to nearest 1 (0,3), other huge
s390 ulps to nearest 1 (0,3), other huge
i686 VC8 errors increase to /-inf, bug for negative args
x86_64 VC6 errors increase to /-inf
sparc ulps 1, 2, 2, 2(8,45,45,9)
sin( 8.50270011698847838 ) 7.99976837364664238e
22
sin( -1.76788240868979430e-31 )
0.0980171403295605760
23
Future Development

• Complete set of hard points for some function
  (with 40 additional bits in double and with
  51 bits in extended double precision)
• Complete list of POSIX functions
• Multiple variable functions

24
Conclusion

• No adequate standards for math librariesSeveral
  standards, sometimes contradictive, highly
  incomplete
• Correct rounding is suitable concept for
  standardization
• It gives interoperability and almost all nice
  properties of implementations
• It is feasible crlibm, INRIA
• Need in methods for hard points calculation
  Current methods give only partial solution

25

• Contact
• E-mail kuliamin_at_ispras.ru
• Web www.ispras.ru/kuliamin
• Thank you!
• Questions?