Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	build-coverage/deps/include/symfpu/core/operations.h	Lines:	182	187	97.3 %
Date:	2021-05-22	Branches:	331	677	48.9 %


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

/*
** operations.h
**
** Martin Brain
** martin.brain@cs.ox.ac.uk
** 15/11/14
**
** A number of compound operations on bit-vectors.  These are default
** implementations to reduce the amount of code that is needed for a
** back-end (and the risk of making a mistake).  Back-ends can provide
** their own implementations of these if they can handle them in a
** smarter way than the default.
**
*/

#include <cassert>
#include <map>

#include "../utils/common.h"


#ifndef SYMFPU_OPERATIONS
#define SYMFPU_OPERATIONS

namespace symfpu {

  /*** Size matching ***/
  template <class t, class bv>
  std::pair<bv, bv> sameWidth(const bv &op1, const bv &op2) {
    typename t::bwt w1(op1.getWidth());
    typename t::bwt w2(op2.getWidth());

    if (w1 == w2) {
      return std::make_pair(op1, op2);
    } else if (w1 < w2) {
      return std::make_pair(op1.matchWidth(op2), op2);
    } else {
      INVARIANT(w1 > w2);
      return std::make_pair(op1.matchWidth(op2), op2);
    }
  }



  /*** Expanding operations ***/
  template <class t, class bv>
  bv expandingAdd (const bv &op1, const bv &op2) {
    PRECONDITION(op1.getWidth() == op2.getWidth());

    bv x(op1.extend(1));
    bv y(op2.extend(1));

    return x + y;
  }

  template <class t, class bv, class prop>
    bv expandingAddWithCarryIn (const bv &op1, const bv &op2, const prop &cin) {
    PRECONDITION(op1.getWidth() == op2.getWidth());

    bv x(op1.extend(1));
    bv y(op2.extend(1));

    bv sum(x + y);

    typename t::bwt w(sum.getWidth());
    bv carry(ITE(cin, bv::one(w), bv::zero(w)));
    bv res(sum.modularAdd(carry)); // Modular is safe due to the extension
                                   // (2^n - 1) + (2^n - 1) + 1 == 2^(n+1) - 1
                                   // -(2^n) + -(2^n) + 1 > 2^(n+1)
    return res;
  }

  template <class t, class bv>
  bv expandingSubtract (const bv &op1, const bv &op2) {
    PRECONDITION(op1.getWidth() == op2.getWidth());

    bv x(op1.extend(1));
    bv y(op2.extend(1));

    return x - y;
  }

  template <class t, class bv>
  bv expandingMultiply (const bv &op1, const bv &op2) {
    typename t::bwt width = op1.getWidth();
    PRECONDITION(width == op2.getWidth());

    bv x(op1.extend(width));
    bv y(op2.extend(width));

    return x * y;
  }


  /*** Conditional Operations ***/
  template <class t, class bv, class prop>
  bv conditionalIncrement (const prop &p, const bv &b) {
    PRECONDITION(IMPLIES(p, b <  bv::maxValue(b.getWidth())));

    typename t::bwt w(b.getWidth());
    bv inc(ITE(p, bv::one(w), bv::zero(w)));

    return b + inc;
  }

  template <class t, class bv, class prop>
  bv conditionalDecrement (const prop &p, const bv &b) {
    PRECONDITION(IMPLIES(p, bv::minValue(b.getWidth()) < b));

    typename t::bwt w(b.getWidth());
    bv inc(ITE(p, bv::one(w), bv::zero(w)));

    return b - inc;
  }

  template <class t, class bv, class prop>
  bv conditionalLeftShiftOne (const prop &p, const bv &b) {
    typename t::bwt w(b.getWidth());
    PRECONDITION(IMPLIES(p, (b.extract(w - 1, w - 1).isAllZeros())));

    bv shifted(b.modularLeftShift(bv::one(w)));
    return bv(ITE(p, shifted, b));
  }

