Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	build-coverage/deps/include/symfpu/core/remainder.h	Lines:	54	54	100.0 %
Date:	2021-09-13	Branches:	102	173	59.0 %


/*
** Copyright (C) 2018 Martin Brain
**
** See the file LICENSE for licensing information.
*/

/*
** remainder.h
**
** Martin Brain
** martin.brain@cs.ox.ac.uk
** 14/12/16
**
** Computing the IEEE-754 remainder of arbitrary precision floats
**
*/

#include "symfpu/core/unpackedFloat.h"
#include "symfpu/core/ite.h"
#include "symfpu/core/rounder.h"
#include "symfpu/core/operations.h"
#include "symfpu/core/add.h"
#include "symfpu/core/sign.h"


#ifndef SYMFPU_REMAINDER
#define SYMFPU_REMAINDER

namespace symfpu {

template <class t>
  unpackedFloat<t> addRemainderSpecialCases (const typename t::fpt &format,
					    const unpackedFloat<t> &left,
					    const unpackedFloat<t> &right,
					    const unpackedFloat<t> &remainderResult) {
  typedef typename t::prop prop;

  prop eitherArgumentNan(left.getNaN() || right.getNaN());
  prop generateNan(left.getInf() || right.getZero());
  prop isNan(eitherArgumentNan || generateNan);

  prop passThrough((!(left.getInf() || left.getNaN()) && right.getInf()) ||
		   left.getZero());

  return ITE(isNan,
	     unpackedFloat<t>::makeNaN(format),
	     ITE(passThrough,
		 left,
		 remainderResult));
 }

/* Let left = x*2^e, right = y*2^f, x \in [1,2), y \in [1,2)
 * x/y \in (0.5,2)   x > y  x/y \in (1,2)   x < y (0.5,1)
 *
 *  rem =  x*2^e     - (y*2^f * int((x*2^e) / (y*2^f)))
 *      =  x*2^e     - (y*2^f * int((x/y) * 2^{e-f}))
 *      = (x*2^{e-f} - (y     * int((x/y) * 2^{e-f}))) * 2^f
 *
 *
 * If e - f >  0
 *      = (x*2^{e-f} - (y     * int((x*2^{e-f})/y)) * 2^f
 *
 *
 * If e - f == 0
 *      = (x         - (y     * int((x/y)          ))) * 2^f
 *      = ITTE(x ?= y,
 *             (x - (y * int[guard=1,sticky=1])) * 2^f
 *             (x -  y) * 2^f,
 *             ...)
 *      = ITTE(x ?= y,
 *             (x - (y * int[guard=1,sticky=1])) * 2^f
 *             left - right,
 *             ...)
 *
 *
 * If e - f == -1
 *      = (x*2^{-1}  - (y     * int((x/y) * 2^{-1 }))) * 2^f
 *      = ITTE(x ?= y,
 *             (x*2^{-1}  - (y * int[guard=0,sticky=1])) * 2^f
 *             (x*2^{-1}  - (y * int[guard=1,sticky=0])) * 2^f
 *             (x*2^{-1}  - (y * int[guard=1,sticky=1])) * 2^f
 *
 * If e - f <= -2
 *      = (x*2^{e-f}  - (y     * int[guard=0,sticky=1])) * 2^f
 *      = ITE(int[guard=0,sticky=1],
 *            x*2^e  - y*2^f,
 *            left)
 *      = ITE(int[guard=0,sticky=1],
 *            left  - right,
 *            left)
 *
 */

// Divide for max(e - f, 0) cycles
// The equal case, if you divide you use to extract the even bit of n, also save the rem.
// Then one more cycle for the guard bit
// Use that remainder to work out the sticky bit
// Round and either subtract or not from saved rem
// Output at 2^f

template <class t>
  unpackedFloat<t> arithmeticRemainder (const typename t::fpt &format,
					const typename t::rm &roundingMode,
					const unpackedFloat<t> &left,
					const unpackedFloat<t> &right) {
  typedef typename t::bwt bwt;
  typedef typename t::prop prop;
  typedef typename t::ubv ubv;
  typedef typename t::sbv sbv;
  //typedef typename t::fpt fpt;

  PRECONDITION(left.valid(format));
  PRECONDITION(right.valid(format));

  // Compute sign
  prop remainderSign(left.getSign());


  // Compute exponent difference
  sbv exponentDifference(expandingSubtract<t>(left.getExponent(), right.getExponent()));
  bwt edWidth(exponentDifference.getWidth());

  // Extend for divide steps
  ubv lsig(left.getSignificand().extend(1));
  ubv rsig(right.getSignificand().extend(1));


  ubv first(divideStep<t>(lsig,rsig).result);
  ubv *running = new ubv(first); // To avoid running out of stack space loop with a pointer

  bwt maxDifference = unpackedFloat<t>::maximumExponentDifference(format);
  for (bwt i = maxDifference - 1; i > 0; i--) {
    prop needPrevious(exponentDifference > sbv(edWidth, i));
    probabilityAnnotation<t>(needPrevious, (i > (maxDifference / 2)) ? VERYUNLIKELY : UNLIKELY);

    ubv r(ITE(needPrevious, *running, lsig));
    delete running;  // We assume the value / reference has been transfered to ITE
    running = new ubv(divideStep<t>(r, rsig).result);
  }

  // The zero exponent difference case is a little different
  // as we need the result bit for the even flag
  // and the actual result for the final
  prop lsbRoundActive(exponentDifference > -sbv::one(edWidth));  // i.e. >= 0

  prop needPrevious(exponentDifference > sbv::zero(edWidth));
  probabilityAnnotation<t>(needPrevious, UNLIKELY);

  ubv r0(ITE(needPrevious, *running, lsig));
  delete running;
  resultWithRemainderBit<t> dsr(divideStep<t>(r0, rsig));

  prop integerEven(!lsbRoundActive || !dsr.remainderBit);  // Note negation of guardBit


  // The same to get the guard flag
  prop guardRoundActive(exponentDifference > -sbv(edWidth,2));  // i.e. >= -1

  ubv rm1(ITE(lsbRoundActive, dsr.result, lsig));
  resultWithRemainderBit<t> dsrg(divideStep<t>(rm1, rsig));

  prop guardBit(guardRoundActive && dsrg.remainderBit);

  prop stickyBit(!ITE(guardRoundActive,
		      dsrg.result,
		      lsig).isAllZeros());


  // The base result if lsbRoundActive
  unpackedFloat<t> reconstruct(remainderSign,
			       right.getExponent(),
			       dsr.result.extract(lsig.getWidth() - 1,1)); // dsr shifts right as last action so is safe


  probabilityAnnotation<t>(needPrevious, UNLIKELY);   // Perhaps stretching it a bit but good for approximation
  unpackedFloat<t> candidateResult(ITE(lsbRoundActive,
				       reconstruct.normaliseUpDetectZero(format),
				       left));

  // The final subtract is a little different as previous ones were
  // guaranteed to be positive
  // TODO : This could be improved as these don't need special cases, etc.

  // From the rounding of the big integer multiple
  prop bonusSubtract(roundingDecision<t>(roundingMode,
					 remainderSign,
					 integerEven,
					 guardBit,
					 stickyBit,
					 prop(false)));
  probabilityAnnotation<t>(bonusSubtract, UNLIKELY); // Again, more like 50/50

  // The big integer has sign left.getSign() ^ right.getSign() so we subtract something of left.getSign().
  // For the integer part we handle this by working with absolutes (ignoring the sign) and
  // adding it back in at the end.
  // However for the correction for the rounded part we need to take it into account
  unpackedFloat<t> signCorrectedRight(right, left.getSign());
  unpackedFloat<t> remainderResult(ITE(bonusSubtract,
				       add<t>(format,
					      roundingMode,
					      candidateResult,
					      signCorrectedRight,
					      false),
				       candidateResult));

  // TODO : fast path if first.isAllZeros()

  POSTCONDITION(remainderResult.valid(format));

  return remainderResult;
 }


// Put it all together...
template <class t>
  unpackedFloat<t> remainderWithRounding (const typename t::fpt &format,
					  const typename t::rm &roundingMode,
					  const unpackedFloat<t> &left,
					  const unpackedFloat<t> &right) {
  //typedef typename t::bwt bwt;
  //typedef typename t::prop prop;
  //typedef typename t::ubv ubv;
  //typedef typename t::sbv sbv;

  PRECONDITION(left.valid(format));
  PRECONDITION(right.valid(format));

  unpackedFloat<t> remainderResult(arithmeticRemainder(format, roundingMode, left, right));

  //unpackedFloat<t> result(addRemainderSpecialCases(format, left, right, roundedRemainderResult));
  unpackedFloat<t> result(addRemainderSpecialCases(format, left, right, remainderResult));

  POSTCONDITION(result.valid(format));

  return result;
 }

// IEEE-754 remainder always uses round to nearest, ties to even
template <class t>
  unpackedFloat<t> remainder (const typename t::fpt &format,
			      const unpackedFloat<t> &left,
			      const unpackedFloat<t> &right) {

  return remainderWithRounding<t>(format, t::RNE(), left, right);
 }


}

