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 |
fpt extendedFormat(format.exponentWidth(), format.significandWidth() + 2); |
101 |
|
// format.exponentWidth() - 1 should also be true but requires shrinking the exponent and |
102 |
|
// then increasing it in the rounder |
103 |
27 |
unpackedFloat<t> sqrtResult(extendedFormat, sqrtSign, exponentHalved, finishedSignificand); |
104 |
|
|
105 |
27 |
POSTCONDITION(sqrtResult.valid(extendedFormat)); |
106 |
|
|
107 |
54 |
return sqrtResult; |
108 |
|
} |
109 |
|
|
110 |
|
|
111 |
|
// Put it all together... |
112 |
|
template <class t> |
113 |
27 |
unpackedFloat<t> sqrt (const typename t::fpt &format, |
114 |
|
const typename t::rm &roundingMode, |
115 |
|
const unpackedFloat<t> &uf) { |
116 |
|
//typedef typename t::bwt bwt; |
117 |
|
typedef typename t::prop prop; |
118 |
|
//typedef typename t::ubv ubv; |
119 |
|
//typedef typename t::sbv sbv; |
120 |
|
|
121 |
27 |
PRECONDITION(uf.valid(format)); |
122 |
|
|
123 |
54 |
unpackedFloat<t> sqrtResult(arithmeticSqrt(format, uf)); |
124 |
|
|
125 |
|
// Exponent is divided by two, thus it can't overflow, underflow or generate a subnormal number. |
126 |
|
// The last one is quite subtle but you can show that the largest number generatable |
127 |
|
// by arithmeticSqrt is 111...111:0:1 with the last two as the guard and sticky bits. |
128 |
|
// Round up (when the sign is positive) and round down (when the sign is negative -- |
129 |
|
// the result will be computed but then discarded) are the only cases when this can increment the significand. |
130 |
27 |
customRounderInfo<t> cri(prop(true), prop(true), prop(false), prop(true), |
131 |
100 |
!((roundingMode == t::RTP() && !sqrtResult.getSign()) || |
132 |
29 |
(roundingMode == t::RTN() && sqrtResult.getSign()))); |
133 |
54 |
unpackedFloat<t> roundedSqrtResult(customRounder(format, roundingMode, sqrtResult, cri)); |
134 |
|
|
135 |
27 |
unpackedFloat<t> result(addSqrtSpecialCases(format, uf, roundedSqrtResult.getSign(), roundedSqrtResult)); |
136 |
|
|
137 |
27 |
POSTCONDITION(result.valid(format)); |
138 |
|
|
139 |
54 |
return result; |
140 |
|
} |
141 |
|
|
142 |
|
|
143 |
|
} |
144 |
|
|
145 |
|
#endif |