1 /*
2 zn_mod.c: functions operating on zn_mod_t objects
3
4 Copyright (C) 2007, 2008, David Harvey
5
6 This file is part of the zn_poly library (version 0.9).
7
8 This program is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 2 of the License, or
11 (at your option) version 3 of the License.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program. If not, see <http://www.gnu.org/licenses/>.
20
21 */
22
23 #include "zn_poly_internal.h"
24
25
26 void
zn_mod_init(zn_mod_t mod,ulong m)27 zn_mod_init (zn_mod_t mod, ulong m)
28 {
29 ZNP_ASSERT (m >= 2);
30
31 mod->m = m;
32 mod->bits = ceil_lg (m);
33
34 mpz_t x, y;
35 mpz_init (x);
36 mpz_init (y);
37
38 // compute B and B^2 mod m
39 mpz_set_ui (x, 1);
40 mpz_mul_2exp (x, x, ULONG_BITS);
41 mpz_mod_ui (x, x, m);
42 mod->B = mpz_get_ui (x);
43
44 mpz_set_ui (x, 1);
45 mpz_mul_2exp (x, x, 2*ULONG_BITS);
46 mpz_mod_ui (x, x, m);
47 mod->B2 = mpz_get_ui (x);
48
49 // compute sh1 and inv1
50 mod->sh1 = ceil_lg (m) - 1;
51 mpz_set_ui (x, 1);
52 mpz_mul_2exp (x, x, mod->sh1 + 1);
53 mpz_sub_ui (x, x, m);
54 mpz_mul_2exp (x, x, ULONG_BITS);
55 mpz_fdiv_q_ui (x, x, m);
56 mpz_add_ui (x, x, 1);
57 mod->inv1 = mpz_get_ui (x);
58
59 // compute sh2, sh3, inv2, m_norm
60 unsigned ell = floor_lg (m) + 1;
61 mod->sh2 = ULONG_BITS - ell;
62 mod->sh3 = ell - 1;
63 mod->m_norm = m << mod->sh2;
64 mpz_set_ui (x, 1);
65 mpz_mul_2exp (x, x, ell);
66 mpz_sub_ui (x, x, m);
67 mpz_mul_2exp (x, x, ULONG_BITS);
68 mpz_sub_ui (x, x, 1);
69 mpz_fdiv_q_ui (x, x, m);
70 mod->inv2 = mpz_get_ui (x);
71
72 // compute inv3, if m is odd
73 if (m & 1)
74 {
75 // m^(-1) = m mod 8
76 ulong minv = m;
77 // lift 2-adically
78 int i;
79 for (i = 3; i < ULONG_BITS; i <<= 1)
80 minv = 2 * minv - m * minv * minv;
81 mod->inv3 = minv;
82 }
83
84 mpz_clear (y);
85 mpz_clear (x);
86 }
87
88
89 void
zn_mod_clear(zn_mod_t mod)90 zn_mod_clear (zn_mod_t mod)
91 {
92 // nothing to do yet, but maybe one day there will be
93 }
94
95
96
97 ulong
zn_mod_pow2(int k,const zn_mod_t mod)98 zn_mod_pow2 (int k, const zn_mod_t mod)
99 {
100 ZNP_ASSERT (mod->m & 1);
101 ZNP_ASSERT (k > -ULONG_BITS && k < ULONG_BITS);
102
103 if (k == 0)
104 return 1;
105
106 if (k > 0)
107 return zn_mod_reduce (1UL << k, mod);
108
109 return zn_mod_pow (zn_mod_divby2 (1, mod), -k, mod);
110 }
111
112
113
114 ulong
zn_mod_pow(ulong x,long k,const zn_mod_t mod)115 zn_mod_pow (ulong x, long k, const zn_mod_t mod)
116 {
117 ZNP_ASSERT (k >= 0);
118
119 // repeated squaring
120 ulong prod = 1;
121 ulong x_pow = x;
122 for (; k; k >>= 1)
123 {
124 if (k & 1)
125 prod = zn_mod_mul (prod, x_pow, mod);
126 x_pow = zn_mod_mul (x_pow, x_pow, mod);
127 }
128 return prod;
129 }
130
131
132
133 ulong
zn_mod_invert(ulong x,const zn_mod_t mod)134 zn_mod_invert (ulong x, const zn_mod_t mod)
135 {
136 ZNP_ASSERT (x < mod->m);
137
138 // for now just use GMP
139 mpz_t a, m;
140 mpz_init (a);
141 mpz_set_ui (a, x);
142 mpz_init (m);
143 mpz_set_ui (m, mod->m);
144 int success = mpz_invert (a, a, m);
145 x = mpz_get_ui (a);
146 mpz_clear (m);
147 mpz_clear (a);
148
149 return success ? x : 0;
150 }
151
152
153 // end of file ****************************************************************
154