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