alpar@1: /* glplib02.c (64-bit arithmetic) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpenv.h"
alpar@1: #include "glplib.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xlset - expand integer to long integer
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  glp_long xlset(int x);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xlset returns x expanded to long integer. */
alpar@1: 
alpar@1: glp_long xlset(int x)
alpar@1: {     glp_long t;
alpar@1:       t.lo = x, t.hi = (x >= 0 ? 0 : -1);
alpar@1:       return t;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xlneg - negate long integer
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  glp_long xlneg(glp_long x);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xlneg returns the difference  0 - x. */
alpar@1: 
alpar@1: glp_long xlneg(glp_long x)
alpar@1: {     if (x.lo)
alpar@1:          x.lo = - x.lo, x.hi = ~x.hi;
alpar@1:       else
alpar@1:          x.hi = - x.hi;
alpar@1:       return x;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xladd - add long integers
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  glp_long xladd(glp_long x, glp_long y);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xladd returns the sum x + y. */
alpar@1: 
alpar@1: glp_long xladd(glp_long x, glp_long y)
alpar@1: {     if ((unsigned int)x.lo <= 0xFFFFFFFF - (unsigned int)y.lo)
alpar@1:          x.lo += y.lo, x.hi += y.hi;
alpar@1:       else
alpar@1:          x.lo += y.lo, x.hi += y.hi + 1;
alpar@1:       return x;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xlsub - subtract long integers
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  glp_long xlsub(glp_long x, glp_long y);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xlsub returns the difference x - y. */
alpar@1: 
alpar@1: glp_long xlsub(glp_long x, glp_long y)
alpar@1: {     return
alpar@1:          xladd(x, xlneg(y));
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xlcmp - compare long integers
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  int xlcmp(glp_long x, glp_long y);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xlcmp returns the sign of the difference x - y. */
alpar@1: 
alpar@1: int xlcmp(glp_long x, glp_long y)
alpar@1: {     if (x.hi >= 0 && y.hi <  0) return +1;
alpar@1:       if (x.hi <  0 && y.hi >= 0) return -1;
alpar@1:       if ((unsigned int)x.hi < (unsigned int)y.hi) return -1;
alpar@1:       if ((unsigned int)x.hi > (unsigned int)y.hi) return +1;
alpar@1:       if ((unsigned int)x.lo < (unsigned int)y.lo) return -1;
alpar@1:       if ((unsigned int)x.lo > (unsigned int)y.lo) return +1;
alpar@1:       return 0;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xlmul - multiply long integers
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  glp_long xlmul(glp_long x, glp_long y);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xlmul returns the product x * y. */
alpar@1: 
alpar@1: glp_long xlmul(glp_long x, glp_long y)
alpar@1: {     unsigned short xx[8], yy[4];
alpar@1:       xx[4] = (unsigned short)x.lo;
alpar@1:       xx[5] = (unsigned short)(x.lo >> 16);
alpar@1:       xx[6] = (unsigned short)x.hi;
alpar@1:       xx[7] = (unsigned short)(x.hi >> 16);
alpar@1:       yy[0] = (unsigned short)y.lo;
alpar@1:       yy[1] = (unsigned short)(y.lo >> 16);
alpar@1:       yy[2] = (unsigned short)y.hi;
alpar@1:       yy[3] = (unsigned short)(y.hi >> 16);
alpar@1:       bigmul(4, 4, xx, yy);
alpar@1:       x.lo = (unsigned int)xx[0] | ((unsigned int)xx[1] << 16);
alpar@1:       x.hi = (unsigned int)xx[2] | ((unsigned int)xx[3] << 16);
alpar@1:       return x;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xldiv - divide long integers
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  glp_ldiv xldiv(glp_long x, glp_long y);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xldiv returns a structure of type glp_ldiv containing
alpar@1: *  members quot (the quotient) and rem (the remainder), both of type
alpar@1: *  glp_long. */
alpar@1: 
alpar@1: glp_ldiv xldiv(glp_long x, glp_long y)
alpar@1: {     glp_ldiv t;
alpar@1:       int m, sx, sy;
alpar@1:       unsigned short xx[8], yy[4];
alpar@1:       /* sx := sign(x) */
alpar@1:       sx = (x.hi < 0);
alpar@1:       /* sy := sign(y) */
alpar@1:       sy = (y.hi < 0);
alpar@1:       /* x := |x| */
alpar@1:       if (sx) x = xlneg(x);
alpar@1:       /* y := |y| */
alpar@1:       if (sy) y = xlneg(y);
alpar@1:       /* compute x div y and x mod y */
alpar@1:       xx[0] = (unsigned short)x.lo;
alpar@1:       xx[1] = (unsigned short)(x.lo >> 16);
alpar@1:       xx[2] = (unsigned short)x.hi;
alpar@1:       xx[3] = (unsigned short)(x.hi >> 16);
alpar@1:       yy[0] = (unsigned short)y.lo;
alpar@1:       yy[1] = (unsigned short)(y.lo >> 16);
alpar@1:       yy[2] = (unsigned short)y.hi;
alpar@1:       yy[3] = (unsigned short)(y.hi >> 16);
alpar@1:       if (yy[3])
alpar@1:          m = 4;
alpar@1:       else if (yy[2])
alpar@1:          m = 3;
alpar@1:       else if (yy[1])
alpar@1:          m = 2;
alpar@1:       else if (yy[0])
alpar@1:          m = 1;
alpar@1:       else
alpar@1:          xerror("xldiv: divide by zero\n");
alpar@1:       bigdiv(4 - m, m, xx, yy);
alpar@1:       /* remainder in x[0], x[1], ..., x[m-1] */
alpar@1:       t.rem.lo = (unsigned int)xx[0], t.rem.hi = 0;
alpar@1:       if (m >= 2) t.rem.lo |= (unsigned int)xx[1] << 16;
alpar@1:       if (m >= 3) t.rem.hi = (unsigned int)xx[2];
alpar@1:       if (m >= 4) t.rem.hi |= (unsigned int)xx[3] << 16;
alpar@1:       if (sx) t.rem = xlneg(t.rem);
alpar@1:       /* quotient in x[m], x[m+1], ..., x[4] */
alpar@1:       t.quot.lo = (unsigned int)xx[m], t.quot.hi = 0;
alpar@1:       if (m <= 3) t.quot.lo |= (unsigned int)xx[m+1] << 16;
alpar@1:       if (m <= 2) t.quot.hi = (unsigned int)xx[m+2];
alpar@1:       if (m <= 1) t.quot.hi |= (unsigned int)xx[m+3] << 16;
alpar@1:       if (sx ^ sy) t.quot = xlneg(t.quot);
alpar@1:       return t;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  xltod - convert long integer to double
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glplib.h"
alpar@1: *  double xltod(glp_long x);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine xltod returns x converted to double. */
alpar@1: 
alpar@1: double xltod(glp_long x)
alpar@1: {     double s, z;
alpar@1:       if (x.hi >= 0)
alpar@1:          s = +1.0;
alpar@1:       else
alpar@1:          s = -1.0, x = xlneg(x);
alpar@1:       if (x.hi >= 0)
alpar@1:          z = 4294967296.0 * (double)x.hi + (double)(unsigned int)x.lo;
alpar@1:       else
alpar@1:       {  xassert(x.hi == 0x80000000 && x.lo == 0x00000000);
alpar@1:          z = 9223372036854775808.0; /* 2^63 */
alpar@1:       }
alpar@1:       return s * z;
alpar@1: }
alpar@1: 
alpar@1: char *xltoa(glp_long x, char *s)
alpar@1: {     /* convert long integer to character string */
alpar@1:       static const char *d = "0123456789";
alpar@1:       glp_ldiv t;
alpar@1:       int neg, len;
alpar@1:       if (x.hi >= 0)
alpar@1:          neg = 0;
alpar@1:       else
alpar@1:          neg = 1, x = xlneg(x);
alpar@1:       if (x.hi >= 0)
alpar@1:       {  len = 0;
alpar@1:          while (!(x.hi == 0 && x.lo == 0))
alpar@1:          {  t = xldiv(x, xlset(10));
alpar@1:             xassert(0 <= t.rem.lo && t.rem.lo <= 9);
alpar@1:             s[len++] = d[t.rem.lo];
alpar@1:             x = t.quot;
alpar@1:          }
alpar@1:          if (len == 0) s[len++] = d[0];
alpar@1:          if (neg) s[len++] = '-';
alpar@1:          s[len] = '\0';
alpar@1:          strrev(s);
alpar@1:       }
alpar@1:       else
alpar@1:          strcpy(s, "-9223372036854775808"); /* -2^63 */
alpar@1:       return s;
alpar@1: }
alpar@1: 
alpar@1: /**********************************************************************/
alpar@1: 
alpar@1: #if 0
alpar@1: #include "glprng.h"
alpar@1: 
alpar@1: #define N_TEST 1000000
alpar@1: /* number of tests */
alpar@1: 
alpar@1: static glp_long myrand(RNG *rand)
alpar@1: {     glp_long x;
alpar@1:       int k;
alpar@1:       k = rng_unif_rand(rand, 4);
alpar@1:       xassert(0 <= k && k <= 3);
alpar@1:       x.lo = rng_unif_rand(rand, 65536);
alpar@1:       if (k == 1 || k == 3)
alpar@1:       {  x.lo <<= 16;
alpar@1:          x.lo += rng_unif_rand(rand, 65536);
alpar@1:       }
alpar@1:       if (k <= 1)
alpar@1:          x.hi = 0;
alpar@1:       else
alpar@1:          x.hi = rng_unif_rand(rand, 65536);
alpar@1:       if (k == 3)
alpar@1:       {  x.hi <<= 16;
alpar@1:          x.hi += rng_unif_rand(rand, 65536);
alpar@1:       }
alpar@1:       if (rng_unif_rand(rand, 2)) x = xlneg(x);
alpar@1:       return x;
alpar@1: }
alpar@1: 
alpar@1: int main(void)
alpar@1: {     RNG *rand;
alpar@1:       glp_long x, y;
alpar@1:       glp_ldiv z;
alpar@1:       int test;
alpar@1:       rand = rng_create_rand();
alpar@1:       for (test = 1; test <= N_TEST; test++)
alpar@1:       {  x = myrand(rand);
alpar@1:          y = myrand(rand);
alpar@1:          if (y.lo == 0 && y.hi == 0) y.lo = 1;
alpar@1:          /* z.quot := x div y, z.rem := x mod y */
alpar@1:          z = xldiv(x, y);
alpar@1:          /* x must be equal to y * z.quot + z.rem */
alpar@1:          xassert(xlcmp(x, xladd(xlmul(y, z.quot), z.rem)) == 0);
alpar@1:       }
alpar@1:       xprintf("%d tests successfully passed\n", N_TEST);
alpar@1:       rng_delete_rand(rand);
alpar@1:       return 0;
alpar@1: }
alpar@1: #endif
alpar@1: 
alpar@1: /* eof */