
/*
 * Computer Algebra Kit (c) 1993,98 by Comp.Alg.Objects.  All Rights Reserved.
 * $Id: intmodp.m,v 1.7 1998/11/13 09:45:30 stes Exp $
 */

#include "cakit.h"

static modp_args 
modp_getargs(id self)
{
  modp_args args;

  args.p = [self modulus];
  return args;
}

static int 
modp_size(modp_args args)
{
  return 1;			/* pointer increment */
}

static int 
modp_bsize(modp_args args)
{
  return sizeof(modp_c);	/* bytes */
}

static modp_t 
modp_id2t(id object)
{
  return [object modp_value];	/* read-only */
}

static modp_t 
modp_id2ref(id object, modp_c * ref)
{
  return [object modp_value];	/* read-only, scratch space "ref" */
}

static modp_c *
modp_id2c(id object)
{
  return [object modp_reference];	/* write for inplace arithmetic */
}

static id 
modp_t2id(id proto, modp_t f)
{
  return [proto modp_value:f];	/* new object containing copy of "f" */
}

static id 
modp_c2id(id proto, modp_c * c)
{
  return [proto modp_reference:c];	/* new object containing "c" */
}

static modp_t 
modp_fromintvalue(int i, modp_args args)
{
  modp_t red;

  if (i == 0)
    return 0;
  if (i < 0) {
    red = -i;
    red %= args.p;		/* it's dangerous to % int & unsigned */
    red = (red) ? args.p - red : 0;
  } else {
    red = i;
    red %= args.p;		/* dito */
  }

  return red;
}

static int 
modp_check(modp_t self, modp_args args)
{
  assert(self < args.p);
  return 1;
}

static modp_t 
modp_zero(modp_args args)
{
  return 0;
}

static modp_t 
modp_one(modp_args args)
{
  return 1;
}

static unsigned 
modp_hash(modp_t a, modp_args args)
{
  return a;
}

static int 
modp_iseq(modp_t a, modp_t b, modp_args args)
{
  assert(modp_check(a, args) && modp_check(b, args));
  return (a == b);
}

static int 
modp_isop(modp_t a, modp_t b, modp_args args)
{
  assert(modp_check(a, args) && modp_check(b, args));
  return (a == 0) ? (b == 0) : (a + b == args.p);
}

static int 
modp_iszero(modp_t a, modp_args args)
{
  return a == 0;
}

static int 
modp_isone(modp_t a, modp_args args)
{
  return a == 1;
}

static int 
modp_isminusone(modp_t a, modp_args args)
{
  return (a + 1) == args.p;
}

static modp_t 
modp_neg(modp_t a, modp_args args)
{
  assert(modp_check(a, args));
  return (a) ? args.p - a : 0;
}

static modp_t 
modp_dbl(modp_t a, modp_args args)
{
  modp_t c = a + a;

  assert(args.p != 2 && modp_check(a, args));
  if (c)
    return (c < args.p) ? c : c - args.p;
  else
    return 0;
}

static modp_t 
modp_add(modp_t a, modp_t b, modp_args args)
{
  modp_t c = a + b;

  assert(modp_check(a, args) && modp_check(b, args));
  return ((c >= args.p) ? c - args.p : c);	/* compiles better than c<p?c:c-p */
}

static modp_t 
modp_sub(modp_t a, modp_t b, modp_args args)
{
  assert(modp_check(a, args) && modp_check(b, args));
  return (a >= b) ? (a - b) : (args.p - (b - a));
}

static modp_t 
modp_mul(modp_t a, modp_t b, modp_args args)
{
  unsigned c = a * b;		/* limitation to 16bits to be able to do this

				 * russian peasant 50% slower
				 */
  assert(sizeof(unsigned) == 2 * sizeof(modp_t));

  return c % args.p;
}

static modp_t 
modp_sqr(modp_t a, modp_args args)
{
  unsigned c = a * a;
  assert(sizeof(unsigned) == 2 * sizeof(modp_t));

  return c % args.p;
}

