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
14 void
mpf_mat_gso(mpf_mat_t B,const mpf_mat_t A)15 mpf_mat_gso(mpf_mat_t B, const mpf_mat_t A)
16 {
17 slong i, j, k;
18 int flag;
19 mpf_t t, s, tmp, eps;
20 flint_bitcnt_t exp;
21
22 if (B->r != A->r || B->c != A->c)
23 {
24 flint_printf("Exception (mpf_mat_gso). Incompatible dimensions.\n");
25 flint_abort();
26 }
27
28 if (B == A)
29 {
30 mpf_mat_t T;
31 mpf_mat_init(T, A->r, A->c, B->prec);
32 mpf_mat_gso(T, A);
33 mpf_mat_swap_entrywise(B, T);
34 mpf_mat_clear(T);
35 return;
36 }
37
38 if (A->r == 0)
39 {
40 return;
41 }
42
43 mpf_init2(t, B->prec);
44 mpf_init2(s, B->prec);
45 mpf_init2(tmp, B->prec);
46 mpf_init2(eps, B->prec);
47 exp = ceil((A->prec) / 64.0) * 64;
48 flint_mpf_set_ui(eps, 1);
49 mpf_div_2exp(eps, eps, exp);
50
51 for (k = 0; k < A->c; k++)
52 {
53 for (j = 0; j < A->r; j++)
54 {
55 mpf_set(mpf_mat_entry(B, j, k), mpf_mat_entry(A, j, k));
56 }
57 flag = 1;
58 while (flag)
59 {
60 flint_mpf_set_ui(t, 0);
61 for (i = 0; i < k; i++)
62 {
63 flint_mpf_set_ui(s, 0);
64 for (j = 0; j < A->r; j++)
65 {
66 mpf_mul(tmp, mpf_mat_entry(B, j, i),
67 mpf_mat_entry(B, j, k));
68 mpf_add(s, s, tmp);
69 }
70 mpf_mul(tmp, s, s);
71 mpf_add(t, t, tmp);
72 for (j = 0; j < A->r; j++)
73 {
74 mpf_mul(tmp, s, mpf_mat_entry(B, j, i));
75 mpf_sub(mpf_mat_entry(B, j, k), mpf_mat_entry(B, j, k),
76 tmp);
77 }
78 }
79 flint_mpf_set_ui(s, 0);
80 for (j = 0; j < A->r; j++)
81 {
82 mpf_mul(tmp, mpf_mat_entry(B, j, k), mpf_mat_entry(B, j, k));
83 mpf_add(s, s, tmp);
84 }
85 mpf_add(t, t, s);
86 flag = 0;
87 if (mpf_cmp(s, t) < 0)
88 {
89 if (mpf_cmp(s, eps) < 0)
90 flint_mpf_set_ui(s, 0);
91 else
92 flag = 1;
93 }
94 }
95 mpf_sqrt(s, s);
96 if (flint_mpf_cmp_ui(s, 0) != 0)
97 mpf_ui_div(s, 1, s);
98 for (j = 0; j < A->r; j++)
99 {
100 mpf_mul(mpf_mat_entry(B, j, k), mpf_mat_entry(B, j, k), s);
101 }
102 }
103 mpf_clears(t, s, tmp, eps, NULL);
104 }
105