chiark - git - mdw - catacomb/blob - calc/gfx.cal

   1 /* -*-apcalc-*-
   2  *
   3  * Testbed for %$\gf{2}$% poltnomial arithmetic
   4  *
   5  * (c) 2000 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Object types ------------------------------------------------------*/
  29
  30 obj gf { x };
  31
  32 /*----- Static variables --------------------------------------------------*/
  33
  34 static obj gf example_gf_object;
  35
  36 /*----- Main code ---------------------------------------------------------*/
  37
  38 dummy = config("lib_debug", -1);
  39
  40 define gf(x)
  41 {
  42   local obj gf g;
  43   g.x = x;
  44   return (g);
  45 }
  46
  47 define gfint(x)
  48 {
  49   if (istype(x, example_gf_object))
  50     return (x.x);
  51   else
  52     return (x);
  53 }
  54
  55 define gf_add(x, y) = gf(xor(gfint(x), gfint(y)));
  56 define gf_sub(x, y) = gf(xor(gfint(x), gfint(y)));
  57 define gf_neg(x) = x;
  58
  59 define gf_mul(x, y)
  60 {
  61   local a = gfint(x), b = gfint(y), z = 0, i, bits = highbit(a);
  62   for (i = 0; i <= bits; i++) {
  63     if (bit(a, i))
  64       z = xor(z, b << i);
  65   }
  66   return gf(z);
  67 }
  68
  69 define gfx_div(rx, dx)
  70 {
  71   local r = gfint(rx), d = gfint(dx), i;
  72   local q = 0, dbits, rbits;
  73   dbits = highbit(d);
  74   rbits = highbit(r);
  75   for (i = rbits - dbits; i >= 0; i--) {
  76     if (bit(r, i + dbits)) {
  77       r = xor(r, d << i);
  78       q |= (1 << i);
  79     }
  80   }
  81   return list(q, r);
  82 }
  83
  84 define gf_div(x, y)
  85 {
  86   local l;
  87   l = gfx_div(x, y);
  88   return gf(l[[0]]);
  89 }
  90
  91 define gf_mod(x, y)
  92 {
  93   local l;
  94   l = gfx_div(x, y);
  95   return gf(l[[1]]);
  96 }
  97
  98 define gf_gcd(a, b)
  99 {
 100   local swap = 0;
 101   local g, x = 1, X = 0, y = 0, Y = 1, q, r, t;
 102   if (a.x < b.x) {
 103     t = a; a = b; b = t;
 104     swap = 1;
 105   }
 106   if (b == gf(0))
 107     g = a;
 108   else {
 109     while (b != gf(0)) {
 110       q = gf_div(a, b); r = gf_mod(a, b);
 111       t = X * q + x; x = X; X = t;
 112       t = Y * q + y; y = Y; Y = t;
 113       a = b; b = r;
 114     }
 115     g = a;
 116   }
 117   if (swap) {
 118     t = x; x = y; y = t;
 119   }
 120   return list(g, x, y);
 121 }
 122
 123 define gf_inv(a, b)
 124 {
 125   local l = gf_gcd(b, a);
 126   if (l[[0]] != gf(1)) quit "not coprime in gf_inv";
 127   return l[[2]];
 128 }
 129
 130 /*----- That's all, folks -------------------------------------------------*/