static modp_t 
modp_inv(modp_t x, modp_args args)
{
  unsigned q, s;
  int u1, u2, us;		/* intermediate values can be negative */
  modp_t a, b, inv;

  assert(x != 0 && x != 1 && x + 1 != args.p);

  a = x;
  b = args.p;
  u1 = 1;
  u2 = 0;
  while (b) {
    s = a % b;
    q = a / b;
    a = b;
    b = s;
    u1 -= q * u2;
    us = u1;
    u1 = u2;
    u2 = us;
  }

  inv = ((u1 < 0) ? u1 + args.p : u1);	/* keep result positive */

  assert(modp_check(inv, args));
  assert(modp_isone(modp_mul(x, inv, args), args));
  return inv;
}

static modp_t 
modp_div(modp_t a, modp_t b, modp_args args)
{
  assert(b != 0);
  return modp_mul(a, modp_inv(b, args), args);
}

static modp_t 
modp_pow(modp_t a, int n, modp_args args)
{
  modp_t o;

  if (n == 0)
    return modp_one(args);
  if (n == 1)
    return a;
  if (n == 2)
    return modp_sqr(a, args);

  o = modp_one(args);
  while (n > 1) {
    if (n & 1)
      o = modp_mul(o, a, args);
    a = modp_sqr(a, args);
    n >>= 1;
  }

  return modp_mul(o, a, args);
}

#if !modp_isvalue
#error "modp_copy_val.fab should not be included"
#endif

static void 
modp_clear(modp_c * c, modp_args args)
{
  /* nothing to clear, but the following might help... */
#ifndef DEBUG
  *c = (modp_t) 0xcafebabe;
#endif
}

static void 
modp_copy(modp_c * c, modp_t a, modp_args args)
{
  *c = a;
}

static void 
modp_deepcopy(modp_c * c, modp_t a, modp_args args)
{
  *c = a;
}

static void 
modp_move(modp_c * c, modp_c * a, modp_args args)
{
  *c = *a;
#ifndef DEBUG
  *a = (modp_t) 0xcafebabe;
#endif
}

#if !modp_isvalue
#error "modp_add_val.fab should not be included"
#endif

static void 
modp_oozero(modp_c * c, modp_args args)
{
  *c = modp_zero(args);
}

static void 
modp_inzero(modp_c * a, modp_args args)
{
  *a = modp_zero(args);
}

static void 
modp_inneg(modp_c * a, modp_args args)
{
  *a = modp_neg(*a, args);
}

static void 
modp_ooneg(modp_c * c, modp_t a, modp_args args)
{
  *c = modp_neg(a, args);
}

static modp_t 
modp_dblv(modp_t a, int v, modp_args args)
{
  modp_t c = modp_dbl(a, args);

  return (v == +1) ? c : modp_neg(c, args);
}

static void 
modp_indbl(modp_c * a, int v, modp_args args)
{
  *a = modp_dblv(*a, v, args);
}

static void 
modp_oodbl(modp_c * c, modp_t a, int v, modp_args args)
{
  *c = modp_dblv(a, v, args);
}

static modp_t 
modp_addvw(modp_t a, int v, modp_t b, int w, modp_args args)
{
  modp_t c;

  assert(v * v == +1 && w * w == 1);

  if (v == +1) {
    if (w == +1) {
      c = modp_add(a, b, args);
    } else {
      c = modp_sub(a, b, args);
    }
  } else {
    if (w == +1) {
      c = modp_sub(b, a, args);
    } else {
      c = modp_neg(modp_add(a, b, args), args);
    }
  }

  return c;
}

static void 
modp_inadd(modp_c * a, int v, modp_t b, int w, modp_args args)
{
  *a = modp_addvw(*a, v, b, w, args);
}

static void 
modp_ooadd(modp_c * c, modp_t a, int v, modp_t b, int w, modp_args args)
{
  *c = modp_addvw(a, v, b, w, args);
}

#if !modp_isvalue
#error "modp_mul_val.fab should not be included"
#endif

static void 
modp_insqr(modp_c * a, modp_args args)
{
  *a = modp_sqr(*a, args);
}

static void 
modp_oosqr(modp_c * c, modp_t a, modp_args args)
{
  *c = modp_sqr(a, args);
}

