5 #include <mLib/alloc.h>
14 /*----- Polynomials over finite fields ------------------------------------*/
21 #define POLY_INIT { 0, 0, 0 }
23 #define POLY_ZEROP(p) (!(p)->len)
24 #define POLY_CONSTANTP(p) ((p)->len == 1)
25 #define POLY_DEGREE(p) ((p)->len - 1)
27 void poly_free(poly
*p
)
31 for (i
= 0; i
< p
->sz
; i
++)
33 xfree(p
->v
); p
->v
= 0;
37 void poly_ensure(field
*f
, poly
*p
, unsigned sz
)
45 n
= p
->sz
; if (!n
) n
= 2; do n
<<= 1; while (n
< sz
);
46 x
= xmalloc(n
* sizeof(mp
*));
47 for (i
= 0; i
< p
->sz
; i
++)
50 x
[i
] = MP_COPY(f
->zero
);
56 void poly_fix(field
*f
, poly
*p
, unsigned n
)
58 while (n
&& F_ZEROP(f
, p
->v
[n
- 1]))
63 #define MAX(a, b) ((a) > (b) ? (a) : (b))
65 static void putmphex(field
*f
, const char *name
, mp
*x
)
69 printf("%s = 0x", name
);
70 z
= F_OUT(f
, MP_NEW
, x
);
71 mp_writefile(z
, stdout
, 16);
76 void poly_dump(field
*f
, const char *name
, poly
*p
)
81 printf("%s = (", name
);
82 for (i
= p
->len
; i
--; ) {
84 z
= F_OUT(f
, z
, p
->v
[i
]);
85 mp_writefile(z
, stdout
, 16);
92 void poly_add(field
*f
, poly
*d
, poly
*x
, poly
*y
)
95 unsigned n
= MAX(x
->len
, y
->len
);
97 /* poly_dump(f, "add x", x); */
98 /* poly_dump(f, "add y", y); */
100 for (i
= 0; i
< n
; i
++) {
101 mp
*a
= (i
< x
->len
) ? x
->v
[i
] : f
->zero
;
102 mp
*b
= (i
< y
->len
) ? y
->v
[i
] : f
->zero
;
103 d
->v
[i
] = F_ADD(f
, d
->v
[i
], a
, b
);
106 /* poly_dump(f, "add d", d); */
109 void poly_sub(field
*f
, poly
*d
, poly
*x
, poly
*y
)
112 unsigned n
= MAX(x
->len
, y
->len
);
114 /* poly_dump(f, "sub x", x); */
115 /* poly_dump(f, "sub y", y); */
116 poly_ensure(f
, d
, n
);
117 for (i
= 0; i
< n
; i
++) {
118 mp
*a
= (i
< x
->len
) ? x
->v
[i
] : f
->zero
;
119 mp
*b
= (i
< y
->len
) ? y
->v
[i
] : f
->zero
;
120 d
->v
[i
] = F_SUB(f
, d
->v
[i
], a
, b
);
123 /* poly_dump(f, "sub d", d); */
126 static void omuladd(field
*f
, poly
*x
, poly
*y
, mp
*z
, unsigned o
)
129 unsigned n
= MAX(x
->len
, y
->len
+ o
);
132 /* poly_dump(f, "omuladd x", x); */
133 /* poly_dump(f, "omuladd y", y); */
134 /* putmphex(f, "omuladd z", z); */
135 /* printf("omuladd o = %u\n", o); */
136 poly_ensure(f
, x
, n
);
137 for (i
= o
; i
< n
; i
++) {
138 if (i
< y
->len
+ o
) {
139 mp
*a
= (i
< x
->len
) ? x
->v
[i
] : f
->zero
;
140 t
= F_MUL(f
, t
, y
->v
[i
- o
], z
);
141 x
->v
[i
] = F_ADD(f
, x
->v
[i
], a
, t
);
146 /* poly_dump(f, "omuladd d", x); */
149 static void omulsub(field
*f
, poly
*x
, poly
*y
, mp
*z
, unsigned o
)
152 unsigned n
= MAX(x
->len
, y
->len
+ o
);
155 /* poly_dump(f, "omulsub x", x); */
156 /* poly_dump(f, "omulsub y", y); */
157 /* putmphex(f, "omulsub z", z); */
158 /* printf("omulsub o = %u\n", o); */
159 poly_ensure(f
, x
, n
);
160 for (i
= o
; i
< n
; i
++) {
161 if (i
< y
->len
+ o
) {
162 mp
*a
= (i
< x
->len
) ? x
->v
[i
] : f
->zero
;
163 t
= F_MUL(f
, t
, y
->v
[i
- o
], z
);
164 x
->v
[i
] = F_SUB(f
, x
->v
[i
], a
, t
);
169 /* poly_dump(f, "omulsub d", x); */
172 void poly_mul(field
*f
, poly
*d
, poly
*x
, poly
*y
)
177 /* poly_dump(f, "mul x", x); */
178 /* poly_dump(f, "mul y", y); */
179 poly_ensure(f
, &dd
, x
->len
+ y
->len
);
180 for (i
= 0; i
< y
->len
; i
++) {
181 if (!F_ZEROP(f
, y
->v
[i
]))
182 omuladd(f
, &dd
, x
, y
->v
[i
], i
);
186 /* poly_dump(f, "mul d", d); */
189 void poly_copy(field
*f
, poly
*d
, poly
*p
)
193 poly_ensure(f
, d
, p
->len
);
194 for (i
= 0; i
< p
->len
; i
++) {
196 d
->v
[i
] = MP_COPY(p
->v
[i
]);
201 void poly_div(field
*f
, poly
*q
, poly
*r
, poly
*x
, poly
*y
)
203 poly qq
= POLY_INIT
, rr
= POLY_INIT
;
209 /* poly_dump(f, "div x", x); */
210 /* poly_dump(f, "div y", y); */
212 i
= F_INV(f
, MP_NEW
, y
->v
[y
->len
- 1]);
213 poly_copy(f
, &rr
, x
);
214 if (rr
.len
>= y
->len
) {
215 o
= oo
= rr
.len
- y
->len
+ 1;
217 if (q
) poly_ensure(f
, &qq
, o
);
220 if (!F_ZEROP(f
, rr
.v
[j
])) {
221 z
= F_MUL(f
, z
, rr
.v
[j
], i
);
222 if (q
) qq
.v
[o
] = MP_COPY(z
);
223 omulsub(f
, &rr
, y
, z
, o
);
228 poly_fix(f
, &qq
, oo
);
231 /* poly_dump(f, "div q", q); */
233 /* poly_dump(f, "div r", &rr); */
244 void poly_monic(field
*f
, poly
*d
, poly
*p
)
251 /* poly_dump(f, "monic p", p); */
252 poly_ensure(f
, d
, n
);
253 z
= F_INV(f
, MP_NEW
, p
->v
[n
- 1]);
254 for (i
= 0; i
< n
; i
++)
255 d
->v
[i
] = F_MUL(f
, d
->v
[i
], p
->v
[i
], z
);
257 /* poly_dump(f, "monic d", d); */
261 void poly_gcd(field
*f
, poly
*g
, poly
*x
, poly
*y
)
263 poly uu
= POLY_INIT
, vv
= POLY_INIT
, rr
= POLY_INIT
;
264 poly
*u
= &uu
, *v
= &vv
, *r
= &rr
;
267 /* poly_dump(f, "gcd x", x); */
268 /* poly_dump(f, "gcd y", y); */
269 poly_copy(f
, u
, x
); poly_copy(f
, v
, y
);
270 if (u
->len
< v
->len
) { t
= u
; u
= v
; v
= t
; }
271 while (!POLY_ZEROP(v
)) {
272 poly_div(f
, 0, r
, u
, v
);
273 t
= u
; u
= v
; v
= r
; r
= t
;
279 /* poly_dump(f, "gcd g", g); */
282 mp
*poly_solve(field
*f
, mp
*d
, mp
*p
, grand
*r
)
289 unsigned m
= f
->nbits
;
292 poly_ensure(f
, &g
, m
+ 1);
293 for (i
= 0; i
<= m
; i
++)
294 g
.v
[i
] = mp_testbit(p
, i
) ?
MP_COPY(f
->one
) : MP_COPY(f
->zero
);
295 poly_fix(f
, &g
, m
+ 1);
296 poly_ensure(f
, &ut
, 2);
297 while (POLY_DEGREE(&g
) > 1) {
298 ut
.v
[1] = mprand(ut
.v
[1], m
, r
, 0);
300 poly_copy(f
, &c
, &ut
);
301 /* poly_dump(f, "c-in", &c); */
303 for (i
= 1; i
< m
; i
++) {
304 poly_mul(f
, &c
, &c
, &c
);
305 poly_add(f
, &c
, &c
, &ut
);
306 poly_div(f
, 0, &c
, &c
, &g
);
307 /* putchar('.'); fflush(stdout); */
309 /* poly_dump(f, "c-out", &c); */
310 poly_gcd(f
, &h
, &c
, &g
);
311 /* poly_dump(f, "h", &h); */
312 if (POLY_CONSTANTP(&h
) || POLY_DEGREE(&h
) == POLY_DEGREE(&g
))
314 if (2 * POLY_DEGREE(&h
) > POLY_DEGREE(&g
))
315 poly_div(f
, &g
, 0, &g
, &h
);
317 poly_copy(f
, &g
, &h
);
328 /*----- Other stuff -------------------------------------------------------*/
330 static int primep(unsigned x
)
334 for (i
= 2; i
*i
< x
; i
++) {
341 static unsigned gcd(unsigned u
, unsigned v
)
345 if (u
< v
) { r
= u
; u
= v
; v
= r
; }
353 static unsigned order(unsigned x
, unsigned p
)
357 if (!x
|| x
== 1) return (0);
358 for (y
= x
, k
= 1; y
!= 1; y
= (y
*x
)%p
, k
++);
362 static int onbtype(unsigned m
)
383 #define PI 3.1415926535897932384626433832795028841971693993751
385 static mp
*fieldpoly(unsigned m
, int t
, grand
*rr
)
393 for (i
= 0; i
<= m
; i
++)
394 p
= mp_setbit(p
, p
, i
);
400 for (i
= 1; i
< m
; i
++) {
403 z
= r
; r
= q
; q
= p
; p
= z
;
410 unsigned pp
= t
*m
+ 1;
413 struct cplx
{ double r
, i
; };
417 do uu
= GR_RANGE(rr
, pp
); while (order(uu
, pp
) != t
);
418 f
= xmalloc(sizeof(struct cplx
) * (m
+ 1));
419 for (i
= 0; i
<= m
; i
++) f
[i
].r
= f
[i
].i
= 0;
421 printf("poly init; type = %u\n", t
);
422 for (i
= m
+ 1; i
--; )
423 printf("%3u: %g + %g i\n", i
, f
[i
].r
, f
[i
].i
);
425 for (i
= 1; i
<= m
; i
++) {
427 for (j
= 0; j
< t
; j
++) {
428 double z
= (pow(2, i
) * pow(uu
, j
) * PI
)/pp
;
429 e
.r
-= cos(z
); e
.i
-= sin(z
);
431 printf("new factor: %g + %g i\n", e
.r
, e
.i
);
433 f
[j
+ 1].r
+= f
[j
].r
;
434 f
[j
+ 1].i
+= f
[j
].i
;
435 n
.r
= f
[j
].r
* e
.r
- f
[j
].i
* e
.i
;
436 n
.i
= f
[j
].r
* e
.i
+ f
[j
].i
* e
.r
;
439 printf("poly after %u iters\n", i
);
440 for (j
= m
+ 1; j
--; )
441 printf("%3u: %g + %g i\n", j
, f
[j
].r
, f
[j
].i
);
454 static int dofip(unsigned m
, mp
**p
, unsigned n
, unsigned x
)
458 for (i
= n
+ 1; i
< x
; i
++) {
459 *p
= mp_setbit(*p
, *p
, i
);
460 if (n ?
dofip(m
, p
, n
- 1, i
) : gf_irreduciblep(*p
))
462 *p
= mp_clearbit(*p
, *p
, i
);
467 static mp
*fip(unsigned m
)
472 p
= mp_setbit(p
, p
, m
);
474 while (!dofip(m
, &p
, n
, m
))
479 static mp
*fnb(mp
*p
)
483 grand
*r
= fibrand_create(0);
485 unsigned m
= mp_bits(p
) - 1;
487 if ((t
= onbtype(m
)) == 0 || t
> 2)
489 f
= field_binpoly(p
);
490 q
= fieldpoly(m
, t
, r
);
491 x
= poly_solve(f
, MP_NEW
, q
, r
);
497 int main(int argc
, char *argv
[])
499 int m
= atoi(argv
[1]);
508 p
= mp_setbit(p
, p
, m
);
509 for (i
= 2; i
< argc
; i
++)
510 p
= mp_setbit(p
, p
, atoi(argv
[i
]));
512 if (!gf_irreduciblep(p
)) {
513 printf("not irreducible\n");
518 for (i
= m
+ 1; i
--; ) {
519 if (mp_testbit(p
, i
))
528 MP_PRINTX("beta", beta
);
532 assert(mparena_count(MPARENA_GLOBAL
) == 0);