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 "gmp.h"
13 #include "flint.h"
14 #include "fft.h"
15
fft_adjust_sqrt2(mp_limb_t * r,mp_limb_t * i1,mp_size_t i,mp_size_t limbs,flint_bitcnt_t w,mp_limb_t * temp)16 void fft_adjust_sqrt2(mp_limb_t * r, mp_limb_t * i1,
17 mp_size_t i, mp_size_t limbs, flint_bitcnt_t w, mp_limb_t * temp)
18 {
19 flint_bitcnt_t wn = limbs*FLINT_BITS;
20 mp_limb_t cy;
21 mp_size_t j = i/2, k = w/2;
22 mp_size_t y;
23 flint_bitcnt_t b1;
24 int negate = 0;
25
26 b1 = j + wn/4 + i*k;
27 if (b1 >= wn)
28 {
29 negate = 1;
30 b1 -= wn;
31 }
32 y = b1/FLINT_BITS;
33 b1 = b1%FLINT_BITS;
34
35 /* multiply by 2^{j + wn/4 + i*k} */
36 if (y)
37 {
38 flint_mpn_copyi(temp + y, i1, limbs - y);
39 cy = mpn_neg_n(temp, i1 + limbs - y, y);
40 temp[limbs] = 0;
41 mpn_addmod_2expp1_1(temp + y, limbs - y, -i1[limbs]);
42 mpn_sub_1(temp + y, temp + y, limbs - y + 1, cy);
43 mpn_mul_2expmod_2expp1(r, temp, limbs, b1);
44 } else
45 mpn_mul_2expmod_2expp1(r, i1, limbs, b1);
46
47 /* multiply by 2^{wn/2} */
48 y = limbs/2;
49 cy = 0;
50
51 flint_mpn_copyi(temp + y, r, limbs - y);
52 temp[limbs] = 0;
53 if (y) cy = mpn_neg_n(temp, r + limbs - y, y);
54 mpn_addmod_2expp1_1(temp + y, limbs - y, -r[limbs]);
55 mpn_sub_1(temp + y, temp + y, limbs - y + 1, cy);
56
57 /* shift by an additional half limb (rare) */
58 if (limbs & 1)
59 mpn_mul_2expmod_2expp1(temp, temp, limbs, FLINT_BITS/2);
60
61 /* subtract */
62 if (negate)
63 mpn_sub_n(r, r, temp, limbs + 1);
64 else
65 mpn_sub_n(r, temp, r, limbs + 1);
66 }
67
68