static void 
modp_inmul(modp_c * a, modp_t b, modp_args args)
{
  *a = modp_mul(*a, b, args);
}

static void 
modp_oomul(modp_c * c, modp_t a, modp_t b, modp_args args)
{
  *c = modp_mul(a, b, args);
}

#if !modp_isvalue
#error "modp_inv_val.fab should not be included"
#endif

static void 
modp_ininv(modp_c * a, modp_args args)
{
  *a = modp_inv(*a, args);
}

static void 
modp_ooinv(modp_c * c, modp_t a, modp_args args)
{
  *c = modp_inv(a, args);
}

#if !modp_isvalue
#error "modp_div_val.fab should not be included"
#endif

static void 
modp_indiv(modp_c * a, modp_t b, modp_args args)
{
  *a = modp_div(*a, b, args);
}

static BOOL 
modp_oodiv(modp_c * c, modp_t a, modp_t b, modp_args args)
{
  *c = modp_div(a, b, args);
  return SUCCESS;
}
#if !modp_isvalue
#error "modp_add_val.fab should not be included"
#endif

static modp_t 
modp_mulv(modp_t a, modp_t b, int v, modp_args args)
{
  modp_t c = modp_mul(a, b, args);

  return (v == +1) ? c : modp_neg(c, args);
}

static void 
modp_inmulv(modp_c * a, modp_t b, int v, modp_args args)
{
  *a = modp_mulv(*a, b, v, args);
}

static void 
modp_oomulv(modp_c * c, modp_t a, modp_t b, int v, modp_args args)
{
  *c = modp_mulv(a, b, v, args);
}

static modp_t 
modp_addmul(modp_t a, int v, modp_t b, modp_t s, int w, modp_args args)
{
  return modp_addvw(a, v, modp_mul(b, s, args), w, args);
}

static modp_t 
modp_addsqr(modp_t a, int v, modp_t b, int w, modp_args args)
{
  return modp_addvw(a, v, modp_sqr(b, args), w, args);
}

static modp_t 
modp_muladd(modp_t a, modp_t s, int v, modp_t b, int w, modp_args args)
{
  return modp_addvw(modp_mul(a, s, args), v, b, w, args);
}

static modp_t 
modp_muladdmul(modp_t a, modp_t s, int v, modp_t b, modp_t t, int w, modp_args args)
{
  return modp_addvw(modp_mul(a, s, args), v, modp_mul(b, t, args), w, args);
}

static void 
modp_inaddmul(modp_c * a, int v, modp_t b, modp_t t, int w, modp_args args)
{
  *a = modp_addmul(*a, v, b, t, w, args);
}

static void 
modp_ooaddmul(modp_c * c, modp_t a, int v, modp_t b, modp_t t, int w, modp_args args)
{
  *c = modp_addmul(a, v, b, t, w, args);
}

static void 
modp_inaddsqr(modp_c * a, int v, modp_t b, int w, modp_args args)
{
  *a = modp_addsqr(*a, v, b, w, args);
}

static void 
modp_ooaddsqr(modp_c * c, modp_t a, int v, modp_t b, int w, modp_args args)
{
  *c = modp_addsqr(a, v, b, w, args);
}

static void 
modp_inmuladd(modp_c * a, modp_t s, int v, modp_t b, int w, modp_args args)
{
  *a = modp_muladd(*a, s, v, b, w, args);
}

static void 
modp_oomuladd(modp_c * c, modp_t a, modp_t s, int v, modp_t b, int w, modp_args args)
{
  *c = modp_muladd(a, s, v, b, w, args);
}

static void 
modp_inmuladdmul(modp_c * a, modp_t s, int v, modp_t b, modp_t t, int w, modp_args args)
{
  *a = modp_muladdmul(*a, s, v, b, t, w, args);
}

static void 
modp_oomuladdmul(modp_c * c, modp_t a, modp_t s, int v, modp_t b, modp_t t, int w, modp_args args)
{
  *c = modp_muladdmul(a, s, v, b, t, w, args);
}
static void 
modp_args_fileOut(id aFiler, modp_args * args)
{
}
static void 
modp_args_fileIn(id aFiler, modp_args * args)
{
}

