1) QuadDouble -> replaced by Quad using templates. (c5b09e1e) · Commits · TNL / tnl-dev

src/TNL/Experimental/Arithmetics/CMakeLists.txt

+3 −3

Original line number	Diff line number	Diff line
		ADD_SUBDIRECTORY( UnitTests )

		set( headers FlopsCounter.h
		QuadDouble.h # Edited by D. Simon
		Quad.h # Modified by D. Simon
		Quad_impl.h # Modified by D. Simon
		Real.h
		MultiPrecision.h #Added by D. Simon
		Multiple.h #Added by D. Simon
		)

		SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/TNL/Experimental/Arithemtics )
		@@ -12,7 +12,7 @@ SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/TNL/Experimental/Arithemtics )
		set( tnl_experimental_arithmetics_SOURCES
		${CURRENT_DIR}/FlopsCounter.cpp
		${CURRENT_DIR}/MultiPrecision.cpp #Added by D.Simon
		${CURRENT_DIR}/QuadDouble.cpp #Added by D.Simon
		)

		set( tnl_experimental_arithmetics_CUDA__SOURCES )

src/TNL/Experimental/Arithmetics/Quad.h

0 → 100644

+104 −0

Original line number	Diff line number	Diff line
		//
		// Created by daniel on 25.11.17.
		/*
		* Quad.h
		*
		* Created by Matěj Novotný on 27.10.2010
		* Modified by Daniel Simon on 26.11.2017
		*
		* INFO: Quad and Quadcpp merged in this single class.
		* Quad edited to use templates.
		*/

		template <class T>
		class Quad
		{
		public:
		/INIT/
		T data[4];

		Quad();
		explicit Quad(const T&);
		explicit Quad(int);
		Quad(const Quad<T>&);

		/OVERLOADED OPERATORS/
		T& operator[](int);
		const T& operator[](int) const;
		Quad<T>& operator =(const Quad<T>&);
		Quad<T>& operator +=(const Quad<T>&);
		Quad<T>& operator -=(const Quad<T>&);
		Quad<T>& operator *=(const Quad<T>&);
		Quad<T>& operator /=(const Quad<T>&);
		Quad<T>& operator =(const T&);
		Quad<T>& operator +=(const T&);
		Quad<T>& operator -=(const T&);
		Quad<T>& operator *=(const T&);
		Quad<T>& operator /=(const T&);
		Quad<T> operator +(const Quad<T>&) const;
		Quad<T> operator -(const Quad<T>&) const;
		Quad<T> operator *(const Quad<T>&) const;
		Quad<T> operator /(const Quad<T>&) const;
		Quad<T> operator +(const T&) const;
		Quad<T> operator -(const T&) const;
		Quad<T> operator *(const T&) const;
		Quad<T> operator /(const T&) const;
		Quad<T> operator +();
		Quad<T> operator -();
		Quad<T> operator +() const;
		Quad<T> operator -() const;
		bool operator ==(const Quad<T>&) const;
		bool operator !=(const Quad<T>&) const;
		bool operator <(const Quad<T>&) const;
		bool operator >(const Quad<T>&) const;
		bool operator >=(const Quad<T>&) const;
		bool operator <=(const Quad<T>&) const;
		explicit operator T() const;
		};


		template <typename T>
		Quad<T> operator +(const T&, const Quad<T>&);
		template <typename T>
		Quad<T> operator -(const T&, const Quad<T>&);
		template <typename T>
		Quad<T> operator *(const T&, const Quad<T>&);
		template <typename T>
		Quad<T> operator /(const T&, const Quad<T>&);

		template <typename T>
		Quad<T> abs(const Quad<T>&);
		template <typename T>
		Quad<T> sqrt(const Quad<T>&);

		template <typename T>
		void quickTwoSum(T a, T b, T s, T e); // Addition of two doubles
		template <typename T>
		void twoSum(T a, T b, T s, T e); // Addition of two doubles
		template <typename T>
		void split(T a, T a_hi, T a_lo); // Split double into two 26 bits parts
		template <typename T>
		void twoProd(T a, T b, T p, T e); // Multiplication of two doubles
		template <typename T>
		void renormalize(T a, T b); // Normalization of number a
		template <typename T>
		void doublePlusQuad(T b, const T a, T s); // Addition of double and quad-double
		template <typename T>
		void doubleTimesQuad(T b, const T a, T s); // Multiplication of double and quad-double
		template <typename T>
		void quadDivDouble(const T a, T b, T s); // Division of two doubles
		template <typename T>
		void quadAdd(const T a, const T b, T *s); // Addition of two quad-doubles
		template <typename T>
		void quadAddAccurate(const T a, const T b, T *s); // Addition of two quad-doubles ! slower algorhitm
		template <typename T>
		void quadMul(const T a, const T b, T *s); // Multiplication of two quad-doubles
		template <typename T>
		void quadMulQuick(const T a, const T b, T *s); // Multiplication of two quad-doubles ! faster algorithm
		template <typename T>
		void quadDiv(const T a, const T b, T *s); // Division of two quad-doubles
		template <typename T>
		void zeroQuad(T *a); // Reset quad-double
		template <typename T>
		void printQuad(T *a); // Print of quad-double

