GCC Code Coverage Report
Directory: . Exec Total Coverage
File: deps/install/include/symfpu/core/divide.h Lines: 43 43 100.0 %
Date: 2021-03-23 Branches: 100 163 61.3 %

Line Exec Source
1
/*
2
** Copyright (C) 2018 Martin Brain
3
**
4
** See the file LICENSE for licensing information.
5
*/
6
7
/*
8
** divide.h
9
**
10
** Martin Brain
11
** martin.brain@cs.ox.ac.uk
12
** 04/02/16
13
**
14
** Division 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
23
#ifndef SYMFPU_DIVIDE
24
#define SYMFPU_DIVIDE
25
26
namespace symfpu {
27
28
template <class t>
29
194
  unpackedFloat<t> addDivideSpecialCases (const typename t::fpt &format,
30
					  const unpackedFloat<t> &left,
31
					  const unpackedFloat<t> &right,
32
					  const typename t::prop &sign,
33
					  const unpackedFloat<t> &divideResult) {
34
  typedef typename t::prop prop;
35
36
200
  prop eitherArgumentNaN(left.getNaN() || right.getNaN());
37
413
  prop generateNaN((left.getInf() && right.getInf()) ||
38
241
		   (left.getZero() && right.getZero()));
39
40
200
  prop isNaN(eitherArgumentNaN || generateNaN);
41
42
413
  prop isInf((!left.getZero() && right.getZero()) ||
43
199
	     (left.getInf() && !right.getInf()));
44
45
413
  prop isZero((!left.getInf() && right.getInf()) ||
46
200
	      (left.getZero() && !right.getZero()));
47
48
  return ITE(isNaN,
49
	     unpackedFloat<t>::makeNaN(format),
50
	     ITE(isInf,
51
		 unpackedFloat<t>::makeInf(format, sign),
52
		 ITE(isZero,
53
		     unpackedFloat<t>::makeZero(format, sign),
54
200
		     divideResult)));
55
 }
56
57
58
 template <class t>
59
194
  unpackedFloat<t> arithmeticDivide (const typename t::fpt &format,
60
				       const unpackedFloat<t> &left,
61
				       const unpackedFloat<t> &right) {
62
  typedef typename t::bwt bwt;
63
  typedef typename t::prop prop;
64
  typedef typename t::ubv ubv;
65
  typedef typename t::sbv sbv;
66
  typedef typename t::fpt fpt;
67
68
194
  PRECONDITION(left.valid(format));
69
194
  PRECONDITION(right.valid(format));
70
71
  // Compute sign
72
200
  prop divideSign(left.getSign() ^ right.getSign());
73
74
  // Subtract up exponents
75
388
  sbv exponentDiff(expandingSubtract<t>(left.getExponent(),right.getExponent()));
76
  // Optimisation : do this late and use the increment as a carry in
77
78
388
  sbv min(unpackedFloat<t>::minSubnormalExponent(format));
79
388
  sbv max(unpackedFloat<t>::maxNormalExponent(format));
80
194
  INVARIANT(expandingSubtract<t>(min,max) <= exponentDiff);
81
194
  INVARIANT(exponentDiff <= expandingSubtract<t>(max, min));
82
  // Optimisation : use the if-then-lazy-else to avoid dividing for underflow and overflow
83
  //                subnormal / greater-than-2^sigwidth does not need to be evaluated
84
85
86
  // Divide the significands
87
  // We need significandWidth() + 1 bits in the result but the top one may cancel, so add two bits
88
388
  ubv extendedNumerator(left.getSignificand().append(ubv::zero(2)));
89
388
  ubv extendedDenominator(right.getSignificand().append(ubv::zero(2)));
90
91
388
  resultWithRemainderBit<t> divided(fixedPointDivide<t>(extendedNumerator, extendedDenominator));
92
93
94
194
  bwt resWidth(divided.result.getWidth());
95
388
  ubv topBit(divided.result.extract(resWidth - 1, resWidth - 1));
96
388
  ubv nextBit(divided.result.extract(resWidth - 2, resWidth - 2));
97
98
  // Alignment of inputs means at least one of the two MSB is 1
99
  //  i.e. [1,2) / [1,2) = [0.5,2)
100
  // Top bit is set by the first round of the divide and thus is 50/50 1 or 0
101
200
  prop topBitSet(topBit.isAllOnes());
102
194
  INVARIANT(topBitSet || nextBit.isAllOnes());
103
194
  INVARIANT(topBitSet == (left.getSignificand() >= right.getSignificand()));
104
105
106
  // Re-align
107
388
  sbv alignedExponent(conditionalDecrement<t>(!topBitSet, exponentDiff)); // Will not overflow as previously expanded
108
388
  ubv alignedSignificand(conditionalLeftShiftOne<t>(!topBitSet, divided.result)); // Will not loose information
109
110
  // Create the sticky bit, it is important that this is after alignment
111
388
  ubv finishedSignificand(alignedSignificand | ubv(divided.remainderBit).extend(resWidth - 1));
112
113
  // Put back together
114
194
  unpackedFloat<t> divideResult(divideSign, alignedExponent.extend(1), finishedSignificand);
115
116
  // A brief word about formats.
117
  // You might think that the extend above is unnecessary : it is from a overflow point of view.
118
  // It's needed so that it is a valid number with exponentWidth() + 2.
119
  // +1 is sufficient in almost all cases.  However:
120
  //    very large normal / very small subnormal
121
  // can have an exponent greater than very large normal * 2 ( + 1)
122
  // because the exponent range is asymmetric with more subnormal than normal.
123
124
194
  fpt extendedFormat(format.exponentWidth() + 2, format.significandWidth() + 2);
125
194
  POSTCONDITION(divideResult.valid(extendedFormat));
126
127
388
  return divideResult;
128
 }
129
130
131
// Put it all together...
132
template <class t>
133
194
  unpackedFloat<t> divide (const typename t::fpt &format,
134
			   const typename t::rm &roundingMode,
135
			   const unpackedFloat<t> &left,
136
			   const unpackedFloat<t> &right) {
137
  //typedef typename t::bwt bwt;
138
  //typedef typename t::prop prop;
139
  //typedef typename t::ubv ubv;
140
  //typedef typename t::sbv sbv;
141
142
194
  PRECONDITION(left.valid(format));
143
194
  PRECONDITION(right.valid(format));
144
145
388
  unpackedFloat<t> divideResult(arithmeticDivide(format, left, right));
146
147
388
  unpackedFloat<t> roundedDivideResult(rounder(format, roundingMode, divideResult));
148
149
194
  unpackedFloat<t> result(addDivideSpecialCases(format, left, right, roundedDivideResult.getSign(), roundedDivideResult));
150
151
194
  POSTCONDITION(result.valid(format));
152
153
388
  return result;
154
 }
155
156
157
}
158
159
#endif