static void 
modp_fileOut(id aFiler, modp_c * c, modp_args args)
{
  [aFiler fileOut:c type:'S'];
}

static void 
modp_fileIn(id aFiler, modp_c * c, modp_args args)
{
  [aFiler fileIn:c type:'S'];
}

@implementation IntegerModp

+ new
{
  return [self error:"You can't create a IntegerModp without modulus"];
}

+ modulus:(unsigned short)aPrime
{
  return [self int:0 modulus:aPrime];
}

- _setUpIntValue:(int)intValue modulus:(unsigned short)aPrime
{
  modulus = aPrime;
  value = modp_fromintvalue(intValue, modp_getargs(self));
  return self;
}

- int:(int)intValue modulus:(unsigned short)aPrime
{
  return [self notImplemented:_cmd];
}

+ int:(int)intValue modulus:(unsigned short)aPrime
{
  unsigned i, m = aPrime;

  if (m == 2)
    [self error:"Modulus Should Not Be Equal To 2"];
  if (m > USHRT_MAX_PRIME)
    [self error:"Modulus Too Big"];

  if (m != USHRT_MAX_PRIME) {
    for (i = 2; i <= m / 2; i++)
      if (m % i == 0)
	return [self error:"Modulus Not Prime"];
  }
  return [[super new] _setUpIntValue:intValue modulus:aPrime];
}

- class_vector
{
  return (id) [modp_vector class];
}

- class_vardnsrecdegsps_polynomial
{
  return (id) [modp_vardnsrecdegsps_polynomial class];
}

- class_vardnsrecdegdns_polynomial
{
  return (id) [modp_vardnsrecdegdns_polynomial class];
}

- class_vardnsexpdegsps_polynomial
{
  return (id) [modp_vardnsexpdegsps_polynomial class];
}

- clone
{
  self = [super clone];
  value = 0;
  assert(modulus);
  return self;
}

- (modp_t) modp_value
{
  return value;
}

- modp_value:(modp_t)aValue
{
  self = [self clone];
  value = aValue;
  return self;
}

- (modp_t *) modp_reference
{
  return &value;
}

- modp_reference:(modp_t *)aValue
{
  self = [self clone];
  value = *aValue;
  return self;
}

- (BOOL)sameClass:b			 {
  return isa == [b class] && modulus == [b modulus];
}
- (BOOL)inAdditiveSemiGroup		 {
  return YES;
}
- (BOOL)inAdditiveMonoid		 {
  return YES;
}
- (BOOL)inAdditiveGroup			 {
  return YES;
}

- (BOOL)inSemiGroup			 {
  return YES;
}
- (BOOL)inMonoid			 {
  return YES;
}
- (BOOL)inGroup				 {
  return YES;
}

- (BOOL)inEuclideanDomain		 {
  return YES;
}
- (BOOL)inIntegralDomain		 {
  return YES;
}
- (BOOL)inField				 {
  return YES;
}

- (BOOL)inOrderedSet			 {
  return NO;
}



- (unsigned) hash
{
  return modp_hash(modp_id2t(self), modp_getargs(self));
}

- (BOOL) isEqual:b
{
  [self checkSameClass:b];
  if (self == b) {
    return YES;
  } else {
    return modp_iseq(modp_id2t(self), modp_id2t(b), modp_getargs(self));
  }
}

- (BOOL) notEqual:b
{
  [self checkSameClass:b];
  if (self == b) {
    return NO;
  } else {
    return !modp_iseq(modp_id2t(self), modp_id2t(b), modp_getargs(self));
  }
}

- (unsigned short) modulus
{
  return modulus;
}

- (int) intValue
{
  return value;
}

- intValue:(int)aValue
{
  modp_t a = modp_fromintvalue(aValue, modp_getargs(self));

  return [self modp_value:a];
}


- zero
{
  return [self intValue:0];
}

- (BOOL) isZero
{
  assert([self check]);
  return modp_iszero(modp_id2t(self), modp_getargs(self));
}

- (BOOL) notZero
{
  assert([self check]);
  return NO == modp_iszero(modp_id2t(self), modp_getargs(self));
}

