/*

** Copyright (C) 2018 Martin Brain

**

** See the file LICENSE for licensing information.

*/

/*

** add.h

**

** Martin Brain

** martin.brain@cs.ox.ac.uk

** 01/09/14

**

** Addition of arbitrary precision floats

**

** The current design is based on a two-path adder but it may be useful to use

** more paths.  There are five cases that are of interest:

**  1. effective add / very far

**     -> set the sticky bit only

**  2. effective add / far or near

**     -> align and add, realign down if needed

**  3. effective sub / very far

**     -> decrement, re-normalise and set sticky bits or (dependent on rounding-mode, skip entirely)

**  4. effective sub / far

**     -> align and subtract, realign up if needed

**  5. effective sub / near

**     -> align, subtract and normalise up

**

*/

/*

** Ideas

**  Optimisation : Collar the exponent difference, convert add to twice the width and thus unify the paths and simplify the shifting.

**  Enable the absolute max underapproximation

*/

#include "symfpu/core/unpackedFloat.h"

#include "symfpu/core/ite.h"

#include "symfpu/core/rounder.h"

#include "symfpu/core/sign.h"

#include "symfpu/core/operations.h"

#ifndef SYMFPU_ADD

#define SYMFPU_ADD

namespace symfpu {

  // There are a number of variants on how this should be done.

  // This is the implementation that handles all of them.

  // Below are restricted versions for particular special cases.

template <class t>

						    const typename t::rm &roundingMode,

						    const unpackedFloat<t> &left,

						    const unpackedFloat<t> &leftID,

						    const unpackedFloat<t> &right,

						    const typename t::prop &returnLeft,

						    const typename t::prop &returnRight,

						    const unpackedFloat<t> &additionResult,

