1 /*
2 Copyright (C) 2009, 2010 William Hart
3 Copyright (C) 2009, 2010 Andy Novocin
4 Copyright (C) 2014 Abhinav Baid
5
6 This file is part of FLINT.
7
8 FLINT is free software: you can redistribute it and/or modify it under
9 the terms of the GNU Lesser General Public License (LGPL) as published
10 by the Free Software Foundation; either version 2.1 of the License, or
11 (at your option) any later version. See <http://www.gnu.org/licenses/>.
12 */
13
14 #include "fmpz_lll.h"
15
16 void
fmpz_lll_storjohann_ulll(fmpz_mat_t FM,slong new_size,const fmpz_lll_t fl)17 fmpz_lll_storjohann_ulll(fmpz_mat_t FM, slong new_size, const fmpz_lll_t fl)
18 {
19 if (fl->rt == Z_BASIS)
20 {
21 slong r, c, mbits, prev_mbits, i, j;
22 int full_prec = 1, done = 0, is_U_I;
23 fmpz_mat_t U, big_td, trunc_data;
24 mpq_t deltax, etax;
25 fmpq_t delta, eta;
26
27 r = FM->r;
28 c = FM->c;
29 mbits = FLINT_ABS(fmpz_mat_max_bits(FM));
30 prev_mbits = mbits;
31
32 fmpz_mat_init(big_td, r, c + r);
33 fmpz_mat_init(trunc_data, r, c);
34
35 mpq_init(deltax);
36 mpq_init(etax);
37
38 fmpq_init(delta);
39 fmpq_init(eta);
40
41 mpq_set_d(deltax, fl->delta);
42 mpq_set_d(etax, fl->eta);
43 fmpq_set_mpq(delta, deltax);
44 fmpq_set_mpq(eta, etax);
45 mpq_clears(deltax, etax, NULL);
46
47 if (mbits > new_size)
48 {
49 full_prec = 0;
50
51 /* do some truncating */
52 fmpz_mat_scalar_tdiv_q_2exp(trunc_data, FM,
53 (ulong) (mbits - new_size));
54
55 /* make a large lattice which has identity in one corner and trunc_data in the other */
56 for (i = 0; i < r; i++)
57 {
58 fmpz_one(fmpz_mat_entry(big_td, i, i));
59
60 for (j = r; j < r + c; j++)
61 fmpz_set(fmpz_mat_entry(big_td, i, j),
62 fmpz_mat_entry(trunc_data, i, j - r));
63 }
64 }
65 else
66 {
67 full_prec = 1;
68 }
69
70 while (done == 0)
71 {
72 if (full_prec == 0)
73 {
74 fmpz_mat_lll_storjohann(big_td, delta, eta);
75 }
76 else
77 {
78 fmpz_mat_lll_storjohann(FM, delta, eta);
79 }
80
81 if (full_prec == 1)
82 done = 1;
83 else
84 {
85 /* get U and compare it to the identity */
86 fmpz_mat_window_init(U, big_td, 0, 0, r, r);
87 is_U_I = fmpz_mat_is_one(U);
88
89 if (is_U_I == 0)
90 {
91 fmpz_mat_mul(FM, U, FM);
92 }
93
94 mbits = FLINT_ABS(fmpz_mat_max_bits(FM));
95 /* make this condition better? */
96 if ((mbits - new_size > 0)
97 && (mbits <= prev_mbits - (slong) (new_size / 4))
98 && is_U_I == 0)
99 {
100 /* do some truncating */
101 fmpz_mat_scalar_tdiv_q_2exp(trunc_data, FM,
102 (ulong) (mbits - new_size));
103
104 /* keep with the big_td concept */
105 for (i = 0; i < r; i++)
106 {
107 for (j = 0; j < i; j++)
108 fmpz_zero(fmpz_mat_entry(big_td, i, j));
109 fmpz_one(fmpz_mat_entry(big_td, i, i));
110
111 for (j = i + 1; j < r; j++)
112 fmpz_zero(fmpz_mat_entry(big_td, i, j));
113
114 for (j = r; j < r + c; j++)
115 fmpz_set(fmpz_mat_entry
116 (big_td, i, j),
117 fmpz_mat_entry(trunc_data, i, j - r));
118 }
119 }
120 else
121 {
122 /* can switch to FM, no need for a new identity */
123 full_prec = 1;
124 }
125
126 prev_mbits = mbits;
127 fmpz_mat_window_clear(U);
128 }
129 }
130
131 fmpz_mat_clear(trunc_data);
132 fmpz_mat_clear(big_td);
133 fmpq_clear(delta);
134 fmpq_clear(eta);
135 }
136 else
137 {
138 fmpz_lll_with_removal_ulll(FM, NULL, new_size, NULL, fl);
139 }
140 }
141