/* Vladimir Kotlyar October 1995
 *
 * Integer manipulation routines
 *
 */

#include <string.h>
#include <assert.h>

#include "alpha_int.h"

int	
a_abs(int x)
{
  return (x > 0) ? x : -x;
}

int
a_max (int x, int y)
{
  return (x > y) ? x : y;
}


int
a_floor_div (int x, int y)
{
  int f;

  /* anyone knows a better way to handle all the different cases
   * for the signs of x, y?
   */

  if (y < 0) {
    x = -x;
    y = -y;
  }

  f = x / y;

  if (y * f > x) f--;

  return f;
}

int
a_gcd (int x, int y)
{
  int r;

  assert (x >= 0);
  assert (y >= 0);

  while (0 != y) {
    r = x % y;
    x = y;
    y = r;
  }

  return x;
}

/*
 * Extended GCD as in [Cohen] algorithm 1.3.6
 *
 * must have a >= 0, b >= 0
 *
 * finds u, v such that a*u + b*v == d == (a,b)
 * also |u|,|v| <= max{ |a|,|b| }. This is necessary for
 * HNF and SNF computations
 */
static int
xa_egcd_pos (int a, int b,
	     int* pu, int* pv)
{
  int v1, v3, t1, t3, q;
  int u, v, d;

  assert (a >= 0);
  assert (b >= 0);

  u = 1;
  d = a;

  if (0 == b) { 
    v = 0;
  }
  else {
    v1 = 0;
    v3 = b;

    while (0 != v3) { 
      q = d / v3;
      t3 = d % v3;
      t1 = u - q * v1;

      u = v1;
      d = v3;
      v1 = t1;
      v3 = t3;
    }

    v = (d - a*u)/b;
  }

  /* Now reduce u, v if b != 0*/
  /* unew = u - floor(u/b)*b */
  /* vnew = v + floor(u/b)*a */  

  if (0 != a) {
    /*
     *pu = u - a_floor_div(u,b) * b;
     *pv = v + a_floor_div(u,b) * a;
     */
    
    *pv = v - a_floor_div(v,a)*a;
    *pu = u + a_floor_div(v,a)*b;
  }
  else {
    *pu = u;
    *pv = v;
  }

  return d;
}

/*
 * This one works fine for all signs of a, b
 */
int
a_egcd (int a, int b, 
	int* pu, int* pv)
{
  int absa, absb;
  int signa, signb;
  int u, v, d;

  signa = (0 > a) ? -1 : 1;
  signb = (0 > b) ? -1 : 1;
  absa = a_abs(a);
  absb = a_abs(b);

  d = xa_egcd_pos (absa, absb, &u, &v);
  if (signa == -1) {
    u = -u;
  }
  if (signb == -1) {
    v = -v;
  }

  *pu = u;
  *pv = v;

  return d;
}

/*
 * Least common multiple: z := lcm(x, y)
 * x and y must be positive!!
 */
int
a_lcm (int x, int y)
{
  int 	g;
 
  assert (0 < x);
  assert (0 < y);

  g = a_gcd(x, y);
  x = x / g;
  y = y / g;

  return x * y * g;
}


#if 0
#include <stdio.h>

main ()
{
  int a, b, u, v, d;

  printf ("%d  %d  %d  %d  %d\n----\n%d %d %d %d %d\n"
	  "-----\n%d %d %d %d %d %d %d\n",
	  a_floor_div (10, 3),
	  a_floor_div (-10, 3),
	  a_floor_div (10, -3),
	  a_floor_div (-10, -3),
	  a_floor_div (9, 3),

	  10 % 3, (-10)%3, 10%(-3), (-10)%(-3), 9%3,

	  a_gcd(1275, 135),
	  a_gcd(135, 1275),
	  a_gcd(100, 0),
	  a_gcd(0, 100),
	  a_gcd(100, 100),
	  a_gcd(100, 1),
	  a_gcd(1, 100));


  while (2 == scanf ("%d%d", &a, &b)) {
    d = a_egcd (a, b, &u, &v);
    printf ("%d %d + %d %d == %d\n", a, u, b, v, d);
    printf ("lcm(|%d|,|%d|) == %d\n", a, b, a_lcm(a_abs(a), a_abs(b)));
  }

  return 0;
}
#endif