						    const typename t::prop &isAdd) {

  typedef typename t::prop prop;

  // NaN

  //prop compatableSigns(ITE(isAdd, signsMatch, !signsMatch));

  // Inf

  // Zero

  // At most one of idLeft, idRight, generatesNaN, generatesInf and bothZero is true.

  // If used in addition additionResult is guaranteed to not be NaN.

  // Subtle trick : as the input to this will have been rounded it will have

  // an ITE with the default values "on top", thus doing the special cases

  // first (inner) rather than last (outer) allows them to be compacted better

	     ITE(isAdd,

		 right,

		 negate(format, right)),

		 leftID,

		 ITE(generatesNaN,

		     unpackedFloat<t>::makeNaN(format),

		     ITE(generatesInf,

			 unpackedFloat<t>::makeInf(format, signOfInf),

			 ITE(bothZero,

			     unpackedFloat<t>::makeZero(format, signOfZero),

 }

  // leftID is the value returned in the idLeft case (i.e. when left is not a

  // special number and right is zero).  This is needed by FMA as the flags

  // for left and leftID are computed differently and need to be handled differently.

template <class t>

						  const typename t::rm &roundingMode,

						  const unpackedFloat<t> &left,

						  const unpackedFloat<t> &leftID,

						  const unpackedFloat<t> &right,

						  const unpackedFloat<t> &additionResult,

						  const typename t::prop &isAdd) {

  return addAdditionSpecialCasesComplete<t>(format, roundingMode, left, leftID,

					    right, typename t::prop(false), typename t::prop(false),

  }

  // This is the usual case; use this one!

  template <class t>

					    const typename t::rm &roundingMode,

					    const unpackedFloat<t> &left,

					    const unpackedFloat<t> &right,

					    const unpackedFloat<t> &additionResult,

					    const typename t::prop &isAdd) {

    return addAdditionSpecialCasesComplete<t>(format, roundingMode, left, left,

					      right, typename t::prop(false), typename t::prop(false),

  }

  // As above but allows the (very) far path to be accelerated

  template <class t>

  unpackedFloat<t> addAdditionSpecialCasesWithBypass (const typename t::fpt &format,

						      const typename t::rm &roundingMode,

						      const unpackedFloat<t> &left,

						      const unpackedFloat<t> &right,

						      const typename t::prop &returnLeft,

						      const typename t::prop &returnRight,

						      const unpackedFloat<t> &additionResult,

						      const typename t::prop &isAdd) {

    return addAdditionSpecialCasesComplete<t>(format, roundingMode, left, left,

					      right, returnLeft, returnRight,

					      additionResult, isAdd);

  }

  /* Computes the normal / subnormal case only.

   * This allows multiple versions of the first phase to be used and

   * the first phase to be used for other things (e.g. FMA).

   */

template <class t>

  typedef typename t::sbv sbv;

  typedef typename t::prop prop;

  prop leftIsMax;

  sbv maxExponent;

  sbv absoluteExponentDifference;

  prop diffIsZero;

  prop diffIsOne;

  prop diffIsGreaterThanPrecision;

  prop diffIsTwoToPrecision;

  prop diffIsGreaterThanPrecisionPlusOne;

		      const prop &diz, const prop &dio, const prop &digtp, const prop &dittp, const prop &disgtppo) :

    leftIsMax(lil), maxExponent(me), absoluteExponentDifference(aed),

  exponentCompareInfo(const exponentCompareInfo<t> &old) :

    leftIsMax(old.leftIsMax), maxExponent(old.maxExponent), absoluteExponentDifference(old.absoluteExponentDifference),

    diffIsZero(old.diffIsZero), diffIsOne(old.diffIsOne),

    diffIsGreaterThanPrecision(old.diffIsGreaterThanPrecision),

    diffIsTwoToPrecision(old.diffIsTwoToPrecision),

    diffIsGreaterThanPrecisionPlusOne(old.diffIsGreaterThanPrecisionPlusOne) {}

};

template <class t>

					    const typename t::bwt significandWidth,

					    const typename t::sbv &leftExponent,

					    const typename t::sbv &rightExponent,

					    const typename t::prop &knownInCorrectOrder) {

   typedef typename t::prop prop;

   typedef typename t::sbv sbv;

#if 0

   // The obvious implementation

   // Compute exponent distance

   sbv maxExponent(max<t,sbv>(leftExponent.extend(1), rightExponent.extend(1)));

   sbv minExponent(min<t,sbv>(leftExponent.extend(1), rightExponent.extend(1)));

   sbv absoluteExponentDifference(maxExponent - minExponent);

   prop leftIsMax(knownInCorrectOrder || leftExponent.extend(1) == maxExponent);

#else

   // A better implementation(?)

#endif

   // Optimisation : compact these comparisons at the bit-level

   return exponentCompareInfo<t>(leftIsMax,

				 maxExponent, absoluteExponentDifference,

 }

// Note that the arithmetic part of add needs the rounding mode.

// This is an oddity due to the way that the sign of zero is generated.

 template <class t>

   unpackedFloat<t> uf;

   customRounderInfo<t> known;

 floatWithCustomRounderInfo(const floatWithCustomRounderInfo<t> &old) : uf(old.uf), known(old.known) {}

 };

 template <class t>

						const typename t::rm &roundingMode,

						const unpackedFloat<t> &left,

						const unpackedFloat<t> &right,

						const typename t::prop &isAdd,

						const typename t::prop &knownInCorrectOrder,

						const exponentCompareInfo<t> &ec) {

   typedef typename t::bwt bwt;

   typedef typename t::fpt fpt;

   typedef typename t::prop prop;

   typedef typename t::ubv ubv;

   typedef typename t::sbv sbv;

   // Work out if an effective subtraction

   /* Exponent difference and effective add implies a large amount about the output exponent and flags

   ** R denotes that this is possible via rounding up and incrementing the exponent

   **

   ** Case       A. max(l,r) + 1,    B. max(l,r)    C. max(l,r) - 1      D. max(l,r) - k       E. zero

   ** Eff. Add      Y                   Y

   **  diff = 0     Y, sticky 0

   **  diff = 1     Y, sticky 0, R      Y, sticky 0

   **  diff : [2,p] decreasing prob., R Y

   **  diff > p     R                   Y

   **

   ** Eff. Sub                          Y              Y                    Y, exact              Y, exact

   **  diff = 0                                        Y, exact             prob. drop with k     low prob.

   **  diff = 1                         Y, sticky 0    Y, exact             prob. drop with k

   **  diff : [2,p]                     Y, R           decreasing prob.

   **  diff > p                         Y, R           low prob.

   **

   */

   // Optimisation : add 'bypass' invariants that link the final exponent to the input using this.

   // Rounder flags

   //prop exact(); // Need to see if it cancels, see below

   // Work out ordering

			prop(true),

   // Optimisation : can we avoid this comparison completely and allow the result to be negative?

   // This may be hard as a compare is cheaper than a negate after, particularly as there has to be an ITE here

   // Extend the significands to give room for carry plus guard and sticky bits

		       left.getSign(),

   // Extended so no info lost, negate before shift so that sign-extension works

		   .resize(negatedSmaller.getWidth()));  // Safe as long as the significand has more bits than the exponent

   // Shift the smaller significand

				 ITE(effectiveAdd,

				 shifted.signExtendedResult));

			    shifted.stickyBit));  // Have to separate otherwise align up may convert it to the guard bit

   // Sum and re-align

						      conditionalRightShiftOne<t,ubv,prop>(overflow, sum)));

				  -sbv::one(exponentWidth),

				  ITE(overflow,

				      sbv::one(exponentWidth),

				      sbv::zero(exponentWidth))));

   // Watch closely...

   // We return something in an extended format

   //  *. One extra exponent bit to deal with the 'overflow' case

   //  *. Two extra significand bits for the guard and sticky bits

   // Put it back together

   // Deal with the major cancellation case

   // It would be nice to use normaliseUpDetectZero but the sign

   // of the zero depends on the rounding mode.

				       ITE(majorCancel,

					   sumResult.normaliseUp(extendedFormat),

					   sumResult)));

   // Some thought is required here to convince yourself that

   // there will be no subnormal values that violate this.

   // See 'all subnormals generated by addition are exact'

   // and the extended exponent.

 }

 template <class t>

   unpackedFloat<t> dualPathArithmeticAdd (const typename t::fpt &format,

					   const typename t::rm &roundingMode,

					   const unpackedFloat<t> &left,

					   const unpackedFloat<t> &right,

					   const typename t::prop &isAdd) {

   typedef typename t::bwt bwt;

   typedef typename t::fpt fpt;

   typedef typename t::prop prop;

   typedef typename t::ubv ubv;

   typedef typename t::sbv sbv;

   PRECONDITION(left.valid(format));

   PRECONDITION(right.valid(format));

   // We return something in an extended format

   //  *. One extra exponent bit to deal with the 'overflow' case

   //  *. Two extra significand bits for the guard and sticky bits

   fpt extendedFormat(format.exponentWidth() + 1, format.significandWidth() + 2);

   // Compute exponent difference and swap the two arguments if needed

   sbv initialExponentDifference(expandingSubtract<t>(left.getExponent(), right.getExponent()));

   bwt edWidth(initialExponentDifference.getWidth());

   sbv edWidthZero(sbv::zero(edWidth));

   prop orderingCorrect( initialExponentDifference >  edWidthZero ||

			(initialExponentDifference == edWidthZero &&

			 left.getSignificand() >= right.getSignificand()));

   unpackedFloat<t> larger(ITE(orderingCorrect, left, right));

   unpackedFloat<t> smaller(ITE(orderingCorrect, right, left));

   sbv exponentDifference(ITE(orderingCorrect,

			      initialExponentDifference,

			      -initialExponentDifference));

   prop resultSign(ITE(orderingCorrect,

		       left.getSign(),

		       prop(!isAdd ^ right.getSign())));

   // Work out if an effective subtraction

   prop effectiveAdd(larger.getSign() ^ smaller.getSign() ^ isAdd);

   // Extend the significands to give room for carry plus guard and sticky bits

   ubv lsig( larger.getSignificand().extend(1).append(ubv::zero(2)));

   ubv ssig(smaller.getSignificand().extend(1).append(ubv::zero(2)));

   // This is a two-path adder, so determine which of the two paths to use

   // The near path is only needed for things that can cancel more than one bit

   prop farPath(exponentDifference > sbv::one(edWidth) || effectiveAdd);

   // Far path : Align

   ubv negatedSmaller(ITE(effectiveAdd, ssig, ssig.modularNegate())); // Extended so no info lost

                                                                      // Negate before shift so that sign-extension works

   sbv significandWidth(edWidth, lsig.getWidth());

   prop noOverlap(exponentDifference > significandWidth);

   ubv shiftAmount(exponentDifference.toUnsigned() // Safe as >= 0

		   .resize(ssig.getWidth()));      // This looses information but the case in which it does is handles by noOverlap

   ubv negatedAlignedSmaller(negatedSmaller.signExtendRightShift(shiftAmount));

   ubv shiftedStickyBit(rightShiftStickyBit<t>(negatedSmaller, shiftAmount));  // Have to separate otherwise align up may convert it to the guard bit

   // Far path : Sum and re-align

   ubv sum(lsig.modularAdd(negatedAlignedSmaller));

   bwt sumWidth(sum.getWidth());

   ubv topBit(sum.extract(sumWidth - 1, sumWidth - 1));

   ubv centerBit(sum.extract(sumWidth - 2, sumWidth - 2));

   prop noOverflow(topBit.isAllZeros()); // Only correct if effectiveAdd is set

   prop noCancel(centerBit.isAllOnes());

   // TODO : Add invariants

   ubv alignedSum(ITE(effectiveAdd,

		      ITE(noOverflow,

			  sum,

			  (sum >> ubv::one(sumWidth)) | (sum & ubv::one(sumWidth))),  // Cheap sticky right shift

		      ITE(noCancel,

			  sum,

			  sum.modularLeftShift(ubv::one(sumWidth))))); // In the case when this looses data, the result is not used

   sbv extendedLargerExponent(larger.getExponent().extend(1));  // So that increment and decrement don't overflow

   sbv correctedExponent(ITE(effectiveAdd,

			     ITE(noOverflow,

				 extendedLargerExponent,

				 extendedLargerExponent.increment()),

			     ITE(noCancel,

				 extendedLargerExponent,

				 extendedLargerExponent.decrement())));

   // Far path : Construct result

   unpackedFloat<t> farPathResult(extendedFormat, resultSign, correctedExponent, (alignedSum | shiftedStickyBit).contract(1));

   // Near path : Align

   prop exponentDifferenceAllZeros(exponentDifference.isAllZeros());

   ubv nearAlignedSmaller(ITE(exponentDifferenceAllZeros, ssig, ssig >> ubv::one(ssig.getWidth())));

   // Near path : Sum and realign

   ubv nearSum(lsig - nearAlignedSmaller);

   // Optimisation : the two paths can be merged up to here to give a pseudo-two path encoding

   prop fullCancel(nearSum.isAllZeros());

   prop nearNoCancel(nearSum.extract(sumWidth - 2, sumWidth - 2).isAllOnes());

   ubv choppedNearSum(nearSum.extract(sumWidth - 3,1)); // In the case this is used, cut bits are all 0

   unpackedFloat<t> cancellation(extendedFormat,

				 resultSign,

				 larger.getExponent().decrement(),

				 choppedNearSum);

   // Near path : Construct result

   unpackedFloat<t> nearPathResult(extendedFormat, resultSign, extendedLargerExponent, nearSum.contract(1));

   // Bring the paths together

   // Optimisation : fix the noOverlap / very far path for directed rounding modes

   unpackedFloat<t> additionResult(ITE(farPath,

				       /* ITE(noOverlap,

					   ITE((isAdd || orderingCorrect),

					       larger,

					       negate(format, larger)),

					       farPathResult), */

				       farPathResult,

				       ITE(fullCancel,

					   unpackedFloat<t>::makeZero(extendedFormat, roundingMode == t::RTN()),

					   ITE(nearNoCancel,

					       nearPathResult,

					       cancellation.normaliseUp(format).extend(1,2)))));

   // Some thought is required here to convince yourself that

   // there will be no subnormal values that violate this.

   // See 'all subnormals generated by addition are exact'

   // and the extended exponent.

   POSTCONDITION(additionResult.valid(extendedFormat));

   return additionResult;

 }

 template <class t>

   unpackedFloat<t> dualPathAdd (const typename t::fpt &format,

				 const typename t::rm &roundingMode,

				 const unpackedFloat<t> &left,

				 const unpackedFloat<t> &right,

				 const typename t::prop &isAdd) {

   PRECONDITION(left.valid(format));

   PRECONDITION(right.valid(format));

   unpackedFloat<t> additionResult(dualPathArithmeticAdd(format, roundingMode, left, right, isAdd));

   unpackedFloat<t> roundedAdditionResult(rounder(format, roundingMode, additionResult));

   unpackedFloat<t> result(addAdditionSpecialCases(format, roundingMode, left, right, roundedAdditionResult, isAdd));

   POSTCONDITION(result.valid(format));

   return result;

 }