- (BOOL) isOpposite:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self isZero] || [self isCharacteristicTwo];
  } else {
    return modp_isop(modp_id2t(self), modp_id2t(b), modp_getargs(self));
  }
}

- (BOOL) notOpposite:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self notZero] && [self notCharacteristicTwo];
  } else {
    return NO == modp_isop(modp_id2t(self), modp_id2t(b), modp_getargs(self));
  }
}

- negate
{
  if ([self isZero] || [self isCharacteristicTwo]) {
    return [self cheapCopy];
  } else {
    modp_c c;

    assert([self check]);
    modp_ooneg(&c, modp_id2t(self), modp_getargs(self));
    return modp_c2id(self, &c);
  }
}

- negateSelf
{
  if (refCount > 1)
    return idReplace(self, [self negate]);

  if ([self isZero] || [self isCharacteristicTwo]) {
    return self;
  } else {
    assert([self check]);
    modp_inneg(modp_id2c(self), modp_getargs(self));
    return [self invalidate];
  }
}

- _double:(int)v
{
  if ([self isZero] || [self isCharacteristicTwo]) {
    return [self zero];
  } else {
    modp_c c;

    assert([self check]);
    modp_oodbl(&c, modp_id2t(self), v, modp_getargs(self));
    return modp_c2id(self, &c);
  }
}

- _doubleSelf:(int)v
{
  if (refCount > 1)
    return idReplace(self, [self _double:v]);

  if ([self isZero] || [self isCharacteristicTwo]) {
    return idReplace(self, [self zero]);
  } else {
    assert([self check]);
    modp_indbl(modp_id2c(self), v, modp_getargs(self));
    return [self invalidate];
  }
}

- double
{
  if ([self isZero] || [self isCharacteristicTwo]) {
    return [self zero];
  } else {
    modp_c c;

    assert([self check]);
    modp_oodbl(&c, modp_id2t(self), +1, modp_getargs(self));
    return modp_c2id(self, &c);
  }
}

- doubleSelf
{
  if (refCount > 1)
    return idReplace(self, [self double]);

  if ([self isZero] || [self isCharacteristicTwo]) {
    return idReplace(self, [self zero]);
  } else {
    assert([self check]);
    modp_indbl(modp_id2c(self), +1, modp_getargs(self));
    return [self invalidate];
  }
}

- _add:(int)v:b:(int)w
{
  assert([self sameClass:b] && v * v == 1 && w * w == 1);

  if (self == b) {
    return (v == w) ? [self _double : v]:[self zero];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return [b _negate:w];
      if ([b isZero])
	return [self _negate:v];
    }
    if (1) {
      modp_c c;

      modp_ooadd(&c, modp_id2t(self), v, modp_id2t(b), w, modp_getargs(self));
      return modp_c2id(self, &c);
    }
  }

  return nil;			/* to quiet some compilers, e.g. WATCOM */
}

- _addSelf:(int)v:b:(int)w
{
  assert([self sameClass:b] && v * v == 1 && w * w == 1);
  if (refCount > 1)
    return idReplace(self, [self _add:v:b:w]);

  if (self == b) {
    return (v == w) ? [self _doubleSelf : v]:idReplace(self, [self zero]);
  } else {
    if ([self inAdditiveMonoid]) {
      if ([b isZero])
	return [self _negateSelf:v];
      if ([self isZero])
	return idReplace(self, [b _negate:w]);
    }
    modp_inadd(modp_id2c(self), v, modp_id2t(b), w, modp_getargs(self));
    return [self invalidate];
  }

  return nil;
}

- add:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self double];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return [b cheapCopy];
      if ([b isZero])
	return [self cheapCopy];
    }
    if (1) {
      modp_c c;

      modp_ooadd(&c, modp_id2t(self), +1, modp_id2t(b), +1, modp_getargs(self));
      return modp_c2id(self, &c);
    }
  }

  return nil;
}

