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 <http://www.gnu.org/licenses/>.
10 */
11
12 #include "fmpz_lll.h"
13 #if HAVE_FENV
14 #include <fenv.h>
15 #endif
16
17 int
fmpz_lll_is_reduced_d(const fmpz_mat_t B,const fmpz_lll_t fl)18 fmpz_lll_is_reduced_d(const fmpz_mat_t B, const fmpz_lll_t fl)
19 {
20 #if HAVE_FENV
21 if (fl->rt == Z_BASIS)
22 {
23 /* NOTE: this algorithm should *not* be changed */
24 slong i, j, k, m, n;
25 d_mat_t A, Q, R, V, Wu, Wd, bound, bound2, bound3, boundt, mm, rm, mn,
26 rn, absR;
27 double *du, *dd;
28 double s, norm = 0, ti, tj;
29 int rounding_direction = fegetround();
30
31 if (B->r == 0 || B->r == 1)
32 return 1;
33
34 m = B->c;
35 n = B->r;
36
37 d_mat_init(A, m, n);
38 d_mat_init(Q, m, n);
39 d_mat_init(R, n, n);
40 d_mat_init(V, n, n);
41 d_mat_zero(R);
42 d_mat_zero(V);
43
44 if (fmpz_mat_get_d_mat_transpose(A, B) == -1)
45 {
46 d_mat_clear(A);
47 d_mat_clear(Q);
48 d_mat_clear(R);
49 d_mat_clear(V);
50 return 0;
51 }
52
53 for (k = 0; k < n; k++)
54 {
55 for (j = 0; j < m; j++)
56 {
57 d_mat_entry(Q, j, k) = d_mat_entry(A, j, k);
58 }
59 for (i = 0; i < k; i++)
60 {
61 s = 0;
62 for (j = 0; j < m; j++)
63 {
64 s += d_mat_entry(Q, j, i) * d_mat_entry(Q, j, k);
65 }
66 d_mat_entry(R, i, k) = s;
67 for (j = 0; j < m; j++)
68 {
69 d_mat_entry(Q, j, k) -= s * d_mat_entry(Q, j, i);
70 }
71 }
72 s = 0;
73 for (j = 0; j < m; j++)
74 {
75 s += d_mat_entry(Q, j, k) * d_mat_entry(Q, j, k);
76 }
77 d_mat_entry(R, k, k) = s = sqrt(s);
78 if (s != 0)
79 {
80 s = 1 / s;
81 for (j = 0; j < m; j++)
82 {
83 d_mat_entry(Q, j, k) *= s;
84 }
85 }
86 }
87 d_mat_clear(Q);
88
89 for (j = n - 1; j >= 0; j--)
90 {
91 d_mat_entry(V, j, j) = 1.0 / d_mat_entry(R, j, j);
92 for (i = j + 1; i < n; i++)
93 {
94 for (k = j + 1; k < n; k++)
95 {
96 d_mat_entry(V, j, i) +=
97 d_mat_entry(V, k, i) * d_mat_entry(R, j, k);
98 }
99 d_mat_entry(V, j, i) *= -d_mat_entry(V, j, j);
100 }
101 }
102
103 d_mat_init(Wu, n, n);
104 d_mat_init(Wd, n, n);
105 du = _d_vec_init(n);
106 dd = _d_vec_init(n);
107
108 fesetround(FE_DOWNWARD);
109 d_mat_mul_classical(Wd, R, V);
110 for (i = 0; i < n; i++)
111 {
112 dd[i] = d_mat_entry(Wd, i, i) - 1;
113 }
114 fesetround(FE_UPWARD);
115 d_mat_mul_classical(Wu, R, V);
116 for (i = 0; i < n; i++)
117 {
118 du[i] = d_mat_entry(Wu, i, i) - 1;
119 }
120 for (i = 0; i < n; i++)
121 {
122 s = 0;
123 for (j = 0; j < n; j++)
124 {
125 if (i != j)
126 s += FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
127 fabs(d_mat_entry(Wu, i, j)));
128 else
129 s += FLINT_MAX(fabs(dd[i]), fabs(du[i]));
130 }
131 norm = FLINT_MAX(norm, s);
132 }
133 if (norm >= 1)
134 {
135 d_mat_clear(A);
136 d_mat_clear(R);
137 d_mat_clear(V);
138 d_mat_clear(Wu);
139 d_mat_clear(Wd);
140 _d_vec_clear(du);
141 _d_vec_clear(dd);
142 fesetround(rounding_direction);
143 return 0;
144 }
145
146 d_mat_init(bound, n, n);
147
148 fesetround(FE_DOWNWARD);
149 for (i = 0; i < n; i++)
150 {
151 dd[i] = d_mat_entry(Wd, i, i) - 2;
152 }
153 fesetround(FE_UPWARD);
154 for (i = 0; i < n; i++)
155 {
156 du[i] = d_mat_entry(Wu, i, i) - 2;
157 }
158 for (i = 0; i < n; i++)
159 {
160 for (j = 0; j < n; j++)
161 {
162 if (j > i)
163 {
164 d_mat_entry(bound, i, j) =
165 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
166 fabs(d_mat_entry(Wu, i, j))) +
167 norm * norm / (1.0 - norm);
168 }
169 else if (j < i)
170 {
171 d_mat_entry(bound, i, j) =
172 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
173 fabs(d_mat_entry(Wu, i, j)));
174 }
175 else
176 {
177 d_mat_entry(bound, i, j) =
178 FLINT_MAX(fabs(dd[i]),
179 fabs(du[i])) + norm * norm / (1.0 - norm);
180 }
181 }
182 }
183 _d_vec_clear(dd);
184 _d_vec_clear(du);
185
186 d_mat_init(mm, n, n);
187 d_mat_init(rm, n, n);
188 d_mat_init(mn, n, n);
189 d_mat_init(rn, n, n);
190 d_mat_init(bound2, n, n);
191
192 for (i = 0; i < n; i++)
193 {
194 for (j = 0; j < n; j++)
195 {
196 d_mat_entry(mm, j, i) =
197 (d_mat_entry(Wu, i, j) + d_mat_entry(Wd, i, j)) / 2;
198 d_mat_entry(rm, j, i) =
199 d_mat_entry(mm, j, i) - d_mat_entry(Wd, i, j);
200 d_mat_entry(mn, i, j) =
201 (d_mat_entry(Wu, i, j) + d_mat_entry(Wd, i, j)) / 2;
202 d_mat_entry(rn, i, j) =
203 d_mat_entry(mn, i, j) - d_mat_entry(Wd, i, j);
204 }
205 }
206 fesetround(FE_DOWNWARD);
207 d_mat_mul_classical(Wd, mm, mn);
208 for (i = 0; i < n; i++)
209 {
210 d_mat_entry(Wd, i, i) -= 1;
211 }
212 fesetround(FE_UPWARD);
213 d_mat_mul_classical(Wu, mm, mn);
214 for (i = 0; i < n; i++)
215 {
216 d_mat_entry(Wu, i, i) -= 1;
217 for (j = 0; j < n; j++)
218 {
219 d_mat_entry(Wu, i, j) =
220 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
221 fabs(d_mat_entry(Wu, i, j)));
222 d_mat_entry(mm, i, j) = fabs(d_mat_entry(mm, i, j));
223 d_mat_entry(mn, i, j) = fabs(d_mat_entry(mn, i, j));
224 }
225 }
226 for (i = 0; i < n; i++)
227 {
228 for (j = 0; j < n; j++)
229 {
230 d_mat_entry(bound2, i, j) =
231 d_mat_entry(mn, i, j) + d_mat_entry(rn, i, j);
232 }
233 }
234 d_mat_mul_classical(bound2, rm, bound2);
235 for (i = 0; i < n; i++)
236 {
237 for (j = 0; j < n; j++)
238 {
239 d_mat_entry(bound2, i, j) += d_mat_entry(Wu, i, j);
240 }
241 }
242 d_mat_mul_classical(Wu, mm, rn);
243 for (i = 0; i < n; i++)
244 {
245 for (j = 0; j < n; j++)
246 {
247 d_mat_entry(bound2, i, j) += d_mat_entry(Wu, i, j);
248 }
249 }
250
251 d_mat_clear(Wu);
252 d_mat_clear(Wd);
253 d_mat_clear(mm);
254 d_mat_clear(mn);
255 d_mat_clear(rm);
256 d_mat_clear(rn);
257
258 d_mat_init(Wu, m, n);
259 d_mat_init(Wd, m, n);
260 d_mat_init(mm, n, m);
261 d_mat_init(mn, m, n);
262 d_mat_init(rm, n, m);
263 d_mat_init(rn, m, n);
264
265 fesetround(FE_DOWNWARD);
266 d_mat_mul_classical(Wd, A, V);
267 fesetround(FE_UPWARD);
268 d_mat_mul_classical(Wu, A, V);
269
270 d_mat_clear(A);
271 d_mat_clear(V);
272
273 d_mat_init(bound3, n, n);
274
275 for (i = 0; i < m; i++)
276 {
277 for (j = 0; j < n; j++)
278 {
279 d_mat_entry(mm, j, i) =
280 (d_mat_entry(Wu, i, j) + d_mat_entry(Wd, i, j)) / 2;
281 d_mat_entry(rm, j, i) =
282 d_mat_entry(mm, j, i) - d_mat_entry(Wd, i, j);
283 d_mat_entry(mn, i, j) =
284 (d_mat_entry(Wu, i, j) + d_mat_entry(Wd, i, j)) / 2;
285 d_mat_entry(rn, i, j) =
286 d_mat_entry(mn, i, j) - d_mat_entry(Wd, i, j);
287 }
288 }
289
290 d_mat_clear(Wd);
291 d_mat_clear(Wu);
292
293 d_mat_init(Wd, n, n);
294 d_mat_init(Wu, n, n);
295
296 fesetround(FE_DOWNWARD);
297 d_mat_mul_classical(Wd, mm, mn);
298 for (i = 0; i < n; i++)
299 {
300 d_mat_entry(Wd, i, i) -= 1;
301 }
302 fesetround(FE_UPWARD);
303 d_mat_mul_classical(Wu, mm, mn);
304 for (i = 0; i < n; i++)
305 {
306 d_mat_entry(Wu, i, i) -= 1;
307 for (j = 0; j < m; j++)
308 {
309 if (j < n)
310 {
311 d_mat_entry(Wu, i, j) =
312 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
313 fabs(d_mat_entry(Wu, i, j)));
314 }
315 d_mat_entry(mm, i, j) = fabs(d_mat_entry(mm, i, j));
316 }
317 }
318
319 d_mat_clear(Wd);
320 d_mat_init(Wd, m, n);
321
322 for (i = 0; i < m; i++)
323 {
324 for (j = 0; j < n; j++)
325 {
326 d_mat_entry(Wd, i, j) =
327 fabs(d_mat_entry(mn, i, j)) + d_mat_entry(rn, i, j);
328 }
329 }
330 d_mat_mul_classical(bound3, rm, Wd);
331 for (i = 0; i < n; i++)
332 {
333 for (j = 0; j < n; j++)
334 {
335 d_mat_entry(bound3, i, j) += d_mat_entry(Wu, i, j);
336 }
337 }
338 d_mat_mul_classical(Wu, mm, rn);
339 for (i = 0; i < n; i++)
340 {
341 for (j = 0; j < n; j++)
342 {
343 d_mat_entry(bound3, i, j) += d_mat_entry(Wu, i, j);
344 }
345 }
346
347 d_mat_clear(Wu);
348 d_mat_clear(Wd);
349 d_mat_clear(mm);
350 d_mat_clear(mn);
351 d_mat_clear(rm);
352 d_mat_clear(rn);
353
354 d_mat_init(boundt, n, n);
355
356 d_mat_transpose(boundt, bound);
357 for (i = 0; i < n; i++)
358 {
359 for (j = 0; j < n; j++)
360 {
361 d_mat_entry(bound2, i, j) =
362 fabs(d_mat_entry(bound2, i, j)) +
363 fabs(d_mat_entry(bound3, i, j));
364 }
365 }
366 d_mat_mul_classical(bound, bound2, bound);
367 d_mat_mul_classical(bound, boundt, bound);
368
369 d_mat_clear(bound2);
370 d_mat_clear(bound3);
371 d_mat_clear(boundt);
372
373 norm = 0;
374 for (i = 0; i < n; i++)
375 {
376 s = 0;
377 for (j = 0; j < n; j++)
378 {
379 s += fabs(d_mat_entry(bound, i, j));
380 }
381 norm = FLINT_MAX(norm, s);
382 }
383 if (norm >= 1)
384 {
385 d_mat_clear(R);
386 d_mat_clear(bound);
387 fesetround(rounding_direction);
388 return 0;
389 }
390
391 d_mat_init(absR, n, n);
392 for (i = 0; i < n; i++)
393 {
394 for (j = 0; j < n; j++)
395 {
396 if (j >= i)
397 {
398 d_mat_entry(bound, i, j) += norm * norm / (1.0 - norm);
399 }
400 else
401 {
402 d_mat_entry(bound, i, j) = 0;
403 }
404 d_mat_entry(absR, i, j) = fabs(d_mat_entry(R, i, j));
405 }
406 }
407 d_mat_mul_classical(bound, bound, absR);
408
409 d_mat_clear(absR);
410
411 for (i = 0; i < n - 1; i++)
412 {
413 fesetround(FE_DOWNWARD);
414 ti = (d_mat_entry(R, i, i) - d_mat_entry(bound, i, i)) * fl->eta;
415 fesetround(FE_UPWARD);
416 for (j = i + 1; j < n; j++)
417 {
418 tj = fabs(d_mat_entry(R, i, j)) + d_mat_entry(bound, i, j);
419 if (tj > ti)
420 {
421 d_mat_clear(R);
422 d_mat_clear(bound);
423 fesetround(rounding_direction);
424 return 0;
425 }
426 }
427 ti = d_mat_entry(R, i, i) + d_mat_entry(bound, i, i);
428 fesetround(FE_DOWNWARD);
429 tj = d_mat_entry(R, i + 1, i + 1) - d_mat_entry(bound, i + 1,
430 i + 1);
431 s = ((fabs(d_mat_entry(R, i, i + 1)) -
432 d_mat_entry(bound, i,
433 i + 1)) / ti) * ((fabs(d_mat_entry(R, i,
434 i + 1)) -
435 d_mat_entry(bound, i,
436 i + 1)) / ti) -
437 fl->delta;
438 s = -s;
439 fesetround(FE_UPWARD);
440 s = sqrt(s) * ti;
441 if (s > tj)
442 {
443 d_mat_clear(R);
444 d_mat_clear(bound);
445 fesetround(rounding_direction);
446 return 0;
447 }
448 }
449
450 d_mat_clear(R);
451 d_mat_clear(bound);
452 fesetround(rounding_direction);
453 }
454 else
455 {
456 /* NOTE: this algorithm should *not* be changed */
457 slong i, j, k, m, n;
458 d_mat_t A, R, V, Wu, Wd, bound, bound2, bound3, boundt, mm, rm, mn,
459 rn, absR;
460 double *du, *dd;
461 double s, norm = 0, ti, tj;
462 int rounding_direction = fegetround();
463
464 if (B->r == 0 || B->r == 1)
465 return 1;
466
467 m = B->c;
468 n = B->r;
469
470 d_mat_init(A, m, n);
471 d_mat_init(R, n, n);
472 d_mat_init(V, n, n);
473 d_mat_zero(R);
474 d_mat_zero(V);
475
476 if (fmpz_mat_get_d_mat_transpose(A, B) == -1)
477 {
478 d_mat_clear(A);
479 d_mat_clear(R);
480 d_mat_clear(V);
481 return 0;
482 }
483
484 for (j = 0; j < n; j++)
485 {
486 d_mat_entry(R, j, j) = d_mat_entry(A, j, j);
487 for (i = 0; i <= j - 1; i++)
488 {
489 d_mat_entry(R, i, j) = d_mat_entry(A, j, i);
490 for (k = 0; k <= i - 1; k++)
491 {
492 d_mat_entry(R, i, j) -=
493 d_mat_entry(R, k, i) * d_mat_entry(R, k, j);
494 }
495 if (d_mat_entry(R, i, i) != 0)
496 {
497 d_mat_entry(R, i, j) /= d_mat_entry(R, i, i);
498 d_mat_entry(R, j, j) -=
499 d_mat_entry(R, i, j) * d_mat_entry(R, i, j);
500 }
501 }
502 d_mat_entry(R, j, j) = sqrt(d_mat_entry(R, j, j));
503 }
504
505 for (j = n - 1; j >= 0; j--)
506 {
507 d_mat_entry(V, j, j) = 1.0 / d_mat_entry(R, j, j);
508 for (i = j + 1; i < n; i++)
509 {
510 for (k = j + 1; k < n; k++)
511 {
512 d_mat_entry(V, j, i) +=
513 d_mat_entry(V, k, i) * d_mat_entry(R, j, k);
514 }
515 d_mat_entry(V, j, i) *= -d_mat_entry(V, j, j);
516 }
517 }
518
519 d_mat_init(Wu, n, n);
520 d_mat_init(Wd, n, n);
521 du = _d_vec_init(n);
522 dd = _d_vec_init(n);
523
524 fesetround(FE_DOWNWARD);
525 d_mat_mul_classical(Wd, R, V);
526 for (i = 0; i < n; i++)
527 {
528 dd[i] = d_mat_entry(Wd, i, i) - 1;
529 }
530 fesetround(FE_UPWARD);
531 d_mat_mul_classical(Wu, R, V);
532 for (i = 0; i < n; i++)
533 {
534 du[i] = d_mat_entry(Wu, i, i) - 1;
535 }
536 for (i = 0; i < n; i++)
537 {
538 s = 0;
539 for (j = 0; j < n; j++)
540 {
541 if (i != j)
542 s += FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
543 fabs(d_mat_entry(Wu, i, j)));
544 else
545 s += FLINT_MAX(fabs(dd[i]), fabs(du[i]));
546 }
547 norm = FLINT_MAX(norm, s);
548 }
549 if (norm >= 1)
550 {
551 d_mat_clear(A);
552 d_mat_clear(R);
553 d_mat_clear(V);
554 d_mat_clear(Wu);
555 d_mat_clear(Wd);
556 _d_vec_clear(du);
557 _d_vec_clear(dd);
558 fesetround(rounding_direction);
559 return 0;
560 }
561
562 d_mat_init(bound, n, n);
563
564 fesetround(FE_DOWNWARD);
565 for (i = 0; i < n; i++)
566 {
567 dd[i] = d_mat_entry(Wd, i, i) - 2;
568 }
569 fesetround(FE_UPWARD);
570 for (i = 0; i < n; i++)
571 {
572 du[i] = d_mat_entry(Wu, i, i) - 2;
573 }
574 for (i = 0; i < n; i++)
575 {
576 for (j = 0; j < n; j++)
577 {
578 if (j > i)
579 {
580 d_mat_entry(bound, i, j) =
581 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
582 fabs(d_mat_entry(Wu, i, j))) +
583 norm * norm / (1.0 - norm);
584 }
585 else if (j < i)
586 {
587 d_mat_entry(bound, i, j) =
588 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
589 fabs(d_mat_entry(Wu, i, j)));
590 }
591 else
592 {
593 d_mat_entry(bound, i, j) =
594 FLINT_MAX(fabs(dd[i]),
595 fabs(du[i])) + norm * norm / (1.0 - norm);
596 }
597 }
598 }
599 _d_vec_clear(dd);
600 _d_vec_clear(du);
601
602 d_mat_init(mm, n, n);
603 d_mat_init(rm, n, n);
604 d_mat_init(mn, n, n);
605 d_mat_init(rn, n, n);
606 d_mat_init(bound2, n, n);
607
608 for (i = 0; i < n; i++)
609 {
610 for (j = 0; j < n; j++)
611 {
612 d_mat_entry(mm, j, i) =
613 (d_mat_entry(Wu, i, j) + d_mat_entry(Wd, i, j)) / 2;
614 d_mat_entry(rm, j, i) =
615 d_mat_entry(mm, j, i) - d_mat_entry(Wd, i, j);
616 d_mat_entry(mn, i, j) =
617 (d_mat_entry(Wu, i, j) + d_mat_entry(Wd, i, j)) / 2;
618 d_mat_entry(rn, i, j) =
619 d_mat_entry(mn, i, j) - d_mat_entry(Wd, i, j);
620 }
621 }
622 fesetround(FE_DOWNWARD);
623 d_mat_mul_classical(Wd, mm, mn);
624 for (i = 0; i < n; i++)
625 {
626 d_mat_entry(Wd, i, i) -= 1;
627 }
628 fesetround(FE_UPWARD);
629 d_mat_mul_classical(Wu, mm, mn);
630 for (i = 0; i < n; i++)
631 {
632 d_mat_entry(Wu, i, i) -= 1;
633 for (j = 0; j < n; j++)
634 {
635 d_mat_entry(Wu, i, j) =
636 FLINT_MAX(fabs(d_mat_entry(Wd, i, j)),
637 fabs(d_mat_entry(Wu, i, j)));
638 d_mat_entry(mm, i, j) = fabs(d_mat_entry(mm, i, j));
639 d_mat_entry(mn, i, j) = fabs(d_mat_entry(mn, i, j));
640 }
641 }
642 for (i = 0; i < n; i++)
643 {
644 for (j = 0; j < n; j++)
645 {
646 d_mat_entry(bound2, i, j) =
647 d_mat_entry(mn, i, j) + d_mat_entry(rn, i, j);
648 }
649 }
650 d_mat_mul_classical(bound2, rm, bound2);
651 for (i = 0; i < n; i++)
652 {
653 for (j = 0; j < n; j++)
654 {
655 d_mat_entry(bound2, i, j) += d_mat_entry(Wu, i, j);
656 }
657 }
658 d_mat_mul_classical(Wu, mm, rn);
659 for (i = 0; i < n; i++)
660 {
661 for (j = 0; j < n; j++)
662 {
663 d_mat_entry(bound2, i, j) += d_mat_entry(Wu, i, j);
664 }
665 }
666
667 d_mat_clear(Wu);
668 d_mat_clear(Wd);
669 d_mat_clear(mm);
670 d_mat_clear(mn);
671 d_mat_clear(rm);
672 d_mat_clear(rn);
673
674 d_mat_init(Wu, m, n);
675 d_mat_init(Wd, m, n);
676 d_mat_init(mm, n, m);
677 d_mat_init(mn, m, n);
678 d_mat_init(rm, n, m);
679 d_mat_init(rn, m, n);
680
681 d_mat_transpose(mm, V);
682 fesetround(FE_DOWNWARD);
683 d_mat_mul_classical(Wd, mm, A);
684 fesetround(FE_UPWARD);
685 d_mat_mul_classical(Wu, mm, A);
686
687 d_mat_clear(A);
688
689 d_mat_init(bound3, n, n);
690
691 fesetround(FE_DOWNWARD);
692 d_mat_mul_classical(mm, Wd, V);
693 for (i = 0; i < n; i++)
694 {
695 d_mat_entry(mm, i, i) -= 1;
696 }
697 fesetround(FE_UPWARD);
698 d_mat_mul_classical(rm, Wd, V);
699 for (i = 0; i < n; i++)
700 {
701 d_mat_entry(rm, i, i) -= 1;
702 }
703
704 fesetround(FE_DOWNWARD);
705 d_mat_mul_classical(mn, Wu, V);
706 for (i = 0; i < n; i++)
707 {
708 d_mat_entry(mn, i, i) -= 1;
709 }
710 fesetround(FE_UPWARD);
711 d_mat_mul_classical(rn, Wu, V);
712 for (i = 0; i < n; i++)
713 {
714 d_mat_entry(rn, i, i) -= 1;
715 }
716
717 d_mat_clear(Wd);
718 d_mat_clear(Wu);
719 d_mat_clear(V);
720
721 for (i = 0; i < n; i++)
722 {
723 for (j = 0; j < n; j++)
724 {
725 d_mat_entry(bound3, i, j) =
726 FLINT_MAX(fabs(d_mat_entry(mm, i, j)),
727 fabs(d_mat_entry(mn, i, j)));
728 d_mat_entry(bound3, i, j) =
729 FLINT_MAX(fabs(d_mat_entry(bound3, i, j)),
730 fabs(d_mat_entry(rm, i, j)));
731 d_mat_entry(bound3, i, j) =
732 FLINT_MAX(fabs(d_mat_entry(bound3, i, j)),
733 fabs(d_mat_entry(rn, i, j)));
734 }
735 }
736
737 d_mat_clear(mm);
738 d_mat_clear(mn);
739 d_mat_clear(rm);
740 d_mat_clear(rn);
741
742 d_mat_init(boundt, n, n);
743
744 d_mat_transpose(boundt, bound);
745 for (i = 0; i < n; i++)
746 {
747 for (j = 0; j < n; j++)
748 {
749 d_mat_entry(bound2, i, j) =
750 fabs(d_mat_entry(bound2, i, j)) +
751 fabs(d_mat_entry(bound3, i, j));
752 }
753 }
754 d_mat_mul_classical(bound, bound2, bound);
755 d_mat_mul_classical(bound, boundt, bound);
756
757 d_mat_clear(bound2);
758 d_mat_clear(bound3);
759 d_mat_clear(boundt);
760
761 norm = 0;
762 for (i = 0; i < n; i++)
763 {
764 s = 0;
765 for (j = 0; j < n; j++)
766 {
767 s += fabs(d_mat_entry(bound, i, j));
768 }
769 norm = FLINT_MAX(norm, s);
770 }
771 if (norm >= 1)
772 {
773 d_mat_clear(R);
774 d_mat_clear(bound);
775 fesetround(rounding_direction);
776 return 0;
777 }
778
779 d_mat_init(absR, n, n);
780 for (i = 0; i < n; i++)
781 {
782 for (j = 0; j < n; j++)
783 {
784 if (j >= i)
785 {
786 d_mat_entry(bound, i, j) += norm * norm / (1.0 - norm);
787 }
788 else
789 {
790 d_mat_entry(bound, i, j) = 0;
791 }
792 d_mat_entry(absR, i, j) = fabs(d_mat_entry(R, i, j));
793 }
794 }
795 d_mat_mul_classical(bound, bound, absR);
796
797 d_mat_clear(absR);
798
799 for (i = 0; i < n - 1; i++)
800 {
801 fesetround(FE_DOWNWARD);
802 ti = (d_mat_entry(R, i, i) - d_mat_entry(bound, i, i)) * fl->eta;
803 fesetround(FE_UPWARD);
804 for (j = i + 1; j < n; j++)
805 {
806 tj = fabs(d_mat_entry(R, i, j)) + d_mat_entry(bound, i, j);
807 if (tj > ti)
808 {
809 d_mat_clear(R);
810 d_mat_clear(bound);
811 fesetround(rounding_direction);
812 return 0;
813 }
814 }
815 ti = d_mat_entry(R, i, i) + d_mat_entry(bound, i, i);
816 fesetround(FE_DOWNWARD);
817 tj = d_mat_entry(R, i + 1, i + 1) - d_mat_entry(bound, i + 1,
818 i + 1);
819 s = ((fabs(d_mat_entry(R, i, i + 1)) -
820 d_mat_entry(bound, i,
821 i + 1)) / ti) * ((fabs(d_mat_entry(R, i,
822 i + 1)) -
823 d_mat_entry(bound, i,
824 i + 1)) / ti) -
825 fl->delta;
826 s = -s;
827 fesetround(FE_UPWARD);
828 s = sqrt(s) * ti;
829 if (s > tj)
830 {
831 d_mat_clear(R);
832 d_mat_clear(bound);
833 fesetround(rounding_direction);
834 return 0;
835 }
836 }
837
838 d_mat_clear(R);
839 d_mat_clear(bound);
840 fesetround(rounding_direction);
841 }
842 return 1;
843 #else
844 return 0;
845 #endif
846 }
847