1 /*
2 Copyright (C) 2009, 2011 William Hart
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 "stdlib.h"
13 #include "gmp.h"
14 #include "flint.h"
15 #include "fft.h"
16
fft_combine_limbs(mp_limb_t * res,mp_limb_t ** poly,slong length,mp_size_t coeff_limbs,mp_size_t output_limbs,mp_size_t total_limbs)17 void fft_combine_limbs(mp_limb_t * res, mp_limb_t ** poly, slong length,
18 mp_size_t coeff_limbs, mp_size_t output_limbs, mp_size_t total_limbs)
19 {
20 mp_size_t skip, i;
21
22 for (skip = 0, i = 0; i < length && skip + output_limbs + 1 <= total_limbs; i++, skip += coeff_limbs)
23 mpn_add(res + skip, res + skip, output_limbs + 1, poly[i], output_limbs);
24
25 while ((skip < total_limbs) && (i < length))
26 {
27 mpn_add(res + skip, res + skip, total_limbs - skip, poly[i], FLINT_MIN(total_limbs - skip, output_limbs));
28
29 i++;
30
31 skip += coeff_limbs;
32 }
33 }
34
fft_combine_bits(mp_limb_t * res,mp_limb_t ** poly,slong length,flint_bitcnt_t bits,mp_size_t output_limbs,mp_size_t total_limbs)35 void fft_combine_bits(mp_limb_t * res, mp_limb_t ** poly, slong length,
36 flint_bitcnt_t bits, mp_size_t output_limbs, mp_size_t total_limbs)
37 {
38 flint_bitcnt_t shift_bits, top_bits = ((FLINT_BITS - 1) & bits);
39 mp_size_t coeff_limbs, i;
40 mp_limb_t * temp, * limb_ptr, * end;
41
42 if (top_bits == 0)
43 {
44 fft_combine_limbs(res, poly, length, bits/FLINT_BITS, output_limbs, total_limbs);
45 return;
46 }
47
48 coeff_limbs = (bits/FLINT_BITS) + 1;
49 temp = flint_malloc((output_limbs + 1)*sizeof(mp_limb_t));
50 shift_bits = 0;
51 limb_ptr = res;
52 end = res + total_limbs;
53
54 for (i = 0; i < length && limb_ptr + output_limbs + 1 < end; i++)
55 {
56 if (shift_bits)
57 {
58 mpn_lshift(temp, poly[i], output_limbs + 1, shift_bits);
59 mpn_add_n(limb_ptr, limb_ptr, temp, output_limbs + 1);
60 } else
61 mpn_add(limb_ptr, limb_ptr, output_limbs + 1, poly[i], output_limbs);
62
63 shift_bits += top_bits;
64 limb_ptr += (coeff_limbs - 1);
65
66 if (shift_bits >= FLINT_BITS)
67 {
68 limb_ptr++;
69 shift_bits -= FLINT_BITS;
70 }
71 }
72
73 while (limb_ptr < end && i < length)
74 {
75 if (shift_bits)
76 {
77 mpn_lshift(temp, poly[i], output_limbs + 1, shift_bits);
78 mpn_add_n(limb_ptr, limb_ptr, temp, end - limb_ptr);
79 } else
80 mpn_add_n(limb_ptr, limb_ptr, poly[i], end - limb_ptr);
81
82 shift_bits += top_bits;
83 limb_ptr += (coeff_limbs - 1);
84
85 if (shift_bits >= FLINT_BITS)
86 {
87 limb_ptr++;
88 shift_bits -= FLINT_BITS;
89 }
90
91 i++;
92 }
93
94 flint_free(temp);
95 }
96