1function gbtest76
2%GBTEST76 test trig and other functions
3
4% SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
5% SPDX-License-Identifier: GPL-3.0-or-later
6
7fprintf ('\ngbtest76: testing trig and special functions\n') ;
8
9rng ('default') ;
10
11types = { 'single', 'double', ...
12    'single complex', 'double complex', 'int32', 'uint32' } ;
13
14for k = 1:length (types)
15    type = types {k} ;
16
17    fprintf ('%s\n', type) ;
18
19    switch (type)
20        case { 'single' }
21            A = single (rand (4)) ;
22            B = single (rand (4)) ;
23            tol = 1e-5 ;
24            G = GrB (A) ;
25            H = GrB (B) ;
26            gbtest76b (A, B, G, H, tol) ;
27        case { 'double' }
28            A = rand (4) ;
29            B = rand (4) ;
30            tol = 1e-10 ;
31            G = GrB (A) ;
32            H = GrB (B) ;
33            gbtest76b (A, B, G, H, tol) ;
34        case { 'single complex' }
35            A = single (rand (4) + 1i* rand (4)) ;
36            B = single (rand (4) + 1i* rand (4)) ;
37            tol = 1e-5 ;
38            G = GrB (A) ;
39            H = GrB (B) ;
40            gbtest76b (A, B, G, H, tol) ;
41        case { 'double complex' }
42            A = rand (4) + 1i* rand (4) ;
43            B = rand (4) + 1i* rand (4) ;
44            tol = 1e-10 ;
45            G = GrB (A) ;
46            H = GrB (B) ;
47            gbtest76b (A, B, G, H, tol) ;
48        case { 'int32' }
49            A = int32 (magic (4)) ;
50            B = int32 (rand (4) * 4) ;
51            tol = 1e-10 ;
52            G = GrB (A) ;
53            H = GrB (B) ;
54            gbtest76b (double (A), double (B), G, H, tol) ;
55        case { 'uint32' }
56            A = uint32 (magic (4)) ;
57            B = uint32 (rand (4) * 4) ;
58            tol = 1e-10 ;
59            G = GrB (A) ;
60            H = GrB (B) ;
61            gbtest76b (double (A), double (B), G, H, tol) ;
62    end
63end
64
65GrB.finalize ;
66fprintf ('\ngbtest76: all tests passed\n') ;
67end
68
69function gbtest76b (A, B, G, H, tol)
70
71    A (1,1) = 0 ;
72    G (1,1) = 0 ;
73    G = GrB.prune (G) ;
74    S = spones (G) ;
75
76    C1 = hypot (A, B) ;
77    C2 = hypot (G, H) ;
78    err = norm (C1-C2, 1) ;
79    assert (err < tol) ;
80
81    C1 = atan2 (real (A), real (B)) ;
82    C2 = atan2 (real (G), real (H)) ;
83    err = norm (C1-C2, 1) ;
84    assert (err < tol) ;
85
86    C1 = atan2 (pi, real (B)) ;
87    C2 = atan2 (pi, real (H)) ;
88    err = norm (C1-C2, 1) ;
89    assert (err < tol) ;
90
91    C1 = atan2 (real (A), 3) ;
92    C2 = atan2 (real (G), 3) ;
93    err = norm (C1-C2, 1) ;
94    assert (err < tol) ;
95
96    C1 = atan2 (pi, 3) ;
97    C2 = atan2 (GrB (pi), 3) ;
98    err = norm (C1-C2, 1) ;
99    assert (err < tol) ;
100
101    C1 = atan2 (0, 3) ;
102    C2 = atan2 (real (GrB (1,1,GrB.type(G))), 3) ;
103    err = norm (C1-C2, 1) ;
104    assert (err < tol) ;
105
106    C1 = complex (real (A), real (B)) ;
107    C2 = complex (real (G), real (H)) ;
108    err = norm (double (C1) - double (C2), 1) ;
109    assert (err < tol) ;
110
111    C1 = complex (pi, real (B)) ;
112    C2 = complex (pi, real (H)) ;
113    err = norm (C1-C2, 1) ;
114    assert (err < tol) ;
115
116    C1 = complex (real (A), 3) ;
117    C2 = complex (real (G), 3) ;
118    err = norm (C1-C2, 1) ;
119    assert (err < tol) ;
120
121    C1 = complex (real (A), 0) ;
122    C2 = complex (real (G), 0) ;
123    err = norm (double (C1) - double (C2), 1) ;
124    assert (err < tol) ;
125
126    C1 = complex (0, real (A)) ;
127    C2 = complex (0, real (G)) ;
128    err = norm (double (C1) - double (C2), 1) ;
129    assert (err < tol) ;
130
131    C1 = complex (pi, 3) ;
132    C2 = complex (GrB (pi), 3) ;
133    err = norm (C1-C2, 1) ;
134    assert (err < tol) ;
135
136    C1 = complex (0, 3) ;
137    C2 = complex (real (GrB (1,1,GrB.type(G))), 3) ;
138    err = norm (C1-C2, 1) ;
139    assert (err < tol) ;
140
141    C1 = acosh (A) ;
142    C2 = acosh (G) ;
143    err = norm (C1-C2, 1) ;
144    assert (err < tol) ;
145
146    C1 = acos (A) ;
147    C2 = acos (G) ;
148    err = norm (S .* (cos (C1) - cos (C2)), 1) ;
149    assert (err < tol) ;
150
151    C1 = acoth (A) ;
152    C2 = acoth (G) ;
153    err = norm (C1-C2, 1) ;
154    assert (err < tol) ;
155
156    C1 = acot (A) ;
157    C2 = acot (G) ;
158    err = norm (C1-C2, 1) ;
159    assert (err < tol) ;
160
161    C1 = acsch (A) ;
162    C2 = acsch (G) ;
163    err = norm (C1-C2, 1) ;
164    assert (err < tol) ;
165
166    C1 = acsc (A) ;
167    if (isreal (C1) || ~ispc)
168        % Windows casinf and casin are broken
169        C2 = acsc (G) ;
170        err = norm (csc (C1) - csc (C2),1) ;
171        assert (err < tol) ;
172    end
173
174    C1 = angle (A) ;
175    C2 = angle (G) ;
176    err = norm (C1-C2, 1) ;
177    assert (err < tol) ;
178
179    C1 = asech (A) ;
180    C2 = asech (G) ;
181    err = norm (C1-C2, 1) ;
182    assert (err < tol) ;
183
184    C1 = asec (A) ;
185    C2 = asec (G) ;
186    err = norm (sec (C1) - sec (C2),1) ;
187    assert (err < tol) ;
188
189    C1 = asinh (A) ;
190    C2 = asinh (G) ;
191    err = norm (C1-C2, 1) ;
192    assert (err < tol) ;
193
194    C1 = asin (A) ;
195    if (isreal (C1) || ~ispc)
196        % Windows casinf and casin are broken
197        C2 = asin (G) ;
198        err = norm (sin (C1) - sin (C2), 1) ;
199        assert (err < tol) ;
200    end
201
202    C1 = asin (2*A) ;
203    if (isreal (C1) || ~ispc)
204        % Windows casinf and casin are broken
205        C2 = asin (2*G) ;
206        err = norm (sin (C1) - sin (C2), 1) ;
207        assert (err < tol) ;
208    end
209
210    C1 = atanh (A) ;
211    C2 = atanh (G) ;
212    err = norm (C1-C2, 1) ;
213    assert (err < tol) ;
214
215    C1 = atan (A) ;
216    C2 = atan (G) ;
217    err = norm (C1-C2, 1) ;
218    assert (err < tol) ;
219
220    C1 = conj (A) ;
221    C2 = conj (G) ;
222    err = norm (C1-C2, 1) ;
223    assert (err < tol) ;
224
225    C1 = cosh (A) ;
226    C2 = cosh (G) ;
227    err = norm (C1-C2, 1) ;
228    assert (err < tol) ;
229
230    C1 = cos (A) ;
231    C2 = cos (G) ;
232    err = norm (C1-C2, 1) ;
233    assert (err < tol) ;
234
235    if (contains (GrB.type (G), 'double'))
236        % MATLAB cannot compute spfun for single;
237        % MATLAB and GraphBLAS can't do this for int*
238        C1 = spfun ('cos', A) ;
239        C2 = spfun ('cos', G) ;
240        err = norm (C1-C2, 1) ;
241        assert (err < tol) ;
242
243        C1 = spfun (@cos, A) ;
244        C2 = spfun (@cos, G) ;
245        err = norm (C1-C2, 1) ;
246        assert (err < tol) ;
247
248        % GraphBLAS doesn't have cosd, so it falls
249        % through to use feval('cosd',x).
250        C1 = spfun ('cosd', A) ;
251        C2 = spfun ('cosd', G) ;
252        err = norm (C1-C2, 1) ;
253        assert (err < tol) ;
254    end
255
256    C1 = coth (A) ;
257    C2 = coth (G) ;
258    err = norm (C1-C2, 1) ;
259    assert (err < tol) ;
260
261    C1 = cot (A) ;
262    C2 = cot (G) ;
263    err = norm (C1-C2, 1) ;
264    assert (err < tol) ;
265
266    C1 = csch (A) ;
267    C2 = csch (G) ;
268    err = norm (C1-C2, 1) ;
269    assert (err < tol) ;
270
271    C1 = csc (A) ;
272    C2 = csc (G) ;
273    err = norm (C1-C2, 1) ;
274    assert (err < tol) ;
275
276    C1 = erfc (real (full (A))) ;
277    C2 = erfc (real (G)) ;
278    err = norm (C1-C2, 1) ;
279    assert (err < tol) ;
280
281    C1 = erf (real (full (A))) ;
282    C2 = erf (real (G)) ;
283    err = norm (C1-C2, 1) ;
284    assert (err < tol) ;
285
286    C1 = exp (A) ;
287    C2 = exp (G) ;
288    err = norm (C1-C2, 1) ;
289    assert (err < tol) ;
290
291    C1 = expm1 (A) ;
292    C2 = expm1 (G) ;
293    err = norm (C1-C2, 1) ;
294    assert (err < tol) ;
295
296    C1 = gammaln (real (full (A))) ;
297    C2 = gammaln (real (G)) ;
298    err = norm (C1-C2, 1) ;
299    assert (err < tol) ;
300
301    C1 = gamma (real (full (A))) ;
302    C2 = gamma (real (G)) ;
303    assert (isinf (C1 (1,1)) == isinf (C2 (1,1)))
304    C1 (1,1) = 0 ;
305    C2 (1,1) = 0 ;
306    err = norm (C1-C2, 1) / norm (C1, 1) ;
307    assert (err < tol) ;
308
309    C1 = imag (A) ;
310    C2 = imag (G) ;
311    err = norm (C1-C2, 1) ;
312    assert (err < tol) ;
313
314    C1 = log10 (A) ;
315    C2 = log10 (G) ;
316    err = norm (C1-C2, 1) ;
317    assert (err < tol) ;
318
319    C1 = log1p (A) ;
320    C2 = log1p (G) ;
321    err = norm (C1-C2, 1) ;
322    assert (err < tol) ;
323
324    C1 = log2 (A) ;
325    C2 = log2 (G) ;
326    err = norm (C1-C2, 1) ;
327    assert (err < tol) ;
328
329    [F1,E1] = log2 (real (A)) ;
330    [F2,E2] = log2 (G) ;
331    C1 = F1 .* (2 .^ E1) ;
332    C2 = F2 .* (2 .^ E2) ;
333    err = norm (C1-C2, 1) ;
334    assert (err < tol) ;
335
336    C1 = pow2 (F1, E1) ;
337    C2 = pow2 (F2, E2) ;
338    err = norm (C1-C2, 1) ;
339    assert (err < tol) ;
340
341    C1 = log (A) ;
342    C2 = log (G) ;
343    err = norm (C1-C2, 1) ;
344    assert (err < tol) ;
345
346    C1 = pow2 (A) ;
347    C2 = pow2 (G) ;
348    err = norm (C1-C2, 1) ;
349    assert (err < tol) ;
350
351    C1 = sech (A) ;
352    C2 = sech (G) ;
353    err = norm (C1-C2, 1) ;
354    assert (err < tol) ;
355
356    C1 = sec (A) ;
357    C2 = sec (G) ;
358    err = norm (C1-C2, 1) ;
359    assert (err < tol) ;
360
361    C1 = sinh (A) ;
362    C2 = sinh (G) ;
363    err = norm (C1-C2, 1) ;
364    assert (err < tol) ;
365
366    C1 = sin (A) ;
367    C2 = sin (G) ;
368    err = norm (C1-C2, 1) ;
369    assert (err < tol) ;
370
371    C1 = tanh (A) ;
372    C2 = tanh (G) ;
373    err = norm (C1-C2, 1) ;
374    assert (err < tol) ;
375
376    C1 = tan (A) ;
377    C2 = tan (G) ;
378    err = norm (C1-C2, 1) ;
379    assert (err < tol) ;
380
381    C1 = A.' ;
382    C2 = G.' ;
383    err = norm (C1-C2, 1) ;
384    assert (err == 0) ;
385
386    C1 = A' ;
387    C2 = G' ;
388    err = norm (C1-C2, 1) ;
389    assert (err == 0) ;
390
391    if (isfloat (G))
392        C1 = eps (A) ;
393        C2 = eps (G) ;
394        err = norm (C1-C2, 1) ;
395        assert (err < tol) ;
396    end
397
398    C1 = A+0 ;
399    C2 = G+0 ;
400    err = norm (C1-C2, 1) ;
401    assert (err == 0) ;
402
403    C1 = 0+A ;
404    C2 = 0+G ;
405    err = norm (C1-C2, 1) ;
406    assert (err == 0) ;
407
408    I = GrB ([1 2 3]) ;
409    J = GrB ([3 1 2]) ;
410
411    C1 = A (I,J) ;
412    C2 = G (I,J) ;
413    err = norm (C1-C2, 1) ;
414    assert (err == 0) ;
415
416    C1 = A.^1 ;
417    C2 = G.^1 ;
418    err = norm (C1-C2, 1) ;
419    assert (err == 0) ;
420
421end
422
423
424