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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "gmp.h"
13 #include "flint.h"
14 #include "fft.h"
15
butterfly_rshB(mp_limb_t * t,mp_limb_t * u,mp_limb_t * i1,mp_limb_t * i2,mp_size_t limbs,mp_size_t x,mp_size_t y)16 void butterfly_rshB(mp_limb_t * t, mp_limb_t * u, mp_limb_t * i1,
17 mp_limb_t * i2, mp_size_t limbs, mp_size_t x, mp_size_t y)
18 {
19 mp_limb_t cy, cy1, cy2, cy3;
20
21 if (x == 0)
22 {
23 if (y == 0)
24 {
25 cy = fft_sumdiff(t, u, i1, i2, limbs + 1);
26 } else /* y != 0 */
27 {
28 cy = fft_sumdiff(t, u, i1, i2 + y, limbs - y);
29 cy1 = (cy>>1);
30 cy2 = -(cy&1);
31 cy = fft_sumdiff(u + limbs - y, t + limbs - y, i1 + limbs - y, i2, y);
32 u[limbs] = (cy>>1) + i1[limbs];
33 t[limbs] = i1[limbs] - (cy&1);
34 mpn_addmod_2expp1_1(t + limbs - y, y, cy1 + i2[limbs]);
35 mpn_addmod_2expp1_1(u + limbs - y, y, cy2 - i2[limbs]);
36 }
37 } else if (y == 0) /* x != 0 */
38 {
39 cy = fft_sumdiff(t, u, i1 + x, i2, limbs - x);
40 cy1 = (cy>>1);
41 cy2 = -(cy&1);
42 cy3 = mpn_neg_n(i1, i1, x);
43 cy = fft_sumdiff(t + limbs - x, u + limbs - x, i1, i2 + limbs - x, x);
44 u[limbs] = -cy3 - (cy&1) - i2[limbs];
45 t[limbs] = -cy3 + i2[limbs] + (cy>>1);
46 mpn_addmod_2expp1_1(t + limbs - x, x, cy1 + i1[limbs]);
47 mpn_addmod_2expp1_1(u + limbs - x, x, cy2 + i1[limbs]);
48 } else if (x == y)
49 {
50 cy = fft_sumdiff(t, u, i1 + x, i2 + x, limbs - x);
51 cy1 = (cy>>1);
52 cy2 = -(cy&1);
53 cy = fft_sumdiff(t + limbs - x, u + limbs - x, i2, i1, x);
54 cy3 = mpn_neg_n(t + limbs - x, t + limbs - x, x);
55 u[limbs] = -(cy&1);
56 t[limbs] = -(cy>>1) - cy3;
57 mpn_addmod_2expp1_1(t + limbs - x, x, cy1 + i1[limbs] + i2[limbs]);
58 mpn_addmod_2expp1_1(u + limbs - x, x, cy2 + i1[limbs] - i2[limbs]);
59 } else if (x > y)
60 {
61 cy = fft_sumdiff(t + limbs - y, u + limbs - y, i2, i1 + x - y, y);
62 cy3 = mpn_neg_n(t + limbs - y, t + limbs - y, y);
63 t[limbs] = -(cy>>1) - cy3;
64 u[limbs] = -(cy&1);
65 cy3 = mpn_neg_n(i1, i1, x - y);
66 cy = fft_sumdiff(t + limbs - x, u + limbs - x, i1, i2 + limbs - x + y, x - y);
67 mpn_addmod_2expp1_1(t + limbs - y, y, (cy>>1) + i2[limbs] - cy3);
68 mpn_addmod_2expp1_1(u + limbs - y, y, -(cy&1) - i2[limbs] - cy3);
69 cy = fft_sumdiff(t, u, i1 + x, i2 + y, limbs - x);
70 mpn_addmod_2expp1_1(t + limbs - x, x, (cy>>1) + i1[limbs]);
71 mpn_addmod_2expp1_1(u + limbs - x, x, -(cy&1) + i1[limbs]);
72 } else /* x < y */
73 {
74 cy = fft_sumdiff(t + limbs - x, u + limbs - x, i2 + y - x, i1, x);
75 cy3 = mpn_neg_n(t + limbs - x, t + limbs - x, x);
76 t[limbs] = -(cy>>1) - cy3;
77 u[limbs] = -(cy&1);
78 cy3 = mpn_neg_n(i2, i2, y - x);
79 cy = fft_sumdiff(t + limbs - y, u + limbs - y, i1 + limbs - y + x, i2, y - x);
80 mpn_addmod_2expp1_1(t + limbs - x, x, (cy>>1) + i1[limbs] - cy3);
81 mpn_addmod_2expp1_1(u + limbs - x, x, -(cy&1) + i1[limbs] + cy3);
82 cy = fft_sumdiff(t, u, i1 + x, i2 + y, limbs - y);
83 mpn_addmod_2expp1_1(t + limbs - y, y, (cy>>1) + i2[limbs]);
84 mpn_addmod_2expp1_1(u + limbs - y, y, -(cy&1) - i2[limbs]);
85 }
86 }
87