  template <class t, class bv, class prop>
  bv conditionalRightShiftOne (const prop &p, const bv &b) {
    typename t::bwt w(b.getWidth());
    // PRECONDITION(IMPLIES(p, (b.extract(0, 0).isAllZeros())));  // Adder uses and compensates for this case.

    bv shifted(b.modularRightShift(bv::one(w)));
    return bv(ITE(p, shifted, b));
  }

  template <class t, class bv, class prop>
  bv conditionalNegate (const prop &p, const bv &b) {
    typename t::bwt w(b.getWidth());
    PRECONDITION(w >= 2);
    PRECONDITION(IMPLIES(p, !(b.extract(w - 1, w - 1).isAllOnes() &&
			      b.extract(w - 2,     0).isAllZeros())));

    return bv(ITE(p, -b, b));
  }

  template <class t, class bv>
  bv abs(const bv &b) {
    return conditionalNegate<t, bv, typename t::prop>(b < bv::zero(b.getWidth()), b);
  }



  /*** Probability Annotations ***/
  enum probability {
    VERYLIKELY = 100,
    LIKELY = 50,
    NEUTRAL = 0,
    UNLIKELY = -50,
    VERYUNLIKELY = -100
  };

  template <class t, class prop>
  void probabilityAnnotation (const prop &, const probability &) {
    // Back-ends can make use of this information if they want
    return;
  }


  /*** Max and min ***/
  template <class t, class bv>
  bv max (const bv &op1, const bv &op2) {
    return ITE(op1 <= op2, op2, op1);
  }

  template <class t, class bv>
  bv min (const bv &op1, const bv &op2) {
    return ITE(op1 <= op2, op1, op2);
  }

  template <class t, class bv>
  bv collar(const bv &op, const bv &lower, const bv &upper) {
    return ITE(op < lower,
	       lower,
	       ITE(upper < op,
		   upper,
		   op));
  }


  /*** Unary/Binary operations ***/
  template <class t, class bv, class prop, class bwt>
  bv countLeadingZerosRec (const bv &op, const bwt position, const prop &allPreceedingZeros) {
    typename t::bwt w(op.getWidth());

    PRECONDITION(0 <= position && position < w);

    bv bit(op.extract(position, position));

    prop isLeadingOne(allPreceedingZeros && (bit.isAllOnes()));
    prop continuingZero(allPreceedingZeros && (bit.isAllZeros()));

    if (position == 0) {
      return ITE(isLeadingOne, bv(w, w - 1), bv(w, w));
    } else {
      return ITE(isLeadingOne,
		 bv(w, w - (position + 1)),
		 countLeadingZerosRec<t>(op, position - 1, continuingZero));
    }
  }

  template <class t, class bv>
  bv countLeadingZeros (const bv &op) {
    typedef typename t::bwt bwt;
    typedef typename t::prop prop;
    bwt w(op.getWidth());

    return countLeadingZerosRec<t>(op, w - 1, prop(true));
  }

  // This is sort of the opposite of count trailing 1's (a.k.a. clz(reverse(not(x))) )
  template <class t, class bv>
  bv orderEncode (const bv &op) {
    typename t::bwt w(op.getWidth());

    //PRECONDITION(bv::zero(w) <= op && op <= bv(w, w)); // Not needed as using modular shift

    bv tmp((bv::one(w + 1).modularLeftShift(op.resize(w + 1))).modularDecrement().extract(w-1,0));
    return tmp;
  }


  // The comparison we need to do is
  //   op.extract(relevantBits - 1, 0) == bv(relevantBits, position + 1)
  // bits can be shared between instances of this.
  // The following is a dynamic programming style solution.
  // HOWEVER : it may not actually give a net saving depending on your application
  template <class t, class bv>
  struct fragmentMap {
    typedef typename t::bwt bwt;
    typedef typename t::prop prop;

    typedef std::pair<bwt, bwt> fragment;
    typedef std::map<fragment, prop> map;

  protected :
    const bv &op;
    map m;

