summaryrefslogtreecommitdiffstats
path: root/scilab/modules/elementary_functions/tests/unit_tests/cos.tst
blob: be5e91e9192fd0a108fad43007aa030bcf6b8c8e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
// =============================================================================
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2008 - INRIA - Pierre MARECHAL <pierre.marechal@inria.fr>
// Copyright (C) 2010 - DIGITEO - Michael Baudin
//
//  This file is distributed under the same license as the Scilab package.
// =============================================================================

// <-- JVM NOT MANDATORY -->

// unit tests for cos() function (element wise cosine)
// =============================================================================

// TODO : complex arithmetic

// 1. Interface
// ============

if execstr("cos()"   ,"errcatch")           == 0 then pause, end
if execstr("cos(1,2)","errcatch")           == 0 then pause, end
if execstr("cos(''my string'')","errcatch") == 0 then pause, end

// 2. Singular Values
// ==================

if ( cos(0) <> 1 ) then pause, end
if ( cos(-0) <> 1 ) then pause, end

// The following tests check cos for a small number of x in the range [0,pi].
// The variable x contains correctly rounded values of the exact value of x. 
// The variable v contains correctly rounded values of the exact value of cos 
// on the exact double representing the various values of x.
// For example, the floating point number which is closest to %pi/2 is 
// 7074237752028440 * 2^-52, which can be written with the decimal string 
// 1.570796326794896558D+00. 
// We have cos(7074237752028440 * 2^-52) = 6.123233995736765886...*10^-17
// exactly.
// If Scilab had support for hex, we would have used it. 
// The exact values are computed in Wolfram Alpha.
// We use more that 17 digits, which, if the decimal to binary conversion is 
// correct, and if the rounding is round-to-nearest, must exactly produce 
// values in this table.
// We avoid using values such as 2*%pi/3, which introduce one multiplication
// and one addition (the test fail is the multiplication or division 
// is not accurate, while the current test is not sensitive to this, i.e.
// we test "cos", not "*"). 
// Failing this test may be caused by:
// * a bad decimal to binary conversion,
// * a wrong implementation of cos.
x = [
  5.235987755982988157D-01 // %pi/6 
  7.853981633974482790D-01 // %pi/4
  1.047197551196597631D+00 // %pi/3
  1.570796326794896558D+00 // %pi/2
  2.094395102393195263D+00 // 2*%pi/3
  2.356194490192344837D+00 // 3*%pi/4
  2.617993877991494411D+00 // 5*%pi/6
  3.141592653589793116D+00 // %pi
];
v = [
    8.660254037844386755D-01  
    7.071067811865475727D-01  
    5.000000000000000994D-01
    6.123233995736765886D-17
   -4.999999999999998011D-01
   -7.071067811865474595D-01
   -8.660254037844386698D-01
   -1.000000000000000000D+00  
];

C = cos(x);
rtol = ceil(abs(C-v)./abs(v)/%eps);
// Our tolerance is : get the exact floating point number, 
// or its neighbour.
ulptol = 1;
if ( or(rtol>ulptol) ) then pause, end
//
// Check symetry on these points
x = -x;
C = cos(x);
rtol = ceil(abs(C-v)./abs(v)/%eps);
// Our tolerance is : get the exact floating point number, 
// or its neighbour.
ulptol = 1;
if ( or(rtol>ulptol) ) then pause, end

//
// cos(x) == 1 for |x| < sqrt(2)*2^-27
// Reference:
// "Worst Cases for Correct Rounding of the Elementary Functions in Double Precision"
// Lefevre, Muller, 2003
// Table 2: Some results for small values in double precision, assuming
// rounding to the nearest.
//
// Check positive and negative normal numbers.
x = [2^(-1022:-27) -2^(-1022:-27)];
// Check that the values are close to 1.
// Our tolerance is : get the exact floating point number, 
// or its neighbour.
C = cos(x);
rtol = ceil(abs(C-1)/%eps);
ulptol = 1;
if ( or(rtol>ulptol) ) then pause, end
// Check that the values are lower or equal to 1.
// No matter how bad our library is, we must have abs(cos(x))<= 1.
// If this test fails, the math library is to be absolutely rejected.
if ( or(abs(C)>1) ) then pause, end
// Compute the number of floats for which cos(x)<>1.
// An excellent library should produce s=0.
// This failure happens for x=2^n and n usually close to -27.
notexact = sum(C<>1);
rtol = 30;
if ( notexact>rtol ) then pause, end

