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 "flint/flint.h"
20 #include "flint/ulong_extras.h"
21 #include "flint/fmpz.h"
22 #include "qfb.h"
23
qfb_pow(qfb_t r,qfb_t f,fmpz_t D,fmpz_t e)24 void qfb_pow(qfb_t r, qfb_t f, fmpz_t D, fmpz_t e)
25 {
26 fmpz_t L, exp;
27 qfb_t pow;
28
29 if (fmpz_is_zero(e))
30 {
31 qfb_principal_form(r, D);
32 return;
33 }
34
35 if (fmpz_is_one(e))
36 {
37 qfb_set(r, f);
38 return;
39 }
40
41 fmpz_init(L);
42 fmpz_init(exp);
43
44 fmpz_set(exp, e);
45
46 fmpz_abs(L, D);
47 fmpz_root(L, L, 4);
48
49 qfb_init(pow);
50
51 qfb_set(pow, f);
52 while (fmpz_is_even(exp))
53 {
54 qfb_nudupl(pow, pow, D, L);
55 qfb_reduce(pow, pow, D);
56 fmpz_fdiv_q_2exp(exp, exp, 1);
57 }
58
59 qfb_set(r, pow);
60 fmpz_fdiv_q_2exp(exp, exp, 1);
61
62 while (!fmpz_is_zero(exp))
63 {
64 qfb_nudupl(pow, pow, D, L);
65 qfb_reduce(pow, pow, D);
66 if (fmpz_is_odd(exp))
67 {
68 qfb_nucomp(r, r, pow, D, L);
69 qfb_reduce(r, r, D);
70 }
71 fmpz_fdiv_q_2exp(exp, exp, 1);
72 }
73
74 qfb_clear(pow);
75 fmpz_clear(L);
76 fmpz_clear(exp);
77 }
78
qfb_pow_with_root(qfb_t r,qfb_t f,fmpz_t D,fmpz_t e,fmpz_t L)79 void qfb_pow_with_root(qfb_t r, qfb_t f, fmpz_t D, fmpz_t e, fmpz_t L)
80 {
81 fmpz_t exp;
82 qfb_t pow;
83
84 if (fmpz_is_zero(e))
85 {
86 qfb_principal_form(r, D);
87 return;
88 }
89
90 if (fmpz_is_one(e))
91 {
92 qfb_set(r, f);
93 return;
94 }
95
96 fmpz_init(exp);
97
98 fmpz_set(exp, e);
99
100 qfb_init(pow);
101
102 qfb_set(pow, f);
103 while (fmpz_is_even(exp))
104 {
105 qfb_nudupl(pow, pow, D, L);
106 qfb_reduce(pow, pow, D);
107 fmpz_fdiv_q_2exp(exp, exp, 1);
108 }
109
110 qfb_set(r, pow);
111 fmpz_fdiv_q_2exp(exp, exp, 1);
112
113 while (!fmpz_is_zero(exp))
114 {
115 qfb_nudupl(pow, pow, D, L);
116 qfb_reduce(pow, pow, D);
117 if (fmpz_is_odd(exp))
118 {
119 qfb_nucomp(r, r, pow, D, L);
120 qfb_reduce(r, r, D);
121 }
122 fmpz_fdiv_q_2exp(exp, exp, 1);
123 }
124
125 qfb_clear(pow);
126 fmpz_clear(exp);
127 }
128