GCC Code Coverage Report
Directory: . Exec Total Coverage
File: deps/install/include/symfpu/core/sqrt.h Lines: 35 35 100.0 %
Date: 2021-03-22 Branches: 65 153 42.5 %

Line Exec Source
1
/*
2
** Copyright (C) 2018 Martin Brain
3
**
4
** See the file LICENSE for licensing information.
5
*/
6
7
/*
8
** sqrt.h
9
**
10
** Martin Brain
11
** martin.brain@cs.ox.ac.uk
12
** 05/02/16
13
**
14
** Square root 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_SQRT
24
#define SYMFPU_SQRT
25
26
namespace symfpu {
27
28
template <class t>
29
27
  unpackedFloat<t> addSqrtSpecialCases (const typename t::fpt &format,
30
					  const unpackedFloat<t> &uf,
31
					  const typename t::prop &sign,
32
					  const unpackedFloat<t> &sqrtResult) {
33
  typedef typename t::prop prop;
34
35
27
  prop generateNaN(uf.getSign() && !uf.getZero());
36
27
  prop isNaN(uf.getNaN() || generateNaN);
37
38
27
  prop isInf(uf.getInf() && !uf.getSign());
39
40
27
  prop isZero(uf.getZero());
41
42
  return ITE(isNaN,
43
	     unpackedFloat<t>::makeNaN(format),
44
	     ITE(isInf,
45
		 unpackedFloat<t>::makeInf(format, prop(false)),
46
		 ITE(isZero,
47
		     unpackedFloat<t>::makeZero(format, sign),
48
27
		     sqrtResult)));
49
 }
50
51
52
 template <class t>
53
27
  unpackedFloat<t> arithmeticSqrt (const typename t::fpt &format,
54
				       const unpackedFloat<t> &uf) {
55
  typedef typename t::bwt bwt;
56
  typedef typename t::prop prop;
57
  typedef typename t::ubv ubv;
58
  typedef typename t::sbv sbv;
59
  typedef typename t::fpt fpt;
60
61
27
  PRECONDITION(uf.valid(format));
62
63
  // Compute sign
64
27
  prop sqrtSign(uf.getSign());
65
66
  // Divide the exponent by 2
67
54
  sbv exponent(uf.getExponent());
68
27
  bwt exponentWidth(exponent.getWidth());
69
27
  prop exponentEven((exponent & sbv::one(exponentWidth)).isAllZeros());
70
  #if 0
71
  sbv exponentHalved(conditionalDecrement<t,sbv,prop>((exponent < sbv::zero(exponentWidth)) && !exponentEven,
72
						      exponent.signExtendRightShift(sbv::one(exponentWidth))));
73
  #endif
74
54
  sbv exponentHalved(exponent.signExtendRightShift(sbv::one(exponentWidth)));
75
  // Right shift rounds down for positive, and away for negative  (-5 >>> 1 == -3)
76
  //  sqrt(1.s * 2^{-(2n + 1)}) = sqrt(1.s * 2^{-2n - 2 + 1)})
77
  //                            = sqrt(1.s * 2^{-2(n + 1)} * 2)
78
  //                            = sqrt(1.s * 2) * 2^{-(n + 1)}
79
  // Optimisation : improve the encoding of this operation
80
81
  // Sqrt the significands
82
  //  extend to allow alignment, pad so result has a guard bit
83
54
  ubv alignedSignificand(conditionalLeftShiftOne<t,ubv,prop>(!exponentEven, uf.getSignificand().extend(1).append(ubv::zero(1))));
84
85
54
  resultWithRemainderBit<t> sqrtd(fixedPointSqrt<t>(alignedSignificand));
86
87
88
27
  bwt resWidth(sqrtd.result.getWidth());
89
54
  ubv topBit(sqrtd.result.extract(resWidth - 1, resWidth - 1));
90
54
  ubv guardBit(sqrtd.result.extract(0,0));
91
92
  // Alignment of inputs means it is the range [1,4) so the result is in [1,2)
93
  // Also, the square root cannot be exactly between two numbers
94
27
  INVARIANT(topBit.isAllOnes());
95
27
  INVARIANT(IMPLIES(guardBit.isAllOnes(), sqrtd.remainderBit));
96
  // This also implies that no alignment of the exponent is needed
97
98
54
  ubv finishedSignificand(sqrtd.result.append(ubv(sqrtd.remainderBit)));
99
100
27
  unpackedFloat<t> sqrtResult(sqrtSign, exponentHalved, finishedSignificand);
101
102
103
27
  fpt extendedFormat(format.exponentWidth(), format.significandWidth() + 2);
104
  // format.exponentWidth() - 1 should also be true but requires shrinking the exponent and
105
  // then increasing it in the rounder
106
27
  POSTCONDITION(sqrtResult.valid(extendedFormat));
107
108
54
  return sqrtResult;
109
 }
110
111
112
// Put it all together...
113
template <class t>
114
27
  unpackedFloat<t> sqrt (const typename t::fpt &format,
115
			   const typename t::rm &roundingMode,
116
			   const unpackedFloat<t> &uf) {
117
  //typedef typename t::bwt bwt;
118
    typedef typename t::prop prop;
119
  //typedef typename t::ubv ubv;
120
  //typedef typename t::sbv sbv;
121
122
27
  PRECONDITION(uf.valid(format));
123
124
54
  unpackedFloat<t> sqrtResult(arithmeticSqrt(format, uf));
125
126
  // Exponent is divided by two, thus it can't overflow, underflow or generate a subnormal number.
127
  // The last one is quite subtle but you can show that the largest number generatable
128
  // by arithmeticSqrt is 111...111:0:1 with the last two as the guard and sticky bits.
129
  // Round up (when the sign is positive) and round down (when the sign is negative --
130
  // the result will be computed but then discarded) are the only cases when this can increment the significand.
131
27
  customRounderInfo<t> cri(prop(true), prop(true), prop(false), prop(true),
132
100
			   !((roundingMode == t::RTP() && !sqrtResult.getSign()) ||
133
29
			     (roundingMode == t::RTN() &&  sqrtResult.getSign())));
134
54
  unpackedFloat<t> roundedSqrtResult(customRounder(format, roundingMode, sqrtResult, cri));
135
136
27
  unpackedFloat<t> result(addSqrtSpecialCases(format, uf, roundedSqrtResult.getSign(), roundedSqrtResult));
137
138
27
  POSTCONDITION(result.valid(format));
139
140
54
  return result;
141
 }
142
143
144
}
145
146
#endif