src/TNL/Experimental/Arithmetics/QuadDouble.h

deleted100644 → 0

+0 −81

Original line number	Diff line number	Diff line
		/*
		* QuadDouble.h
		*
		* Created by Matěj Novotný on 27.10.2010
		* Modified by Daniel Simon on 28.10.2017
		*
		* INFO: Quad and Quadcpp merged in this single class.
		*/


		class QuadDouble
		{
		public:
		/INIT/
		double data[4];
		QuadDouble();
		explicit QuadDouble(const double&);
		explicit QuadDouble(int);
		QuadDouble(const QuadDouble&);

		/OVERLOADED OPERATORS/
		double& operator[](int);
		const double& operator[](int) const;
		QuadDouble& operator =(const QuadDouble&);
		QuadDouble& operator +=(const QuadDouble&);
		QuadDouble& operator -=(const QuadDouble&);
		QuadDouble& operator *=(const QuadDouble&);
		QuadDouble& operator /=(const QuadDouble&);
		QuadDouble& operator =(const double&);
		QuadDouble& operator +=(const double&);
		QuadDouble& operator -=(const double&);
		QuadDouble& operator *=(const double&);
		QuadDouble& operator /=(const double&);
		QuadDouble operator +(const QuadDouble&) const;
		QuadDouble operator -(const QuadDouble&) const;
		QuadDouble operator *(const QuadDouble&) const;
		QuadDouble operator /(const QuadDouble&) const;
		QuadDouble operator +(const double&) const;
		QuadDouble operator -(const double&) const;
		QuadDouble operator *(const double&) const;
		QuadDouble operator /(const double&) const;
		QuadDouble operator +();
		QuadDouble operator -();
		QuadDouble operator +() const;
		QuadDouble operator -() const;
		bool operator ==(const QuadDouble&) const;
		bool operator !=(const QuadDouble&) const;
		bool operator <(const QuadDouble&) const;
		bool operator >(const QuadDouble&) const;
		bool operator >=(const QuadDouble&) const;
		bool operator <=(const QuadDouble&) const;
		operator double() const;
		};


		/NON MEMBER OPERATORS/
		QuadDouble operator +(const double&, const QuadDouble&);
		QuadDouble operator -(const double&, const QuadDouble&);
		QuadDouble operator *(const double&, const QuadDouble&);
		QuadDouble operator /(const double&, const QuadDouble&);

		/NON MEMBER CLASS MATH FUNCTIONS/
		QuadDouble abs(const QuadDouble&);
		QuadDouble sqrt(const QuadDouble&);

		/NON MEMBER HELP FUNCTIONS FOR QUAD/
		void quickTwoSum(double a, double b, double s, double e); // Addition of two doubles
		void twoSum(double a, double b, double s, double e); // Addition of two doubles
		void split(double a, double a_hi, double a_lo); // Split double into two 26 bits parts
		void twoProd(double a, double b, double p, double e); // Multiplication of two doubles
		void renormalize(double a, double b); // Normalization of number a
		void doublePlusQuad(double b, const double a, double s); // Addition of double and quad-double
		void doubleTimesQuad(double b, const double a, double s); // Multiplication of double and quad-double
		void quadDivDouble(const double a, double b, double s); // Division of two doubles
		void quadAdd(const double a, const double b, double *s); // Addition of two quad-doubles
		void quadAddAccurate(const double a, const double b, double *s); // Addition of two quad-doubles ! slower algorhitm
		void quadMul(const double a, const double b, double *s); // Multiplication of two quad-doubles
		void quadMulQuick(const double a, const double b, double *s); // Multiplication of two quad-doubles ! faster algorithm
		void quadDiv(const double a, const double b, double *s); // Division of two quad-doubles
		void zeroQuad(double *a); // Reset quad-double
		void printQuad(double *a); // Print of quad-double

