1 /*
2 Copyright (C) 2012 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "arb.h"
13
14 void
arb_fma_naive(arb_t res,const arb_t x,const arb_t y,const arb_t z,slong prec)15 arb_fma_naive(arb_t res, const arb_t x, const arb_t y, const arb_t z, slong prec)
16 {
17 arb_t t;
18 arb_init(t);
19 arb_set(t, z);
20 arb_addmul(t, x, y, prec);
21 arb_set(res, t);
22 arb_clear(t);
23 }
24
main()25 int main()
26 {
27 slong iter;
28 flint_rand_t state;
29
30 flint_printf("fma....");
31 fflush(stdout);
32
33 flint_randinit(state);
34
35 for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
36 {
37 arb_t x, y, z, res1, res2;
38 slong prec;
39 int aliasing;
40
41 arb_init(x);
42 arb_init(y);
43 arb_init(z);
44 arb_init(res1);
45 arb_init(res2);
46
47 prec = 2 + n_randint(state, 200);
48
49 arb_randtest_special(x, state, 200, 100);
50 arb_randtest_special(y, state, 200, 100);
51 arb_randtest_special(z, state, 200, 100);
52 arb_randtest_special(res1, state, 200, 100);
53 arb_randtest_special(res2, state, 200, 100);
54
55 if (n_randint(state, 10) == 0 &&
56 fmpz_bits(ARF_EXPREF(arb_midref(x))) < 10 &&
57 fmpz_bits(ARF_EXPREF(arb_midref(y))) < 10 &&
58 fmpz_bits(ARF_EXPREF(arb_midref(z))) < 10)
59 {
60 prec = ARF_PREC_EXACT;
61 }
62
63 aliasing = n_randint(state, 7);
64
65 switch (aliasing)
66 {
67 case 0:
68 arb_fma(res1, x, y, z, prec);
69 arb_fma_naive(res2, x, y, z, prec);
70 break;
71 case 1:
72 arb_set(res1, z);
73 arb_fma(res1, x, y, res1, prec);
74 arb_fma_naive(res2, x, y, z, prec);
75 break;
76 case 2:
77 arb_set(res1, x);
78 arb_fma(res1, res1, y, z, prec);
79 arb_fma_naive(res2, x, y, z, prec);
80 break;
81 case 3:
82 arb_set(res1, y);
83 arb_fma(res1, x, res1, z, prec);
84 arb_fma_naive(res2, x, y, z, prec);
85 break;
86 case 4:
87 arb_fma(res1, x, x, z, prec);
88 arb_fma_naive(res2, x, x, z, prec);
89 break;
90 case 5:
91 arb_set(res1, x);
92 arb_fma(res1, res1, res1, z, prec);
93 arb_fma_naive(res2, x, x, z, prec);
94 break;
95 default:
96 arb_set(res1, x);
97 arb_fma(res1, res1, res1, res1, prec);
98 arb_fma_naive(res2, x, x, x, prec);
99 break;
100 }
101
102 if (!arb_equal(res1, res2))
103 {
104 flint_printf("FAIL!\n");
105 flint_printf("prec = %wd, aliasing = %d\n\n", prec, aliasing);
106 flint_printf("x = "); arb_printd(x, 30); flint_printf("\n\n");
107 flint_printf("y = "); arb_printd(y, 30); flint_printf("\n\n");
108 flint_printf("z = "); arb_printd(z, 30); flint_printf("\n\n");
109 flint_printf("res1 = "); arb_printd(res1, 30); flint_printf("\n\n");
110 flint_printf("res2 = "); arb_printd(res2, 30); flint_printf("\n\n");
111 flint_abort();
112 }
113
114 arb_clear(x);
115 arb_clear(y);
116 arb_clear(z);
117 arb_clear(res1);
118 arb_clear(res2);
119 }
120
121 flint_randclear(state);
122 flint_cleanup();
123 flint_printf("PASS\n");
124 return EXIT_SUCCESS;
125 }
126