#endif



Generated by: GCOVR (Version 3.2)

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	119	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	129	prop eitherArgumentNan(left.getNaN() \|\| right.getNaN());
39	129	prop generateNan(left.getInf() \|\| right.getZero());
40	129	prop isNan(eitherArgumentNan \|\| generateNan);
41
42	129	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	129	remainderResult));
50		}
51
52		/* Let left = x2^e, right = y2^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 = x2^e - (y2^f * int((x2^e) / (y2^f)))
56		* = x2^e - (y2^f * int((x/y) * 2^{e-f}))
57		* = (x2^{e-f} - (y int((x/y) * 2^{e-f}))) * 2^f
58		*
59		*
60		* If e - f > 0
61		* = (x2^{e-f} - (y int((x2^{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		* = (x2^{-1} - (y int((x/y) * 2^{-1 }))) * 2^f
78		* = ITTE(x ?= y,
79		* (x2^{-1} - (y int[guard=0,sticky=1])) * 2^f
80		* (x2^{-1} - (y int[guard=1,sticky=0])) * 2^f
81		* (x2^{-1} - (y int[guard=1,sticky=1])) * 2^f
82		*
83		* If e - f <= -2
84		* = (x2^{e-f} - (y int[guard=0,sticky=1])) * 2^f
85		* = ITE(int[guard=0,sticky=1],
86		* x2^e - y2^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	119	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	119	PRECONDITION(left.valid(format));
113	119	PRECONDITION(right.valid(format));
114
115		// Compute sign
116	129	prop remainderSign(left.getSign());
117
118
119		// Compute exponent difference
120	238	sbv exponentDifference(expandingSubtract<t>(left.getExponent(), right.getExponent()));
121	119	bwt edWidth(exponentDifference.getWidth());
122
123		// Extend for divide steps
124	238	ubv lsig(left.getSignificand().extend(1));
125	238	ubv rsig(right.getSignificand().extend(1));
126
127
128	238	ubv first(divideStep<t>(lsig,rsig).result);
129	119	ubv *running = new ubv(first); // To avoid running out of stack space loop with a pointer
130
131	119	bwt maxDifference = unpackedFloat<t>::maximumExponentDifference(format);
132	4878	for (bwt i = maxDifference - 1; i > 0; i--) {
133	5376	prop needPrevious(exponentDifference > sbv(edWidth, i));
134	4759	probabilityAnnotation<t>(needPrevious, (i > (maxDifference / 2)) ? VERYUNLIKELY : UNLIKELY);
135
136	9518	ubv r(ITE(needPrevious, *running, lsig));
137	4759	delete running; // We assume the value / reference has been transfered to ITE
138	4759	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	129	prop lsbRoundActive(exponentDifference > -sbv::one(edWidth)); // i.e. >= 0
145
146	129	prop needPrevious(exponentDifference > sbv::zero(edWidth));
147	119	probabilityAnnotation<t>(needPrevious, UNLIKELY);
148
149	238	ubv r0(ITE(needPrevious, *running, lsig));
150	119	delete running;
151	238	resultWithRemainderBit<t> dsr(divideStep<t>(r0, rsig));
152
153	129	prop integerEven(!lsbRoundActive \|\| !dsr.remainderBit); // Note negation of guardBit
154
155
156		// The same to get the guard flag
157	129	prop guardRoundActive(exponentDifference > -sbv(edWidth,2)); // i.e. >= -1
158
159	238	ubv rm1(ITE(lsbRoundActive, dsr.result, lsig));
160	238	resultWithRemainderBit<t> dsrg(divideStep<t>(rm1, rsig));
161
162	129	prop guardBit(guardRoundActive && dsrg.remainderBit);
163
164	129	prop stickyBit(!ITE(guardRoundActive,
165		dsrg.result,
166		lsig).isAllZeros());
167
168
169		// The base result if lsbRoundActive
170	238	unpackedFloat<t> reconstruct(remainderSign,
171		right.getExponent(),
172	119	dsr.result.extract(lsig.getWidth() - 1,1)); // dsr shifts right as last action so is safe
173
174
175	119	probabilityAnnotation<t>(needPrevious, UNLIKELY); // Perhaps stretching it a bit but good for approximation
176	238	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	129	prop bonusSubtract(roundingDecision<t>(roundingMode,
186		remainderSign,
187		integerEven,
188		guardBit,
189		stickyBit,
190		prop(false)));
191	119	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	238	unpackedFloat<t> signCorrectedRight(right, left.getSign());
198	119	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	119	POSTCONDITION(remainderResult.valid(format));
209
210	238	return remainderResult;
211		}
212
213
214		// Put it all together...
215		template <class t>
216	119	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	119	PRECONDITION(left.valid(format));
226	119	PRECONDITION(right.valid(format));
227
228	238	unpackedFloat<t> remainderResult(arithmeticRemainder(format, roundingMode, left, right));
229
230		//unpackedFloat<t> result(addRemainderSpecialCases(format, left, right, roundedRemainderResult));
231	119	unpackedFloat<t> result(addRemainderSpecialCases(format, left, right, remainderResult));
232
233	119	POSTCONDITION(result.valid(format));
234
235	238	return result;
236		}
237
238		// IEEE-754 remainder always uses round to nearest, ties to even
239		template <class t>
240	119	unpackedFloat<t> remainder (const typename t::fpt &format,
241		const unpackedFloat<t> &left,
242		const unpackedFloat<t> &right) {
243
244	119	return remainderWithRounding<t>(format, t::RNE(), left, right);
245		}
246
247
248		}
249
250		#endif
251