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

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