GCC Code Coverage Report
Directory: . Exec Total Coverage
File: deps/install/include/symfpu/core/operations.h Lines: 182 187 97.3 %
Date: 2021-03-22 Branches: 329 677 48.6 %

Line Exec Source
1
/*
2
** Copyright (C) 2018 Martin Brain
3
**
4
** See the file LICENSE for licensing information.
5
*/
6
7
/*
8
** operations.h
9
**
10
** Martin Brain
11
** martin.brain@cs.ox.ac.uk
12
** 15/11/14
13
**
14
** A number of compound operations on bit-vectors.  These are default
15
** implementations to reduce the amount of code that is needed for a
16
** back-end (and the risk of making a mistake).  Back-ends can provide
17
** their own implementations of these if they can handle them in a
18
** smarter way than the default.
19
**
20
*/
21
22
#include <cassert>
23
#include <map>
24
25
#include "../utils/common.h"
26
27
28
#ifndef SYMFPU_OPERATIONS
29
#define SYMFPU_OPERATIONS
30
31
namespace symfpu {
32
33
  /*** Expanding operations ***/
34
  template <class t, class bv>
35
  bv expandingAdd (const bv &op1, const bv &op2) {
36
    PRECONDITION(op1.getWidth() == op2.getWidth());
37
38
    bv x(op1.extend(1));
39
    bv y(op2.extend(1));
40
41
    return x + y;
42
  }
43
44
  template <class t, class bv, class prop>
45
141
    bv expandingAddWithCarryIn (const bv &op1, const bv &op2, const prop &cin) {
46
141
    PRECONDITION(op1.getWidth() == op2.getWidth());
47
48
282
    bv x(op1.extend(1));
49
282
    bv y(op2.extend(1));
50
51
282
    bv sum(x + y);
52
53
141
    typename t::bwt w(sum.getWidth());
54
282
    bv carry(ITE(cin, bv::one(w), bv::zero(w)));
55
141
    bv res(sum.modularAdd(carry)); // Modular is safe due to the extension
56
                                   // (2^n - 1) + (2^n - 1) + 1 == 2^(n+1) - 1
57
                                   // -(2^n) + -(2^n) + 1 > 2^(n+1)
58
282
    return res;
59
  }
60
61
  template <class t, class bv>
62
1986
  bv expandingSubtract (const bv &op1, const bv &op2) {
63
1986
    PRECONDITION(op1.getWidth() == op2.getWidth());
64
65
3972
    bv x(op1.extend(1));
66
3972
    bv y(op2.extend(1));
67
68
3972
    return x - y;
69
  }
70
71
  template <class t, class bv>
72
465
  bv expandingMultiply (const bv &op1, const bv &op2) {
73
465
    typename t::bwt width = op1.getWidth();
74
465
    PRECONDITION(width == op2.getWidth());
75
76
930
    bv x(op1.extend(width));
77
930
    bv y(op2.extend(width));
78
79
930
    return x * y;
80
  }
81
82
83
  /*** Conditional Operations ***/
84
  template <class t, class bv, class prop>
85
1299
  bv conditionalIncrement (const prop &p, const bv &b) {
86
1299
    PRECONDITION(IMPLIES(p, b <  bv::maxValue(b.getWidth())));
87
88
1299
    typename t::bwt w(b.getWidth());
89
2598
    bv inc(ITE(p, bv::one(w), bv::zero(w)));
90
91
2598
    return b + inc;
92
  }
93
94
  template <class t, class bv, class prop>
95
194
  bv conditionalDecrement (const prop &p, const bv &b) {
96
194
    PRECONDITION(IMPLIES(p, bv::minValue(b.getWidth()) < b));
97
98
194
    typename t::bwt w(b.getWidth());
99
388
    bv inc(ITE(p, bv::one(w), bv::zero(w)));
100
101
388
    return b - inc;
102
  }
103
104
  template <class t, class bv, class prop>
105
1123
  bv conditionalLeftShiftOne (const prop &p, const bv &b) {
106
1123
    typename t::bwt w(b.getWidth());
107
1123
    PRECONDITION(IMPLIES(p, (b.extract(w - 1, w - 1).isAllZeros())));
108
109
2246
    bv shifted(b.modularLeftShift(bv::one(w)));
110
2246
    return bv(ITE(p, shifted, b));
111
  }
112
113
  template <class t, class bv, class prop>
114
761
  bv conditionalRightShiftOne (const prop &p, const bv &b) {
115
761
    typename t::bwt w(b.getWidth());
116
    // PRECONDITION(IMPLIES(p, (b.extract(0, 0).isAllZeros())));  // Adder uses and compensates for this case.
117
118
1522
    bv shifted(b.modularRightShift(bv::one(w)));
119
1522
    return bv(ITE(p, shifted, b));
120
  }
121
122
  template <class t, class bv, class prop>
123
765
  bv conditionalNegate (const prop &p, const bv &b) {
124
765
    typename t::bwt w(b.getWidth());
125
765
    PRECONDITION(w >= 2);
126
765
    PRECONDITION(IMPLIES(p, !(b.extract(w - 1, w - 1).isAllOnes() &&
127
			      b.extract(w - 2,     0).isAllZeros())));
128
129
765
    return bv(ITE(p, -b, b));
130
  }
131
132
  template <class t, class bv>
133
4
  bv abs(const bv &b) {
134
4
    return conditionalNegate<t, bv, typename t::prop>(b < bv::zero(b.getWidth()), b);
135
  }
136
137
138
139
  /*** Probability Annotations ***/
140
  enum probability {
141
    VERYLIKELY = 100,
142
    LIKELY = 50,
143
    NEUTRAL = 0,
144
    UNLIKELY = -50,
145
    VERYUNLIKELY = -100
146
  };
147
148
  template <class t, class prop>
149
36868
  void probabilityAnnotation (const prop &, const probability &) {
150
    // Back-ends can make use of this information if they want
151
36868
    return;
152
  }
153
154
155
  /*** Max and min ***/
156
  template <class t, class bv>
157
32709
  bv max (const bv &op1, const bv &op2) {
158
32709
    return ITE(op1 <= op2, op2, op1);
159
  }
160
161
  template <class t, class bv>
162
  bv min (const bv &op1, const bv &op2) {
163
    return ITE(op1 <= op2, op1, op2);
164
  }
165
166
  template <class t, class bv>
167
1299
  bv collar(const bv &op, const bv &lower, const bv &upper) {
168
2476
    return ITE(op < lower,
169
	       lower,
170
2476
	       ITE(upper < op,
171
		   upper,
172
5013
		   op));
173
  }
174
175
176
  /*** Unary/Binary operations ***/
177
  template <class t, class bv, class prop, class bwt>
178
  bv countLeadingZerosRec (const bv &op, const bwt position, const prop &allPreceedingZeros) {
179
    typename t::bwt w(op.getWidth());
180
181
    PRECONDITION(0 <= position && position < w);
182
183
    bv bit(op.extract(position, position));
184
185
    prop isLeadingOne(allPreceedingZeros && (bit.isAllOnes()));
186
    prop continuingZero(allPreceedingZeros && (bit.isAllZeros()));
187
188
    if (position == 0) {
189
      return ITE(isLeadingOne, bv(w, w - 1), bv(w, w));
190
    } else {
191
      return ITE(isLeadingOne,
192
		 bv(w, w - (position + 1)),
193
		 countLeadingZerosRec<t>(op, position - 1, continuingZero));
194
    }
195
  }
196
197
  template <class t, class bv>
198
  bv countLeadingZeros (const bv &op) {
199
    typedef typename t::bwt bwt;
200
    typedef typename t::prop prop;
201
    bwt w(op.getWidth());
202
203
    return countLeadingZerosRec<t>(op, w - 1, prop(true));
204
  }
205
206
  // This is sort of the opposite of count trailing 1's (a.k.a. clz(reverse(not(x))) )
207
  template <class t, class bv>
208
22687
  bv orderEncode (const bv &op) {
209
22687
    typename t::bwt w(op.getWidth());
210
211
    //PRECONDITION(bv::zero(w) <= op && op <= bv(w, w)); // Not needed as using modular shift
212
213
22687
    bv tmp((bv::one(w + 1).modularLeftShift(op.resize(w + 1))).modularDecrement().extract(w-1,0));
214
22687
    return tmp;
215
  }
216
217
218
  // The comparison we need to do is
219
  //   op.extract(relevantBits - 1, 0) == bv(relevantBits, position + 1)
220
  // bits can be shared between instances of this.
221
  // The following is a dynamic programming style solution.
222
  // HOWEVER : it may not actually give a net saving depending on your application
223
  template <class t, class bv>
224
5075
  struct fragmentMap {
225
    typedef typename t::bwt bwt;
226
    typedef typename t::prop prop;
227
228
    typedef std::pair<bwt, bwt> fragment;
229
    typedef std::map<fragment, prop> map;
230
231
  protected :
232
    const bv &op;
233
    map m;
234
235
    prop getComparitorRec(bwt length, bwt value)
236
    {
237
      PRECONDITION(length > 0);
238
      PRECONDITION(bitsToRepresent(value) <= length);
239
      typename map::const_iterator it = m.find(std::make_pair(length, value));
240
241
      if (it != m.end()) {
242
	return it->second;
243
      } else {
244
	bwt leadingBit = bwt(1) << (length - 1);
245
	prop leadingBitIsOne(op.extract(length - 1, length - 1).isAllOnes());
246
	prop correctComparison((value & leadingBit) ?
247
			       leadingBitIsOne : !leadingBitIsOne);
248
	prop *step = NULL;
249
250
	if (length == 1) {
251
	  step = new prop(correctComparison);
252
	} else {
253
	  prop rec(getComparitorRec(length - 1, value & (~leadingBit)));
254
	  step = new prop(correctComparison && rec);
255
	}
256
257
	prop res(*step);
258
	delete step;
259
260
	m.insert(std::make_pair(std::make_pair(length, value), res));
261
	return res;
262
      }
263
    }
264
265
  public :
266
5075
    fragmentMap(const bv &_op) : op(_op) {}
267
268
    prop getComparitor(bwt length, bwt value)
269
    {
270
      PRECONDITION(length > 0);
271
      PRECONDITION(bitsToRepresent(value) <= length);
272
273
      prop res(getComparitorRec(length, value));
274
275
      POSTCONDITION(bv(res) == (op.extract(length - 1, 0) == bv(length, value)));
276
      return res;
277
    }
278
  };
279
280
281
  // A more compact, bitwise implementation of orderEncode for SAT encoding
282
  // Intended to be used in specialisation of orderEncode
283
  template <class t, class bv>
284
5075
  bv orderEncodeBitwise (const bv &op) {
285
    typedef typename t::bwt bwt;
286
5075
    bwt w(op.getWidth());
287
288
10150
    fragmentMap<t,bv> m(op);
289
290
    // If op is too large, then set everything to 1
291
10150
    bv outOfRange(op >= bv(w, w));
292
293
    // SAND to fill in the remaining bits
294
5075
    bv * working = new bv(outOfRange);
295
58060
    for (bwt i = w; i > 0; --i) {
296
52985
      bwt position = i - 1;    // Position in the output bitvectors
297
52985
      bwt relevantBits = bitsToRepresent(position + 1);
298
52985
      INVARIANT(relevantBits > 0);
299
300
      //bv activateBit(m.getComparitor(relevantBits, position + 1));  // No more compact and slower
301
105970
      bv activateBit(op.extract(relevantBits - 1, 0) == bv(relevantBits, position + 1));
302
105970
      bv nextBit(working->extract(0,0) | activateBit);
303
304
52985
      bv * tmp = working;
305
52985
      working = new bv(working->append(nextBit));
306
52985
      delete tmp;
307
    }
308
309
5075
    bv output(working->extract(w - 1,0));
310
5075
    delete working;
311
312
5075
    POSTCONDITION(output == (bv::one(w + 1).modularLeftShift(op.resize(w + 1))).modularDecrement().extract(w-1,0));
313
314
10150
    return output;
315
  }
316
317
318
  /*** Custom shifts ***/
319
  // 1 if and only if the right shift moves at least one 1 out of the word
320
  template <class t, class bv>
321
761
  bv rightShiftStickyBit (const bv &op, const bv &shift) {
322
2283
    bv stickyBit(ITE((orderEncode<t>(shift) & op).isAllZeros(),
323
1522
		     bv::zero(op.getWidth()),
324
1522
		     bv::one(op.getWidth())));
325
326
761
    return stickyBit;
327
  }
328
329
330
  // It is easier to compute along with the shift
331
  template <class t>
332
761
  struct stickyRightShiftResult {
333
    typedef typename t::ubv ubv;
334
335
    ubv signExtendedResult;
336
    ubv stickyBit;
337
338
761
    stickyRightShiftResult(const ubv &ser, const ubv &sb) : signExtendedResult(ser), stickyBit(sb) {}
339
  stickyRightShiftResult(const stickyRightShiftResult &old) : signExtendedResult(old.signExtendedResult), stickyBit(old.stickyBit) {}
340
  };
341
342
343
  template <class t>
344
735
  stickyRightShiftResult<t> stickyRightShift (const typename t::ubv &input, const typename t::ubv &shiftAmount) {
345
735
    stickyRightShiftResult<t> res(input.signExtendRightShift(shiftAmount), rightShiftStickyBit<t>(input, shiftAmount));
346
347
735
    return res;
348
  }
349
350
351
  template <class t>
352
26
  stickyRightShiftResult<t> stickyRightShiftBitwise (const typename t::ubv &input, const typename t::ubv &shiftAmount) {
353
    typedef typename t::bwt bwt;
354
    typedef typename t::prop prop;
355
    typedef typename t::ubv ubv;
356
357
26
    bwt width(input.getWidth());
358
26
    bwt startingPosition(positionOfLeadingOne(width));
359
26
    INVARIANT(0 < startingPosition && startingPosition < width);
360
361
    // Catch the out of bounds case
362
26
    PRECONDITION(shiftAmount.getWidth() == width);
363
52
    prop fullShift(shiftAmount >= ubv(width, width));
364
    // Note the shiftAmount is treated as unsigned...
365
366
26
    ubv *working = new ubv(input);
367
26
    prop *stickyBit = new prop(ITE(fullShift, !input.isAllZeros(), prop(false)));
368
369
147
    for (bwt i = startingPosition + 1; i > 0; --i)
370
    {
371
121
      bwt shiftAmountPosition = i - 1;
372
373
242
      prop shiftEnabled = fullShift || shiftAmount.extract(shiftAmountPosition, shiftAmountPosition).isAllOnes();
374
375
242
      prop stickyAccumulate(shiftEnabled && !(working->extract((1ULL << shiftAmountPosition) - 1, 0).isAllZeros())); // Would this loose data?
376
377
121
      prop * tmp = stickyBit;
378
121
      stickyBit = new prop(*stickyBit || stickyAccumulate);
379
121
      delete tmp;
380
381
382
      // Note the slightly unexpected sign extension
383
242
      ubv shifted(working->signExtendRightShift(ubv::one(width) << ubv(width, shiftAmountPosition)));
384
385
121
      ubv * workingTmp = working;
386
121
      working = new ubv(ITE(shiftEnabled, shifted, *working));
387
121
      delete workingTmp;
388
    }
389
390
26
    stickyRightShiftResult<t> res(*working, ubv(*stickyBit).extend(width - 1));
391
392
26
    delete working;
393
26
    delete stickyBit;
394
395
26
    POSTCONDITION(res.signExtendedResult == input.signExtendRightShift(shiftAmount));
396
26
    POSTCONDITION(res.stickyBit == rightShiftStickyBit<t>(input, shiftAmount));
397
398
52
    return res;
399
  }
400
401
402
403
  template <class t>
404
1197
  struct normaliseShiftResult {
405
    typedef typename t::ubv ubv;
406
    typedef typename t::prop prop;
407
408
    ubv normalised;
409
    ubv shiftAmount;
410
    prop isZero;
411
412
1197
  normaliseShiftResult(const ubv &n, const ubv &s, const prop &z) : normalised(n), shiftAmount(s), isZero(z) {}
413
  normaliseShiftResult(const normaliseShiftResult<t> &old) :  normalised(old.normalised), shiftAmount(old.shiftAmount), isZero(old.isZero) {}
414
  };
415
416
  template <class t>
417
1197
  normaliseShiftResult<t> normaliseShift (const typename t::ubv input) {
418
    typedef typename t::bwt bwt;
419
    typedef typename t::prop prop;
420
    typedef typename t::ubv ubv;
421
422
1197
    bwt width(input.getWidth());
423
1197
    bwt startingMask(previousPowerOfTwo(width));
424
1197
    INVARIANT(startingMask < width);
425
426
    // Catch the zero case
427
1247
    prop zeroCase(input.isAllZeros());
428
429
1197
    ubv *working = new ubv(input);
430
1197
    ubv *shiftAmount = NULL;
431
1197
    prop *deactivateShifts = new prop(zeroCase);
432
433
5662
    for (bwt i = startingMask; i > 0; i >>= 1) {
434
4709
      prop newDeactivateShifts = *deactivateShifts || working->extract(width-1,width-1).isAllOnes();
435
4465
      delete deactivateShifts;
436
4465
      deactivateShifts = new prop(newDeactivateShifts);
437
438
8930
      ubv mask(ubv::allOnes(i).append(ubv::zero(width - i)));
439
4709
      prop shiftNeeded(!(*deactivateShifts) && (mask & *working).isAllZeros());
440
441
      // Modular is safe because of the mask comparison
442
8930
      ubv shifted(ITE(shiftNeeded, working->modularLeftShift(ubv(width, i)), *working));
443
4465
      delete working;
444
4465
      working = new ubv(shifted);
445
446
4465
      if (shiftAmount == NULL) {
447
1197
	shiftAmount = new ubv(shiftNeeded);
448
      } else {
449
6536
	ubv newShiftAmount = shiftAmount->append(ubv(shiftNeeded));
450
3268
	delete shiftAmount;
451
3268
	shiftAmount = new ubv(newShiftAmount);
452
      }
453
    }
454
455
1197
    normaliseShiftResult<t> res(*working, *shiftAmount, zeroCase);
456
457
1197
    delete deactivateShifts;
458
1197
    delete working;
459
1197
    delete shiftAmount;
460
461
1197
    POSTCONDITION(res.normalised.extract(width-1,width-1).isAllZeros() == res.isZero);
462
1197
    POSTCONDITION(IMPLIES(res.isZero, res.shiftAmount.isAllZeros()));
463
464
1197
    bwt shiftAmountWidth(res.shiftAmount.getWidth());
465
1197
    bwt widthBits(bitsToRepresent(width));
466
1197
    POSTCONDITION(shiftAmountWidth == widthBits ||
467
		  shiftAmountWidth == widthBits - 1); // If width is an exact power of 2
468
2394
    ubv widthBV(widthBits, width);
469
1197
    POSTCONDITION(res.shiftAmount.matchWidth(widthBV) < widthBV);
470
471
2394
    return res;
472
  }
473
474
475
  /*** Dividers ***/
476
  template <class t>
477
4763
  struct resultWithRemainderBit {
478
    typedef typename t::ubv ubv;
479
    typedef typename t::prop prop;
480
481
    ubv result;
482
    prop remainderBit;
483
484
4763
  resultWithRemainderBit(const ubv &o, const prop &r) : result(o), remainderBit(r) {}
485
  resultWithRemainderBit(const resultWithRemainderBit<t> &old) : result(old.result), remainderBit(old.remainderBit) {}
486
  };
487
488
  // x and y are fixed-point numbers in the range [1,2)
489
  // Compute o \in [0.5,2), r \in [0,\delta) such that:  x = o*y + r
490
  // Return (o, r != 0)
491
  template <class t>
492
194
  resultWithRemainderBit<t> fixedPointDivide (const typename t::ubv &x, const typename t::ubv &y) {
493
194
    typename t::bwt w(x.getWidth());
494
495
    // Same width and both have MSB ones
496
194
    PRECONDITION(y.getWidth() == w);
497
194
    PRECONDITION(x.extract(w - 1, w - 1).isAllOnes());
498
194
    PRECONDITION(y.extract(w - 1, w - 1).isAllOnes());
499
500
    typedef typename t::ubv ubv;
501
502
    // Not the best way of doing this but pretty universal
503
388
    ubv ex(x.append(ubv::zero(w - 1)));
504
388
    ubv ey(y.extend(w - 1));
505
506
388
    ubv div(ex / ey);
507
388
    ubv rem(ex % ey);
508
509
388
    return resultWithRemainderBit<t>(div.extract(w - 1, 0), !(rem.isAllZeros()));
510
  }
511
512
513
  // x is a fixed-point number in the range [1,4) with 2/p bits
514
  // Compute o \in [1,sqrt(2)), r \in [0,o*2 + 1) such that x = o*o + r with 1/p bits
515
  // Return (o, r != 0)
516
  template <class t>
517
27
  resultWithRemainderBit<t> fixedPointSqrt (const typename t::ubv &x) {
518
    typedef typename t::bwt bwt;
519
    typedef typename t::ubv ubv;
520
    typedef typename t::prop prop;
521
522
    // The default algorithm given here isn't a great one
523
    // However it is simple and has a very simple termination criteria.
524
    // Plus most symbolic back ends will prefer
525
    // o = nondet(), r = nondet(), assert(r < o*2 + 1), assert(x = o*o + r)
526
527
27
    bwt inputWidth(x.getWidth());
528
27
    bwt outputWidth(inputWidth - 1);
529
530
    // To compare against, we need to pad x to 2/2p
531
54
    ubv xcomp(x.append(ubv::zero(inputWidth - 2)));
532
533
    // Start at 1
534
54
    ubv working(ubv::one(outputWidth) << ubv(outputWidth, outputWidth - 1));
535
536
    bwt location;
537
324
    for (location = outputWidth - 1; location > 0; --location) { // Offset by 1 for easy termination
538
594
      ubv shift(ubv(outputWidth, location - 1));
539
540
594
      ubv candidate(working | (ubv::one(outputWidth) << shift));
541
542
297
      prop addBit(expandingMultiply<t, ubv>(candidate, candidate) <= xcomp);
543
544
297
      working = working | (ubv(addBit).extend(outputWidth - 1) << shift);
545
    }
546
547
54
    return resultWithRemainderBit<t>(working, !(expandingMultiply<t, ubv>(working, working) == xcomp));
548
  }
549
550
  // One step of a divider
551
  // Here the "remainder bit" is actual the result bit and
552
  // The result is the remainder
553
  template <class t>
554
4542
  resultWithRemainderBit<t> divideStep (const typename t::ubv &x, const typename t::ubv &y) {
555
    typedef typename t::bwt bwt;
556
    typedef typename t::ubv ubv;
557
    typedef typename t::prop prop;
558
559
4542
    bwt xWidth(x.getWidth());
560
4542
    bwt yWidth(y.getWidth());
561
562
4542
    PRECONDITION(xWidth == yWidth);
563
4542
    PRECONDITION(yWidth >= 2);
564
4542
    PRECONDITION(y.extract(yWidth - 2, yWidth - 2).isAllOnes());  // Assume y is aligned
565
566
5189
    prop canSubtract(x >= y);
567
9084
    ubv sub(x.modularAdd(y.modularNegate())); // TODO : modular subtract or better
568
9084
    ubv step(ITE(canSubtract, sub, x));
569
570
9084
    return resultWithRemainderBit<t>(step << ubv::one(xWidth), canSubtract);
571
  }
572
}
573
574
#endif