- addSelf:b
{
  assert([self sameClass:b]);
  if (refCount > 1)
    return idReplace(self, [self add:b]);

  if (self == b) {
    return [self doubleSelf];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([b isZero])
	return self;
      if ([self isZero])
	return idReplace(self, [b cheapCopy]);
    }
    modp_inadd(modp_id2c(self), +1, modp_id2t(b), +1, modp_getargs(self));
    return [self invalidate];
  }

  return nil;
}

- subtract:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self zero];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return [b negate];
      if ([b isZero])
	return [self cheapCopy];
    }
    if (1) {
      modp_c c;

      modp_ooadd(&c, modp_id2t(self), +1, modp_id2t(b), -1, modp_getargs(self));
      return modp_c2id(self, &c);
    }
  }

  return nil;
}

- subtractSelf:b
{
  assert([self sameClass:b]);
  if (refCount > 1)
    return idReplace(self, [self subtract:b]);

  if (self == b) {
    return idReplace(self, [self zero]);
  } else {
    if ([self inAdditiveMonoid]) {
      if ([b isZero])
	return self;
      if ([self isZero])
	return idReplace(self, [b negate]);
    }
    modp_inadd(modp_id2c(self), +1, modp_id2t(b), -1, modp_getargs(self));
    return [self invalidate];
  }

  return nil;
}


- one
{
  return [self intValue:+1];
}

- minusOne
{
  return [self intValue:-1];
}

- (BOOL) isOne
{
  assert([self check]);
  return modp_isone(modp_id2t(self), modp_getargs(self));
}

- (BOOL) notOne
{
  assert([self check]);
  return NO == modp_isone(modp_id2t(self), modp_getargs(self));
}

- (BOOL) isMinusOne
{
  assert([self check]);
  return modp_isminusone(modp_id2t(self), modp_getargs(self));
}

- (BOOL) notMinusOne
{
  assert([self check]);
  return NO == modp_isminusone(modp_id2t(self), modp_getargs(self));
}

- square
{
  if ([self inMonoid] && [self isOne])
    return [self cheapCopy];
  if ([self inAdditiveMonoid]) {
    if ([self isMinusOne])
      return [self one];
    if ([self isZero])
      return [self cheapCopy];
  }
  if (1) {
    modp_c c;

    assert([self check]);
    modp_oosqr(&c, modp_id2t(self), modp_getargs(self));
    return modp_c2id(self, &c);
  }
  return nil;
}

- squareSelf
{
  if (refCount > 1)
    return idReplace(self, [self square]);

  if ([self inMonoid] && [self isOne])
    return self;
  if ([self inAdditiveMonoid]) {
    if ([self isMinusOne])
      return idReplace(self, [self one]);
    if ([self isZero])
      return self;
  }
  assert([self check]);
  modp_insqr(modp_id2c(self), modp_getargs(self));
  return [self invalidate];
}

- multiply:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self square];
  } else {
    if ([self inMonoid]) {
      if ([self isOne])
	return [b cheapCopy];
      if ([b isOne])
	return [self cheapCopy];
    }
    if ([self inAdditiveMonoid]) {
      if ([self isZero] || [b isZero])
	return [self zero];
      if ([self isMinusOne])
	return [b negate];
      if ([b isMinusOne])
	return [self negate];
    }
    if (1) {
      modp_c c;

      modp_oomul(&c, modp_id2t(self), modp_id2t(b), modp_getargs(self));
      return modp_c2id(self, &c);
    }
  }

  return nil;
}

- multiplySelf:b
{
  assert([self sameClass:b]);
  if (refCount > 1)
    return idReplace(self, [self multiply:b]);

  if (self == b) {
    return [self squareSelf];
  } else {
    if ([self inMonoid]) {
      if ([b isOne])
	return self;
      if ([self isOne])
	return idReplace(self, [b cheapCopy]);
    }
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return self;
      if ([b isZero])
	return idReplace(self, [b cheapCopy]);
      if ([b isMinusOne])
	return [self negateSelf];
      if ([self isMinusOne])
	return idReplace(self, [b negate]);
    }
    modp_inmul(modp_id2c(self), modp_id2t(b), modp_getargs(self));
    return [self invalidate];
  }

  return nil;
}