    prop getComparitorRec(bwt length, bwt value)
    {
      PRECONDITION(length > 0);
      PRECONDITION(bitsToRepresent(value) <= length);
      typename map::const_iterator it = m.find(std::make_pair(length, value));

      if (it != m.end()) {
	return it->second;
      } else {
	bwt leadingBit = bwt(1) << (length - 1);
	prop leadingBitIsOne(op.extract(length - 1, length - 1).isAllOnes());
	prop correctComparison((value & leadingBit) ?
			       leadingBitIsOne : !leadingBitIsOne);
	prop *step = NULL;

	if (length == 1) {
	  step = new prop(correctComparison);
	} else {
	  prop rec(getComparitorRec(length - 1, value & (~leadingBit)));
	  step = new prop(correctComparison && rec);
	}

	prop res(*step);
	delete step;

	m.insert(std::make_pair(std::make_pair(length, value), res));
	return res;
      }
    }

  public :
    fragmentMap(const bv &_op) : op(_op) {}

    prop getComparitor(bwt length, bwt value)
    {
      PRECONDITION(length > 0);
      PRECONDITION(bitsToRepresent(value) <= length);

      prop res(getComparitorRec(length, value));

      POSTCONDITION(bv(res) == (op.extract(length - 1, 0) == bv(length, value)));
      return res;
    }
  };


  // A more compact, bitwise implementation of orderEncode for SAT encoding
  // Intended to be used in specialisation of orderEncode
  template <class t, class bv>
  bv orderEncodeBitwise (const bv &op) {
    typedef typename t::bwt bwt;
    bwt w(op.getWidth());

    fragmentMap<t,bv> m(op);

    // If op is too large, then set everything to 1
    bv outOfRange(op >= bv(w, w));

    // SAND to fill in the remaining bits
    bv * working = new bv(outOfRange);
    for (bwt i = w; i > 0; --i) {
      bwt position = i - 1;    // Position in the output bitvectors
      bwt relevantBits = bitsToRepresent(position + 1);
      INVARIANT(relevantBits > 0);

      //bv activateBit(m.getComparitor(relevantBits, position + 1));  // No more compact and slower
      bv activateBit(op.extract(relevantBits - 1, 0) == bv(relevantBits, position + 1));
      bv nextBit(working->extract(0,0) | activateBit);

      bv * tmp = working;
      working = new bv(working->append(nextBit));
      delete tmp;
    }

    bv output(working->extract(w - 1,0));
    delete working;

    POSTCONDITION(output == (bv::one(w + 1).modularLeftShift(op.resize(w + 1))).modularDecrement().extract(w-1,0));

    return output;
  }


  /*** Custom shifts ***/
  // 1 if and only if the right shift moves at least one 1 out of the word
  template <class t, class bv>
  bv rightShiftStickyBit (const bv &op, const bv &shift) {
    bv stickyBit(ITE((orderEncode<t>(shift) & op).isAllZeros(),
		     bv::zero(op.getWidth()),
		     bv::one(op.getWidth())));

    return stickyBit;
  }


  // It is easier to compute along with the shift
  template <class t>
  struct stickyRightShiftResult {
    typedef typename t::ubv ubv;

    ubv signExtendedResult;
    ubv stickyBit;

    stickyRightShiftResult(const ubv &ser, const ubv &sb) : signExtendedResult(ser), stickyBit(sb) {}
  stickyRightShiftResult(const stickyRightShiftResult &old) : signExtendedResult(old.signExtendedResult), stickyBit(old.stickyBit) {}
  };


  template <class t>
  stickyRightShiftResult<t> stickyRightShift (const typename t::ubv &input, const typename t::ubv &shiftAmount) {
    stickyRightShiftResult<t> res(input.signExtendRightShift(shiftAmount), rightShiftStickyBit<t>(input, shiftAmount));

    return res;
  }


