1 /****************************************************************/
2 /* file aritz.c
3
4 ARIBAS interpreter for Arithmetic
5 Copyright (C) 1996-2007 O.Forster
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20
21 Address of the author
22
23 Otto Forster
24 Math. Institut der LMU
25 Theresienstr. 39
26 D-80333 Muenchen, Germany
27
28 Email forster@mathematik.uni-muenchen.de
29 WWW http://www.mathematik.uni-muenchen.de
30 */
31 /****************************************************************/
32 /*
33 ** aritz.c
34 ** functions for polynomials
35 ** and for arithmetic in GF(2**n)
36 **
37 ** date of last change
38 **
39 ** 2002-04-07: created
40 ** 2003-03-05: gf2n functions
41 ** 2007-04-04: gf2X functions
42 ** 2007-09-08: bugfix in gf2X_mod
43 ** 2007-09-18: gf2polsquare
44 ** 2007-10-01: gf2X_divide
45 */
46 #include "common.h"
47
48 /*-----------------------------------------------------------------*/
49 /* field extension Fp[sqrt(D)] */
50 typedef struct {
51 word2 *pp;
52 int plen;
53 word2 *D;
54 int dlen;
55 } FP2D;
56
57 typedef struct {
58 word2 *xx;
59 int xlen;
60 word2 *yy;
61 int ylen;
62 } PAIRXY;
63
64 /*-----------------------------------------------------------------*/
65 /* setbit and testbit suppose that vv is an array of word2 */
66 #define setbit(vv,i) vv[(i)>>4] |= (1 << ((i)&0xF))
67 #define testbit(vv,i) (vv[(i)>>4] & (1 << ((i)&0xF)))
68 /*-----------------------------------------------------------------*/
69 #define MULTFLAG 0
70 #define DIVFLAG 1
71 #define MODFLAG 2
72 #define MODNFLAG 4
73 #define DDIVFLAG (DIVFLAG | MODFLAG)
74 /*-----------------------------------------------------------------*/
75
76 PUBLIC void iniaritz _((void));
77
78 PUBLIC truc gf2nzero, gf2none;
79 PUBLIC truc gf2nintsym, gf2n_sym;
80 PUBLIC truc znXmultsym, znXsqsym, znXddivsym, znXdivsym, znXmodsym;
81 PUBLIC truc polmultsym, polNmultsym;
82 PUBLIC truc polmodsym, polNmodsym, poldivsym, polNdivsym;
83 PUBLIC truc addgf2ns _((truc *ptr));
84 PUBLIC truc multgf2ns _((truc *ptr));
85 PUBLIC truc exptgf2n _((truc *ptr));
86 PUBLIC truc divgf2ns _((truc *ptr));
87 PUBLIC int fpSqrt _((word2 *pp, int plen, word2 *aa, int alen,
88 word2 *zz, word2 *hilf));
89 PUBLIC int fp2Sqrt _((word2 *pp, int plen, word2 *aa, int alen,
90 word2 *zz, word2 *hilf));
91 PUBLIC unsigned fp_sqrt _((unsigned p, unsigned a));
92
93
94 PRIVATE truc Fgf2nint _((void));
95 PRIVATE truc FznXmult _((void));
96 PRIVATE truc FznXsquare _((void));
97 PRIVATE truc FznXdivide _((void));
98 PRIVATE truc FznXdiv _((void));
99 PRIVATE truc FznXmod _((void));
100 PRIVATE truc Fpolmult _((void));
101 PRIVATE truc FpolNmult _((void));
102 PRIVATE truc Fpolmod _((void));
103 PRIVATE truc FpolNmod _((void));
104 PRIVATE truc Fpoldiv _((void));
105 PRIVATE truc FpolNdiv _((void));
106 PRIVATE truc znXdivmodsymb _((int mode));
107 PRIVATE int chkznXmultargs _((truc sym, truc *argptr));
108 PRIVATE int chkpolmultargs _((truc sym, truc *argptr));
109 PRIVATE int chkpoldivargs _((truc sym, truc *argptr));
110 PRIVATE truc znXpolmult _((truc *mptr, truc *argptr));
111 PRIVATE truc znXpolsquare _((truc *mptr, truc *argptr));
112 PRIVATE truc znXpoldiv _((truc *mptr, truc *argptr, int mode));
113 PRIVATE truc znXpolmod _((truc *mptr, truc *argptr));
114 PRIVATE truc multintpols _((truc *argptr, int mode));
115 PRIVATE truc modintpols _((truc *argptr, int mode));
116
117 PRIVATE int gf2polmod _((word2 *x, int n, word2 *y, int m));
118 PRIVATE int gf2polmult _((word2 *x, int n, word2 *y, int m, word2 *z));
119 PRIVATE int gf2poldivide _((word2 *x, int n, word2 *y, int m, word2 *z, int *rlenptr));
120 PRIVATE int gf2poldiv _((word2 *x, int n, word2 *y, int m, word2 *z));
121 PRIVATE int gf2ninverse _((word2 *x, int n, word2 *z, word2 *uu));
122 PRIVATE int gf2npower _((word2 *x, int n, word2 *y, int m,
123 word2 *z, word2 *hilf));
124 PRIVATE int gf2polsquare _((word2 *x, int n, word2 *z));
125 PRIVATE int gf2polmodpow _((word2 *x, int n, word2 *y, int m,
126 word2 *f, int flen, word2 *z, word2 *hilf));
127 PRIVATE int gf2polgcd _((word2 *x, int n, word2 *y, int m));
128 PRIVATE int gf2polgcdx _((word2 *x, int n, word2 *y, int m,
129 word2 *z, int *zlenptr, word2 *hilf));
130 PRIVATE int gf2polirred _((word2 *x, int n, word2 *y, word2 *hilf));
131 PRIVATE int gf2polirred1 _((word2 *x, int n, word2 *y, word2 *hilf));
132 PRIVATE int gf2ntrace _((word2 *x, int n));
133 PRIVATE int shiftleft1 _((word2 *x, int n));
134 PRIVATE int bitxorshift _((word2 *x, int n, word2 *y, int m, int s));
135 PRIVATE unsigned gf2polfindirr _((int n));
136
137 PRIVATE truc gf2ndegsym, gf2npolsym, gf2ninisym, gf2ntrsym, maxgf2nsym;
138 PRIVATE truc Fgf2ninit _((void));
139 PRIVATE truc Fgf2ndegree _((void));
140 PRIVATE truc Fgf2nfieldpol _((void));
141 PRIVATE truc Fgf2ntrace _((void));
142 PRIVATE truc Fmaxgf2n _((void));
143 PRIVATE int gf2nmod _((word2 *x, int n));
144
145 PRIVATE truc gf2Xmulsym, gf2Xsqsym, gf2Xddivsym, gf2Xdivsym, gf2Xmodsym, gf2Xgcdsym;
146 PRIVATE truc gf2Xmpowsym, gf2Xprimsym;
147 PRIVATE truc Fgf2Xmult _((void));
148 PRIVATE truc Fgf2Xsquare _((void));
149 PRIVATE truc Fgf2Xdivide _((void));
150 PRIVATE truc Fgf2Xdiv _((void));
151 PRIVATE truc Fgf2Xmod _((void));
152 PRIVATE truc Fgf2Xgcd _((void));
153 PRIVATE truc Fgf2Xmodpow _((void));
154 PRIVATE truc Fgf2Xprimtest _((void));
155
156 /*----------------------------------------------------*/
157 PRIVATE void fp2Dmult _((FP2D *pfp2D, PAIRXY *pZ1, PAIRXY *pZ2, word2 *hilf));
158 PRIVATE void fp2Dsquare _((FP2D *pfp2D, PAIRXY *pZ, word2 *hilf));
159 PRIVATE void fp2Dpower _((FP2D *pfp2D, PAIRXY *pZ, word2 *ex, int exlen,
160 word2 *hilf));
161
162 PRIVATE long nonresdisc _((word2 *pp, int plen, word2 *aa, int alen,
163 word2 *hilf));
164
165 PRIVATE int fpSqrt58 _((word2 *pp, int plen, word2 *aa, int alen,
166 word2 *zz, word2 *hilf));
167 PRIVATE int fpSqrt14 _((word2 *pp, int plen, word2 *aa, int alen,
168 word2 *zz, word2 *hilf));
169 PRIVATE unsigned fp_sqrt14 _((unsigned p, unsigned a));
170 PRIVATE void fp2pow _((unsigned p, unsigned D, unsigned *uu, unsigned n));
171
172 PRIVATE truc Ffpsqrt _((void));
173
174 PRIVATE truc fpsqrtsym;
175
176 /* #define NAUSKOMM */
177 #ifdef NAUSKOMM
178 PRIVATE truc gggsym, gg1sym, gg2sym;
179 PRIVATE truc Fggg _((void));
180 PRIVATE truc Fgg1 _((void));
181 PRIVATE truc Fgg2 _((void));
182 #endif
183
184 #define GF2XARITH
185 /* #define ZnXARITH */
186 /*------------------------------------------------------------------*/
iniaritz()187 PUBLIC void iniaritz()
188 {
189 word2 x[1];
190
191 gf2nzero = mk0gf2n(x,0);
192 x[0] = 1;
193 gf2none = mk0gf2n(x,1);
194
195 gf2nintsym = newsym ("gf2nint", sTYPESPEC, gf2nzero);
196 gf2n_sym = new0symsig("gf2nint", sFBINARY, (wtruc)Fgf2nint, s_1);
197 gf2ntrsym = newsymsig("gf2n_trace",sFBINARY, (wtruc)Fgf2ntrace, s_1);
198 gf2ninisym = newsymsig("gf2n_init", sFBINARY, (wtruc)Fgf2ninit, s_1);
199 gf2ndegsym = newsymsig("gf2n_degree", sFBINARY, (wtruc)Fgf2ndegree, s_0);
200 gf2npolsym = newsymsig("gf2n_fieldpol",sFBINARY,(wtruc)Fgf2nfieldpol,s_0);
201 maxgf2nsym = newsymsig("max_gf2nsize",sFBINARY, (wtruc)Fmaxgf2n, s_0);
202
203 #ifdef GF2XARITH
204 gf2Xmulsym = newsymsig("gf2X_mult", sFBINARY, (wtruc)Fgf2Xmult, s_2);
205 gf2Xsqsym = newsymsig("gf2X_square", sFBINARY, (wtruc)Fgf2Xsquare, s_1);
206 gf2Xddivsym= newsymsig("gf2X_divide", sFBINARY, (wtruc)Fgf2Xdivide, s_2);
207 gf2Xdivsym = newsymsig("gf2X_div", sFBINARY, (wtruc)Fgf2Xdiv, s_2);
208 gf2Xmodsym = newsymsig("gf2X_mod", sFBINARY, (wtruc)Fgf2Xmod, s_2);
209 gf2Xgcdsym = newsymsig("gf2X_gcd", sFBINARY, (wtruc)Fgf2Xgcd, s_2);
210 gf2Xmpowsym = newsymsig("gf2X_modpower", sFBINARY, (wtruc)Fgf2Xmodpow, s_3);
211 gf2Xprimsym = newsymsig("gf2X_primetest",sFBINARY,(wtruc)Fgf2Xprimtest, s_1);
212 #endif /* GF2XARITH */
213
214 #ifdef ZnXARITH
215 znXmultsym = newsymsig("ZnX_mult", sFBINARY, (wtruc)FznXmult, s_3);
216 znXsqsym = newsymsig("ZnX_square",sFBINARY, (wtruc)FznXsquare, s_2);
217 znXddivsym = newsymsig("ZnX_divide",sFBINARY, (wtruc)FznXdivide, s_3);
218 znXdivsym = newsymsig("ZnX_div", sFBINARY, (wtruc)FznXdiv, s_3);
219 znXmodsym = newsymsig("ZnX_mod", sFBINARY, (wtruc)FznXmod, s_3);
220 #endif
221 #ifdef POLYARITH
222 polmultsym = newsymsig("pol_mult", sFBINARY, (wtruc)Fpolmult, s_2);
223 polNmultsym= newintsym("pol_mult() mod", sFBINARY,(wtruc)FpolNmult);
224 polmodsym = newsymsig("pol_mod", sFBINARY, (wtruc)Fpolmod, s_2);
225 polNmodsym = newintsym("pol_mod() mod", sFBINARY,(wtruc)FpolNmod);
226 poldivsym = newsymsig("pol_div", sFBINARY, (wtruc)Fpoldiv, s_2);
227 polNdivsym = newintsym("pol_div() mod", sFBINARY,(wtruc)FpolNdiv);
228 #endif /* POLYARITH */
229
230 fpsqrtsym = newsymsig("gfp_sqrt", sFBINARY, (wtruc)Ffpsqrt, s_2);
231
232 #ifdef NAUSKOMM
233 gggsym = newsymsig("ggg", sFBINARY, (wtruc)Fggg, s_2);
234 gg1sym = newsymsig("gg1", sFBINARY, (wtruc)Fgg1, s_2);
235 gg2sym = newsymsig("gg2", sFBINARY, (wtruc)Fgg2, s_3);
236 #endif
237 }
238 /*-------------------------------------------------------------------*/
Ffpsqrt()239 PRIVATE truc Ffpsqrt()
240 {
241 word2 *pp, *aa, *hilf;
242 int plen, alen, len, sign;
243
244 if(chkints(fpsqrtsym,argStkPtr-1,2) == aERROR)
245 return(brkerr());
246 plen = bigref(argStkPtr-1,&pp,&sign);
247 if(plen == 0 || (pp[0] & 1) == 0 || (plen == 1 && pp[0] <= 2)) {
248 error(fpsqrtsym,err_oddprim,argStkPtr[-1]);
249 return brkerr();
250 }
251 else if(plen > scrbufSize/16) {
252 error(fpsqrtsym,err_ovfl,voidsym);
253 return brkerr();
254 }
255 aa = AriScratch;
256 alen = bigretr(argStkPtr,aa,&sign);
257 if(sign)
258 alen = modnegbig(aa,alen,pp,plen,AriBuf);
259 else if(alen >= plen)
260 alen = modbig(aa,alen,pp,plen,AriBuf);
261
262 hilf = AriScratch + plen + 2;
263 len = fpSqrt(pp,plen,aa,alen,AriBuf,hilf);
264 if(len < 0) {
265 error(scratch("gfp_sqrt(p,a)"),
266 "p not prime or a not a square mod p",voidsym);
267 return brkerr();
268 }
269 return mkint(0,AriBuf,len);
270 }
271 /*-------------------------------------------------------------------*/
272 #ifdef NAUSKOMM
273 /*-------------------------------------------------------------------*/
Fggg()274 PRIVATE truc Fggg()
275 {
276 #if 1
277 word2 *pp, *aa, *hilf;
278 size_t N;
279 int plen, alen, len, sign;
280
281 if(chkints(gggsym,argStkPtr-1,2) == aERROR)
282 return(brkerr());
283 aa = AriScratch;
284 plen = bigref(argStkPtr-1,&pp,&sign);
285 alen = bigretr(argStkPtr,aa,&sign);
286 N = 2*plen + alen;
287 if(N > scrbufSize/16) {
288 error(gggsym,err_ovfl,voidsym);
289 return brkerr();
290 }
291 hilf = AriScratch + N;
292
293 len = fp2Sqrt(pp,plen,aa,alen,AriBuf,hilf);
294 if(len < 0) {
295 error(gggsym,"sqrt mod p*p failed",voidsym);
296 return brkerr();
297 }
298 return mkint(0,AriBuf,len);
299 #else
300 return zero;
301 #endif
302 }
303 /*-------------------------------------------------------------------*/
Fgg1()304 PRIVATE truc Fgg1()
305 {
306 word2 *pp, *aa, *hilf;
307 size_t N;
308 int plen, alen, len, sign;
309
310 if(chkints(gg1sym,argStkPtr-1,2) == aERROR)
311 return(brkerr());
312 aa = AriScratch;
313 plen = bigref(argStkPtr-1,&pp,&sign);
314 alen = bigretr(argStkPtr,aa,&sign);
315 N = plen + alen;
316 if(N > scrbufSize/12) {
317 error(gg1sym,err_ovfl,voidsym);
318 return brkerr();
319 }
320 hilf = AriScratch + N;
321 if(sign)
322 alen = modnegbig(aa,alen,pp,plen,hilf);
323
324 len = fpSqrt14(pp,plen,aa,alen,AriBuf,hilf);
325 if(len < 0) {
326 error(gg1sym,"sqrt mod p failed",voidsym);
327 return brkerr();
328 }
329 return mkint(0,AriBuf,len);
330 }
331 /*-------------------------------------------------------------------*/
332 /*
333 ** gg2(p,D,c) calculates (c + sqrt D)**(p+1)/2 in Fp[sqrt D]
334 ** and returns its x-coordinate
335 */
Fgg2()336 PRIVATE truc Fgg2()
337 {
338 word2 *pp, *D, *ex, *xx, *yy, *hilf;
339 int N, plen, dlen, exlen, xlen, sign;
340 FP2D fp2D;
341 PAIRXY Z;
342
343 if(chkints(gg1sym,argStkPtr-2,3) == aERROR)
344 return(brkerr());
345 D = hilf = AriScratch;
346 plen = bigref(argStkPtr-2,&pp,&sign);
347 dlen = bigretr(argStkPtr-1,D,&sign);
348 fp2D.pp = pp;
349 fp2D.plen = plen;
350 fp2D.D = D;
351 fp2D.dlen = dlen;
352
353 xx = AriBuf;
354 N = plen + dlen + 2;
355 ex = hilf + 2*N;
356 yy = hilf + 4*N;
357 hilf += 5*N;
358
359 xlen = bigretr(argStkPtr,xx,&sign);
360 Z.xx = xx;
361 Z.xlen = xlen;
362 Z.yy = yy; yy[0] = 1;
363 Z.ylen = 1;
364 cpyarr(pp,plen,ex);
365 exlen = incarr(ex,plen,1);
366 exlen = shiftarr(ex,exlen,-1);
367 fp2Dpower(&fp2D,&Z,ex,exlen,hilf);
368 xlen = Z.xlen;
369
370 return mkint(0,xx,xlen);
371 }
372 /*-------------------------------------------------------------------*/
373 #endif /* NAUSKOMM */
374 /*-----------------------------------------------------------------*/
375 /*
376 ** Hypothesis: (pp,plen) an odd prime, (aa, alen) a QR mod (pp,plen)
377 ** Function calculates a square root (zz,zlen) of (aa,alen)
378 ** (aa, alen) is reduced mod (pp,plen)
379 ** space in zz must be at least 2*plen
380 ** Return value is zlen
381 ** In case of error, -1 is returned
382 */
fpSqrt(pp,plen,aa,alen,zz,hilf)383 PUBLIC int fpSqrt(pp,plen,aa,alen,zz,hilf)
384 word2 *pp, *aa, *zz, *hilf;
385 int plen, alen;
386 {
387 word2 *ex, *uu, *vv;
388 int m8, exlen, zlen, ulen, vlen;
389
390 if((plen < 1) || (pp[0] & 1) != 1)
391 return -1;
392
393 if(cmparr(pp,plen,aa,alen) <= 0) {
394 alen = modbig(aa,alen,pp,plen,hilf);
395 }
396 if(!alen)
397 return 0;
398
399 m8 = (pp[0] & 0x7);
400 if(m8 == 5) {
401 return fpSqrt58(pp,plen,aa,alen,zz,hilf);
402 }
403 else if(m8 == 1) {
404 return fpSqrt14(pp,plen,aa,alen,zz,hilf);
405 }
406
407 /******** else p = 3 mod 4 ********/
408 ex = hilf;
409 uu = hilf + plen + 2;
410 vv = hilf + 3*plen + 4;
411 hilf += 5*plen + 6;
412
413 cpyarr(pp,plen,ex);
414 exlen = shiftarr(ex,plen,-2); /* ex = (p-3)/4 */
415 ulen = modpower(aa,alen,ex,exlen,pp,plen,uu,hilf); /* a**((p-3)/4) */
416 zlen = modmultbig(uu,ulen,aa,alen,pp,plen,zz,hilf); /* a**((p+1)/4 */
417 vlen = modmultbig(uu,ulen,zz,zlen,pp,plen,vv,hilf); /* a**((p-1)/2 */
418 if(vv[0] != 1 || vlen != 1)
419 return -1;
420 return zlen;
421 }
422 /*-----------------------------------------------------------------*/
423 /*
424 ** Hypothesis: (pp,plen) an odd prime = 5 mod 8,
425 ** (aa, alen) is a QR mod (pp,plen)
426 ** Function calculates a square root (zz,zlen) of (aa,alen)
427 ** Return value is zlen
428 ** In case of error, -1 is returned
429 */
fpSqrt58(pp,plen,aa,alen,zz,hilf)430 PRIVATE int fpSqrt58(pp,plen,aa,alen,zz,hilf)
431 word2 *pp, *aa, *zz, *hilf;
432 int plen, alen;
433 {
434 word2 *ex, *uu, *vv, *ww;
435 int zlen, exlen, ulen, vlen, wlen, cmp;
436
437 ex = hilf;
438 uu = hilf + plen + 2;
439 vv = hilf + 3*plen + 4;
440 ww = hilf + 5*plen + 6;
441 hilf += 7*plen + 8;
442
443 cpyarr(pp,plen,ex);
444 exlen = shiftarr(ex,plen,-3); /* ex = (p - 5)/8 */
445 vlen = modpower(aa,alen,ex,exlen,pp,plen,vv,hilf); /* v = a**((p-5)/8) */
446 ulen = modmultbig(vv,vlen,aa,alen,pp,plen,uu,hilf); /* u = a**((p+3)/8) */
447 wlen = modmultbig(uu,ulen,vv,vlen,pp,plen,ww,hilf); /* w = a**((p-1)/4) */
448 if(ww[0] == 1 && wlen == 1) {
449 cpyarr(uu,ulen,zz);
450 zlen = ulen;
451 }
452 else {
453 wlen = incarr(ww,wlen,1);
454 cmp = cmparr(ww,wlen,pp,plen);
455 if(cmp != 0)
456 return -1;
457 exlen = shiftarr(ex,exlen,1);
458 exlen = incarr(ex,exlen,1); /* ex = (p - 1)/4 */
459 ww[0] = 2; wlen = 1;
460 vlen = modpower(ww,wlen,ex,exlen,pp,plen,vv,hilf);
461 zlen = modmultbig(uu,ulen,vv,vlen,pp,plen,ww,hilf);
462 cpyarr(ww,zlen,zz);
463 }
464 if((zz[0] & 1) == 1)
465 zlen = sub1arr(zz,zlen,pp,plen);
466 return zlen;
467 }
468 /*-----------------------------------------------------------------*/
469 /*
470 ** Hypothesis: (pp,plen) an odd prime = 1 mod 4,
471 ** (aa, alen) is a QR mod (pp,plen)
472 ** Function calculates a square root (zz,zlen) of (aa,alen)
473 ** Return value is zlen
474 ** In case of error, -1 is returned
475 */
fpSqrt14(pp,plen,aa,alen,zz,hilf)476 PRIVATE int fpSqrt14(pp,plen,aa,alen,zz,hilf)
477 word2 *pp, *aa, *zz, *hilf;
478 int plen, alen;
479 {
480 long c;
481 word2 *ex, *D, *xx, *yy;
482 int exlen, dlen, xlen, zlen, cmp;
483 FP2D fp2D;
484 PAIRXY Z;
485
486 c = nonresdisc(pp,plen,aa,alen,hilf);
487 if(c < 0)
488 return -1;
489 xx = zz;
490 yy = hilf;
491 D = hilf + 2*plen + 2;
492 ex = hilf + 4*plen + 4;
493 hilf += 5*plen + 6;
494
495 xlen = long2big(c,xx);
496 Z.xx = xx;
497 Z.xlen = xlen;
498 yy[0] = 1;
499 Z.yy = yy;
500 Z.ylen = 1;
501
502 dlen = multbig(xx,xlen,xx,xlen,D,hilf); /* c**2 */
503 cmp = cmparr(D,dlen,aa,alen);
504 if(cmp < 0) {
505 dlen = sub1arr(D,dlen,aa,alen);
506 dlen = modnegbig(D,dlen,pp,plen,hilf);
507 }
508 else {
509 dlen = subarr(D,dlen,aa,alen);
510 dlen = modbig(D,dlen,pp,plen,hilf);
511 }
512 fp2D.pp = pp;
513 fp2D.plen = plen;
514 fp2D.D = D;
515 fp2D.dlen = dlen;
516 cpyarr(pp,plen,ex);
517 exlen = incarr(ex,plen,1);
518 exlen = shiftarr(ex,exlen,-1);
519 fp2Dpower(&fp2D,&Z,ex,exlen,hilf);
520 if(Z.ylen != 0)
521 return -1;
522 else {
523 zlen = Z.xlen;
524 if((zz[0] & 1) == 1)
525 zlen = sub1arr(zz,zlen,pp,plen);
526 return zlen;
527 }
528 }
529 /*---------------------------------------------------------------*/
530 /*
531 ** returns square root of a mod p, where p is a prime
532 ** Hypothesis: jac(a,p) = 1.
533 ** p and a should be 16-bit numbers.
534 */
fp_sqrt(p,a)535 PUBLIC unsigned fp_sqrt(p,a)
536 unsigned p,a;
537 {
538 if((p & 3) == 3)
539 return(modpow(a,(p+1)/4,p));
540 else
541 return(fp_sqrt14(p,a));
542 }
543 /*---------------------------------------------------------------*/
544 /*
545 ** returns square root of a mod p, where p is a prime = 1 mod 4.
546 ** Hypothesis: jac(a,p) = 1.
547 ** p and a should be 16-bit numbers.
548 */
fp_sqrt14(p,a)549 PRIVATE unsigned fp_sqrt14(p,a)
550 unsigned p,a;
551 {
552 word4 c;
553 unsigned D, u;
554 unsigned uu[2];
555
556 a = a % p;
557 a = p - a;
558 for(c=1; c < p; c++) {
559 D = (unsigned)((c*c + a) % p);
560 if(jac(D,p) == -1)
561 break;
562 }
563 uu[0] = (unsigned)c;
564 uu[1] = 1;
565 fp2pow(p,D,uu,(p+1)/2);
566 u = uu[0];
567 if(u & 1)
568 u = p-u;
569 return u;
570 }
571 /*---------------------------------------------------------------*/
572 /*
573 ** calculates uu**n in the field Fp(sqrt(D))
574 ** Hypothesis: jac(D,p) = -1,
575 ** p 16-bit prime
576 ** (uu[0],uu[1]) is destructively replaced by result
577 */
fp2pow(p,D,uu,n)578 PRIVATE void fp2pow(p,D,uu,n)
579 unsigned p,D;
580 unsigned *uu;
581 unsigned n;
582 {
583 word4 x,x0,y,y0,X,Y;
584
585 if(n == 0) {
586 uu[0] = 1;
587 uu[1] = 0;
588 return;
589 }
590 x = uu[0]; y = uu[1];
591 X = 1; Y = 0;
592 while(n > 1) {
593 if(n & 1) {
594 x0 = X;
595 y0 = (Y*y) % p;
596 /*
597 ** X = (X*x + D*y0) % p;
598 ** Y = (x0*y + Y*x) % p;
599 */
600 X *= x; X %= p;
601 X += D*y0; X %= p;
602 Y *= x; X %= p;
603 Y += x0*y; Y %= p;
604 }
605 x0 = x;
606 y0 = (y*y) % p;
607 /*
608 ** x = (x*x + D*y0) % p;
609 ** y = (2*x0*y) % p;
610 */
611 x *= x; x %= p;
612 x += D*y0; x %= p;
613 y *= x0; y %= p;
614 y += y; y %= p;
615 n >>= 1;
616 }
617 x0 = X;
618 y0 = (Y*y) % p;
619 /*
620 ** uu[0] = (X*x + D*y0) % p;
621 ** uu[1] = (X*y + Y*x) % p;
622 */
623 X *= x; X %= p;
624 X += D*y0;
625 uu[0] = X % p;
626 Y *= x;
627 Y += x0*y;
628 uu[1] = Y % p;
629 }
630 /*-----------------------------------------------------------------*/
631 /*
632 ** Calculates a square root of (aa,alen) mod (pp,plen)**2
633 ** Hypothesis: pp odd prime, jacobi(aq,pq) = 1
634 */
635 /*
636 z := fp_sqrt(p,a);
637 xi := (z*z - a) div p;
638 eta := mod_inverse(2*z,p);
639 delta := xi*eta mod p;
640 z := z - delta*p;
641 return (z mod p**2);
642 */
fp2Sqrt(pp,plen,aa,alen,zz,hilf)643 PUBLIC int fp2Sqrt(pp,plen,aa,alen,zz,hilf)
644 word2 *pp, *aa, *zz, *hilf;
645 int plen, alen;
646 {
647 word2 *xi, *eta, *delta, *ww;
648 int m, xilen, elen, dlen, zlen, rlen, wlen, cmp, sign;
649
650 m = 2*plen + 2;
651 xi = hilf;
652 eta = xi + m;
653 delta = eta + m;
654 ww = delta + m;
655 hilf = ww + (alen >= m ? alen + 2: m);
656
657 cpyarr(aa,alen,ww);
658 wlen = modbig(ww,alen,pp,plen,hilf);
659 zlen = fpSqrt(pp,plen,ww,wlen,zz,hilf);
660 if(zlen < 0) /* error */
661 return zlen;
662 xilen = multbig(zz,zlen,zz,zlen,xi,hilf);
663 cmp = cmparr(xi,xilen,aa,alen);
664 if(cmp >= 0) {
665 xilen = subarr(xi,xilen,aa,alen);
666 sign = 0;
667 }
668 else {
669 xilen = sub1arr(xi,xilen,aa,alen);
670 sign = MINUSBYTE;
671 }
672 xilen = divbig(xi,xilen,pp,plen,ww,&rlen,hilf);
673 if(sign) {
674 xilen = modnegbig(ww,xilen,pp,plen,hilf);
675 }
676 cpyarr(ww,xilen,xi);
677 /* xi := (z*z - a) div p */
678 cpyarr(zz,zlen,ww);
679 elen = shiftarr(ww,zlen,1);
680 elen = modinverse(ww,elen,pp,plen,eta,hilf);
681 /* eta = mod_inverse(2*z,p) */
682 dlen = modmultbig(xi,xilen,eta,elen,pp,plen,delta,hilf);
683 /* delta := xi*eta mod p */
684 wlen = multbig(delta,dlen,pp,plen,ww,hilf);
685 cmp = cmparr(zz,zlen,ww,wlen);
686 if(cmp >= 0) {
687 zlen = subarr(zz,zlen,ww,wlen);
688 sign = 0;
689 }
690 else {
691 zlen = sub1arr(zz,zlen,ww,wlen);
692 sign = MINUSBYTE;
693 }
694 /* z := z - delta*p, sign! */
695 wlen = multbig(pp,plen,pp,plen,ww,hilf);
696 if(sign == 0)
697 zlen = modbig(zz,zlen,ww,wlen,hilf);
698 else {
699 zlen = modnegbig(zz,zlen,ww,wlen,hilf);
700 }
701 return zlen;
702 }
703 /*-----------------------------------------------------------------*/
704 /*
705 ** returns a number c such that jacobi(c*c - (aa,alen),(pp,plen)) = -1
706 ** In case of error (possibly (pp,plen) a square) returns -1
707 */
nonresdisc(pp,plen,aa,alen,hilf)708 PRIVATE long nonresdisc(pp,plen,aa,alen,hilf)
709 word2 *pp, *aa, *hilf;
710 int plen, alen;
711 {
712 word4 c,v;
713 word2 *xx, *yy;
714 int k, vlen, xlen, cmp, sign, res;
715 unsigned u;
716
717 if(alen > plen)
718 alen = modbig(aa,alen,pp,plen,hilf);
719 if(!alen)
720 return -1;
721
722 xx = hilf;
723 hilf += (plen >= alen ? plen : alen) + 2;
724 yy = hilf;
725 hilf += plen + 2;
726
727 u = 0x1000;
728 if(plen == 1 && pp[0] < u) {
729 u = pp[0];
730 }
731 c = random2(u);
732 for(k=1; k<=60000; k++,c++) {
733 v = c*c;
734 vlen = long2big(v,xx);
735 cmp = cmparr(aa,alen,xx,vlen);
736 if(cmp > 0) {
737 xlen = sub1arr(xx,vlen,aa,alen);
738 sign = MINUSBYTE;
739 }
740 else if(cmp < 0) {
741 xlen = subarr(xx,vlen,aa,alen);
742 sign = 0;
743 }
744 else
745 continue;
746 cpyarr(pp,plen,yy);
747 res = jacobi(sign,xx,xlen,yy,plen,hilf);
748 if(res < 0)
749 return c;
750 if((k & 0xFF) == 0) {
751 if(!rabtest(pp,plen,hilf))
752 break;
753 }
754 }
755 return -1;
756 }
757 /*-----------------------------------------------------------------*/
758 /*
759 ** Destructively calculates *pZ1 := (*pZ1) * (*pZ2)
760 ** in the field given by *pfp2D
761 */
fp2Dmult(pfp2D,pZ1,pZ2,hilf)762 PRIVATE void fp2Dmult(pfp2D,pZ1,pZ2,hilf)
763 FP2D *pfp2D;
764 PAIRXY *pZ1, *pZ2;
765 word2 *hilf;
766 {
767 word2 *x0, *zz, *ww;
768 int zlen, plen, wlen, x0len;
769
770 plen = pfp2D->plen;
771 x0 = hilf;
772 zz = hilf + plen + 2;
773 ww = hilf + 3*plen + 4;
774 hilf += 5*plen + 6;
775
776 /* x0 := x1 */
777 x0len = pZ1->xlen;
778 cpyarr(pZ1->xx,x0len,x0);
779 /* x1 := x1*x2 + y1*y2*D */
780 zlen = multbig(pZ1->yy,pZ1->ylen,pZ2->yy,pZ2->ylen,zz,hilf);
781 zlen = modbig(zz,zlen,pfp2D->pp,plen,hilf);
782 wlen = multbig(zz,zlen,pfp2D->D,pfp2D->dlen,ww,hilf);
783 zlen = multbig(pZ1->xx,pZ1->xlen,pZ2->xx,pZ2->xlen,zz,hilf);
784 wlen = addarr(ww,wlen,zz,zlen);
785 wlen = modbig(ww,wlen,pfp2D->pp,plen,hilf);
786 cpyarr(ww,wlen,pZ1->xx);
787 pZ1->xlen = wlen;
788 /* y1 := x0*y2 + y1*x2 */
789 wlen = multbig(x0,x0len,pZ2->yy,pZ2->ylen,ww,hilf);
790 zlen = multbig(pZ1->yy,pZ1->ylen,pZ2->xx,pZ2->xlen,zz,hilf);
791 wlen = addarr(ww,wlen,zz,zlen);
792 wlen = modbig(ww,wlen,pfp2D->pp,plen,hilf);
793 cpyarr(ww,wlen,pZ1->yy);
794 pZ1->ylen = wlen;
795 }
796 /*-----------------------------------------------------------------*/
797 /*
798 ** Destructively calculates *pZ := (*pZ)**2
799 ** in the field given by *pfp2D
800 */
fp2Dsquare(pfp2D,pZ,hilf)801 PRIVATE void fp2Dsquare(pfp2D,pZ,hilf)
802 FP2D *pfp2D;
803 PAIRXY *pZ;
804 word2 *hilf;
805 {
806 word2 *x0, *zz, *ww;
807 int zlen, plen, wlen, x0len;
808
809 plen = pfp2D->plen;
810 x0 = hilf;
811 zz = hilf + plen + 2;
812 ww = hilf + 3*plen + 4;
813 hilf += 5*plen + 6;
814
815 /* x0 := x */
816 x0len = pZ->xlen;
817 cpyarr(pZ->xx,x0len,x0);
818 /* x := x*x + y*y*D */
819 zlen = multbig(pZ->yy,pZ->ylen,pZ->yy,pZ->ylen,zz,hilf);
820 zlen = modbig(zz,zlen,pfp2D->pp,plen,hilf);
821 wlen = multbig(zz,zlen,pfp2D->D,pfp2D->dlen,ww,hilf);
822 zlen = multbig(pZ->xx,pZ->xlen,pZ->xx,pZ->xlen,zz,hilf);
823 wlen = addarr(ww,wlen,zz,zlen);
824 wlen = modbig(ww,wlen,pfp2D->pp,plen,hilf);
825 cpyarr(ww,wlen,pZ->xx);
826 pZ->xlen = wlen;
827 /* y := 2*x0*y */
828 zlen = multbig(x0,x0len,pZ->yy,pZ->ylen,zz,hilf);
829 zlen = shiftarr(zz,zlen,1);
830 zlen = modbig(zz,zlen,pfp2D->pp,plen,hilf);
831 cpyarr(zz,zlen,pZ->yy);
832 pZ->ylen = zlen;
833 }
834 /*-----------------------------------------------------------------*/
835 /*
836 ** Destructively calculates *pZ := (*pZ)**(ex,exlen)
837 ** in the field given by *pfp2D
838 */
fp2Dpower(pfp2D,pZ,ex,exlen,hilf)839 PRIVATE void fp2Dpower(pfp2D,pZ,ex,exlen,hilf)
840 FP2D *pfp2D;
841 PAIRXY *pZ;
842 word2 *ex, *hilf;
843 int exlen;
844 {
845 PAIRXY Z0;
846 word2 *xx, *yy;
847 int plen, bitl, k;
848 int allowintr;
849
850 if(exlen == 0) {
851 pZ->xx[0] = 1;
852 pZ->xlen = 1;
853 pZ->ylen = 0;
854 return;
855 }
856
857 plen = pfp2D->plen;
858 xx = hilf;
859 yy = hilf + 2*plen + 2;
860 hilf += 4*plen + 4;
861
862 /* z0 := z */
863 cpyarr(pZ->xx,pZ->xlen,xx);
864 cpyarr(pZ->yy,pZ->ylen,yy);
865 Z0.xx = xx; Z0.yy = yy;
866 Z0.xlen = pZ->xlen;
867 Z0.ylen = pZ->ylen;
868 bitl = (exlen-1)*16 + bitlen(ex[exlen-1]);
869 allowintr = (plen >= 16 && (exlen + plen >= 256) ? 1 : 0);
870 for(k=bitl-2; k>=0; k--) {
871 fp2Dsquare(pfp2D,pZ,hilf);
872 if(testbit(ex,k))
873 fp2Dmult(pfp2D,pZ,&Z0,hilf);
874 if(allowintr && INTERRUPT) {
875 setinterrupt(0);
876 reset(err_intr);
877 }
878 }
879 return;
880 }
881 /*-----------------------------------------------------------------*/
Fpolmult()882 PRIVATE truc Fpolmult()
883 {
884 int type;
885
886 type = chkpolmultargs(polmultsym,argStkPtr-1);
887 if(type == aERROR)
888 return brkerr();
889 if(type <= fBIGNUM) {
890 return multintpols(argStkPtr-1,MULTFLAG);
891 }
892 else {
893 error(polmultsym,err_imp,voidsym);
894 return mkvect0(0);
895 }
896 }
897 /*------------------------------------------------------------------*/
FznXmult()898 PRIVATE truc FznXmult()
899 {
900 int type;
901
902 if(chkintnz(znXmultsym,argStkPtr-2) == aERROR)
903 return brkerr();
904 type = chkznXmultargs(znXmultsym,argStkPtr-1);
905 if(type == aERROR)
906 return brkerr();
907 return znXpolmult(argStkPtr-2,argStkPtr-1);
908 }
909 /*------------------------------------------------------------------*/
FznXsquare()910 PRIVATE truc FznXsquare()
911 {
912 int type;
913
914 if(chkintnz(znXsqsym,argStkPtr-1) == aERROR)
915 return brkerr();
916 type = *FLAGPTR(argStkPtr);
917 if(type != fVECTOR) {
918 error(znXsqsym,err_vect,*argStkPtr);
919 return brkerr();
920 }
921 if(chkintvec(znXsqsym,argStkPtr) == aERROR)
922 return brkerr();
923
924 return znXpolsquare(argStkPtr-1,argStkPtr);
925 }
926 /*------------------------------------------------------------------*/
FznXdivide()927 PRIVATE truc FznXdivide()
928 {
929 int type, mode;
930
931 if(chkintnz(znXdivsym,argStkPtr-2) == aERROR)
932 return brkerr();
933 type = chkznXmultargs(znXdivsym,argStkPtr-1);
934 if(type == aERROR)
935 return brkerr();
936 mode = DDIVFLAG;
937 return znXpoldiv(argStkPtr-2,argStkPtr-1,mode);
938 }
939 /*------------------------------------------------------------------*/
FznXdiv()940 PRIVATE truc FznXdiv()
941 {
942 int type, mode;
943
944 if(chkintnz(znXdivsym,argStkPtr-2) == aERROR)
945 return brkerr();
946 type = chkznXmultargs(znXdivsym,argStkPtr-1);
947 if(type == aERROR)
948 return brkerr();
949 mode = DIVFLAG;
950 return znXpoldiv(argStkPtr-2,argStkPtr-1,mode);
951 }
952 /*------------------------------------------------------------------*/
FznXmod()953 PRIVATE truc FznXmod()
954 {
955 int type, mode;
956
957 if(chkintnz(znXmodsym,argStkPtr-2) == aERROR)
958 return brkerr();
959 type = chkznXmultargs(znXdivsym,argStkPtr-1);
960 if(type == aERROR)
961 return brkerr();
962 mode = MODFLAG;
963 return znXpoldiv(argStkPtr-2,argStkPtr-1,mode);
964 }
965 /*------------------------------------------------------------------*/
FpolNmult()966 PRIVATE truc FpolNmult()
967 {
968 int type;
969
970 type = chkpolmultargs(polNmultsym,argStkPtr-2);
971 if(type == aERROR || chkintnz(polNmultsym,argStkPtr) == aERROR)
972 return brkerr();
973 return multintpols(argStkPtr-2,MODNFLAG);
974 }
975 /*------------------------------------------------------------------*/
znXpolsquare(mptr,argptr)976 PRIVATE truc znXpolsquare(mptr,argptr)
977 truc *mptr;
978 truc *argptr;
979 {
980 truc *ptr1;
981 truc *workptr, *ptr, *wptr;
982 struct vector *vecptr;
983 truc obj;
984 word2 *x, *y, *zz, *aa, *hilf;
985 int len, len1, k, i, j0, j1;
986 int n1, n2, n3, n, m, sign, sign1, sign2, sign3;
987 unsigned mlen, offshilf;
988
989 len1 = *VECLENPTR(argptr);
990 ptr1 = argptr;
991 if(len1 == 0)
992 return mkvect0(0);
993
994 /* now len1 >= 1 */
995 workptr = workStkPtr+1;
996 if(!WORKspace(len1)) {
997 error(znXsqsym,err_memev,voidsym);
998 return brkerr();
999 }
1000 mlen = aribufSize/3 - 6;
1001 ptr = VECTORPTR(ptr1);
1002 wptr = workptr;
1003 for(i=0; i<len1; i++) { /* vector stored in workptr[] */
1004 if(*FLAGPTR(ptr) == fBIGNUM && *BIGLENPTR(ptr) >= mlen)
1005 goto ovflexit;
1006 *wptr++ = *ptr++;
1007 }
1008 len = 2*len1 - 1;
1009 obj = mkvect0(len);
1010 WORKpush(obj);
1011
1012 zz = AriBuf;
1013 n3 = bigretr(mptr,zz,&sign3);
1014 if(n3 >= mlen)
1015 goto ovflexit;
1016 aa = AriBuf + ((mlen + 6) & 0xFFFE);
1017 offshilf = (scrbufSize/2) & 0xFFFE;
1018 for(k=0; k<len; k++) {
1019 hilf = AriScratch + offshilf;
1020 n = 0; sign = 0;
1021 j0 = (k < len1 ? 0 : k+1-len1);
1022 j1 = (k+1)/2;
1023 for(i=j0; i<j1; i++) {
1024 n1 = bigref(workptr+i,&x,&sign1);
1025 n2 = bigref(workptr+(k-i),&y,&sign2);
1026 m = multbig(x,n1,y,n2,AriScratch,hilf);
1027 sign1 = (sign1 == sign2 ? 0 : MINUSBYTE);
1028 n = addsarr(aa,n,sign,AriScratch,m,sign1,&sign);
1029 }
1030 n = shlarr(aa,n,1);
1031 if ((k&1) == 0) {
1032 n1 = bigref(workptr+(k/2),&x,&sign1);
1033 m = multbig(x,n1,x,n1,AriScratch,hilf);
1034 sign1 = 0;
1035 n = addsarr(aa,n,sign,AriScratch,m,sign1,&sign);
1036 }
1037 n = modbig(aa,n,zz,n3,hilf);
1038 if(n && (sign != sign3)) {
1039 n = sub1arr(aa,n,zz,n3);
1040 sign = sign3;
1041 }
1042 obj = mkint(sign,aa,n);
1043 *(VECTORPTR(workStkPtr) + k) = obj;
1044 }
1045 vecptr = VECSTRUCTPTR(workStkPtr);
1046 ptr = &(vecptr->ele0) + len - 1;
1047 while(len > 0 && *ptr-- == zero)
1048 len--;
1049 vecptr->len = len;
1050 obj = WORKretr();
1051 workStkPtr = workptr-1;
1052 return obj;
1053 ovflexit:
1054 workStkPtr = workptr-1;
1055 error(znXsqsym,err_ovfl,voidsym);
1056 return(brkerr());
1057 }
1058 /*------------------------------------------------------------------*/
1059 /*
1060 ** multiplies two integer polynomials given by argptr[0] and argptr[1]
1061 ** and mods them out modulo *mptr
1062 */
znXpolmult(mptr,argptr)1063 PRIVATE truc znXpolmult(mptr,argptr)
1064 truc *mptr;
1065 truc *argptr;
1066 {
1067 truc *ptr1, *ptr2;
1068 truc *workptr, *w2ptr, *ptr, *wptr;
1069 struct vector *vecptr;
1070 truc obj;
1071 word2 *x, *y, *zz, *aa, *hilf;
1072 int len, len1, len2, k, i, j0, j1;
1073 int n1, n2, n3, n, m, sign, sign1, sign2, sign3;
1074 unsigned mlen, offshilf;
1075
1076 len = *VECLENPTR(argptr);
1077 len2 = *VECLENPTR(argptr+1);
1078 if(len >= len2) {
1079 len1 = len;
1080 ptr1 = argptr;
1081 ptr2 = argptr + 1;
1082 }
1083 else {
1084 len1 = len2;
1085 len2 = len;
1086 ptr1 = argptr + 1;
1087 ptr2 = argptr;
1088 }
1089 /* now len1 >= len2, lenk = length of vector *ptrk */
1090 if(len2 == 0)
1091 return mkvect0(0);
1092
1093 /* now len2 >= 1 */
1094 workptr = workStkPtr+1;
1095 if(!WORKspace(len1+len2)) {
1096 error(znXmultsym,err_memev,voidsym);
1097 return brkerr();
1098 }
1099 mlen = aribufSize/3 - 6;
1100 ptr = VECTORPTR(ptr1);
1101 wptr = workptr;
1102 for(i=0; i<len1; i++) { /* first vector stored in workptr[] */
1103 if(*FLAGPTR(ptr) == fBIGNUM && *BIGLENPTR(ptr) >= mlen)
1104 goto ovflexit;
1105 *wptr++ = *ptr++;
1106 }
1107 ptr = VECTORPTR(ptr2);
1108 wptr = w2ptr = workptr + len1;
1109 for(k=0; k<len2; k++) { /* store second polynomial */
1110 if(*FLAGPTR(ptr) == fBIGNUM && *BIGLENPTR(ptr) >= mlen)
1111 goto ovflexit;
1112 *wptr++ = *ptr++;
1113 }
1114 len = len1 + len2 - 1;
1115 obj = mkvect0(len);
1116 WORKpush(obj);
1117
1118 zz = AriBuf;
1119 n3 = bigretr(mptr,zz,&sign3);
1120 if(n3 >= mlen)
1121 goto ovflexit;
1122 aa = AriBuf + ((mlen + 6) & 0xFFFE);
1123 offshilf = (scrbufSize/2) & 0xFFFE;
1124 for(k=0; k<len; k++) {
1125 hilf = AriScratch + offshilf;
1126 n = 0; sign = 0;
1127 j0 = (k < len2 ? 0 : k+1-len2);
1128 j1 = (k < len1 ? k : len1-1);
1129 for(i=j0; i<=j1; i++) {
1130 n1 = bigref(workptr+i,&x,&sign1);
1131 n2 = bigref(w2ptr+(k-i),&y,&sign2);
1132 m = multbig(x,n1,y,n2,AriScratch,hilf);
1133 sign1 = (sign1 == sign2 ? 0 : MINUSBYTE);
1134 n = addsarr(aa,n,sign,AriScratch,m,sign1,&sign);
1135 }
1136 n = modbig(aa,n,zz,n3,hilf);
1137 if(n && (sign != sign3)) {
1138 n = sub1arr(aa,n,zz,n3);
1139 sign = sign3;
1140 }
1141 obj = mkint(sign,aa,n);
1142 *(VECTORPTR(workStkPtr) + k) = obj;
1143 }
1144 vecptr = VECSTRUCTPTR(workStkPtr);
1145 ptr = &(vecptr->ele0) + len - 1;
1146 while(len > 0 && *ptr-- == zero)
1147 len--;
1148 vecptr->len = len;
1149 obj = WORKretr();
1150 workStkPtr = workptr-1;
1151 return obj;
1152 ovflexit:
1153 workStkPtr = workptr-1;
1154 error(znXmultsym,err_ovfl,voidsym);
1155 return(brkerr());
1156 }
1157 /*------------------------------------------------------------------*/
znXdivmodsymb(mode)1158 PRIVATE truc znXdivmodsymb(mode)
1159 int mode;
1160 {
1161 if (mode == MODFLAG)
1162 return znXmodsym;
1163 else if (mode == DIVFLAG)
1164 return znXdivsym;
1165 else if (mode == DDIVFLAG)
1166 return znXddivsym;
1167 else
1168 return voidsym;
1169 }
1170 /*------------------------------------------------------------------*/
1171 /*
1172 ** division of integer polynomial given by argptr[0] by argptr[1],
1173 ** modded out modulo *mptr
1174 */
znXpoldiv(mptr,argptr,mode)1175 PRIVATE truc znXpoldiv(mptr,argptr,mode)
1176 truc *mptr;
1177 truc *argptr;
1178 int mode;
1179 {
1180 truc *workptr, *w1ptr, *w2ptr, *leadptr, *ptr1, *ptr2;
1181 truc *vptr;
1182 truc obj;
1183 word2 *yy, *zz, *aa, *bb, *hilf;
1184 unsigned mlen, offsbb, offshilf;
1185 int sign, sign1, sign2, sign3, signL;
1186 int j, k, m, m1, m2, n, len1, len2, len3, len4;
1187 int leadcorr = 0;
1188 int mode1;
1189
1190 len1 = *VECLENPTR(argptr);
1191 len2 = *VECLENPTR(argptr+1);
1192 if(!len2) {
1193 error(znXdivmodsymb(mode),err_div,voidsym);
1194 return brkerr();
1195 }
1196
1197 if(len2 > len1) {
1198 if (mode == MODFLAG)
1199 return *argptr;
1200 else
1201 obj = mkvect0(0);
1202 if (mode == DIVFLAG)
1203 return obj;
1204 else { /* mode == DDIVFLAG */
1205 WORKpush(obj);
1206 obj = mkvect0(2);
1207 vptr = VECTOR(obj);
1208 vptr[0] = WORKretr();
1209 vptr[1] = *argptr;
1210 return obj;
1211 }
1212 }
1213
1214 mlen = aribufSize/3 - 6;
1215 if(*FLAGPTR(mptr) == fBIGNUM && *BIGLENPTR(mptr) >= mlen)
1216 goto ovflexit;
1217 if(*SIGNPTR(mptr) == MINUSBYTE || *mptr == zero) {
1218 error(znXdivmodsymb(mode),err_pint,*mptr);
1219 return brkerr();
1220 }
1221
1222 workptr = workStkPtr+1;
1223 if(!WORKspace(len1 + len2)) {
1224 error(znXdivmodsymb(mode),err_memev,voidsym);
1225 return brkerr();
1226 }
1227 w2ptr = workptr + len1;
1228
1229 offsbb = (mlen + 6) & 0xFFFE;
1230 offshilf = (scrbufSize/2) & 0xFFFE;
1231
1232 ptr1 = workptr;
1233 ptr2 = VECTORPTR(argptr);
1234 for(k=0; k<len1; k++) { /* store first polynomial */
1235 if(*FLAGPTR(ptr2) == fBIGNUM && *BIGLENPTR(ptr2) >= mlen)
1236 goto ovflexit;
1237 *ptr1++ = *ptr2++;
1238 }
1239 ptr1 = w2ptr;
1240 ptr2 = VECTORPTR(argptr+1);
1241 for(k=0; k<len2; k++) { /* store second polynomial */
1242 if(*FLAGPTR(ptr2) == fBIGNUM && *BIGLENPTR(ptr2) >= mlen)
1243 goto ovflexit;
1244 *ptr1++ = *ptr2++;
1245 }
1246 aa = AriBuf;
1247 bb = AriBuf + offsbb;
1248
1249 if(w2ptr[len2-1] != constone) {
1250 leadcorr = 1;
1251 leadptr = w2ptr + len2 - 1;
1252 n = bigretr(leadptr,bb,&signL);
1253 m2 = bigref(mptr,&zz,&sign2);
1254 m = modinverse(bb,n,zz,m2,aa,AriScratch);
1255 if (m==0) {
1256 error(znXdivmodsymb(mode),
1257 "leading coeff must be invertible mod n",w2ptr[len2-1]);
1258 return brkerr();
1259 }
1260 obj = mkint(signL,aa,m);
1261 *leadptr = obj;
1262 for(k=0; k<len2-1; k++) {
1263 n = bigref(leadptr,&zz,&signL);
1264 m = bigref(w2ptr+k, &yy, &sign1);
1265 n = multbig(yy,m,zz,n,bb,AriScratch);
1266 sign = (signL == sign1 ? 0 : MINUSBYTE);
1267 m2 = bigref(mptr,&zz,&sign2);
1268 n = modbig(bb,n,zz,m2,AriScratch);
1269 obj = mkint(sign,bb,n);
1270 w2ptr[k] = obj;
1271 }
1272 }
1273
1274 len3 = len1 - len2 + 1;
1275 for(k=len3-1; k>=0; k--) {
1276 n = bigretr(workptr+len2-1+k,aa,&sign);
1277 m2 = bigref(mptr,&zz,&sign2);
1278 n = modbig(aa,n,zz,m2,AriScratch);
1279 if(!n)
1280 continue;
1281 for(j=0; j<=len2-2; j++) {
1282 m = bigref(w2ptr+j,&yy,&sign1);
1283 if(!m)
1284 continue;
1285 m = multbig(aa,n,yy,m,AriScratch,AriScratch+offshilf);
1286 sign1 = (sign == sign1 ? MINUSBYTE : 0);
1287 /* sign of the negative product */
1288 m1 = bigretr(workptr+k+j, bb, &sign2);
1289 m = addsarr(bb,m1,sign2,AriScratch,m,sign1,&sign3);
1290 obj = mkint(sign3,bb,m);
1291 workptr[k+j] = obj;
1292 }
1293 }
1294 if (mode & MODFLAG) {
1295 len4 = len2 - 1;
1296 for(k=0; k<len4; k++) {
1297 n = bigretr(workptr+k,aa,&sign);
1298 m2 = bigref(mptr,&zz,&sign2);
1299 n = modbig(aa,n,zz,m2,AriScratch);
1300 if(n && (sign != sign2)) {
1301 n = sub1arr(aa,n,zz,m2);
1302 sign = sign2;
1303 }
1304 obj = mkint(sign,aa,n);
1305 workptr[k] = obj;
1306 }
1307 while(len4>0 && workptr[len4-1] == zero)
1308 len4--;
1309 obj = mkvect0(len4);
1310 ptr1 = VECTOR(obj);
1311 ptr2 = workptr;
1312 for(k=0; k<len4; k++)
1313 *ptr1++ = *ptr2++;
1314 WORKpush(obj);
1315 }
1316 if (mode & DIVFLAG) {
1317 w1ptr = workptr + len2 - 1;
1318 if (leadcorr) {
1319 m = bigretr(leadptr,aa,&signL);
1320 for(k=0; k<len3; k++) {
1321 n = bigref(w1ptr+k,&yy,&sign);
1322 if (!n)
1323 continue;
1324 n = multbig(yy,n,aa,m,bb,AriScratch);
1325 sign = (signL == sign ? 0 : MINUSBYTE);
1326 m2 = bigref(mptr,&zz,&sign2);
1327 n = modbig(bb,n,zz,m2,AriScratch);
1328 if(n && (sign != sign2)) {
1329 n = sub1arr(bb,n,zz,m2);
1330 sign = sign2;
1331 }
1332 obj = mkint(sign,bb,n);
1333 w1ptr[k] = obj;
1334 }
1335 }
1336 else {
1337 for(k=0; k<len3; k++) {
1338 n = bigretr(w1ptr+k,aa,&sign);
1339 m2 = bigref(mptr,&zz,&sign2);
1340 n = modbig(aa,n,zz,m2,AriScratch);
1341 if(n && (sign != sign2)) {
1342 n = sub1arr(aa,n,zz,m2);
1343 sign = sign2;
1344 }
1345 obj = mkint(sign,aa,n);
1346 w1ptr[k] = obj;
1347 }
1348 }
1349 while(len3>0 && w1ptr[len3-1] == zero)
1350 len3--;
1351 obj = mkvect0(len3);
1352 ptr1 = VECTOR(obj);
1353 ptr2 = w1ptr;
1354 for(k=0; k<len3; k++)
1355 *ptr1++ = *ptr2++;
1356 WORKpush(obj);
1357 }
1358 if (mode == DDIVFLAG) {
1359 obj = mkvect0(2);
1360 vptr = VECTOR(obj);
1361 vptr[0] = WORKretr();
1362 vptr[1] = WORKretr();
1363 }
1364 else {
1365 obj = WORKretr();
1366 }
1367 workStkPtr = workptr-1;
1368 return obj;
1369
1370 ovflexit:
1371 workStkPtr = workptr-1;
1372 error(znXdivmodsymb(mode),err_ovfl,voidsym);
1373 return(brkerr());
1374 }
1375 /*------------------------------------------------------------------*/
1376 /*
1377 ** multiplies two integer polynomials given by argptr[0] and argptr[1]
1378 */
multintpols(argptr,mode)1379 PRIVATE truc multintpols(argptr,mode)
1380 truc *argptr;
1381 int mode;
1382 {
1383 truc *ptr1, *ptr2;
1384 truc *workptr, *w2ptr, *ptr, *wptr;
1385 struct vector *vecptr;
1386 truc obj;
1387 word2 *x, *y, *zz, *aa, *hilf;
1388 int len, len1, len2, k, i, j0, j1;
1389 int n1, n2, n3, n, m, sign, sign1, sign2, sign3;
1390 unsigned mlen, offshilf;
1391
1392 len = *VECLENPTR(argptr);
1393 len2 = *VECLENPTR(argptr+1);
1394 if(len >= len2) {
1395 len1 = len;
1396 ptr1 = argptr;
1397 ptr2 = argptr + 1;
1398 }
1399 else {
1400 len1 = len2;
1401 len2 = len;
1402 ptr1 = argptr + 1;
1403 ptr2 = argptr;
1404 }
1405 /* now len1 >= len2, lenk = length of vector *ptrk */
1406 if(len2 == 0)
1407 return mkvect0(0);
1408
1409 /* now len2 >= 1 */
1410 workptr = workStkPtr+1;
1411 if(!WORKspace(len1+len2)) {
1412 error(polmultsym,err_memev,voidsym);
1413 return brkerr();
1414 }
1415 mlen = aribufSize/3 - 6;
1416 ptr = VECTORPTR(ptr1);
1417 wptr = workptr;
1418 for(i=0; i<len1; i++) { /* first vector stored in workptr[] */
1419 if(*FLAGPTR(ptr) == fBIGNUM && *BIGLENPTR(ptr) >= mlen)
1420 goto ovflexit;
1421 *wptr++ = *ptr++;
1422 }
1423 ptr = VECTORPTR(ptr2);
1424 wptr = w2ptr = workptr + len1;
1425 for(k=0; k<len2; k++) { /* store second polynomial */
1426 if(*FLAGPTR(ptr) == fBIGNUM && *BIGLENPTR(ptr) >= mlen)
1427 goto ovflexit;
1428 *wptr++ = *ptr++;
1429 }
1430 len = len1 + len2 - 1;
1431 obj = mkvect0(len);
1432 WORKpush(obj);
1433 if(mode & MODNFLAG) {
1434 zz = AriBuf;
1435 n3 = bigretr(argStkPtr,zz,&sign3);
1436 if(n3 >= mlen)
1437 goto ovflexit;
1438 }
1439 aa = AriBuf + ((mlen + 6) & 0xFFFE);
1440 offshilf = (scrbufSize/2) & 0xFFFE;
1441 for(k=0; k<len; k++) {
1442 hilf = AriScratch + offshilf;
1443 n = 0; sign = 0;
1444 j0 = (k < len2 ? 0 : k+1-len2);
1445 j1 = (k < len1 ? k : len1-1);
1446 for(i=j0; i<=j1; i++) {
1447 n1 = bigref(workptr+i,&x,&sign1);
1448 n2 = bigref(w2ptr+(k-i),&y,&sign2);
1449 m = multbig(x,n1,y,n2,AriScratch,hilf);
1450 sign1 = (sign1 == sign2 ? 0 : MINUSBYTE);
1451 n = addsarr(aa,n,sign,AriScratch,m,sign1,&sign);
1452 }
1453 if(mode & MODNFLAG) {
1454 n = modbig(aa,n,zz,n3,hilf);
1455 if(n && (sign != sign3)) {
1456 n = sub1arr(aa,n,zz,n3);
1457 sign = sign3;
1458 }
1459 }
1460 obj = mkint(sign,aa,n);
1461 *(VECTORPTR(workStkPtr) + k) = obj;
1462 }
1463 vecptr = VECSTRUCTPTR(workStkPtr);
1464 ptr = &(vecptr->ele0) + len - 1;
1465 while(len > 0 && *ptr-- == zero)
1466 len--;
1467 vecptr->len = len;
1468 obj = WORKretr();
1469 workStkPtr = workptr-1;
1470 return obj;
1471 ovflexit:
1472 workStkPtr = workptr-1;
1473 error(polmultsym,err_ovfl,voidsym);
1474 return(brkerr());
1475 }
1476 /*------------------------------------------------------------------*/
Fpolmod()1477 PRIVATE truc Fpolmod()
1478 {
1479 int type;
1480
1481 type = chkpoldivargs(polmodsym,argStkPtr-1);
1482 if(type == aERROR)
1483 return brkerr();
1484
1485 return modintpols(argStkPtr-1,MODFLAG);
1486 }
1487 /*------------------------------------------------------------------*/
FpolNmod()1488 PRIVATE truc FpolNmod()
1489 {
1490 int type;
1491
1492 type = chkpoldivargs(polNmodsym,argStkPtr-2);
1493 if(type == aERROR)
1494 return brkerr();
1495
1496 return modintpols(argStkPtr-2, MODFLAG | MODNFLAG);
1497 }
1498 /*------------------------------------------------------------------*/
Fpoldiv()1499 PRIVATE truc Fpoldiv()
1500 {
1501 int type;
1502
1503 type = chkpoldivargs(poldivsym,argStkPtr-1);
1504 if(type == aERROR)
1505 return brkerr();
1506
1507 return modintpols(argStkPtr-1,DIVFLAG);
1508 }
1509 /*------------------------------------------------------------------*/
FpolNdiv()1510 PRIVATE truc FpolNdiv()
1511 {
1512 int type;
1513
1514 type = chkpoldivargs(polNdivsym,argStkPtr-2);
1515 if(type == aERROR)
1516 return brkerr();
1517
1518 return modintpols(argStkPtr-2, DIVFLAG | MODNFLAG);
1519 }
1520 /*------------------------------------------------------------------*/
modintpols(argptr,mode)1521 PRIVATE truc modintpols(argptr,mode)
1522 truc *argptr;
1523 int mode;
1524 {
1525 truc *workptr, *w1ptr, *w2ptr, *ptr1, *ptr2;
1526 truc obj;
1527 word2 *yy, *zz, *aa, *bb, *hilf;
1528 unsigned mlen, offsbb, offshilf;
1529 int sign, sign1, sign2, sign3;
1530 int j, k, m, m1, m2, n, len1, len2, len3;
1531 int mode1;
1532
1533 len1 = *VECLENPTR(argptr);
1534 len2 = *VECLENPTR(argptr+1);
1535 if(!len2) {
1536 error(polmodsym,err_div,voidsym);
1537 return brkerr();
1538 }
1539 ptr2 = VECTORPTR(argptr+1);
1540 if(ptr2[len2-1] != constone) {
1541 error(polmodsym,"divisor must have leading coeff = 1",ptr2[len2-1]);
1542 return brkerr();
1543 }
1544 if(len2 > len1)
1545 return(*argptr);
1546 else if(len2 == 1)
1547 return mkvect0(0);
1548 workptr = workStkPtr+1;
1549 if(!WORKspace(len1 + len2)) {
1550 error(polmodsym,err_memev,voidsym);
1551 return brkerr();
1552 }
1553 mlen = aribufSize/3 - 6;
1554 offsbb = (mlen + 6) & 0xFFFE;
1555 offshilf = (scrbufSize/2) & 0xFFFE;
1556 w2ptr = workptr + len1;
1557 ptr1 = workptr;
1558 ptr2 = VECTORPTR(argptr);
1559 for(k=0; k<len1; k++) { /* store first polynomial */
1560 if(*FLAGPTR(ptr2) == fBIGNUM && *BIGLENPTR(ptr2) >= mlen)
1561 goto ovflexit;
1562 *ptr1++ = *ptr2++;
1563 }
1564 ptr1 = w2ptr;
1565 ptr2 = VECTORPTR(argptr+1);
1566 for(k=0; k<len2; k++) { /* store second polynomial */
1567 if(*FLAGPTR(ptr2) == fBIGNUM && *BIGLENPTR(ptr2) >= mlen)
1568 goto ovflexit;
1569 *ptr1++ = *ptr2++;
1570 }
1571 aa = AriBuf;
1572 bb = AriBuf + offsbb;
1573 len3 = len1 - len2;
1574 for(k=len3; k>=0; k--) {
1575 n = bigretr(workptr+len2+k-1,aa,&sign);
1576 if(mode & MODNFLAG) {
1577 m2 = bigref(argStkPtr,&zz,&sign2);
1578 n = modbig(aa,n,zz,m2,AriScratch);
1579 if(n && (sign != sign2)) {
1580 n = sub1arr(aa,n,zz,m2);
1581 sign = sign2;
1582 }
1583 }
1584 if(!n)
1585 continue;
1586 for(j=0; j<=len2-2; j++) {
1587 hilf = AriScratch + offshilf;
1588 m = bigref(w2ptr+j,&yy,&sign1);
1589 if(!m)
1590 continue;
1591 m = multbig(aa,n,yy,m,AriScratch,hilf);
1592 sign1 = (sign == sign1 ? MINUSBYTE : 0);
1593 /* sign of the negative product */
1594 m1 = bigretr(workptr+k+j,bb,&sign2);
1595 m = addsarr(bb,m1,sign2,AriScratch,m,sign1,&sign3);
1596 obj = mkint(sign3,bb,m);
1597 workptr[k+j] = obj;
1598 }
1599 }
1600 mode1 = (mode & DDIVFLAG);
1601 if(mode1 == MODFLAG) {
1602 if(mode & MODNFLAG) {
1603 for(k=0; k<=len2-1; k++) {
1604 n = bigretr(workptr+k,aa,&sign);
1605 m2 = bigref(argStkPtr,&zz,&sign2);
1606 n = modbig(aa,n,zz,m2,AriScratch);
1607 if(n && (sign != sign2)) {
1608 n = sub1arr(aa,n,zz,m2);
1609 sign = sign2;
1610 }
1611 obj = mkint(sign,aa,n);
1612 workptr[k] = obj;
1613 }
1614 }
1615 k = len2-2;
1616 while(k >= 0 && workptr[k] == zero) {
1617 k--; len2--;
1618 }
1619 obj = mkvect0(len2-1);
1620 ptr1 = VECTOR(obj);
1621 ptr2 = workptr;
1622 for(k=0; k<=len2-2; k++)
1623 *ptr1++ = *ptr2++;
1624 }
1625 else if(mode1 == DIVFLAG) {
1626 w1ptr = workptr + len2 - 1;
1627 if(mode & MODNFLAG) {
1628 for(k=0; k<=len3; k++) {
1629 n = bigretr(w1ptr+k,aa,&sign);
1630 m2 = bigref(argStkPtr,&zz,&sign2);
1631 n = modbig(aa,n,zz,m2,AriScratch);
1632 if(n && (sign != sign2)) {
1633 n = sub1arr(aa,n,zz,m2);
1634 sign = sign2;
1635 }
1636 obj = mkint(sign,aa,n);
1637 w1ptr[k] = obj;
1638 }
1639 }
1640 while(len3>=0 && w1ptr[len3] == zero)
1641 len3--;
1642 obj = mkvect0(len3+1);
1643 ptr1 = VECTOR(obj);
1644 ptr2 = w1ptr;
1645 for(k=0; k<=len3; k++)
1646 *ptr1++ = *ptr2++;
1647 }
1648 workStkPtr = workptr-1;
1649 return obj;
1650 ovflexit:
1651 workStkPtr = workptr-1;
1652 error(polmodsym,err_ovfl,voidsym);
1653 return(brkerr());
1654 }
1655 /*------------------------------------------------------------------*/
chkznXmultargs(sym,argptr)1656 PRIVATE int chkznXmultargs(sym,argptr)
1657 truc sym;
1658 truc *argptr;
1659 {
1660 int flg1, flg2;
1661 truc *ptr;
1662
1663 flg1 = *FLAGPTR(argptr);
1664 flg2 = *FLAGPTR(argptr+1);
1665 if(flg1 != fVECTOR || flg2 != fVECTOR) {
1666 ptr = (flg1 == fVECTOR ? argptr+1 : argptr);
1667 return error(sym,err_vect,*ptr);
1668 }
1669 flg1 = chkintvec(sym,argptr);
1670 if(flg1 != aERROR)
1671 flg2 = chkintvec(sym,argptr+1);
1672 if(flg1 == aERROR || flg2 == aERROR)
1673 return aERROR;
1674
1675 return (flg1 >= flg2 ? flg1 : flg2);
1676 }
1677 /*------------------------------------------------------------------*/
chkpolmultargs(sym,argptr)1678 PRIVATE int chkpolmultargs(sym,argptr)
1679 truc sym;
1680 truc *argptr;
1681 {
1682 int flg1, flg2;
1683 truc *ptr;
1684
1685 flg1 = *FLAGPTR(argptr);
1686 flg2 = *FLAGPTR(argptr+1);
1687 if(flg1 != fVECTOR || flg2 != fVECTOR) {
1688 ptr = (flg1 == fVECTOR ? argptr+1 : argptr);
1689 return error(sym,err_vect,*ptr);
1690 }
1691 if(sym == polmultsym) {
1692 flg1 = chknumvec(sym,argptr);
1693 if(flg1 != aERROR)
1694 flg2 = chknumvec(sym,argptr+1);
1695 }
1696 else {
1697 flg1 = chkintvec(sym,argptr);
1698 if(flg1 != aERROR)
1699 flg2 = chkintvec(sym,argptr+1);
1700 }
1701 if(flg1 == aERROR || flg2 == aERROR)
1702 return aERROR;
1703
1704 return (flg1 >= flg2 ? flg1 : flg2);
1705 }
1706 /*------------------------------------------------------------------*/
chkpoldivargs(sym,argptr)1707 PRIVATE int chkpoldivargs(sym,argptr)
1708 truc sym;
1709 truc *argptr;
1710 {
1711 int flg1, flg2;
1712 truc *ptr;
1713
1714 flg1 = *FLAGPTR(argptr);
1715 flg2 = *FLAGPTR(argptr+1);
1716 if(flg1 != fVECTOR || flg2 != fVECTOR) {
1717 ptr = (flg1 == fVECTOR ? argptr+1 : argptr);
1718 return error(sym,err_vect,*ptr);
1719 }
1720 flg1 = chkintvec(sym,argptr);
1721 if(flg1 != aERROR)
1722 flg2 = chkintvec(sym,argptr+1);
1723 if(flg1 == aERROR || flg2 == aERROR)
1724 return aERROR;
1725
1726
1727 return (flg1 >= flg2 ? flg1 : flg2);
1728 }
1729 /*******************************************************************/
1730 typedef struct {
1731 int mode;
1732 unsigned deg;
1733 word2 ftail;
1734 } GF2n_Field;
1735
1736 static GF2n_Field gf2nField = {1, 8, 0x1B};
1737 static int MaxGf2n = 4099;
1738
1739 /*-------------------------------------------------------------*/
1740 /*
1741 ** if k = sum(b_i * 2**i), then
1742 ** spreadbyte[k] = sum(b_i * 4**i).
1743 */
1744 static word2 spreadbyte[256] = {
1745 0x0000, 0x0001, 0x0004, 0x0005, 0x0010, 0x0011, 0x0014, 0x0015,
1746 0x0040, 0x0041, 0x0044, 0x0045, 0x0050, 0x0051, 0x0054, 0x0055,
1747 0x0100, 0x0101, 0x0104, 0x0105, 0x0110, 0x0111, 0x0114, 0x0115,
1748 0x0140, 0x0141, 0x0144, 0x0145, 0x0150, 0x0151, 0x0154, 0x0155,
1749 0x0400, 0x0401, 0x0404, 0x0405, 0x0410, 0x0411, 0x0414, 0x0415,
1750 0x0440, 0x0441, 0x0444, 0x0445, 0x0450, 0x0451, 0x0454, 0x0455,
1751 0x0500, 0x0501, 0x0504, 0x0505, 0x0510, 0x0511, 0x0514, 0x0515,
1752 0x0540, 0x0541, 0x0544, 0x0545, 0x0550, 0x0551, 0x0554, 0x0555,
1753 0x1000, 0x1001, 0x1004, 0x1005, 0x1010, 0x1011, 0x1014, 0x1015,
1754 0x1040, 0x1041, 0x1044, 0x1045, 0x1050, 0x1051, 0x1054, 0x1055,
1755 0x1100, 0x1101, 0x1104, 0x1105, 0x1110, 0x1111, 0x1114, 0x1115,
1756 0x1140, 0x1141, 0x1144, 0x1145, 0x1150, 0x1151, 0x1154, 0x1155,
1757 0x1400, 0x1401, 0x1404, 0x1405, 0x1410, 0x1411, 0x1414, 0x1415,
1758 0x1440, 0x1441, 0x1444, 0x1445, 0x1450, 0x1451, 0x1454, 0x1455,
1759 0x1500, 0x1501, 0x1504, 0x1505, 0x1510, 0x1511, 0x1514, 0x1515,
1760 0x1540, 0x1541, 0x1544, 0x1545, 0x1550, 0x1551, 0x1554, 0x1555,
1761 0x4000, 0x4001, 0x4004, 0x4005, 0x4010, 0x4011, 0x4014, 0x4015,
1762 0x4040, 0x4041, 0x4044, 0x4045, 0x4050, 0x4051, 0x4054, 0x4055,
1763 0x4100, 0x4101, 0x4104, 0x4105, 0x4110, 0x4111, 0x4114, 0x4115,
1764 0x4140, 0x4141, 0x4144, 0x4145, 0x4150, 0x4151, 0x4154, 0x4155,
1765 0x4400, 0x4401, 0x4404, 0x4405, 0x4410, 0x4411, 0x4414, 0x4415,
1766 0x4440, 0x4441, 0x4444, 0x4445, 0x4450, 0x4451, 0x4454, 0x4455,
1767 0x4500, 0x4501, 0x4504, 0x4505, 0x4510, 0x4511, 0x4514, 0x4515,
1768 0x4540, 0x4541, 0x4544, 0x4545, 0x4550, 0x4551, 0x4554, 0x4555,
1769 0x5000, 0x5001, 0x5004, 0x5005, 0x5010, 0x5011, 0x5014, 0x5015,
1770 0x5040, 0x5041, 0x5044, 0x5045, 0x5050, 0x5051, 0x5054, 0x5055,
1771 0x5100, 0x5101, 0x5104, 0x5105, 0x5110, 0x5111, 0x5114, 0x5115,
1772 0x5140, 0x5141, 0x5144, 0x5145, 0x5150, 0x5151, 0x5154, 0x5155,
1773 0x5400, 0x5401, 0x5404, 0x5405, 0x5410, 0x5411, 0x5414, 0x5415,
1774 0x5440, 0x5441, 0x5444, 0x5445, 0x5450, 0x5451, 0x5454, 0x5455,
1775 0x5500, 0x5501, 0x5504, 0x5505, 0x5510, 0x5511, 0x5514, 0x5515,
1776 0x5540, 0x5541, 0x5544, 0x5545, 0x5550, 0x5551, 0x5554, 0x5555};
1777 /*-------------------------------------------------------------------*/
1778 /*
1779 ** Adds two gf2n_int's in ptr[0] and ptr[1]
1780 */
addgf2ns(ptr)1781 PUBLIC truc addgf2ns(ptr)
1782 truc *ptr;
1783 {
1784 word2 *y;
1785 int n, m, deg;
1786 int sign;
1787
1788 n = bigretr(ptr,AriBuf,&sign);
1789 m = bigref(ptr+1,&y,&sign);
1790 deg = gf2nField.deg;
1791 if(deg < bit_length(AriBuf,n) || deg < bit_length(y,m)) {
1792 error(plussym,"gf2nint summand too big",voidsym);
1793 return brkerr();
1794 }
1795 n = xorbitvec(AriBuf,n,y,m);
1796 return(mkgf2n(AriBuf,n));
1797 }
1798 /*-------------------------------------------------------------------*/
1799 /*
1800 ** Multiplies two gf2n_int's in ptr[0] and ptr[1]
1801 */
multgf2ns(ptr)1802 PUBLIC truc multgf2ns(ptr)
1803 truc *ptr;
1804 {
1805 int n, m, sign, deg;
1806 word2 *x, *y;
1807
1808 n = bigref(ptr,&x,&sign);
1809 m = bigref(ptr+1,&y,&sign);
1810 deg = gf2nField.deg;
1811 if(deg < bit_length(x,n) || deg < bit_length(y,m)) {
1812 error(timessym,"gf2nint factor too big",voidsym);
1813 return brkerr();
1814 }
1815 n = gf2polmult(x,n,y,m,AriBuf);
1816 n = gf2nmod(AriBuf,n);
1817 return mkgf2n(AriBuf,n);
1818 }
1819 /*-------------------------------------------------------------------*/
1820 /*
1821 ** Divide gf2nint ptr[0] by gf2nint ptr[1]
1822 */
divgf2ns(ptr)1823 PUBLIC truc divgf2ns(ptr)
1824 truc *ptr;
1825 {
1826 word2 *x, *y, *z;
1827 int n, m, sign, deg;
1828
1829 n = bigref(ptr,&x,&sign);
1830 deg = gf2nField.deg;
1831 if(deg < bit_length(x,n)) {
1832 error(divfsym,"gf2nint argument too big",*ptr);
1833 return brkerr();
1834 }
1835 y = AriBuf;
1836 m = bigretr(ptr+1,y,&sign);
1837 if(deg < bit_length(y,m)) {
1838 error(divfsym,"gf2nint argument too big",ptr[1]);
1839 return brkerr();
1840 }
1841 z = y + m + 1;
1842 m = gf2ninverse(y,m,z,AriScratch);
1843 if(m == 0) {
1844 error(divfsym,err_div,voidsym);
1845 return brkerr();
1846 }
1847 n = gf2polmult(z,m,x,n,AriScratch);
1848 cpyarr(AriScratch,n,AriBuf);
1849 n = gf2nmod(AriBuf,n);
1850 return mkgf2n(AriBuf,n);
1851 }
1852 /*-------------------------------------------------------------------*/
1853 /*
1854 ** gf2nint in ptr[0] is raised to power ptr[1], which may
1855 ** be a positive or negative integer
1856 */
exptgf2n(ptr)1857 PUBLIC truc exptgf2n(ptr)
1858 truc *ptr;
1859 {
1860 word2 *x, *y, *z;
1861 int n, m, N, deg, sign;
1862
1863 n = bigref(ptr,&x,&sign);
1864 deg = gf2nField.deg;
1865 if(deg < bit_length(x,n)) {
1866 error(powersym,"gf2nint argument too big",*ptr);
1867 return brkerr();
1868 }
1869 m = bigref(ptr+1,&y,&sign);
1870 if(sign) {
1871 cpyarr(x,n,AriBuf);
1872 x = AriBuf;
1873 z = AriBuf + n + 1;
1874 n = gf2ninverse(AriBuf,n,z,AriScratch);
1875 if(n == 0) {
1876 error(powersym,err_div,voidsym);
1877 return brkerr();
1878 }
1879 else {
1880 cpyarr(z,n,x);
1881 z = AriBuf + n + 1;
1882 }
1883 if(m == 1 && y[0] == 1) {
1884 return mkgf2n(x,n);
1885 }
1886 }
1887 else if(m == 0) {
1888 return gf2none;
1889 }
1890 else {
1891 z = AriBuf;
1892 }
1893 N = gf2npower(x,n,y,m,z,AriScratch);
1894 return mkgf2n(z,N);
1895 }
1896 /*-------------------------------------------------------------------*/
1897 /*
1898 ** Transforms object in *argStkPtr to data type gf2nint
1899 */
Fgf2nint()1900 PRIVATE truc Fgf2nint()
1901 {
1902 word2 *x;
1903 byte *bpt;
1904 unsigned u;
1905 unsigned len;
1906 int i, n, flg, sign;
1907
1908 flg = *FLAGPTR(argStkPtr);
1909
1910 if(flg == fFIXNUM || flg == fBIGNUM || flg == fGF2NINT) {
1911 n = bigretr(argStkPtr,AriBuf,&sign);
1912 }
1913 else if(flg == fBYTESTRING) {
1914 len = *STRLENPTR(argStkPtr);
1915 if(len >= aribufSize*2 - 2) {
1916 error(gf2n_sym,err_2long,mkfixnum(len));
1917 return(brkerr());
1918 }
1919 bpt = (byte *)STRINGPTR(argStkPtr);
1920 n = len / 2;
1921 x = AriBuf;
1922 for(i=0; i<n; i++) {
1923 u = bpt[1];
1924 u <<= 8;
1925 u += bpt[0];
1926 *x++ = u;
1927 bpt += 2;
1928 }
1929 if(len & 1) {
1930 *x = *bpt;
1931 n++;
1932 }
1933 }
1934 else {
1935 error(gf2n_sym,err_int,voidsym);
1936 return brkerr();
1937 }
1938 n = gf2nmod(AriBuf,n);
1939 return mkgf2n(AriBuf,n);
1940 }
1941 /*-------------------------------------------------------------------*/
Fgf2ninit()1942 PRIVATE truc Fgf2ninit()
1943 {
1944 int n, m;
1945 unsigned u;
1946 word2 *x;
1947
1948 if(*FLAGPTR(argStkPtr) != fFIXNUM) {
1949 error(gf2ninisym,err_pfix,voidsym);
1950 return brkerr();
1951 }
1952 n = *WORD2PTR(argStkPtr);
1953 if(n < 2) {
1954 error(gf2ninisym,"degree >= 2 expected",*argStkPtr);
1955 return brkerr();
1956 }
1957 else if(n > MaxGf2n) {
1958 error(gf2ninisym,"maximal degree is",mkfixnum(MaxGf2n));
1959 return brkerr();
1960 }
1961 u = gf2polfindirr(n);
1962 if(!u) {
1963 error(gf2ninisym, "no irreducible polynomial found", voidsym);
1964 return brkerr();
1965 }
1966 gf2nField.deg = n;
1967 gf2nField.ftail = u;
1968
1969 m = n/16 + 1;
1970 x = AriBuf;
1971 setarr(x,m,0);
1972 x[0] = u;
1973 setbit(x,n);
1974
1975 return mkint(0,x,m);
1976 }
1977 /*-------------------------------------------------------------------*/
Fgf2ndegree()1978 PRIVATE truc Fgf2ndegree()
1979 {
1980 unsigned deg = gf2nField.deg;
1981
1982 return mkfixnum(deg);
1983 }
1984 /*-------------------------------------------------------------------*/
Fgf2nfieldpol()1985 PRIVATE truc Fgf2nfieldpol()
1986 {
1987 int n;
1988 unsigned deg;
1989 word2 *x;
1990
1991 deg = gf2nField.deg;
1992
1993 n = (deg+1)/16 + 1;
1994 x = AriBuf;
1995 setarr(x,n,0);
1996 x[0] = gf2nField.ftail;
1997 setbit(x,deg);
1998 return mkint(0,x,n);
1999 }
2000 /*-------------------------------------------------------------------*/
2001 #if 0
2002 PRIVATE truc Fgf2nparms()
2003 {
2004 int mode, n;
2005 unsigned deg;
2006 word2 *x;
2007 truc fdeg, fpol, vec;
2008 truc *ptr;
2009
2010 mode = gf2nField.mode;
2011 deg = gf2nField.deg;
2012
2013 if(mode == 1) {
2014 n = (deg+1)/16 + 1;
2015 x = AriBuf;
2016 setarr(x,n,0);
2017 x[0] = gf2nField.ftail;
2018 setbit(x,deg);
2019 fpol = mkint(0,x,n);
2020 WORKpush(fpol);
2021 fdeg = mkfixnum(deg);
2022 vec = mkvect0(2);
2023 ptr = VECTOR(vec);
2024 ptr[0] = fdeg;
2025 ptr[1] = WORKretr();
2026 return vec;
2027 }
2028 else {
2029 error(gf2parmsym,err_imp,voidsym);
2030 return brkerr();
2031 }
2032 }
2033 #endif
2034 /*-------------------------------------------------------------------*/
Fmaxgf2n()2035 PRIVATE truc Fmaxgf2n()
2036 {
2037 return mkfixnum(MaxGf2n);
2038 }
2039 /*-------------------------------------------------------------------*/
Fgf2ntrace()2040 PRIVATE truc Fgf2ntrace()
2041 {
2042 int n, t, sign;
2043 word2 *x;
2044
2045 if(*FLAGPTR(argStkPtr) != fGF2NINT) {
2046 error(gf2ntrsym,"gfnint expected",*argStkPtr);
2047 return brkerr();
2048 }
2049 x = AriBuf;
2050 n = bigretr(argStkPtr,x,&sign);
2051 if(n == aERROR)
2052 return brkerr();
2053 if(gf2nField.deg < bit_length(x,n)) {
2054 error(gf2ntrsym,"gf2nint argument too big",voidsym);
2055 return brkerr();
2056 }
2057 t = gf2ntrace(x,n);
2058 return mkfixnum(t);
2059 }
2060 /*------------------------------------------------------------------*/
Fgf2Xmult()2061 PRIVATE truc Fgf2Xmult()
2062 {
2063 int n,m,len,sign;
2064 word2 *x, *y;
2065
2066 n = bigref(argStkPtr-1,&x,&sign);
2067 if (n == aERROR)
2068 goto errexit;
2069 m = bigref(argStkPtr,&y,&sign);
2070 if (m == aERROR)
2071 goto errexit;
2072 if(n + m >= aribufSize) {
2073 error(gf2Xmulsym,err_ovfl,voidsym);
2074 return(brkerr());
2075 }
2076 else if(!n || !m)
2077 return(zero);
2078 len = gf2polmult(x,n,y,m,AriBuf);
2079 return mkint(0,AriBuf,len);
2080
2081 errexit:
2082 error(gf2Xmulsym,err_intt,voidsym);
2083 return brkerr();
2084 }
2085 /*------------------------------------------------------------------*/
Fgf2Xsquare()2086 PRIVATE truc Fgf2Xsquare()
2087 {
2088 int n,len,sign;
2089 word2 *x;
2090
2091 n = bigref(argStkPtr,&x,&sign);
2092 if (n == aERROR) {
2093 error(gf2Xsqsym,err_intt,voidsym);
2094 return brkerr();
2095 }
2096 if (2*n >= aribufSize) {
2097 error(gf2Xmulsym,err_ovfl,voidsym);
2098 return(brkerr());
2099 }
2100 else if (!n)
2101 return zero;
2102
2103 len = gf2polsquare(x,n,AriBuf);
2104 return mkint(0,AriBuf,len);
2105 }
2106 /*------------------------------------------------------------------*/
Fgf2Xdivide()2107 PRIVATE truc Fgf2Xdivide()
2108 {
2109 int n,m,len,sign,rlen;
2110 word2 *x, *y;
2111 truc obj, vec;
2112 truc *vptr;
2113
2114 x = AriScratch;
2115 n = bigretr(argStkPtr-1,x,&sign);
2116 if (n == aERROR)
2117 goto errexit;
2118 else if(n >= aribufSize) {
2119 error(gf2Xddivsym,err_ovfl,voidsym);
2120 return(brkerr());
2121 }
2122 m = bigref(argStkPtr,&y,&sign);
2123 if (m == aERROR)
2124 goto errexit;
2125 else if(m == 0) {
2126 error(gf2Xddivsym,err_div,voidsym);
2127 return brkerr();
2128 }
2129 len = gf2poldivide(x,n,y,m,AriBuf,&rlen);
2130 y = AriBuf + len;
2131 cpyarr(x,rlen,y);
2132
2133 obj = mkint(0,AriBuf,len);
2134 WORKpush(obj);
2135 obj = mkint(0,y,rlen);
2136 WORKpush(obj);
2137 vec = mkvect0(2);
2138 vptr = VECTOR(vec);
2139 vptr[1] = WORKretr();
2140 vptr[0] = WORKretr();
2141 return vec;
2142
2143 errexit:
2144 error(gf2Xddivsym,err_intt,voidsym);
2145 return brkerr();
2146 }
2147 /*------------------------------------------------------------------*/
Fgf2Xdiv()2148 PRIVATE truc Fgf2Xdiv()
2149 {
2150 int n,m,len,sign;
2151 word2 *x, *y;
2152
2153 x = AriScratch;
2154 n = bigretr(argStkPtr-1,x,&sign);
2155 if (n == aERROR)
2156 goto errexit;
2157 m = bigref(argStkPtr,&y,&sign);
2158 if (m == aERROR)
2159 goto errexit;
2160 if (m > n)
2161 return zero;
2162 else if(n >= aribufSize) {
2163 error(gf2Xdivsym,err_ovfl,voidsym);
2164 return(brkerr());
2165 }
2166 else if(!m) {
2167 error(gf2Xdivsym,err_div,voidsym);
2168 return brkerr();
2169 }
2170 len = gf2poldiv(x,n,y,m,AriBuf);
2171 return mkint(0,AriBuf,len);
2172
2173 errexit:
2174 error(gf2Xdivsym,err_intt,voidsym);
2175 return brkerr();
2176 }
2177 /*------------------------------------------------------------------*/
Fgf2Xmod()2178 PRIVATE truc Fgf2Xmod()
2179 {
2180 int n,m,len,sign;
2181 word2 *x, *y;
2182
2183 x = AriBuf;
2184 n = bigretr(argStkPtr-1,x,&sign);
2185 if (n == aERROR)
2186 goto errexit;
2187 m = bigref(argStkPtr,&y,&sign);
2188 if (m == aERROR)
2189 goto errexit;
2190 if(n >= aribufSize) {
2191 error(gf2Xmodsym,err_ovfl,voidsym);
2192 return(brkerr());
2193 }
2194 else if(!m) {
2195 error(gf2Xmodsym,err_div,voidsym);
2196 return brkerr();
2197 }
2198 len = gf2polmod(x,n,y,m);
2199 return mkint(0,AriBuf,len);
2200
2201 errexit:
2202 error(gf2Xmodsym,err_intt,voidsym);
2203 return brkerr();
2204 }
2205 /*------------------------------------------------------------------*/
Fgf2Xgcd()2206 PRIVATE truc Fgf2Xgcd()
2207 {
2208 int n,m,len,sign;
2209 word2 *x, *y;
2210
2211 x = AriBuf;
2212 y = AriScratch;
2213 n = bigretr(argStkPtr-1,x,&sign);
2214 if (n == aERROR)
2215 goto errexit;
2216 m = bigretr(argStkPtr,y,&sign);
2217 if (m == aERROR)
2218 goto errexit;
2219 len = gf2polgcd(x,n,y,m);
2220 return mkint(0,AriBuf,len);
2221
2222 errexit:
2223 error(gf2Xgcdsym,err_intt,voidsym);
2224 return brkerr();
2225 }
2226 /*------------------------------------------------------------------*/
Fgf2Xmodpow()2227 PRIVATE truc Fgf2Xmodpow()
2228 {
2229 word2 *x, *y, *f;
2230 int n,m,flen,len,sign;
2231 int flg;
2232 truc symb;
2233
2234 n = bigref(argStkPtr-2,&x,&sign);
2235 if (n == aERROR) {
2236 symb = argStkPtr[-2];
2237 goto errexit;
2238 }
2239 flg = *FLAGPTR(argStkPtr-1);
2240 if (flg != fFIXNUM && flg != fBIGNUM) {
2241 error(gf2Xmpowsym,err_int,argStkPtr[-1]);
2242 return brkerr();
2243 }
2244 flen = bigref(argStkPtr,&f,&sign);
2245 if (flen == aERROR) {
2246 symb = argStkPtr[0];
2247 goto errexit;
2248 }
2249 else if (flen >= aribufSize/2 - 1) {
2250 error(gf2Xmpowsym,err_ovfl,argStkPtr[0]);
2251 return brkerr();
2252 }
2253 m = bigref(argStkPtr-1,&y,&sign);
2254 if (m == aERROR) {
2255 symb = argStkPtr[-1];
2256 goto errexit;
2257 }
2258 else if (sign) {
2259 /******************************/
2260 /* TODO: negative exponent */
2261 error(gf2Xmpowsym,err_p0int,argStkPtr[-1]);
2262 return brkerr();
2263 }
2264 len = gf2polmodpow(x,n,y,m,f,flen,AriBuf,AriScratch);
2265
2266 return mkint(0,AriBuf,len);
2267
2268 errexit:
2269 error(gf2Xmpowsym,err_intt,symb);
2270 return brkerr();
2271 }
2272 /*------------------------------------------------------------------*/
2273 /*
2274 ** P := (x,n) and Q := (y,m) are considered as polynomials over GF(2)
2275 ** Calculates P mod Q; works destructively on x; returns new length
2276 */
gf2polmod(x,n,y,m)2277 PRIVATE int gf2polmod(x,n,y,m)
2278 word2 *x, *y;
2279 int n,m;
2280 {
2281 int N, M, k, s;
2282
2283 if(m == 0)
2284 return n;
2285
2286 N = bit_length(x,n) - 1;
2287 M = bit_length(y,m) - 1;
2288 if(M > N)
2289 return n;
2290 for(k=N; k>=M; --k) {
2291 if(testbit(x,k)) {
2292 s = k-M;
2293 n = bitxorshift(x,n,y,m,s);
2294 }
2295 }
2296 return n;
2297 }
2298 /*------------------------------------------------------------------*/
2299 /*
2300 ** P := (x,n) and Q := (y,m) are considered as polynomials over GF(2)
2301 ** Calculates P div Q; works destructively on x;
2302 ** Result returned in z; return value=zlen
2303 */
gf2poldiv(x,n,y,m,z)2304 PRIVATE int gf2poldiv(x,n,y,m,z)
2305 word2 *x, *y, *z;
2306 int n,m;
2307 {
2308 int N, M, k, s, zlen;
2309
2310 if(m == 0)
2311 return n;
2312
2313 N = bit_length(x,n) - 1;
2314 M = bit_length(y,m) - 1;
2315 if(M > N)
2316 return 0;
2317 zlen = (N - M)/16 + 1;
2318 setarr(z,zlen,0);
2319 for(k=N; k>=M; --k) {
2320 if(testbit(x,k)) {
2321 s = k-M;
2322 n = bitxorshift(x,n,y,m,s);
2323 setbit(z,s);
2324 }
2325 }
2326 return zlen;
2327 }
2328 /*------------------------------------------------------------------*/
2329 /*
2330 ** P := (x,n) and Q := (y,m) are considered as polynomials over GF(2)
2331 ** Calculates P div Q; works destructively on x;
2332 ** Quotient returned in z; return value=zlen
2333 ** x becomes rest, ist len is returned in *rlentpr
2334 */
gf2poldivide(x,n,y,m,z,rlenptr)2335 PRIVATE int gf2poldivide(x,n,y,m,z,rlenptr)
2336 word2 *x, *y, *z;
2337 int n,m;
2338 int *rlenptr;
2339 {
2340 int N, M, k, s, zlen;
2341
2342 if(m == 0) { /* division by zero */
2343 *rlenptr = n;
2344 return 0;
2345 }
2346 N = bit_length(x,n) - 1;
2347 M = bit_length(y,m) - 1;
2348 if(M > N) {
2349 *rlenptr = n;
2350 return 0;
2351 }
2352 zlen = (N - M)/16 + 1;
2353 setarr(z,zlen,0);
2354 for(k=N; k>=M; --k) {
2355 if(testbit(x,k)) {
2356 s = k-M;
2357 n = bitxorshift(x,n,y,m,s);
2358 setbit(z,s);
2359 }
2360 }
2361 *rlenptr = n;
2362 return zlen;
2363 }
2364 /*------------------------------------------------------------------*/
gf2ntrace(x,n)2365 PRIVATE int gf2ntrace(x,n)
2366 word2 *x;
2367 int n;
2368 {
2369 int deg = gf2nField.deg;
2370 unsigned ftail = gf2nField.ftail;
2371 int k, m, t;
2372
2373 m = (deg/16) + 1;
2374 for(k=n; k<m; k++)
2375 x[k] = 0;
2376 t = (x[0] & 1);
2377 for(k=1; k<deg; k++) {
2378 n = shiftleft1(x,n);
2379 if(n > m)
2380 n = m;
2381 if(testbit(x,deg))
2382 x[0] ^= ftail;
2383 if(testbit(x,k))
2384 t++;
2385 }
2386 return (t & 1);
2387 }
2388 /*------------------------------------------------------------------*/
2389 /*
2390 ** Shift (x,n) to the left by 1 bit
2391 ** Works destructively on x
2392 */
shiftleft1(x,n)2393 PRIVATE int shiftleft1(x,n)
2394 word2 *x;
2395 int n;
2396 {
2397 int k, blen;
2398 unsigned maskhi = 0x8000, u = 1;
2399
2400 if(n == 0)
2401 return n;
2402 blen = bitlen(x[n-1]);
2403 for(k=n-1; k>=1; k--) {
2404 x[k] <<= 1;
2405 if(x[k-1] & maskhi)
2406 x[k] |= u;
2407 }
2408 x[0] <<= 1;
2409 if(blen == 15) {
2410 x[n] = u;
2411 n++;
2412 }
2413 return n;
2414 }
2415 /*------------------------------------------------------------------*/
2416 /*
2417 ** (y,m) is shifted by s >= 0 and xored with (x,n)
2418 ** If bitlength of (x,n) < bitlength of (y,m) plus s,
2419 ** higher entries of x must be 0.
2420 */
bitxorshift(x,n,y,m,s)2421 PRIVATE int bitxorshift(x,n,y,m,s)
2422 word2 *x, *y;
2423 int n,m,s;
2424 {
2425 int s0, s1, t1, k, k1;
2426 unsigned u;
2427
2428 if(!m)
2429 return n;
2430
2431 s0 = (s >> 4);
2432 s1 = (s & 0xF);
2433
2434 if(s1 == 0) {
2435 for(k=0, k1=s0-1; k<m; k++) {
2436 x[++k1] ^= y[k];
2437 }
2438 }
2439 else {
2440 t1 = 16 - s1;
2441 x[s0] ^= (y[0] << s1);
2442 for(k=1, k1=s0; k<m; k++) {
2443 x[++k1] ^= ((y[k] << s1) | (y[k-1] >> t1));
2444 }
2445 u = (y[m-1] >> t1);
2446 if(u) {
2447 x[++k1] ^= u;
2448 }
2449 }
2450 if(k1 >= n)
2451 n = k1+1;
2452 while((n > 0) && (x[n-1] == 0))
2453 n--;
2454 return n;
2455 }
2456 /*------------------------------------------------------------------*/
2457 /*
2458 ** Reduce (x,n) modulo the field polynomial given in gf2nField
2459 ** Works destructively on (x,n)
2460 */
gf2nmod(x,n)2461 PRIVATE int gf2nmod(x,n)
2462 word2 *x;
2463 int n;
2464 {
2465 int N, m, s, s0, s1;
2466 int deg = gf2nField.deg;
2467 unsigned ftail = gf2nField.ftail;
2468 unsigned mask = 0xFFFF;
2469 unsigned t1, t2;
2470
2471 N = bit_length(x,n) - 1;
2472 if(N < deg)
2473 return n;
2474
2475 m = (deg >> 4);
2476 mask >>= (16 - (deg & 0xF));
2477
2478 for(s = N-deg; s >= 0; s--) {
2479 if(testbit(x,deg+s)) {
2480 s0 = (s >> 4);
2481 s1 = (s & 0xF);
2482 t1 = (ftail << s1);
2483 x[s0] ^= t1;
2484 if(s1 && (t2 = (ftail >> (16-s1))))
2485 x[s0+1] ^= t2;
2486 }
2487 }
2488 x[m] &= mask;
2489 while(m >= 0 && (x[m] == 0))
2490 m--;
2491 return m+1;
2492 }
2493 /*------------------------------------------------------------------*/
2494 /*
2495 ** Multiplies gf2pols (x,n) and (y,m);
2496 ** does not alter (x,n) or (y,m)
2497 ** Result returned in z
2498 */
gf2polmult(x,n,y,m,z)2499 PRIVATE int gf2polmult(x,n,y,m,z)
2500 word2 *x, *y, *z;
2501 int n,m;
2502 {
2503 int k, M, n1;
2504
2505 M = bit_length(y,m) - 1;
2506 if(M < 0)
2507 return 0;
2508 cpyarr(x,n,z);
2509 n1 = shiftarr(z,n,M);
2510 for(k=M-1; k>=0; k--) {
2511 if(testbit(y,k)) {
2512 n1 = bitxorshift(z,n1,x,n,k);
2513 }
2514 }
2515 return n1;
2516 }
2517 /*------------------------------------------------------------------*/
2518 /*
2519 ** Squares the gf2pol (x,n)
2520 ** Result returned in z
2521 */
gf2polsquare(x,n,z)2522 PRIVATE int gf2polsquare(x,n,z)
2523 word2 *x, *z;
2524 int n;
2525 {
2526 int k, N;
2527 unsigned u;
2528
2529 if(n == 0)
2530 return 0;
2531 N = 0;
2532 for(k=0; k<n-1; k++) {
2533 u = x[k];
2534 z[N] = spreadbyte[u & 0x00FF];
2535 N++;
2536 z[N] = spreadbyte[(u >> 8) & 0x00FF];
2537 N++;
2538 }
2539 u = x[n-1];
2540 z[N] = spreadbyte[u & 0x00FF];
2541 N++;
2542 u = (u >> 8) & 0x00FF;
2543 if(u) {
2544 z[N] = spreadbyte[u];
2545 N++;
2546 }
2547 return N;
2548 }
2549 /*-------------------------------------------------------------------*/
2550 /*
2551 ** Calculates inverse of (x,n).
2552 ** Result in z; if (x,n) is not invertible, 0 is returned
2553 ** uu is an auxiliary array needed for intermediate calculations
2554 */
gf2ninverse(x,n,z,uu)2555 PRIVATE int gf2ninverse(x,n,z,uu)
2556 word2 *x, *z, *uu;
2557 int n;
2558 {
2559 int deg, m, zlen;
2560 word2 *y;
2561
2562 y = uu;
2563 deg = gf2nField.deg;
2564 m = (deg+1)/16 + 1;
2565 setarr(y,m,0);
2566 y[0] = gf2nField.ftail;
2567 setbit(y,deg);
2568
2569 uu = y + m + 1;
2570 n = gf2polgcdx(x,n,y,m,z,&zlen,uu);
2571 if(x[0] != 1 || n != 1) {
2572 return 0;
2573 }
2574 return zlen;
2575 }
2576 /*------------------------------------------------------------------*/
2577 /*
2578 ** gf2nint (x,n) is raised to power (y,n)
2579 ** Result in z; hilf is an auxiliary array
2580 */
gf2npower(x,n,y,m,z,hilf)2581 PRIVATE int gf2npower(x,n,y,m,z,hilf)
2582 word2 *x,*y,*z,*hilf;
2583 int n, m;
2584 {
2585 int N, k, exlen;
2586
2587 if(m == 0) {
2588 z[0] = 1;
2589 return 1;
2590 }
2591 else if(n == 0)
2592 return 0;
2593 exlen = bit_length(y,m);
2594 cpyarr(x,n,z);
2595 N = n;
2596 for(k=exlen-2; k>=0; k--) {
2597 cpyarr(z,N,hilf);
2598 N = gf2polsquare(hilf,N,z);
2599 N = gf2nmod(z,N);
2600 if(testbit(y,k)) {
2601 cpyarr(z,N,hilf);
2602 N = gf2polmult(hilf,N,x,n,z);
2603 N = gf2nmod(z,N);
2604 }
2605 }
2606 return N;
2607 }
2608 /*------------------------------------------------------------------*/
2609 /*
2610 ** gf2pol (x,n) is raised to power (y,n) modulo (f,flen)
2611 ** Result in z; hilf is an auxiliary array
2612 ** (x,n), (y,m) and (f,flen) are not altered
2613 */
gf2polmodpow(x,n,y,m,f,flen,z,hilf)2614 PRIVATE int gf2polmodpow(x,n,y,m,f,flen,z,hilf)
2615 word2 *x,*y,*z,*f,*hilf;
2616 int n, m, flen;
2617 {
2618 int N, k, exlen;
2619
2620 if(m == 0) {
2621 z[0] = 1;
2622 return 1;
2623 }
2624 else if(n == 0)
2625 return 0;
2626 exlen = bit_length(y,m);
2627 cpyarr(x,n,z);
2628 N = n;
2629 for(k=exlen-2; k>=0; k--) {
2630 cpyarr(z,N,hilf);
2631 N = gf2polsquare(hilf,N,z);
2632 N = gf2polmod(z,N,f,flen);
2633 if(testbit(y,k)) {
2634 cpyarr(z,N,hilf);
2635 N = gf2polmult(hilf,N,x,n,z);
2636 N = gf2polmod(z,N,f,flen);
2637 }
2638 }
2639 return N;
2640 }
2641 /*------------------------------------------------------------------*/
2642 /*
2643 ** Calculates greatest common divisor of gf2pols (x,n) and (y,m);
2644 ** Works destructively on x and y
2645 ** Result is returned in x
2646 */
gf2polgcd(x,n,y,m)2647 PRIVATE int gf2polgcd(x,n,y,m)
2648 word2 *x, *y;
2649 int n,m;
2650 {
2651 while(m > 0) {
2652 n = gf2polmod(x,n,y,m);
2653 if(n == 0) {
2654 cpyarr(y,m,x);
2655 n = m;
2656 break;
2657 }
2658 m = gf2polmod(y,m,x,n);
2659 }
2660 return n;
2661 }
2662 /*-----------------------------------------------------------------*/
2663 /*
2664 ** Calculates the gcd d of (x,n) and (y,m) (considered as polynomials
2665 ** over GF(2)) and calculates a coefficient lambda such that
2666 ** d = lambda*(x,n) mod (y,m).
2667 ** Works destructively on x and y.
2668 ** The gcd is stored in x, its length is the return value;
2669 ** lambda = (z, *zlenptr)
2670 ** uu is an auxiliary array needed for intermediate calculations
2671 */
gf2polgcdx(x,n,y,m,z,zlenptr,uu)2672 PRIVATE int gf2polgcdx(x,n,y,m,z,zlenptr,uu)
2673 word2 *x, *y, *z, *uu;
2674 int n,m;
2675 int *zlenptr;
2676 {
2677 int s, N, M, zlen, ulen, nn, m0;
2678 word2 *y0;
2679
2680 nn = (n >= m ? n : m);
2681 setarr(z,nn,0);
2682 setarr(uu,nn,0);
2683 z[0] = 1;
2684 zlen = 1;
2685 ulen = 0;
2686 m0 = m;
2687 y0 = uu + nn + 1;
2688 cpyarr(y,m0,y0);
2689
2690 N = bit_length(x,n);
2691 M = bit_length(y,m);
2692
2693 while(M > 0) {
2694 /* loop invariants:
2695 (x,n) = z*(x0,n0) mod (y0,m0);
2696 (y,m) = uu*(x0,n0) mod (y0,m0);
2697 where (x0,n0) and (y0,m0) are the initial
2698 values of (x,n) and (y,m)
2699 */
2700 if(N >= M) {
2701 s = N - M;
2702 n = bitxorshift(x,n,y,m,s);
2703 N = bit_length(x,n);
2704 zlen = bitxorshift(z,zlen,uu,ulen,s);
2705 if(zlen > m0)
2706 zlen = gf2polmod(z,zlen,y0,m0);
2707 }
2708 else {
2709 s = M - N;
2710 m = bitxorshift(y,m,x,n,s);
2711 M = bit_length(y,m);
2712 ulen = bitxorshift(uu,ulen,z,zlen,s);
2713 if(ulen > m0)
2714 ulen = gf2polmod(uu,ulen,y0,m0);
2715 }
2716 if(N == 0) {
2717 cpyarr(y,m,x);
2718 n = m;
2719 cpyarr(uu,ulen,z);
2720 zlen = ulen;
2721 break;
2722 }
2723 }
2724 if(bit_length(z,zlen) >= bit_length(y0,m0))
2725 zlen = gf2polmod(z,zlen,y0,m0);
2726 *zlenptr = zlen;
2727 return n;
2728 }
2729 /*-----------------------------------------------------------------*/
Fgf2Xprimtest()2730 PRIVATE truc Fgf2Xprimtest()
2731 {
2732 int n, sign, res;
2733 word2 *x;
2734
2735 n = bigref(argStkPtr,&x,&sign);
2736 if (n == aERROR) {
2737 error(gf2Xprimsym,err_intt,argStkPtr[0]);
2738 return brkerr();
2739 }
2740 res = gf2polirred1(x,n,AriBuf,AriScratch);
2741 return (res ? true : false);
2742 }
2743 /*-----------------------------------------------------------------*/
2744 /*
2745 ** Tests whether the gf2 polynomial given by (x,n) is irreducible.
2746 ** yy and zz are auxiliary arrays needed for intermediate calculations
2747 */
gf2polirred1(x,n,yy,zz)2748 PRIVATE int gf2polirred1(x,n,yy,zz)
2749 word2 *x, *yy, *zz;
2750 int n;
2751 {
2752 int m, m1, N, N0, k, mode;
2753 word2 xi;
2754 word2 *x2k;
2755
2756 if (n == 0)
2757 return 0;
2758 if ((x[0] & 1) == 0) /* polynomial divisible by X */
2759 return 0;
2760
2761 N = bit_length(x,n) - 1;
2762 N0 = 5*intsqrt(N); /* somewhat arbitrary */
2763 if (N0 < N/3) {
2764 mode = 1;
2765 }
2766 else {
2767 mode = 0;
2768 N0 = N/2;
2769 }
2770 x2k = zz + n + 1;
2771 xi = 2;
2772 m = 1;
2773 x2k[0] = xi;
2774 for(k=1; k<=N0; k++) {
2775 m = gf2polsquare(x2k,m,yy);
2776 m = gf2polmod(yy,m,x,n);
2777 if((m == 1 && yy[0] == xi) || m == 0)
2778 return 0;
2779 cpyarr(yy,m,x2k);
2780 yy[0] ^= xi;
2781 cpyarr(x,n,zz);
2782 m1 = gf2polgcd(yy,m,zz,n);
2783 if(m1 != 1 || yy[0] != 1) {
2784 return 0;
2785 }
2786 }
2787 if (mode == 0)
2788 return 1;
2789 for(k=N0+1; k<N; k++) {
2790 m = gf2polsquare(x2k,m,yy);
2791 m = gf2polmod(yy,m,x,n);
2792 if((yy[0] == xi && m == 1) || m == 0)
2793 return 0;
2794 cpyarr(yy,m,x2k);
2795 }
2796 m = gf2polsquare(x2k,m,yy);
2797 m = gf2polmod(yy,m,x,n);
2798 if (yy[0] == xi && m == 1)
2799 return 1;
2800 else
2801 return 0;
2802 }
2803 /*-----------------------------------------------------------------*/
2804 /*
2805 ** Tests whether the gf2 polynomial given by (x,n) is irreducible.
2806 ** yy and zz are auxiliary arrays needed for intermediate calculations
2807 */
gf2polirred(x,n,yy,zz)2808 PRIVATE int gf2polirred(x,n,yy,zz)
2809 word2 *x, *yy, *zz;
2810 int n;
2811 {
2812 int m, m1, N, k;
2813 word2 xi;
2814 word2 *x2k;
2815
2816 if(n == 0)
2817 return 0;
2818 N = bit_length(x,n) - 1;
2819 if(N <= 1)
2820 return N;
2821 x2k = zz + n + 1;
2822 xi = 2;
2823 m = 1;
2824 yy[0] = xi;
2825 N = N/2;
2826 for(k=1; k<=N; k++) {
2827 m = gf2polmult(yy,m,yy,m,x2k);
2828 cpyarr(x,n,zz);
2829 m = gf2polmod(x2k,m,zz,n);
2830 if(m == 1 && x2k[0] == xi)
2831 return 0;
2832 cpyarr(x2k,m,yy);
2833 if (m == 0) {
2834 x2k[0] = xi;
2835 m1 = 1;
2836 }
2837 else {
2838 x2k[0] ^= xi;
2839 m1 = m;
2840 }
2841 m1 = gf2polgcd(x2k,m1,zz,n);
2842 if(m1 != 1 || x2k[0] != 1) {
2843 return 0;
2844 }
2845 }
2846 return 1;
2847 }
2848 /*-----------------------------------------------------------------*/
gf2polfindirr(n)2849 PRIVATE unsigned gf2polfindirr(n)
2850 int n;
2851 {
2852 unsigned u;
2853 int m;
2854 word2 *x, *yy, *zz;
2855
2856 if(n < 2)
2857 return 0;
2858 m = (n+1)/16 + 1;
2859 x = AriBuf;
2860 yy = AriBuf + m + 1;
2861 zz = AriScratch;
2862
2863 for(u=3; u<0xFFF0; u+=2) {
2864 if(bitcount(u) & 1)
2865 continue;
2866 setarr(x,m,0);
2867 x[0] = u;
2868 setbit(x,n);
2869 if(gf2polirred(x,m,yy,zz))
2870 return u;
2871 }
2872 return 0;
2873 }
2874 /*******************************************************************/
2875