General robustification.
[u/mdw/catacomb] / f-prime.c
1 /* -*-c-*-
2 *
3 * $Id: f-prime.c,v 1.11 2004/04/03 03:32:05 mdw Exp $
4 *
5 * Prime fields with Montgomery arithmetic
6 *
7 * (c) 2001 Straylight/Edgeware
8 */
9
10 /*----- Licensing notice --------------------------------------------------*
11 *
12 * This file is part of Catacomb.
13 *
14 * Catacomb is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Library General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * Catacomb is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Library General Public License for more details.
23 *
24 * You should have received a copy of the GNU Library General Public
25 * License along with Catacomb; if not, write to the Free
26 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27 * MA 02111-1307, USA.
28 */
29
30 /*----- Revision history --------------------------------------------------*
31 *
32 * $Log: f-prime.c,v $
33 * Revision 1.11 2004/04/03 03:32:05 mdw
34 * General robustification.
35 *
36 * Revision 1.10 2004/04/02 01:03:49 mdw
37 * Miscellaneous constification.
38 *
39 * Revision 1.9 2004/04/01 21:28:41 mdw
40 * Normal basis support (translates to poly basis internally). Rewrite
41 * EC and prime group table generators in awk, so that they can reuse data
42 * for repeated constants.
43 *
44 * Revision 1.8 2004/04/01 12:50:09 mdw
45 * Add cyclic group abstraction, with test code. Separate off exponentation
46 * functions for better static linking. Fix a buttload of bugs on the way.
47 * Generally ensure that negative exponents do inversion correctly. Add
48 * table of standard prime-field subgroups. (Binary field subgroups are
49 * currently unimplemented but easy to add if anyone ever finds a good one.)
50 *
51 * Revision 1.7 2004/03/27 17:54:11 mdw
52 * Standard curves and curve checking.
53 *
54 * Revision 1.6 2004/03/23 15:19:32 mdw
55 * Test elliptic curves more thoroughly.
56 *
57 * Revision 1.5 2004/03/23 12:08:26 mdw
58 * Random field-element selection.
59 *
60 * Revision 1.4 2004/03/21 22:52:06 mdw
61 * Merge and close elliptic curve branch.
62 *
63 * Revision 1.3.4.3 2004/03/21 22:39:46 mdw
64 * Elliptic curves on binary fields work.
65 *
66 * Revision 1.3.4.2 2004/03/20 00:13:31 mdw
67 * Projective coordinates for prime curves
68 *
69 * Revision 1.3.4.1 2003/06/10 13:43:53 mdw
70 * Simple (non-projective) curves over prime fields now seem to work.
71 *
72 * Revision 1.3 2003/05/15 23:25:59 mdw
73 * Make elliptic curve stuff build.
74 *
75 * Revision 1.2 2002/01/13 13:48:44 mdw
76 * Further progress.
77 *
78 * Revision 1.1 2001/04/29 18:12:33 mdw
79 * Prototype version.
80 *
81 */
82
83 /*----- Header files ------------------------------------------------------*/
84
85 #include <mLib/sub.h>
86
87 #include "field.h"
88 #include "mpmont.h"
89 #include "mprand.h"
90
91 /*----- Main code ---------------------------------------------------------*/
92
93 typedef struct fctx {
94 field f;
95 mpmont mm;
96 } fctx;
97
98 /* --- Field operations --- */
99
100 static void fdestroy(field *ff)
101 { fctx *f = (fctx *)ff; mpmont_destroy(&f->mm); DESTROY(f); }
102
103 static mp *frand(field *ff, mp *d, grand *r)
104 { fctx *f = (fctx *)ff; return (mprand_range(d, f->mm.m, r, 0)); }
105
106 static mp *fin(field *ff, mp *d, mp *x) {
107 fctx *f = (fctx *)ff;
108 mp_div(0, &d, x, f->mm.m);
109 return (mpmont_mul(&f->mm, d, d, f->mm.r2));
110 }
111
112 static mp *fout(field *ff, mp *d, mp *x)
113 { fctx *f = (fctx *)ff; return (mpmont_reduce(&f->mm, d, x)); }
114
115 static int fzerop(field *ff, mp *x) { return (!MP_LEN(x)); }
116
117 static mp *fneg(field *ff, mp *d, mp *x)
118 { fctx *f = (fctx *)ff; return (mp_sub(d, f->mm.m, x)); }
119
120 static mp *fadd(field *ff, mp *d, mp *x, mp *y) {
121 fctx *f = (fctx *)ff; d = mp_add(d, x, y);
122 if (d->f & MP_NEG) d = mp_add(d, d, f->mm.m);
123 else if (MP_CMP(d, >, f->mm.m)) d = mp_sub(d, d, f->mm.m);
124 return (d);
125 }
126
127 static mp *fsub(field *ff, mp *d, mp *x, mp *y) {
128 fctx *f = (fctx *)ff; d = mp_sub(d, x, y);
129 if (d->f & MP_NEG) d = mp_add(d, d, f->mm.m);
130 else if (MP_CMP(d, >, f->mm.m)) d = mp_sub(d, d, f->mm.m);
131 return (d);
132 }
133
134 static mp *fmul(field *ff, mp *d, mp *x, mp *y)
135 { fctx *f = (fctx *)ff; return (mpmont_mul(&f->mm, d, x, y)); }
136
137 static mp *fsqr(field *ff, mp *d, mp *x) {
138 fctx *f = (fctx *)ff; d = mp_sqr(d, x);
139 return (mpmont_reduce(&f->mm, d, d));
140 }
141
142 static mp *finv(field *ff, mp *d, mp *x) {
143 fctx *f = (fctx *)ff; d = mpmont_reduce(&f->mm, d, x);
144 mp_gcd(0, 0, &d, f->mm.m, d); return (mpmont_mul(&f->mm, d, d, f->mm.r2));
145 }
146
147 static mp *freduce(field *ff, mp *d, mp *x)
148 { fctx *f = (fctx *)ff; mp_div(0, &d, x, f->mm.m); return (d); }
149
150 static mp *fsqrt(field *ff, mp *d, mp *x) {
151 fctx *f = (fctx *)ff; d = mpmont_reduce(&f->mm, d, x);
152 d = mp_modsqrt(d, d, f->mm.m); if (!d) return (d);
153 return (mpmont_mul(&f->mm, d, d, f->mm.r2));
154 }
155
156 static mp *fdbl(field *ff, mp *d, mp *x) {
157 fctx *f = (fctx *)ff; d = mp_lsl(d, x, 1);
158 if (MP_CMP(d, >, f->mm.m)) d = mp_sub(d, d, f->mm.m);
159 return (d);
160 }
161
162 static mp *ftpl(field *ff, mp *d, mp *x) {
163 fctx *f = (fctx *)ff; MP_DEST(d, MP_LEN(x) + 1, x->f);
164 MPX_UMULN(d->v, d->vl, x->v, x->vl, 3);
165 while (MP_CMP(d, >, f->mm.m)) d = mp_sub(d, d, f->mm.m);
166 return (d);
167 }
168
169 static mp *fqdl(field *ff, mp *d, mp *x) {
170 fctx *f = (fctx *)ff; d = mp_lsl(d, x, 2);
171 while (MP_CMP(d, >, f->mm.m)) d = mp_sub(d, d, f->mm.m);
172 return (d);
173 }
174
175 static mp *fhlv(field *ff, mp *d, mp *x) {
176 fctx *f = (fctx *)ff;
177 if (!MP_LEN(x)) { MP_COPY(x); MP_DROP(d); return (x); }
178 if (x->v[0] & 1) { d = mp_add(d, x, f->mm.m); x = d; }
179 return (mp_lsr(d, x, 1));
180 }
181
182 /* --- Field operations table --- */
183
184 static const field_ops fops = {
185 FTY_PRIME, "prime",
186 fdestroy, frand, field_stdsamep,
187 fin, fout,
188 fzerop, fneg, fadd, fsub, fmul, fsqr, finv, freduce, fsqrt,
189 0,
190 fdbl, ftpl, fqdl, fhlv
191 };
192
193 /* --- @field_prime@ --- *
194 *
195 * Arguments: @mp *p@ = the characteristic of the field
196 *
197 * Returns: A pointer to the field or null.
198 *
199 * Use: Creates a field structure for a prime field of size %$p$%,
200 * using Montgomery reduction for arithmetic.
201 */
202
203 field *field_prime(mp *p)
204 {
205 fctx *f;
206
207 if (!MP_ISPOS(p) || !MP_ISODD(p))
208 return (0);
209 f = CREATE(fctx);
210 f->f.ops = &fops;
211 mpmont_create(&f->mm, p);
212 f->f.zero = MP_ZERO;
213 f->f.one = f->mm.r;
214 f->f.m = f->mm.m;
215 f->f.nbits = mp_bits(p);
216 f->f.noctets = (f->f.nbits + 7) >> 3;
217 return (&f->f);
218 }
219
220 /*----- That's all, folks -------------------------------------------------*/