  template <class t>
  stickyRightShiftResult<t> stickyRightShiftBitwise (const typename t::ubv &input, const typename t::ubv &shiftAmount) {
    typedef typename t::bwt bwt;
    typedef typename t::prop prop;
    typedef typename t::ubv ubv;

    bwt width(input.getWidth());
    bwt startingPosition(positionOfLeadingOne(width));
    INVARIANT(0 < startingPosition && startingPosition < width);

    // Catch the out of bounds case
    PRECONDITION(shiftAmount.getWidth() == width);
    prop fullShift(shiftAmount >= ubv(width, width));
    // Note the shiftAmount is treated as unsigned...

    ubv *working = new ubv(input);
    prop *stickyBit = new prop(ITE(fullShift, !input.isAllZeros(), prop(false)));

    for (bwt i = startingPosition + 1; i > 0; --i)
    {
      bwt shiftAmountPosition = i - 1;

      prop shiftEnabled = fullShift || shiftAmount.extract(shiftAmountPosition, shiftAmountPosition).isAllOnes();

      prop stickyAccumulate(shiftEnabled && !(working->extract((1ULL << shiftAmountPosition) - 1, 0).isAllZeros())); // Would this loose data?

      prop * tmp = stickyBit;
      stickyBit = new prop(*stickyBit || stickyAccumulate);
      delete tmp;


      // Note the slightly unexpected sign extension
      ubv shifted(working->signExtendRightShift(ubv::one(width) << ubv(width, shiftAmountPosition)));

      ubv * workingTmp = working;
      working = new ubv(ITE(shiftEnabled, shifted, *working));
      delete workingTmp;
    }

    stickyRightShiftResult<t> res(*working, ubv(*stickyBit).extend(width - 1));

    delete working;
    delete stickyBit;

    POSTCONDITION(res.signExtendedResult == input.signExtendRightShift(shiftAmount));
    POSTCONDITION(res.stickyBit == rightShiftStickyBit<t>(input, shiftAmount));

    return res;
  }



  template <class t>
  struct normaliseShiftResult {
    typedef typename t::ubv ubv;
    typedef typename t::prop prop;

    ubv normalised;
    ubv shiftAmount;
    prop isZero;

  normaliseShiftResult(const ubv &n, const ubv &s, const prop &z) : normalised(n), shiftAmount(s), isZero(z) {}
  normaliseShiftResult(const normaliseShiftResult<t> &old) :  normalised(old.normalised), shiftAmount(old.shiftAmount), isZero(old.isZero) {}
  };

  template <class t>
  normaliseShiftResult<t> normaliseShift (const typename t::ubv input) {
    typedef typename t::bwt bwt;
    typedef typename t::prop prop;
    typedef typename t::ubv ubv;

    bwt width(input.getWidth());
    bwt startingMask(previousPowerOfTwo(width));
    INVARIANT(startingMask < width);

    // Catch the zero case
    prop zeroCase(input.isAllZeros());

    ubv *working = new ubv(input);
    ubv *shiftAmount = NULL;
    prop *deactivateShifts = new prop(zeroCase);

    for (bwt i = startingMask; i > 0; i >>= 1) {
      prop newDeactivateShifts = *deactivateShifts || working->extract(width-1,width-1).isAllOnes();
      delete deactivateShifts;
      deactivateShifts = new prop(newDeactivateShifts);

      ubv mask(ubv::allOnes(i).append(ubv::zero(width - i)));
      prop shiftNeeded(!(*deactivateShifts) && (mask & *working).isAllZeros());

      // Modular is safe because of the mask comparison
      ubv shifted(ITE(shiftNeeded, working->modularLeftShift(ubv(width, i)), *working));
      delete working;
      working = new ubv(shifted);

      if (shiftAmount == NULL) {
	shiftAmount = new ubv(shiftNeeded);
      } else {
	ubv newShiftAmount = shiftAmount->append(ubv(shiftNeeded));
	delete shiftAmount;
	shiftAmount = new ubv(newShiftAmount);
      }
    }

    normaliseShiftResult<t> res(*working, *shiftAmount, zeroCase);

    delete deactivateShifts;
    delete working;
    delete shiftAmount;

    POSTCONDITION(res.normalised.extract(width-1,width-1).isAllZeros() == res.isZero);
    POSTCONDITION(IMPLIES(res.isZero, res.shiftAmount.isAllZeros()));

    bwt shiftAmountWidth(res.shiftAmount.getWidth());
    bwt widthBits(bitsToRepresent(width));
    POSTCONDITION(shiftAmountWidth == widthBits ||
		  shiftAmountWidth == widthBits - 1); // If width is an exact power of 2
    ubv widthBV(widthBits, width);
    POSTCONDITION(res.shiftAmount.matchWidth(widthBV) < widthBV);

    return res;
  }