src/TNL/Experimental/Arithmetics/QuadDouble.cpp→src/TNL/Experimental/Arithmetics/Quad_impl.h

+188 −137

Original line number	Diff line number	Diff line
		/*
		* QuadDouble.cpp
		* Quad.cpp
		*
		* Created by Matěj Novotný on 27.10.2010
		* Modified by Daniel Simon on 28.10.2017
		* Modified by Daniel Simon on 23.11.2017
		*
		* INFO: Quad and Quadcpp merged in this single class.
		* Quad edited to use templates.
		*/

		#include "QuadDouble.h"
		#include <math.h>
		#include <stdio.h>
		#include <cmath>
		#include <cstdio>

		#include "Quad.h"

		#define ABS(n) ((n) > 0 ? (n): -(n))

		/CLASS METHODS IMPLEMENTATION/

		QuadDouble::QuadDouble() {
		template <class T>
		Quad<T>::Quad() {
		zeroQuad(data);
		}

		QuadDouble::QuadDouble(const double& value) {
		template <class T>
		Quad<T>::Quad(const T& value) {
		data[0] = value;
		data[1] = 0;
		data[2] = 0;
		data[3] = 0;
		}

		QuadDouble::QuadDouble(int value) {
		data[0] = (double)value;
		template <class T>
		Quad<T>::Quad(int value) {
		data[0] = (T)value;
		data[1] = 0;
		data[2] = 0;
		data[3] = 0;
		}

		QuadDouble::QuadDouble(const QuadDouble& other) {
		template <class T>
		Quad<T>::Quad(const Quad<T>& other) {
		data[0] = other[0];
		data[1] = other[1];
		data[2] = other[2];
		data[3] = other[3];
		}

		double& QuadDouble::operator [](int idx) {
		template <class T>
		T& Quad<T>::operator [](int idx) {
		return data[idx];
		}

		const double& QuadDouble::operator [](int idx) const{
		template <class T>
		const T& Quad<T>::operator [](int idx) const{
		return data[idx];
		}

		QuadDouble& QuadDouble::operator =(const QuadDouble& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator =(const Quad<T>& rhs) {
		data[0] = rhs[0];
		data[1] = rhs[1];
		data[2] = rhs[2];
		@@ -55,27 +64,32 @@ QuadDouble& QuadDouble::operator =(const QuadDouble& rhs) {
		return *this;
		}

		QuadDouble& QuadDouble::operator +=(const QuadDouble& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator +=(const Quad<T>& rhs) {
		quadAddAccurate(data, rhs.data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator -=(const QuadDouble& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator -=(const Quad<T>& rhs) {
		quadAddAccurate(data, (-rhs).data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator *=(const QuadDouble& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator *=(const Quad<T>& rhs) {
		quadMul(data, rhs.data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator /=(const QuadDouble& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator /=(const Quad<T>& rhs) {
		quadDiv(data, rhs.data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator =(const double& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator =(const T& rhs) {
		data[0] = rhs;
		data[1] = 0;
		data[2] = 0;
		@@ -83,84 +97,95 @@ QuadDouble& QuadDouble::operator =(const double& rhs) {
		return *this;
		}

		QuadDouble& QuadDouble::operator +=(const double& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator +=(const T& rhs) {
		doublePlusQuad(rhs, data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator -=(const double& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator -=(const T& rhs) {
		doublePlusQuad(-rhs, data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator *=(const double& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator *=(const T& rhs) {
		doubleTimesQuad(rhs, data, data);
		return *this;
		}

		QuadDouble& QuadDouble::operator /=(const double& rhs) {
		template <class T>
		Quad<T>& Quad<T>::operator /=(const T& rhs) {
		quadDivDouble(data, rhs, data);
		return *this;
		}



		QuadDouble QuadDouble::operator +(const QuadDouble& value) const{
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator +(const Quad<T>& value) const{
		Quad<T> qd(*this);
		qd += value;
		return qd;
		}

		QuadDouble QuadDouble::operator -(const QuadDouble& value) const{
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator -(const Quad<T>& value) const{
		Quad<T> qd(*this);
		qd -= value;
		return qd;
		}

		QuadDouble QuadDouble::operator *(const QuadDouble& value) const{
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator *(const Quad<T>& value) const{
		Quad<T> qd(*this);
		qd *= value;
		return qd;
		}

		QuadDouble QuadDouble::operator /(const QuadDouble& value) const{
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator /(const Quad<T>& value) const{
		Quad<T> qd(*this);
		qd /= value;
		return qd;
		}


		QuadDouble QuadDouble::operator +(const double& value) const {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator +(const T& value) const {
		Quad<T> qd(*this);
		qd += value;
		return qd;
		}

		QuadDouble QuadDouble::operator -(const double& value) const {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator -(const T& value) const {
		Quad<T> qd(*this);
		qd -= value;
		return qd;
		}

		QuadDouble QuadDouble::operator *(const double& value) const {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator *(const T& value) const {
		Quad<T> qd(*this);
		qd *= value;
		return qd;
		}

		QuadDouble QuadDouble::operator /(const double& value) const {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator /(const T& value) const {
		Quad<T> qd(*this);
		qd /= value;
		return qd;
		}

		QuadDouble QuadDouble::operator +() {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator +() {
		Quad<T> qd(*this);
		return qd;
		}

		QuadDouble QuadDouble::operator -() {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator -() {
		Quad<T> qd(*this);
		qd[0] = -qd[0];
		qd[1] = -qd[1];
		qd[2] = -qd[2];
		@@ -168,13 +193,15 @@ QuadDouble QuadDouble::operator -() {
		return qd;
		}

		QuadDouble QuadDouble::operator +() const {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator +() const {
		Quad<T> qd(*this);
		return qd;
		}

		QuadDouble QuadDouble::operator -() const {
		QuadDouble qd(*this);
		template <class T>
		Quad<T> Quad<T>::operator -() const {
		Quad<T> qd(*this);
		qd[0] = -qd[0];
		qd[1] = -qd[1];
		qd[2] = -qd[2];
		@@ -182,20 +209,22 @@ QuadDouble QuadDouble::operator -() const {
		return qd;
		}

		bool QuadDouble::operator ==(const QuadDouble& rhs) const {
		template <class T>
		bool Quad<T>::operator ==(const Quad<T>& rhs) const {
		if (data[0] == rhs[0] && data[1] == rhs[1] && data[2] == rhs[2] && data[3] == rhs[3]) {
		return true;
		}
		return false;
		}

		bool QuadDouble::operator !=(const QuadDouble& rhs) const {
		template <class T>
		bool Quad<T>::operator !=(const Quad<T>& rhs) const {
		return !(*this == rhs);
		}


		bool QuadDouble::operator <(const QuadDouble& rhs) const {
		QuadDouble qd(*this);
		template <class T>
		bool Quad<T>::operator <(const Quad<T>& rhs) const {
		Quad<T> qd(*this);
		qd -= rhs;
		if (qd[0] < 0.) {
		return true;
		@@ -203,8 +232,9 @@ bool QuadDouble::operator <(const QuadDouble& rhs) const {
		return false;
		}

		bool QuadDouble::operator >(const QuadDouble& rhs) const {
		QuadDouble qd(*this);
		template <class T>
		bool Quad<T>::operator >(const Quad<T>& rhs) const {
		Quad<T> qd(*this);
		qd -= rhs;
		if (qd[0] > 0.) {
		return true;
		@@ -212,8 +242,9 @@ bool QuadDouble::operator >(const QuadDouble& rhs) const {
		return false;
		}

		bool QuadDouble::operator >=(const QuadDouble& rhs) const {
		QuadDouble qd(*this);
		template <class T>
		bool Quad<T>::operator >=(const Quad<T>& rhs) const {
		Quad<T> qd(*this);
		qd -= rhs;
		if (qd[0] >= 0.) {
		return true;
		@@ -221,8 +252,9 @@ bool QuadDouble::operator >=(const QuadDouble& rhs) const {
		return false;
		}

		bool QuadDouble::operator <=(const QuadDouble& rhs) const {
		QuadDouble qd(*this);
		template <class T>
		bool Quad<T>::operator <=(const Quad<T>& rhs) const {
		Quad<T> qd(*this);
		qd -= rhs;
		if (qd[0] <= 0.) {
		return true;
		@@ -230,52 +262,53 @@ bool QuadDouble::operator <=(const QuadDouble& rhs) const {
		return false;
		}


		QuadDouble::operator double() const{
		template <class T>
		Quad<T>::operator T() const{
		return data[0];
		}


		/NON MEMBER OPERATORS/

		QuadDouble operator+(const double& v1, const QuadDouble& v2) {
		QuadDouble qd(v1);
		template <typename T>
		Quad<T> operator+(const T& v1, const Quad<T>& v2) {
		Quad<T> qd(v1);
		qd += v2;
		return qd;
		}

		QuadDouble operator-(const double& v1, const QuadDouble& v2) {
		QuadDouble qd(v1);
		template <typename T>
		Quad<T> operator-(const T& v1, const Quad<T>& v2) {
		Quad<T> qd(v1);
		qd -= v2;
		return qd;
		}

		QuadDouble operator*(const double& v1, const QuadDouble& v2) {
		QuadDouble qd(v1);
		template <typename T>
		Quad<T> operator*(const T& v1, const Quad<T>& v2) {
		Quad<T> qd(v1);
		qd *= v2;
		return qd;
		}

		QuadDouble operator/(const double& v1, const QuadDouble& v2) {
		QuadDouble qd(v1);
		template <typename T>
		Quad<T> operator/(const T& v1, const Quad<T>& v2) {
		Quad<T> qd(v1);
		qd /= v2;
		return qd;
		}

		/NON MEMBER CLASS MATH FUNCTIONS/

		QuadDouble abs(const QuadDouble& value) {
		QuadDouble qd(value);
		template <typename T>
		Quad<T> abs(const Quad<T>& value) {
		Quad<T> qd(value);
		if (value[0] < 0) {
		qd = -qd;
		}
		return qd;
		}

		QuadDouble sqrt(const QuadDouble& value) {
		QuadDouble qd(value);
		QuadDouble x(1/sqrt((double)qd));
		QuadDouble step;
		template <typename T>
		Quad<T> sqrt(const Quad<T>& value) {
		Quad<T> qd(value);
		Quad<T> x(1/sqrt((T)qd));
		Quad<T> step;
		//TODO zjednodušit dělení 2
		step = x * (1. - qd * x * x);
		step[0] /= 2;
		@@ -299,23 +332,22 @@ QuadDouble sqrt(const QuadDouble& value) {
		return qd;
		}

		/NON MEMBER HELP FUNCTIONS FOR QUAD/

		#define ABS(n) ((n) > 0 ? (n): -(n))

		void threeThreeSum(double i0, double i1, double i2, double o0, double o1, double *o2) {
		template <typename T>
		void threeThreeSum(T i0, T i1, T i2, T o0, T o1, T *o2) {
		twoSum(i0, i1, &i1, &i0); // 1
		twoSum(i1, i2, o0, &i1); // 2
		twoSum(i0, i1, o1, o2); // 3
		}

		void threeTwoSum(double i0, double i1, double i2, double o0, double o1) {
		template <typename T>
		void threeTwoSum(T i0, T i1, T i2, T o0, T o1) {
		twoSum(i0, i1, &i1, &i0); // 1
		twoSum(i1, i2, o0, &i1); // 2
		*o1 = i1 + i0; // 3
		}

		void fourTwoSum(double i0, double i1, double i2, double i3, double o0, double o1) {
		template <typename T>
		void fourTwoSum(T i0, T i1, T i2, T i3, T o0, T o1) {
		twoSum(i0, i2, &i0, &i2); // 1
		twoSum(i1, i3, &i1, &i3); // 2
		i2 += i1; // 3
		@@ -324,8 +356,9 @@ void fourTwoSum(double i0, double i1, double i2, double i3, double o0, double
		quickTwoSum(i0, i3, o0, o1); // 6
		}

		void sixThreeSum(double i0, double i1, double i2, double i3, double i4, double i5, double o0, double o1,
		double *o2) {
		template <typename T>
		void sixThreeSum(T i0, T i1, T i2, T i3, T i4, T i5, T o0, T o1,
		T *o2) {
		threeThreeSum(i0, i1, i2, &i2, &i0, &i1); // 1
		threeThreeSum(i3, i4, i5, &i5, &i3, &i4); // 2
		twoSum(i2, i5, o0, &i5); // 3
		@@ -334,15 +367,17 @@ void sixThreeSum(double i0, double i1, double i2, double i3, double i4, double i
		*o2 = i1 + i4 + i3 + i5; // 6
		}

		void sixTwoSum(double i0, double i1, double i2, double i3, double i4, double i5, double o0, double o1) {
		template <typename T>
		void sixTwoSum(T i0, T i1, T i2, T i3, T i4, T i5, T o0, T o1) {
		threeTwoSum(i0, i1, i2, &i1, &i0); // 1
		threeTwoSum(i3, i4, i5, &i4, &i3); // 2
		twoSum(i1, i4, o0, &i1); // 3
		*o1 = i0 + i3 + i1; // 4
		}

		void nineTwoSum(double i0, double i1, double i2, double i3, double i4, double i5, double i6, double i7,
		double i8, double o0, double o1) {
		template <typename T>
		void nineTwoSum(T i0, T i1, T i2, T i3, T i4, T i5, T i6, T i7,
		T i8, T o0, T o1) {
		twoSum(i5, i6, &i5, &i6); // 1
		twoSum(i4, i7, &i4, &i7); // 2
		twoSum(i1, i2, &i1, &i2); // 3
		@@ -353,50 +388,56 @@ void nineTwoSum(double i0, double i1, double i2, double i3, double i4, double i5
		threeTwoSum(i3, i0, i8, o0, o1); // 8
		}

		void quickTwoSum(double a, double b, double s, double e) {
		template <typename T>
		void doubleAccumulate(T i0, T i1, T i2, T o0, T o1, T *o2) {
		twoSum(i1, i2, o0, o2);
		twoSum(i0, *o0, o0, o1);
		if (*o1 == 0) {
		o1 = o0;
		*o0 = 0;
		}
		if (*o2 == 0) {
		o2 = o1;
		o1 = o0;
		*o0 = 0;
		}
		}

		template <typename T>
		void quickTwoSum(T a, T b, T s, T e) {
		*s = a + b;
		e = b - (s - a);
		}

		void twoSum(double a, double b, double s, double e) {
		template <typename T>
		void twoSum(T a, T b, T s, T e) {
		*s = a + b;
		double v = *s - a;
		e = (a - (s - v)) + (b - v);
		}

		void split(double a, double a_hi, double a_lo) {
		double t = 134217729 * a;
		template <typename T>
		void split(T a, T a_hi, T a_lo) {
		T t = 134217729 * a;
		*a_hi = t - (t - a);
		a_lo = a - a_hi;
		}

		void twoProd(double a, double b, double p, double e) {
		template <typename T>
		void twoProd(T a, T b, T p, T e) {
		p = a b;
		double a_hi, a_lo, b_hi, b_lo;
		T a_hi, a_lo, b_hi, b_lo;
		split(a, &a_hi, &a_lo);
		split(b, &b_hi, &b_lo);
		e = ((a_hi b_hi - p) + a_hi b_lo + a_lo * b_hi) + a_lo * b_lo;
		}

		void doubleAccumulate(double i0, double i1, double i2, double o0, double o1, double *o2) {
		twoSum(i1, i2, o0, o2);
		twoSum(i0, *o0, o0, o1);
		if (*o1 == 0) {
		o1 = o0;
		*o0 = 0;
		}
		if (*o2 == 0) {
		o2 = o1;
		o1 = o0;
		*o0 = 0;
		}
		}

		void renormalize(double a, double b) {
		double s;
		double t[5];
		template <typename T>
		void renormalize(T a, T b) {
		T s;
		T t[5];
		int k = 0;
		double e;
		T e;
		quickTwoSum(a[3], a[4], &s, t + 4);
		quickTwoSum(a[2], s, &s, t + 3);
		quickTwoSum(a[1], s, &s, t + 2);
		@@ -414,9 +455,10 @@ void renormalize(double a, double b) {
		}
		}

		void doublePlusQuad(double b, const double a, double s) {
		double m[5];
		double e = b;
		template <typename T>
		void doublePlusQuad(T b, const T a, T s) {
		T m[5];
		T e = b;
		int i;
		for (i = 0; i < 4; i++) {
		twoSum(e, a[i], m + i, &e);
		@@ -425,8 +467,9 @@ void doublePlusQuad(double b, const double a, double s) {
		renormalize(m, s);
		}

		void doubleTimesQuad(double b, const double a, double s) {
		double m[7];
		template <typename T>
		void doubleTimesQuad(T b, const T a, T s) {
		T m[7];
		twoProd(b, a[0], m, m + 1); // 1
		twoProd(b, a[1], m + 2, m + 3); // 2
		twoSum(m[1], m[2], m + 1, m + 2); // 3
		@@ -438,10 +481,11 @@ void doubleTimesQuad(double b, const double a, double s) {
		renormalize(m, s);
		}

		void quadDivDouble(const double a, double b, double s) {
		template <typename T>
		void quadDivDouble(const T a, T b, T s) {
		//double b1[] = {b, 0, 0, 0};
		//quadDiv(a, b1, s);
		double m[13];
		T m[13];
		int i; // ten půjde odstranit
		m[5] = a[0];
		m[6] = a[1];
		@@ -458,8 +502,9 @@ void quadDivDouble(const double a, double b, double s) {
		renormalize(m, s);
		}

		void quadAdd(const double a, const double b, double *s) {
		double m[8];
		template <typename T>
		void quadAdd(const T a, const T b, T *s) {
		T m[8];
		twoSum(a[0], b[0], m, m + 1); // 1
		twoSum(a[1], b[1], m + 2, m + 3); // 2
		twoSum(m[2], m[1], m + 1, m + 2); // 3
		@@ -473,8 +518,9 @@ void quadAdd(const double a, const double b, double *s) {
		renormalize(m, s);
		}

		void quadAddAccurate(const double a, const double b, double *s) {
		double m[11];
		template <typename T>
		void quadAddAccurate(const T a, const T b, T *s) {
		T m[11];
		int i = 0;
		int j = 0;
		int k = 0;
		@@ -513,8 +559,9 @@ void quadAddAccurate(const double a, const double b, double *s) {
		renormalize(m, s);
		}

		void quadMul(const double a, const double b, double *s) {
		double m[20];
		template <typename T>
		void quadMul(const T a, const T b, T *s) {
		T m[20];
		twoProd(a[0], b[0], m, m + 1); // 1
		twoProd(a[0], b[1], m + 2, m + 3); // 2
		twoProd(a[1], b[0], m + 4, m + 5); // 3
		@@ -532,8 +579,9 @@ void quadMul(const double a, const double b, double *s) {
		renormalize(m, s);
		}

		void quadMulQuick(const double a, const double b, double *s) {
		double m[12];
		template <typename T>
		void quadMulQuick(const T a, const T b, T *s) {
		T m[12];
		twoProd(a[0], b[0], m, m + 1); // 1
		twoProd(a[0], b[1], m + 2, m + 3); // 2
		twoProd(a[1], b[0], m + 4, m + 5); // 3
		@@ -547,8 +595,9 @@ void quadMulQuick(const double a, const double b, double *s) {
		renormalize(m, s);
		}

		void quadDiv(const double a, const double b, double *s) {
		double m[13];
		template <typename T>
		void quadDiv(const T a, const T b, T *s) {
		T m[13];
		//double n[4];
		//double k[4];
		int i; // ten půjde odstranit
		@@ -564,13 +613,15 @@ void quadDiv(const double a, const double b, double *s) {
		renormalize(m, s);
		}

		void zeroQuad(double *a) {
		template <typename T>
		void zeroQuad(T *a) {
		a[0] = 0;
		a[1] = 0;
		a[2] = 0;
		a[3] = 0;
		}

		void printQuad(double *a) {
		template <typename T>
		void printQuad(T *a) {
		printf("%.15le + %.15le + %.15le + %.15le\n", a[0], a[1], a[2], a[3]);
		}

src/TNL/Experimental/Arithmetics/UnitTests/MultiPrecisionTest.cpp

0 → 100644

+0 −0

Empty file added.