- inverse
{
  if ([self inMonoid] && [self isOne])
    return [self cheapCopy];
  if ([self inAdditiveMonoid]) {
    if ([self isMinusOne])
      return [self cheapCopy];
    if ([self isZero])
      return nil;
  }
  if (1) {
    modp_c c;

    assert([self check]);
    modp_ooinv(&c, modp_id2t(self), modp_getargs(self));
    return modp_c2id(self, &c);
  }
  return nil;
}

- inverseSelf
{
  if (refCount > 1)
    return idReplace(self, [self inverse]);

  if ([self inMonoid] && [self isOne])
    return self;
  if ([self inAdditiveMonoid]) {
    if ([self isMinusOne])
      return self;
    if ([self isZero])
      return nil;
  }
  assert([self check]);
  modp_ininv(modp_id2c(self), modp_getargs(self));
  return [self invalidate];
}

- divide:b
{
  if ([b isZero])
    return nil;
  if (self == b) {
    return [self one];
  } else {
    id tmp = [b inverse];

    self = [self multiply:tmp];
    tmp = [tmp free];
    return self;
  }
}

- _add:(int)v:B multiply:b:(int)w
{
  assert([self sameClass:B] && [self sameClass:b] && v * v == 1 && w * w == 1);

  if ([self isZero])
    return [B _multiply:b:+w];
  if ([B isZero])
    return [self _negate:v];
  if ([b isZero])
    return [self _negate:v];

  if ([b isOne])
    return [self _add:v:B:+w];
  if ([B isOne])
    return [self _add:v:b:+w];
  if ([b isMinusOne])
    return [self _add:v:B:-w];
  if ([B isMinusOne])
    return [self _add:v:b:-w];

  if (1) {
    modp_c c;

    modp_ooaddmul(&c, modp_id2t(self), v, modp_id2t(B), modp_id2t(b), w, modp_getargs(self));
    return modp_c2id(self, &c);
  }
  return nil;
}

- _addSelf:(int)v:B multiply:b:(int)w
{
  if (refCount > 1)
    return idReplace(self, [self _add:v:B multiply:b:w]);
  assert([self sameClass:B] && [self sameClass:b] && v * v == 1 && w * w == 1);

  if ([self isZero])
    return idReplace(self, [B _multiply:b:+w]);
  if ([B isZero])
    return [self _negateSelf:v];
  if ([b isZero])
    return [self _negateSelf:v];

  if ([b isOne])
    return [self _addSelf:v:B:+w];
  if ([B isOne])
    return [self _addSelf:v:b:+w];
  if ([b isMinusOne])
    return [self _addSelf:v:B:-w];
  if ([B isMinusOne])
    return [self _addSelf:v:b:-w];

  modp_inaddmul(modp_id2c(self), v, modp_id2t(B), modp_id2t(b), w, modp_getargs(self));
  return [self invalidate];
}

- _multiply:a:(int)v add:B:(int)w
{
  assert([self sameClass:a] && [self sameClass:B] && v * v == 1 && w * w == 1);

  if ([self isZero])
    return [B _negate:+w];
  if ([a isZero])
    return [B _negate:+w];
  if ([B isZero])
    return [self _multiply:a:+v];

  if ([self isOne])
    return [a _add:+v:B:+w];
  if ([a isOne])
    return [self _add:+v:B:+w];
  if ([self isMinusOne])
    return [a _add:-v:B:+w];
  if ([a isMinusOne])
    return [self _add:-v:B:+w];

  if (1) {
    modp_c c;

    modp_oomuladd(&c, modp_id2t(self), modp_id2t(a), v, modp_id2t(B), w, modp_getargs(self));
    return modp_c2id(self, &c);
  }
  return nil;
}

- _multiplySelf:a:(int)v add:B:(int)w
{
  if (refCount > 1)
    return idReplace(self, [self _multiply:a:v add:B:w]);
  assert([self sameClass:a] && [self sameClass:B] && v * v == 1 && w * w == 1);

  if ([self isZero])
    return idReplace(self, [B _negate:+w]);
  if ([a isZero])
    return idReplace(self, [B _negate:+w]);
  if ([B isZero])
    return [self _multiplySelf:a:+v];

  if ([self isOne])
    return idReplace(self, [a _add:+v:B:+w]);
  if ([a isOne])
    return [self _addSelf:+v:B:+w];
  if ([self isMinusOne])
    return idReplace(self, [a _add:-v:B:+w]);
  if ([a isMinusOne])
    return [self _addSelf:-v:B:+w];

  modp_inmuladd(modp_id2c(self), modp_id2t(a), v, modp_id2t(B), w, modp_getargs(self));
  return [self invalidate];
}

