GCC Code Coverage Report
Directory: . Exec Total Coverage
File: deps/install/include/symfpu/core/remainder.h Lines: 54 54 100.0 %
Date: 2021-03-22 Branches: 102 173 59.0 %

Line Exec Source
1
/*
2
** Copyright (C) 2018 Martin Brain
3
**
4
** See the file LICENSE for licensing information.
5
*/
6
7
/*
8
** remainder.h
9
**
10
** Martin Brain
11
** martin.brain@cs.ox.ac.uk
12
** 14/12/16
13
**
14
** Computing the IEEE-754 remainder of arbitrary precision floats
15
**
16
*/
17
18
#include "symfpu/core/unpackedFloat.h"
19
#include "symfpu/core/ite.h"
20
#include "symfpu/core/rounder.h"
21
#include "symfpu/core/operations.h"
22
#include "symfpu/core/add.h"
23
#include "symfpu/core/sign.h"
24
25
26
#ifndef SYMFPU_REMAINDER
27
#define SYMFPU_REMAINDER
28
29
namespace symfpu {
30
31
template <class t>
32
105
  unpackedFloat<t> addRemainderSpecialCases (const typename t::fpt &format,
33
					    const unpackedFloat<t> &left,
34
					    const unpackedFloat<t> &right,
35
					    const unpackedFloat<t> &remainderResult) {
36
  typedef typename t::prop prop;
37
38
115
  prop eitherArgumentNan(left.getNaN() || right.getNaN());
39
115
  prop generateNan(left.getInf() || right.getZero());
40
115
  prop isNan(eitherArgumentNan || generateNan);
41
42
115
  prop passThrough((!(left.getInf() || left.getNaN()) && right.getInf()) ||
43
		   left.getZero());
44
45
  return ITE(isNan,
46
	     unpackedFloat<t>::makeNaN(format),
47
	     ITE(passThrough,
48
		 left,
49
115
		 remainderResult));
50
 }
51
52
/* Let left = x*2^e, right = y*2^f, x \in [1,2), y \in [1,2)
53
 * x/y \in (0.5,2)   x > y  x/y \in (1,2)   x < y (0.5,1)
54
 *
55
 *  rem =  x*2^e     - (y*2^f * int((x*2^e) / (y*2^f)))
56
 *      =  x*2^e     - (y*2^f * int((x/y) * 2^{e-f}))
57
 *      = (x*2^{e-f} - (y     * int((x/y) * 2^{e-f}))) * 2^f
58
 *
59
 *
60
 * If e - f >  0
61
 *      = (x*2^{e-f} - (y     * int((x*2^{e-f})/y)) * 2^f
62
 *
63
 *
64
 * If e - f == 0
65
 *      = (x         - (y     * int((x/y)          ))) * 2^f
66
 *      = ITTE(x ?= y,
67
 *             (x - (y * int[guard=1,sticky=1])) * 2^f
68
 *             (x -  y) * 2^f,
69
 *             ...)
70
 *      = ITTE(x ?= y,
71
 *             (x - (y * int[guard=1,sticky=1])) * 2^f
72
 *             left - right,
73
 *             ...)
74
 *
75
 *
76
 * If e - f == -1
77
 *      = (x*2^{-1}  - (y     * int((x/y) * 2^{-1 }))) * 2^f
78
 *      = ITTE(x ?= y,
79
 *             (x*2^{-1}  - (y * int[guard=0,sticky=1])) * 2^f
80
 *             (x*2^{-1}  - (y * int[guard=1,sticky=0])) * 2^f
81
 *             (x*2^{-1}  - (y * int[guard=1,sticky=1])) * 2^f
82
 *
83
 * If e - f <= -2
84
 *      = (x*2^{e-f}  - (y     * int[guard=0,sticky=1])) * 2^f
85
 *      = ITE(int[guard=0,sticky=1],
86
 *            x*2^e  - y*2^f,
87
 *            left)
88
 *      = ITE(int[guard=0,sticky=1],
89
 *            left  - right,
90
 *            left)
91
 *
92
 */
93
94
// Divide for max(e - f, 0) cycles
95
// The equal case, if you divide you use to extract the even bit of n, also save the rem.
96
// Then one more cycle for the guard bit
97
// Use that remainder to work out the sticky bit
98
// Round and either subtract or not from saved rem
99
// Output at 2^f
100
101
template <class t>
102
105
  unpackedFloat<t> arithmeticRemainder (const typename t::fpt &format,
103
					const typename t::rm &roundingMode,
104
					const unpackedFloat<t> &left,
105
					const unpackedFloat<t> &right) {
106
  typedef typename t::bwt bwt;
107
  typedef typename t::prop prop;
108
  typedef typename t::ubv ubv;
109
  typedef typename t::sbv sbv;
110
  //typedef typename t::fpt fpt;
111
112
105
  PRECONDITION(left.valid(format));
113
105
  PRECONDITION(right.valid(format));
114
115
  // Compute sign
116
115
  prop remainderSign(left.getSign());
117
118
119
  // Compute exponent difference
120
210
  sbv exponentDifference(expandingSubtract<t>(left.getExponent(), right.getExponent()));
121
105
  bwt edWidth(exponentDifference.getWidth());
122
123
  // Extend for divide steps
124
210
  ubv lsig(left.getSignificand().extend(1));
125
210
  ubv rsig(right.getSignificand().extend(1));
126
127
128
210
  ubv first(divideStep<t>(lsig,rsig).result);
129
105
  ubv *running = new ubv(first); // To avoid running out of stack space loop with a pointer
130
131
105
  bwt maxDifference = unpackedFloat<t>::maximumExponentDifference(format);
132
4332
  for (bwt i = maxDifference - 1; i > 0; i--) {
133
4844
    prop needPrevious(exponentDifference > sbv(edWidth, i));
134
4227
    probabilityAnnotation<t>(needPrevious, (i > (maxDifference / 2)) ? VERYUNLIKELY : UNLIKELY);
135
136
8454
    ubv r(ITE(needPrevious, *running, lsig));
137
4227
    delete running;  // We assume the value / reference has been transfered to ITE
138
4227
    running = new ubv(divideStep<t>(r, rsig).result);
139
  }
140
141
  // The zero exponent difference case is a little different
142
  // as we need the result bit for the even flag
143
  // and the actual result for the final
144
115
  prop lsbRoundActive(exponentDifference > -sbv::one(edWidth));  // i.e. >= 0
145
146
115
  prop needPrevious(exponentDifference > sbv::zero(edWidth));
147
105
  probabilityAnnotation<t>(needPrevious, UNLIKELY);
148
149
210
  ubv r0(ITE(needPrevious, *running, lsig));
150
105
  delete running;
151
210
  resultWithRemainderBit<t> dsr(divideStep<t>(r0, rsig));
152
153
115
  prop integerEven(!lsbRoundActive || !dsr.remainderBit);  // Note negation of guardBit
154
155
156
  // The same to get the guard flag
157
115
  prop guardRoundActive(exponentDifference > -sbv(edWidth,2));  // i.e. >= -1
158
159
210
  ubv rm1(ITE(lsbRoundActive, dsr.result, lsig));
160
210
  resultWithRemainderBit<t> dsrg(divideStep<t>(rm1, rsig));
161
162
115
  prop guardBit(guardRoundActive && dsrg.remainderBit);
163
164
115
  prop stickyBit(!ITE(guardRoundActive,
165
		      dsrg.result,
166
		      lsig).isAllZeros());
167
168
169
  // The base result if lsbRoundActive
170
210
  unpackedFloat<t> reconstruct(remainderSign,
171
			       right.getExponent(),
172
105
			       dsr.result.extract(lsig.getWidth() - 1,1)); // dsr shifts right as last action so is safe
173
174
175
105
  probabilityAnnotation<t>(needPrevious, UNLIKELY);   // Perhaps stretching it a bit but good for approximation
176
210
  unpackedFloat<t> candidateResult(ITE(lsbRoundActive,
177
				       reconstruct.normaliseUpDetectZero(format),
178
				       left));
179
180
  // The final subtract is a little different as previous ones were
181
  // guaranteed to be positive
182
  // TODO : This could be improved as these don't need special cases, etc.
183
184
  // From the rounding of the big integer multiple
185
115
  prop bonusSubtract(roundingDecision<t>(roundingMode,
186
					 remainderSign,
187
					 integerEven,
188
					 guardBit,
189
					 stickyBit,
190
					 prop(false)));
191
105
  probabilityAnnotation<t>(bonusSubtract, UNLIKELY); // Again, more like 50/50
192
193
  // The big integer has sign left.getSign() ^ right.getSign() so we subtract something of left.getSign().
194
  // For the integer part we handle this by working with absolutes (ignoring the sign) and
195
  // adding it back in at the end.
196
  // However for the correction for the rounded part we need to take it into account
197
210
  unpackedFloat<t> signCorrectedRight(right, left.getSign());
198
105
  unpackedFloat<t> remainderResult(ITE(bonusSubtract,
199
				       add<t>(format,
200
					      roundingMode,
201
					      candidateResult,
202
					      signCorrectedRight,
203
					      false),
204
				       candidateResult));
205
206
  // TODO : fast path if first.isAllZeros()
207
208
105
  POSTCONDITION(remainderResult.valid(format));
209
210
210
  return remainderResult;
211
 }
212
213
214
// Put it all together...
215
template <class t>
216
105
  unpackedFloat<t> remainderWithRounding (const typename t::fpt &format,
217
					  const typename t::rm &roundingMode,
218
					  const unpackedFloat<t> &left,
219
					  const unpackedFloat<t> &right) {
220
  //typedef typename t::bwt bwt;
221
  //typedef typename t::prop prop;
222
  //typedef typename t::ubv ubv;
223
  //typedef typename t::sbv sbv;
224
225
105
  PRECONDITION(left.valid(format));
226
105
  PRECONDITION(right.valid(format));
227
228
210
  unpackedFloat<t> remainderResult(arithmeticRemainder(format, roundingMode, left, right));
229
230
  //unpackedFloat<t> result(addRemainderSpecialCases(format, left, right, roundedRemainderResult));
231
105
  unpackedFloat<t> result(addRemainderSpecialCases(format, left, right, remainderResult));
232
233
105
  POSTCONDITION(result.valid(format));
234
235
210
  return result;
236
 }
237
238
// IEEE-754 remainder always uses round to nearest, ties to even
239
template <class t>
240
105
  unpackedFloat<t> remainder (const typename t::fpt &format,
241
			      const unpackedFloat<t> &left,
242
			      const unpackedFloat<t> &right) {
243
244
105
  return remainderWithRounding<t>(format, t::RNE(), left, right);
245
 }
246
247
248
}
249
250
#endif
251