1 /*=============================================================================
2
3 This file is part of Antic.
4
5 Antic is free software: you can redistribute it and/or modify it under
6 the terms of the GNU Lesser General Public License (LGPL) as published
7 by the Free Software Foundation; either version 2.1 of the License, or
8 (at your option) any later version. See <http://www.gnu.org/licenses/>.
9
10 =============================================================================*/
11 /******************************************************************************
12
13 Copyright (C) 2012 William Hart
14
15 ******************************************************************************/
16
17 #include <stdlib.h>
18 #include <gmp.h>
19 #include "qfb.h"
20
qfb_reduce(qfb_t r,qfb_t f,fmpz_t D)21 void qfb_reduce(qfb_t r, qfb_t f, fmpz_t D)
22 {
23 int done = 0;
24 fmpz_t t;
25
26 qfb_set(r, f);
27
28 fmpz_init(t);
29
30 while(!done)
31 {
32 done = 1;
33
34 if (fmpz_cmp(r->c, r->a) < 0)
35 {
36 fmpz_swap(r->a, r->c);
37 fmpz_neg(r->b, r->b);
38
39 done = 0;
40 }
41
42 if (fmpz_cmpabs(r->b, r->a) > 0)
43 {
44 fmpz_add(t, r->a, r->a);
45 fmpz_fdiv_r(r->b, r->b, t);
46 if (fmpz_cmp(r->b, r->a) > 0)
47 fmpz_sub(r->b, r->b, t);
48
49 fmpz_add(t, t, t);
50 fmpz_mul(r->c, r->b, r->b);
51 fmpz_sub(r->c, r->c, D);
52 fmpz_divexact(r->c, r->c, t);
53
54 done = 0;
55 }
56 }
57
58 if (fmpz_cmpabs(r->a, r->b) == 0 || fmpz_cmp(r->a, r->c) == 0)
59 if (fmpz_sgn(r->b) < 0)
60 fmpz_neg(r->b, r->b);
61
62 fmpz_clear(t);
63 }
64