  /*** Dividers ***/
  template <class t>
  struct resultWithRemainderBit {
    typedef typename t::ubv ubv;
    typedef typename t::prop prop;

    ubv result;
    prop remainderBit;

  resultWithRemainderBit(const ubv &o, const prop &r) : result(o), remainderBit(r) {}
  resultWithRemainderBit(const resultWithRemainderBit<t> &old) : result(old.result), remainderBit(old.remainderBit) {}
  };

  // x and y are fixed-point numbers in the range [1,2)
  // Compute o \in [0.5,2), r \in [0,\delta) such that:  x = o*y + r
  // Return (o, r != 0)
  template <class t>
  resultWithRemainderBit<t> fixedPointDivide (const typename t::ubv &x, const typename t::ubv &y) {
    typename t::bwt w(x.getWidth());

    // Same width and both have MSB ones
    PRECONDITION(y.getWidth() == w);
    PRECONDITION(x.extract(w - 1, w - 1).isAllOnes());
    PRECONDITION(y.extract(w - 1, w - 1).isAllOnes());

    typedef typename t::ubv ubv;

    // Not the best way of doing this but pretty universal
    ubv ex(x.append(ubv::zero(w - 1)));
    ubv ey(y.extend(w - 1));

    ubv div(ex / ey);
    ubv rem(ex % ey);

    return resultWithRemainderBit<t>(div.extract(w - 1, 0), !(rem.isAllZeros()));
  }


  // x is a fixed-point number in the range [1,4) with 2/p bits
  // Compute o \in [1,sqrt(2)), r \in [0,o*2 + 1) such that x = o*o + r with 1/p bits
  // Return (o, r != 0)
  template <class t>
  resultWithRemainderBit<t> fixedPointSqrt (const typename t::ubv &x) {
    typedef typename t::bwt bwt;
    typedef typename t::ubv ubv;
    typedef typename t::prop prop;

    // The default algorithm given here isn't a great one
    // However it is simple and has a very simple termination criteria.
    // Plus most symbolic back ends will prefer
    // o = nondet(), r = nondet(), assert(r < o*2 + 1), assert(x = o*o + r)

    bwt inputWidth(x.getWidth());
    bwt outputWidth(inputWidth - 1);

    // To compare against, we need to pad x to 2/2p
    ubv xcomp(x.append(ubv::zero(inputWidth - 2)));

    // Start at 1
    ubv working(ubv::one(outputWidth) << ubv(outputWidth, outputWidth - 1));

    bwt location;
    for (location = outputWidth - 1; location > 0; --location) { // Offset by 1 for easy termination
      ubv shift(ubv(outputWidth, location - 1));

      ubv candidate(working | (ubv::one(outputWidth) << shift));

      prop addBit(expandingMultiply<t, ubv>(candidate, candidate) <= xcomp);

      working = working | (ubv(addBit).extend(outputWidth - 1) << shift);
    }

    return resultWithRemainderBit<t>(working, !(expandingMultiply<t, ubv>(working, working) == xcomp));
  }

  // One step of a divider
  // Here the "remainder bit" is actual the result bit and
  // The result is the remainder
  template <class t>
  resultWithRemainderBit<t> divideStep (const typename t::ubv &x, const typename t::ubv &y) {
    typedef typename t::bwt bwt;
    typedef typename t::ubv ubv;
    typedef typename t::prop prop;

    bwt xWidth(x.getWidth());
    bwt yWidth(y.getWidth());

    PRECONDITION(xWidth == yWidth);
    PRECONDITION(yWidth >= 2);
    PRECONDITION(y.extract(yWidth - 2, yWidth - 2).isAllOnes());  // Assume y is aligned

    prop canSubtract(x >= y);
    ubv sub(x.modularAdd(y.modularNegate())); // TODO : modular subtract or better
    ubv step(ITE(canSubtract, sub, x));

    return resultWithRemainderBit<t>(step << ubv::one(xWidth), canSubtract);
  }
}

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