1 /*
2 Copyright (C) 2014 Abhinav Baid
3
4 This file is part of FLINT.
5
6 FLINT 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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "mpf_mat.h"
13 #include "ulong_extras.h"
14
15 int
main(void)16 main(void)
17 {
18 int i;
19 FLINT_TEST_INIT(state);
20
21 flint_printf("gso....");
22 fflush(stdout);
23
24 for (i = 0; i < 100 * flint_test_multiplier(); i++)
25 {
26 mpf_t dot, tmp;
27 int j, k, l;
28 mpf_mat_t A, B;
29 flint_bitcnt_t prec;
30
31 slong m, n;
32
33 m = n_randint(state, 10);
34 n = n_randint(state, 10);
35 prec = n_randint(state, 200) + 3;
36
37 mpf_mat_init(A, m, n, prec);
38 mpf_mat_init(B, m, n, prec);
39 mpf_init2(dot, prec);
40 mpf_init2(tmp, prec);
41
42 mpf_mat_randtest(A, state, prec);
43
44 mpf_mat_set(B, A);
45
46 mpf_mat_gso(A, A);
47
48 for (j = 0; j < n; j++)
49 {
50 mpf_t norm;
51 mpf_init2(norm, prec);
52 for (l = 0; l < m; l++)
53 {
54 mpf_mul(tmp, mpf_mat_entry(A, l, j), mpf_mat_entry(A, l, j));
55 mpf_add(norm, norm, tmp);
56 }
57
58 mpf_sub_ui(dot, norm, 1);
59 mpf_abs(dot, dot);
60 flint_mpf_set_ui(tmp, 1);
61 mpf_div_2exp(tmp, tmp, prec - 3);
62 if (flint_mpf_cmp_ui(norm, 0) != 0 && mpf_cmp(dot, tmp) > 0)
63 {
64 flint_printf("FAIL:\n");
65 flint_printf("A:\n");
66 mpf_mat_print(A);
67 mpf_out_str(stdout, 10, 0, norm);
68 flint_printf("\n");
69 mpf_out_str(stdout, 10, 0, dot);
70 flint_printf("\n");
71 mpf_out_str(stdout, 10, 0, tmp);
72 flint_printf("\n");
73 flint_printf("%d\n", prec);
74 flint_printf("%d\n", j);
75 abort();
76 }
77 mpf_clear(norm);
78 for (k = j + 1; k < n; k++)
79 {
80 flint_mpf_set_ui(dot, 0);
81 for (l = 0; l < m; l++)
82 {
83 mpf_mul(tmp, mpf_mat_entry(A, l, j),
84 mpf_mat_entry(A, l, k));
85 mpf_add(dot, dot, tmp);
86 }
87
88 flint_mpf_set_ui(tmp, 1);
89 mpf_div_2exp(tmp, tmp, prec);
90 mpf_abs(dot, dot);
91 if (mpf_cmp(dot, tmp) > 0)
92 {
93 flint_printf("FAIL:\n");
94 flint_printf("A:\n");
95 mpf_mat_print(A);
96 flint_printf("B:\n");
97 mpf_mat_print(B);
98 mpf_out_str(stdout, 10, 0, dot);
99 flint_printf("\n");
100 flint_printf("%d %d\n", j, k);
101 abort();
102 }
103 }
104 }
105
106 mpf_mat_clear(A);
107 mpf_mat_clear(B);
108 mpf_clears(dot, tmp, NULL);
109 }
110
111 FLINT_TEST_CLEANUP(state);
112
113 flint_printf("PASS\n");
114 return EXIT_SUCCESS;
115 }
116