// Check that abs(cos(x))<= 1, for large normal floating point numbers x.
// If this test fails, the math library is to be absolutely rejected.
x = [2^(0:1023) -2^(0:1023)];
C = cos(x);
if ( or(abs(C)>1) ) then pause, end

// 3. Not A Number
// ===============

if ~isnan(cos(%nan)) then pause, end
if ~isnan(cos(-%nan)) then pause, end


// 4. Limit values
// ===============

if ~isnan(real(cos(%inf)))    then pause, end
if imag(cos(%inf)) <> 0       then pause, end

if ~isnan(real(cos(-%inf)))   then pause, end
if imag(cos(-%inf)) <> 0      then pause, end



// 5. Properties
// =============

// All the tests below are based on equalities of the form C=0, 
// with C = f(A,B) and A, B matrices.
// We consider the elementwise absolute error abs(C), and 
// focus on the maximum of this error, that is, we compute max(abs(C)).
// This absolute error is compared to atolratio * %eps, but it is 
// expressed as max(abs(C))/%eps.
// Indeed, if this test fail, we first compute max(abs(C))/%eps, see 
// its value (e.g. 0.9, 1.2 or 2.3) and set the tolerance to the minimum
// integer larger than this.

rand("seed",0);
A = rand(100,100);
B = rand(100,100);

// The ratio between the absolute tolerance and %eps.
atolratio = 5;

// cos(-x) = cos(x)
C = cos(-A) - cos(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(%pi - x) = - cos(x)
C = cos(%pi - A) + cos(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(%pi + x) = - cos(x)
C = cos(%pi + A) + cos(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(2a) = 2 cos^2(a) - 1
C = cos(2*A) - 2 * (cos(A)).^2 + 1;
if ( max(abs(C))/%eps > atolratio) then pause, end

//
// At this point, we do not test the accuracy of cos anymore:
// we test the matching between cos and sin.
// Thus, these tests may be put into cos.tst or sin.tst
// 

// cos(%pi/2 - x) = sin(x)
C = cos(%pi/2 - A) - sin(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(%pi/2 + x) = - sin(x)
C = cos(%pi/2 + A) + sin(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos^2(a) + sin^2(a) = 1
C = (cos(A)).^2 + (sin(A)).^2 - 1;
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(a + b) = cos(a) cos(b) - sin(a) sin(b)
C = cos(A + B) - cos(A).*cos(B) + sin(A).*sin(B);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(a -b) = cos(a) cos(b) + sin(a) sin(b)
C = cos(A - B) - cos(A).*cos(B) - sin(A).*sin(B);
if ( max(abs(C))/%eps > atolratio) then pause, end

// sin(a + b) = sin(a) cos(b) + sin(b) cos(a)
C = sin(A + B) - sin(A).*cos(B) - sin(B).*cos(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// sin(a -b) = sin(a) cos(b) - sin(b) cos(a)
C = sin(A - B) - sin(A).*cos(B) + sin(B).*cos(A);
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(2a) = cos^2(a) - sin^2(a) 
C = cos(2*A) - (cos(A)).^2 + (sin(A)).^2;
if ( max(abs(C))/%eps > atolratio) then pause, end

// cos(2a) = 1 - 2 sin^2(a)
C = cos(2*A) - 1 + 2 * (sin(A)).^2;
if ( max(abs(C))/%eps > atolratio) then pause, end