template <class t>

			 const typename t::rm &roundingMode,

			 const unpackedFloat<t> &left,

			 const unpackedFloat<t> &right,

			 const typename t::prop &isAdd) {

   //typedef typename t::bwt bwt;

   typedef typename t::prop prop;

   //typedef typename t::ubv ubv;

   //typedef typename t::sbv sbv;

   // Optimisation : add a flag which assumes that left and right are in the correct order

						   left.getExponent(), right.getExponent(), knownInCorrectOrder));

 }

 template <class t>

   unpackedFloat<t> addWithBypass (const typename t::fpt &format,

				   const typename t::rm &roundingMode,

				   const unpackedFloat<t> &left,

				   const unpackedFloat<t> &right,

				   const typename t::prop &isAdd) {

   //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));

   // Optimisation : add a flag which assumes that left and right are in the correct order

   prop knownInCorrectOrder(false);

   exponentCompareInfo<t> ec(addExponentCompare<t>(left.getExponent().getWidth() + 1, left.getSignificand().getWidth(),

						   left.getExponent(), right.getExponent(), knownInCorrectOrder));

   floatWithCustomRounderInfo<t> additionResult(arithmeticAdd(format, roundingMode, left, right, isAdd, knownInCorrectOrder, ec));

   unpackedFloat<t> roundedAdditionResult(customRounder(format, roundingMode, additionResult.uf, additionResult.known));

   // In the "very far path", i.e. exponent difference is greater than significand length + 1

   // addition becomes max(left,right) or max(left,right)+/-1 ULP.  This is rare in execution

   // (to the point of being a software quality issue) but common in theorem proving.

   // Given we have to have cases for "return left" and "return right" to handle zeros,

   // we might as well make use of these cases to handle when addition behaves like max...

   // Note that this is possible but more complex with just diffIsGreaterThanPrecision.

   prop enableBypass(ec.diffIsGreaterThanPrecisionPlusOne &&

		      !left.getNaN() &&  !left.getInf() &&  !left.getZero() && // Handle this as special cases

		     !right.getNaN() && !right.getInf() && !right.getZero());

   // Duplication but easier to recompute than to pass

   prop effectiveAdd((left.getSign() ^ right.getSign()) ^ isAdd);

   prop resultSign(ITE((knownInCorrectOrder || ec.leftIsMax),  // CAUTION : only true in the enableBypass case!

		       left.getSign(),

		       prop(!isAdd ^ right.getSign())));

   prop significandEven(true);  // This is an optimisation that assumes only RNE uses this bit

                                // This needs to be changed to implement things like roundToOdd

                                // or for the diffIsGreaterThanPrecision case.

   prop farRoundUp(roundingDecision<t>(roundingMode, resultSign, significandEven, !effectiveAdd, prop(true), prop(false)));

   // Returns left or right unchanged if adding and rounded down or subtracting and rounded up

   prop roundInCorrectDirection(effectiveAdd ^ farRoundUp);

   prop  returnLeft(enableBypass &&  ec.leftIsMax && roundInCorrectDirection);

   prop returnRight(enableBypass && !ec.leftIsMax && roundInCorrectDirection);

   unpackedFloat<t> result(addAdditionSpecialCasesWithBypass(format, roundingMode, left, right, returnLeft, returnRight, roundedAdditionResult, isAdd));

   POSTCONDITION(result.valid(format));

   return result;

 }

 // True if and only if adding these would result in a catastrophic cancellation

 // I.E. if the addition cancells out cancelAmount or more MSBs leaving only LSBs

 template <class t>

   typename t::prop isCatastrophicCancellation (const typename t::fpt &format,

						const unpackedFloat<t> &left,

						const unpackedFloat<t> &right,

						const typename t::bwt &cancelAmount,

						const typename t::prop &isAdd) {

   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));

   PRECONDITION(cancelAmount >= 2);   // cancel 0 is not meaningful

                                      // cancel 1 is common on subtract and arguably not an error

   PRECONDITION(cancelAmount <= format.significandWidth());  // can't cancel more than you have

   // 1. It has to be an effective subtraction

   // Duplication but easier to recompute than to pass

   prop effectiveAdd((left.getSign() ^ right.getSign()) ^ isAdd);

   // 2. It must be either normal or subnormal numbers

   prop leftSpecial(left.getNaN() || left.getInf() || left.getZero());

   prop rightSpecial(right.getNaN() || right.getInf() || right.getZero());

   // 3.A. exponents are equal and so are the leading cancelAmount bits

   // 3.B. exponent diff is one and the smaller one is 11111, larger one is 10000

   // Optimisation : add a flag which assumes that left and right are in the correct order

   prop knownInCorrectOrder(false);

   exponentCompareInfo<t> ec(addExponentCompare<t>(left.getExponent().getWidth() + 1, left.getSignificand().getWidth(),

						   left.getExponent(), right.getExponent(), knownInCorrectOrder));

   // Can ignore the MSB of the significand as by invariants this is always 1

   bwt significandWidth(format.significandWidth());

   bwt topBit(significandWidth - 2);

   bwt bottomBit(significandWidth - cancelAmount);

   ubv leftExtract(left.getSignificand.extract(topBit, bottomBit));

   ubv rightExtract(right.getSignificand.extract(topBit, bottomBit));

   prop result(ITE(!effectiveAdd && !leftSpecial && !rightSpecial,

		   ITE(ec.diffIsZero,

		       leftExtract == rightExtract,

		       ITE(ec.diffIsOne,

			   ITE(ec.leftIsMax,

			        leftExtract.isAllZeros() && rightExtract.isAllOnes(),

			       rightExtract.isAllZeros() &&  leftExtract.isAllOnes()),

			   false)),

		   false));

   return result;

 }

}

#endif