- _multiply:a:(int)v add:B multiply:b:(int)w
{
  assert([self sameClass:a] && [self sameClass:B] && [self sameClass:b] && v * v == 1 && w * w == 1);

  if ([self isZero])
    return [B _multiply:b:+w];
  if ([a isZero])
    return [B _multiply:b:+w];
  if ([B isZero])
    return [self _multiply:a:+v];
  if ([b isZero])
    return [self _multiply:a:+v];

  if ([self isOne])
    return [a _add:+v:B multiply:b:+w];
  if ([a isOne])
    return [self _add:+v:B multiply:b:+w];
  if ([b isOne])
    return [self _multiply:a:+v add:B:+w];
  if ([B isOne])
    return [self _multiply:a:+v add:b:+w];

  if ([self isMinusOne])
    return [a _add:-v:B multiply:b:+w];
  if ([a isMinusOne])
    return [self _add:-v:B multiply:b:+w];
  if ([b isMinusOne])
    return [self _multiply:a:+v add:B:-w];
  if ([B isMinusOne])
    return [self _multiply:a:+v add:b:-w];

  if (1) {
    modp_c c;

    modp_oomuladdmul(&c, modp_id2t(self), modp_id2t(a), v, modp_id2t(B), modp_id2t(b), w, modp_getargs(self));
    return modp_c2id(self, &c);
  }
  return nil;
}

- _multiplySelf:a:(int)v add:B multiply:b:(int)w
{
  if (refCount > 1)
    return idReplace(self, [self _multiply:a:v add:B multiply:b:w]);
  assert([self sameClass:a] && [self sameClass:B] && [self sameClass:b] && v * v == 1 && w * w == 1);

  if ([self isZero])
    return idReplace(self, [B _multiply:b:+w]);
  if ([a isZero])
    return idReplace(self, [B _multiply:b:+w]);
  if ([B isZero])
    return [self _multiplySelf:a:+v];
  if ([b isZero])
    return [self _multiplySelf:a:+v];

  if ([self isOne])
    return idReplace(self, [a _add:+v:B multiply:b:+w]);
  if ([a isOne])
    return [self _addSelf:+v:B multiply:b:+w];
  if ([b isOne])
    return [self _multiplySelf:a:+v add:B:+w];
  if ([B isOne])
    return [self _multiplySelf:a:+v add:b:+w];

  if ([self isMinusOne])
    return idReplace(self, [a _add:-v:B multiply:b:+w]);
  if ([a isMinusOne])
    return [self _addSelf:-v:B multiply:b:+w];
  if ([b isMinusOne])
    return [self _multiplySelf:a:+v add:B:-w];
  if ([B isMinusOne])
    return [self _multiplySelf:a:+v add:b:-w];

  modp_inmuladdmul(modp_id2c(self), modp_id2t(a), v, modp_id2t(B), modp_id2t(b), w, modp_getargs(self));
  return [self invalidate];
}



- (int) characteristic
{
  return modulus;
}

- frobenius
{
  return [self cheapCopy];
}

- frobeniusInverse
{
  return [self cheapCopy];
}


- (BOOL) printsLeadingSign
{
  return NO;
}

- printOn:(IOD)aFile
{
  fprintf(aFile, "%i", value);
  return self;
}


- fileOutOn:aFiler
{
  [super fileOutOn:aFiler];
  [aFiler fileOut:&value type:'S'];
  [aFiler fileOut:&modulus type:'S'];
  return self;
}

- fileInFrom:aFiler
{
  [super fileInFrom:aFiler];
  [aFiler fileIn:&value type:'S'];
  [aFiler fileIn:&modulus type:'S'];
  return self;
}

@end

