Return to dblsub.c CVS log

Up to [gforth] / gforth / engine

changed FSF address in copyright messages

    1: /* some routines for double-cell arithmetic
    2:    only used if BUGGY_LONG_LONG
    3: 
    4:    Copyright (C) 1996 Free Software Foundation, Inc.
    5:  * Copyright (C) 1995  Dirk Uwe Zoller
    6:  *
    7:  * This library is free software; you can redistribute it and/or
    8:  * modify it under the terms of the GNU Library General Public
    9:  * License as published by the Free Software Foundation; either
   10:  * version 2 of the License, or (at your option) any later version.
   11:  *
   12:  * This library is distributed in the hope that it will be useful,
   13:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
   14:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   15:  * See the GNU Library General Public License for more details.
   16:  *
   17:  * You should have received a copy of the GNU Library General Public
   18:  * License along with this library; if not, write to the Free
   19:  * Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA.
   20: 
   21:  This has been adapted from pfe-0.9.14
   22:  */
   23: 
   24: #include "config.h"
   25: #include "forth.h"
   26: 
   27: /* !! a bit machine dependent */
   28: #define HALFCELL_BITS	(CELL_BITS/2)
   29: #define UH(x)		(((UCell)(x))>>HALFCELL_BITS)
   30: #define LH(x)		((x)&((~(UCell)0)>>HALFCELL_BITS))
   31: #define L2U(x)		(((UCell)(x))<<HALFCELL_BITS)
   32: #define HIGHBIT(x)	(((UCell)(x))>>(CELL_BITS-1))
   33: #define UD2D(ud)	({UDCell _ud=(ud); (DCell){_ud.hi,_ud.lo};})
   34: #define D2UD(d)		({DCell _d=(d); (UDCell){_d.hi,_d.lo};})
   35: 
   36: DCell dnegate(DCell d1)
   37: {
   38:   DCell res;
   39: 
   40:   res.hi = ~d1.hi + (d1.lo==0);
   41:   res.lo = -d1.lo;
   42:   return res;
   43: }
   44: 
   45: UDCell ummul (UCell a, UCell b)	/* unsigned multiply, mixed precision */
   46: {
   47:   UDCell res;
   48:   UCell m,ul,lu,uu;
   49: 
   50:   res.lo = a*b;
   51: /*ll = LH(a)*LH(b); dead code */
   52:   ul = UH(a)*LH(b);
   53:   lu = LH(a)*UH(b);
   54:   uu = UH(a)*UH(b);
   55:   m = ul+lu;
   56:   res.hi = (uu
   57: 	    + L2U(m<ul) /* the carry of ul+lu */
   58: 	    + UH(m)
   59: 	    + (res.lo<L2U(m)) /* the carry of ll+L2U(m) */
   60: 	    );
   61:   return res;
   62: }
   63: 
   64: DCell mmul (Cell a, Cell b)		/* signed multiply, mixed precision */
   65: {
   66:   DCell res;
   67: 
   68:   res = UD2D(ummul (a, b));
   69:   if (a < 0)
   70:     res.hi -= b;
   71:   if (b < 0)
   72:     res.hi -= a;
   73:   return res;
   74: }
   75: 
   76: UDCell umdiv (UDCell u, UCell v)
   77: /* Divide unsigned double by single precision using shifts and subtracts.
   78:    Return quotient in lo, remainder in hi. */
   79: {
   80:   int i = CELL_BITS, c = 0;
   81:   UCell q = 0, h = u.hi, l = u.lo;
   82:   UDCell res;
   83: 
   84:   for (;;)
   85:     {
   86:       if (c || h >= v)
   87: 	{
   88: 	  q++;
   89: 	  h -= v;
   90: 	}
   91:       if (--i < 0)
   92: 	break;
   93:       c = HIGHBIT (h);
   94:       h <<= 1;
   95:       h += HIGHBIT (l);
   96:       l <<= 1;
   97:       q <<= 1;
   98:     }
   99:   res.hi = h;
  100:   res.lo = q;
  101:   return res;
  102: }
  103: 
  104: DCell smdiv (DCell num, Cell denom)	/* symmetric divide procedure, mixed prec */
  105: {
  106:   DCell res;
  107:   Cell numsign=num.hi;
  108:   Cell denomsign=denom;
  109: 
  110:   if (numsign < 0)
  111:     num = dnegate (num);
  112:   if (denomsign < 0)
  113:     denom = -denom;
  114:   res = UD2D(umdiv (D2UD(num), denom));
  115:   if ((numsign^denomsign)<0)
  116:     res.lo = -res.lo;
  117:   if (numsign<0)
  118:     res.hi = -res.hi;
  119:   return res;
  120: }
  121: 
  122: DCell fmdiv (DCell num, Cell denom)	/* floored divide procedure, mixed prec */
  123: {
  124:   /* I have this technique from Andrew Haley */
  125:   DCell res;
  126:   Cell denomsign=denom;
  127: 
  128:   if (denom < 0) {
  129:     denom = -denom;
  130:     num = dnegate(num);
  131:   }
  132:   if (num.hi < 0)
  133:     num.hi += denom;
  134:   res = UD2D(umdiv(D2UD(num),denom));
  135:   if (denomsign<0)
  136:     res.hi = -res.hi;
  137:   return res;
  138: }