src/poly.c

   1 /*
   2  * transsip - the telephony toolkit
   3  * By Daniel Borkmann <daniel@transsip.org>
   4  * Copyright 2012 Daniel Borkmann <dborkma@tik.ee.ethz.ch>
   5  * Subject to the GPL, version 2.
   6  * Based on Bhaskar Biswas and Nicolas Sendrier McEliece
   7  * implementation, LGPL 2.1.
   8  */
   9
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12 #include <string.h>
  13
  14 #include "gf.h"
  15 #include "poly.h"
  16 #include "xmalloc.h"
  17
  18 struct polynomial *poly_alloc(int d)
  19 {
  20         struct polynomial *p = NULL;
  21
  22         p = xmalloc(sizeof(*p));
  23         p->deg = -1;
  24         p->size = d + 1;
  25         p->coeff = xzmalloc(p->size * sizeof(gf16_t));
  26
  27         return p;
  28 }
  29
  30 struct polynomial *poly_alloc_from_string(int d, const unsigned char *str,
  31                                           size_t len)
  32 {
  33         struct polynomial *p = NULL;
  34
  35         p = xmalloc(sizeof(*p));
  36         p->deg = -1;
  37         p->size = d + 1;
  38         p->coeff = xmemdup(str, len);
  39
  40         bug_on(p->size * sizeof(gf16_t) != len);
  41
  42         return p;
  43 }
  44
  45 struct polynomial *poly_copy(struct polynomial *p)
  46 {
  47         struct polynomial *q = NULL;
  48
  49         q = xmalloc(sizeof(*q));
  50         q->deg = p->deg;
  51         q->size = p->size;
  52         q->coeff = xzmalloc(q->size * sizeof(gf16_t));
  53
  54         memcpy(q->coeff, p->coeff, p->size * sizeof(gf16_t));
  55
  56         return q;
  57 }
  58
  59 void poly_free(struct polynomial *p)
  60 {
  61         xfree(p->coeff);
  62         xfree(p);
  63 }
  64
  65 void poly_set_to_zero(struct polynomial *p)
  66 {
  67         memset(p->coeff, 0, p->size * sizeof(gf16_t));
  68         p->deg = -1;
  69 }
  70
  71 int poly_calcule_deg(struct polynomial *p)
  72 {
  73         int d = p->size - 1;
  74
  75         while ((d >= 0) && (p->coeff[d] == gf_zero()))
  76                 --d;
  77         p->deg = d;
  78
  79         return d;
  80 }
  81
  82 /* copy q in p */
  83 void poly_set(struct polynomial *p, struct polynomial *q)
  84 {
  85         int d = p->size - q->size;
  86
  87         if (d < 0) {
  88                 memcpy(p->coeff, q->coeff, p->size * sizeof(gf16_t));
  89                 poly_calcule_deg(p);
  90         } else {
  91                 memcpy(p->coeff, q->coeff, q->size * sizeof(gf16_t));
  92                 memset(p->coeff + q->size, 0, d * sizeof(gf16_t));
  93                 p->deg = q->deg;
  94         }
  95 }
  96
  97 static gf16_t poly_eval_aux(gf16_t *coeff, gf16_t a, int d)
  98 {
  99         gf16_t b;
 100
 101         b = coeff[d--];
 102         for (; d >= 0; --d) {
 103                 if (b != gf_zero())
 104                         b = gf_add(gf_mul(b, a), coeff[d]);
 105                 else
 106                         b = coeff[d];
 107         }
 108
 109         return b;
 110 }
 111
 112 struct polynomial *poly_mul(struct polynomial *p, struct polynomial *q)
 113 {
 114         int i, j, dp, dq;
 115         struct polynomial *r;
 116
 117         poly_calcule_deg(p);
 118         poly_calcule_deg(q);
 119
 120         dp = poly_deg(p);
 121         dq = poly_deg(q);
 122
 123         r = poly_alloc(dp + dq);
 124
 125         for (i = 0; i <= dp; ++i)
 126                 for (j = 0; j <= dq; ++j)
 127                         poly_addto_coeff(r, i + j, gf_mul(poly_coeff(p, i),
 128                                          poly_coeff(q, j)));
 129         poly_calcule_deg(r);
 130
 131         return r;
 132 }
 133
 134 gf16_t poly_eval(struct polynomial *p, gf16_t a)
 135 {
 136         return poly_eval_aux(p->coeff, a, poly_deg(p));
 137 }
 138
 139 /* p contain it's remainder modulo g */
 140 void poly_rem(struct polynomial *p, struct polynomial *g)
 141 {
 142         int i, j, d;
 143         gf16_t a, b;
 144
 145         d = poly_deg(p) - poly_deg(g);
 146         if (d >= 0) {
 147                 a = gf_inv(poly_tete(g));
 148
 149                 for (i = poly_deg(p); d >= 0; --i, --d) {
 150                         if (poly_coeff(p, i) != gf_zero()) {
 151                                 b = gf_mul_fast(a, poly_coeff(p, i));
 152
 153                                 for (j = 0; j < poly_deg(g); ++j) {
 154                                         poly_addto_coeff(p, j + d,
 155                                                 gf_mul_fast(b, poly_coeff(g, j)));
 156                                 }
 157
 158                                 poly_set_coeff(p, i, gf_zero());
 159                         }
 160                 }
 161
 162                 poly_set_deg(p, poly_deg(g) - 1);
 163
 164                 while ((poly_deg(p) >= 0) &&
 165                        (poly_coeff(p, poly_deg(p)) == gf_zero()))
 166                         poly_set_deg(p, poly_deg(p) - 1);
 167         }
 168 }
 169
 170 void poly_sqmod_init(struct polynomial *g, struct polynomial **sq)
 171 {
 172         int i, d;
 173
 174         d = poly_deg(g);
 175
 176         for (i = 0; i < d / 2; ++i) {
 177                 /* sq[i] = x^(2i) mod g = x^(2i) */
 178                 poly_set_to_zero(sq[i]);
 179                 poly_set_deg(sq[i], 2 * i);
 180                 poly_set_coeff(sq[i], 2 * i, gf_unit());
 181         }
 182
 183         for (; i < d; ++i) {
 184                 /* sq[i] = x^(2i) mod g = (x^2 * sq[i-1]) mod g */
 185                 memset(sq[i]->coeff, 0, 2 * sizeof(gf16_t));
 186                 memcpy(sq[i]->coeff + 2, sq[i - 1]->coeff, d * sizeof(gf16_t));
 187                 poly_set_deg(sq[i], poly_deg(sq[i - 1]) + 2);
 188                 poly_rem(sq[i], g);
 189         }
 190 }
 191
 192 /*
 193  * Modulo p square of a certain polynomial g, sq[] contains the square
 194  * Modulo g of the base canonical polynomials of degree < d, where d is
 195  * the degree of G. The table sq[] will be calculated by poly_sqmod_init
 196  */
 197 void poly_sqmod(struct polynomial *res, struct polynomial *p,
 198                 struct polynomial **sq, int d)
 199 {
 200         int i, j;
 201         gf16_t a;
 202
 203         poly_set_to_zero(res);
 204
 205         /* terms of low degree */
 206         for (i = 0; i < d / 2; ++i)
 207                 poly_set_coeff(res, i * 2, gf_square(poly_coeff(p, i)));
 208
 209         /* terms of high degree */
 210         for (; i < d; ++i) {
 211                 if (poly_coeff(p, i) != gf_zero()) {
 212                         a = gf_square(poly_coeff(p, i));
 213
 214                         for (j = 0; j < d; ++j)
 215                                 poly_addto_coeff(res, j,
 216                                         gf_mul_fast(a, poly_coeff(sq[i], j)));
 217                 }
 218         }
 219
 220         /* update degre */
 221         poly_set_deg(res, d - 1);
 222
 223         while ((poly_deg(res) >= 0) &&
 224                (poly_coeff(res, poly_deg(res)) == gf_zero()))
 225                 poly_set_deg(res, poly_deg(res) - 1);
 226 }
 227
 228 /* destructive */
 229 static struct polynomial *poly_gcd_aux(struct polynomial *p1,
 230                                        struct polynomial *p2)
 231 {
 232         if (poly_deg(p2) == -1)
 233                 return p1;
 234         else {
 235                 poly_rem(p1, p2);
 236                 return poly_gcd_aux(p2, p1);
 237         }
 238 }
 239
 240 struct polynomial *poly_gcd(struct polynomial *p1, struct polynomial *p2)
 241 {
 242         struct polynomial *a, *b, *c;
 243
 244         a = poly_copy(p1);
 245         b = poly_copy(p2);
 246
 247         if (poly_deg(a) < poly_deg(b))
 248                 c = poly_copy(poly_gcd_aux(b, a));
 249         else
 250                 c = poly_copy(poly_gcd_aux(a, b));
 251
 252         poly_free(a);
 253         poly_free(b);
 254
 255         return c;
 256 }
 257
 258 struct polynomial *poly_quo(struct polynomial *p, struct polynomial *d)
 259 {
 260         int i, j, dd, dp;
 261         gf16_t a, b;
 262         struct polynomial *quo, *rem;
 263
 264         dd = poly_calcule_deg(d);
 265         dp = poly_calcule_deg(p);
 266
 267         rem = poly_copy(p);
 268         quo = poly_alloc(dp - dd);
 269
 270         poly_set_deg(quo, dp - dd);
 271
 272         a = gf_inv(poly_coeff(d, dd));
 273
 274         for (i = dp; i >= dd; --i) {
 275                 b = gf_mul_fast(a, poly_coeff(rem, i));
 276
 277                 poly_set_coeff(quo, i - dd, b);
 278                 if (b != gf_zero()) {
 279                         poly_set_coeff(rem, i, gf_zero());
 280
 281                         for (j = i - 1; j >= i - dd; --j)
 282                                 poly_addto_coeff(rem, j,
 283                                         gf_mul_fast(b, poly_coeff(d, dd - i + j)));
 284                 }
 285         }
 286
 287         poly_free(rem);
 288
 289         return quo;
 290 }
 291
 292 /* returns the degree of the smallest factor */
 293 int poly_degppf(struct polynomial *g)
 294 {
 295         int i, d, res;
 296         struct polynomial **u, *p, *r, *s;
 297
 298         d = poly_deg(g);
 299
 300         u = malloc(d * sizeof(*u));
 301         for (i = 0; i < d; ++i)
 302                 u[i] = poly_alloc(d + 1);
 303
 304         poly_sqmod_init(g, u);
 305
 306         p = poly_alloc(d - 1);
 307         poly_set_deg(p, 1);
 308         poly_set_coeff(p, 1, gf_unit());
 309
 310         r = poly_alloc(d - 1);
 311         res = d;
 312
 313         for (i = 1; i <= (d / 2) * gf_extd(); ++i) {
 314                 poly_sqmod(r, p, u, d);
 315
 316                 /* r = x^(2^i) mod g */
 317                 if ((i % gf_extd()) == 0) {
 318                         /* so 2^i = (2^m)^j (m ext. degree) */
 319                         poly_addto_coeff(r, 1, gf_unit());
 320                         poly_calcule_deg(r); /* the degree may change */
 321
 322                         s = poly_gcd(g, r);
 323                         if (poly_deg(s) > 0) {
 324                                 poly_free(s);
 325                                 res = i / gf_extd();
 326                                 break;
 327                         }
 328
 329                         poly_free(s);
 330                         poly_addto_coeff(r, 1, gf_unit());
 331                         poly_calcule_deg(r); /* the degree may change */
 332                 }
 333
 334                 /* no need for the exchange s */
 335                 s = p;
 336                 p = r;
 337                 r = s;
 338         }
 339
 340         poly_free(p);
 341         poly_free(r);
 342
 343         for (i = 0; i < d; ++i)
 344                 poly_free(u[i]);
 345         xfree(u);
 346
 347         return res;
 348 }
 349
 350 /* we suppose deg(g) >= deg(p) */
 351 void poly_eeaux(struct polynomial **u, struct polynomial **v,
 352                 struct polynomial *p, struct polynomial *g, int t)
 353 {
 354         int i, j, dr, du, delta;
 355         gf16_t a;
 356         struct polynomial *aux, *r0, *r1, *u0, *u1;
 357
 358         /* initialisation of the local variables */
 359         /* r0 <- g, r1 <- p, u0 <- 0, u1 <- 1 */
 360         dr = poly_deg(g);
 361
 362         r0 = poly_alloc(dr);
 363         r1 = poly_alloc(dr - 1);
 364         u0 = poly_alloc(dr - 1);
 365         u1 = poly_alloc(dr - 1);
 366
 367         poly_set(r0, g);
 368         poly_set(r1, p);
 369         poly_set_to_zero(u0);
 370         poly_set_to_zero(u1);
 371         poly_set_coeff(u1, 0, gf_unit());
 372         poly_set_deg(u1, 0);
 373
 374         /* invariants:
 375          * r1 = u1 * p + v1 * g
 376          * r0 = u0 * p + v0 * g
 377          * and deg(u1) = deg(g) - deg(r0)
 378          * it stops when deg (r1) <t (deg (r0)> = t)
 379          * and therefore deg (u1) = deg (g) - deg (r0) <deg (g) - t
 380          */
 381
 382         du = 0;
 383         dr = poly_deg(r1);
 384         delta = poly_deg(r0) - dr;
 385
 386         while (dr >= t) {
 387                 for (j = delta; j >= 0; --j) {
 388                         a = gf_div(poly_coeff(r0, dr + j),
 389                                    poly_coeff(r1, dr));
 390                         if (a != gf_zero()) {
 391                                 /* u0(z) <- u0(z) + a * u1(z) * z^j */
 392                                 for (i = 0; i <= du; ++i)
 393                                         poly_addto_coeff(u0, i + j,
 394                                                 gf_mul_fast(a, poly_coeff(u1, i)));
 395
 396                                 /* r0(z) <- r0(z) + a * r1(z) * z^j */
 397                                 for (i = 0; i <= dr; ++i)
 398                                         poly_addto_coeff(r0, i + j,
 399                                                 gf_mul_fast(a, poly_coeff(r1, i)));
 400                         }
 401                 }
 402
 403                 /* exchange */
 404                 aux = r0;
 405                 r0 = r1;
 406                 r1 = aux;
 407
 408                 aux = u0;
 409                 u0 = u1;
 410                 u1 = aux;
 411
 412                 du = du + delta;
 413                 delta = 1;
 414                 while (poly_coeff(r1, dr - delta) == gf_zero())
 415                         delta++;
 416                 dr -= delta;
 417         }
 418
 419         poly_set_deg(u1, du);
 420         poly_set_deg(r1, dr);
 421
 422         /* return u1 and r1 */
 423         *u = u1;
 424         *v = r1;
 425
 426         poly_free(r0);
 427         poly_free(u0);
 428 }
 429
 430 /* the field is already defined
 431  * return a polynomial of degree t irreducible in the field */
 432 struct polynomial *poly_randgen_irred(int t, int (*u8rnd)(void))
 433 {
 434         int i;
 435         struct polynomial *g;
 436
 437         g = poly_alloc(t);
 438         poly_set_deg(g, t);
 439         poly_set_coeff(g, t, gf_unit());
 440
 441         i = 0;
 442         do {
 443                 for (i = 0; i < t; ++i)
 444                         poly_set_coeff(g, i, gf_rand(u8rnd));
 445         } while (poly_degppf(g) < t);
 446
 447         return g;
 448 }
 449
 450 /* p = p * x mod g
 451  * p of degree <= deg(g)-1 */
 452 static void poly_shiftmod(struct polynomial *p, struct polynomial *g)
 453 {
 454         int i, t;
 455         gf16_t a;
 456
 457         t = poly_deg(g);
 458         a = gf_div(p->coeff[t - 1], g->coeff[t]);
 459
 460         for (i = t - 1; i > 0; --i)
 461                 p->coeff[i] = gf_add(p->coeff[i - 1],
 462                                      gf_mul(a, g->coeff[i]));
 463         p->coeff[0] = gf_mul(a, g->coeff[0]);
 464 }
 465
 466 struct polynomial **poly_sqrtmod_init(struct polynomial *g)
 467 {
 468         int i, t;
 469         struct polynomial **sqrt, *aux, *p, *q, **sq_aux;
 470
 471         t = poly_deg(g);
 472
 473         sq_aux = xmalloc(t * sizeof(*sq_aux));
 474         for (i = 0; i < t; ++i)
 475                 sq_aux[i] = poly_alloc(t + 1);
 476         poly_sqmod_init(g, sq_aux);
 477
 478         q = poly_alloc(t - 1);
 479         p = poly_alloc(t - 1);
 480         poly_set_deg(p, 1);
 481
 482         poly_set_coeff(p, 1, gf_unit());
 483         /* q(z) = 0, p(z) = z */
 484         for (i = 0; i < t * gf_extd() - 1; ++i) {
 485                 /* q(z) <- p(z)^2 mod g(z) */
 486                 poly_sqmod(q, p, sq_aux, t);
 487                 /* q(z) <-> p(z) */
 488                 aux = q;
 489                 q = p;
 490                 p = aux;
 491         }
 492         /* p(z) = z^(2^(tm-1)) mod g(z) = sqrt(z) mod g(z) */
 493
 494         sqrt = xmalloc(t * sizeof(*sqrt));
 495         for (i = 0; i < t; ++i)
 496                 sqrt[i] = poly_alloc(t - 1);
 497
 498         poly_set(sqrt[1], p);
 499         poly_calcule_deg(sqrt[1]);
 500
 501         for(i = 3; i < t; i += 2) {
 502                 poly_set(sqrt[i], sqrt[i - 2]);
 503                 poly_shiftmod(sqrt[i], g);
 504                 poly_calcule_deg(sqrt[i]);
 505         }
 506
 507         for (i = 0; i < t; i += 2) {
 508                 poly_set_to_zero(sqrt[i]);
 509                 sqrt[i]->coeff[i / 2] = gf_unit();
 510                 sqrt[i]->deg = i / 2;
 511         }
 512
 513         for (i = 0; i < t; ++i)
 514                 poly_free(sq_aux[i]);
 515         xfree(sq_aux);
 516
 517         poly_free(p);
 518         poly_free(q);
 519
 520         return sqrt;
 521 }
 522
 523 struct polynomial **poly_syndrome_init(struct polynomial *generator,
 524                                        gf16_t *support, int n)
 525 {
 526         int i, j, t;
 527         gf16_t a;
 528         struct polynomial **F;
 529
 530         F = xmalloc(n * sizeof(*F));
 531         t = poly_deg(generator);
 532
 533         /* g(z)=g_t+g_(t-1).z^(t-1)+......+g_1.z+g_0 */
 534         /* f(z)=f_(t-1).z^(t-1)+......+f_1.z+f_0 */
 535
 536         for (j = 0; j < n; ++j) {
 537                 F[j] = poly_alloc(t - 1);
 538                 poly_set_coeff(F[j], t - 1, gf_unit());
 539
 540                 for (i = t - 2; i >= 0; --i) {
 541                         poly_set_coeff(F[j], i,
 542                                        gf_add(poly_coeff(generator, i + 1),
 543                                               gf_mul(support[j],
 544                                                      poly_coeff(F[j], i + 1))));
 545                 }
 546
 547                 a = gf_add(poly_coeff(generator, 0),
 548                            gf_mul(support[j], poly_coeff(F[j], 0)));
 549
 550                 for (i = 0; i < t; ++i) {
 551                         poly_set_coeff(F[j], i, gf_div(poly_coeff(F[j], i), a));
 552                 }
 553         }
 554
 555         return F;
 556 }