Return to dblsub.c CVS log

Up to [gforth] / gforth / engine

Annotation of gforth/engine/dblsub.c, revision 1.7

1.1       anton       1: /* some routines for double-cell arithmetic
                      2:    only used if BUGGY_LONG_LONG
                      3: 
1.4       anton       4:    Copyright (C) 1996,2000,2003 Free Software Foundation, Inc.
1.1       anton       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
1.2       anton      19:  * Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA.
1.1       anton      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: 
1.5       anton      84:   if (v==0)
                     85:     throw(BALL_DIVZERO);
                     86:   if (h>=v)
                     87:     throw(BALL_RESULTRANGE);
1.1       anton      88:   for (;;)
                     89:     {
                     90:       if (c || h >= v)
                     91:        {
                     92:          q++;
                     93:          h -= v;
                     94:        }
                     95:       if (--i < 0)
                     96:        break;
                     97:       c = HIGHBIT (h);
                     98:       h <<= 1;
                     99:       h += HIGHBIT (l);
                    100:       l <<= 1;
                    101:       q <<= 1;
                    102:     }
                    103:   res.hi = h;
                    104:   res.lo = q;
                    105:   return res;
                    106: }
                    107: 
                    108: DCell smdiv (DCell num, Cell denom)    /* symmetric divide procedure, mixed prec */
                    109: {
                    110:   DCell res;
                    111:   Cell numsign=num.hi;
                    112:   Cell denomsign=denom;
                    113: 
                    114:   if (numsign < 0)
                    115:     num = dnegate (num);
                    116:   if (denomsign < 0)
                    117:     denom = -denom;
                    118:   res = UD2D(umdiv (D2UD(num), denom));
1.5       anton     119:   if ((numsign^denomsign)<0) {
1.6       anton     120:     res.lo = -res.lo;
1.7     ! anton     121:     if (((Cell)res.lo) > 0) /* note: == 0 is possible */
        !           122:       throw(BALL_RESULTRANGE);
1.6       anton     123:   } else {
1.7     ! anton     124:     if (((Cell)res.lo) < 0)
1.5       anton     125:       throw(BALL_RESULTRANGE);
                    126:   }
1.1       anton     127:   if (numsign<0)
                    128:     res.hi = -res.hi;
                    129:   return res;
                    130: }
                    131: 
                    132: DCell fmdiv (DCell num, Cell denom)    /* floored divide procedure, mixed prec */
                    133: {
                    134:   /* I have this technique from Andrew Haley */
                    135:   DCell res;
                    136:   Cell denomsign=denom;
1.6       anton     137:   Cell numsign;
1.1       anton     138: 
                    139:   if (denom < 0) {
                    140:     denom = -denom;
                    141:     num = dnegate(num);
                    142:   }
1.6       anton     143:   numsign = num.hi;
                    144:   if (numsign < 0)
1.1       anton     145:     num.hi += denom;
                    146:   res = UD2D(umdiv(D2UD(num),denom));
1.7     ! anton     147:   if ((numsign^((Cell)res.lo)) < 0)
        !           148:     throw(BALL_RESULTRANGE);
1.1       anton     149:   if (denomsign<0)
                    150:     res.hi = -res.hi;
                    151:   return res;
                    152: }