// Copyright 2006. Peter Gottschling, Matthias Troyer, Rolf Bonderer// Software License for MTL// // Copyright (c) 2007 The Trustees of Indiana University. All rights reserved.// Authors: Peter Gottschling and Andrew Lumsdaine// // This file is part of the Matrix Template Library// // See also license.mtl.txt in the distribution.#ifndef LA_CONCEPTS_INCLUDE#define LA_CONCEPTS_INCLUDE#include <boost/config/concept_macros.hpp>#ifdef __GXX_CONCEPTS__#  include <concepts>#else#  ifdef LA_SHOW_WARNINGS#    warning "Concepts are not used"#  endif#endif#include <boost/numeric/linear_algebra/identity.hpp>#include <boost/numeric/linear_algebra/is_invertible.hpp>#include <boost/numeric/linear_algebra/inverse.hpp>#include <boost/numeric/linear_algebra/operators.hpp>#include <boost/numeric/linear_algebra/algebraic_concepts.hpp>#include <complex>// If desired one can disable the default concept maps with LA_NO_CONCEPT_MAPS// We consider to change the namespace from math to numeric// More precisely, the concepts may be moved into namespace numeric and the standard functions stay in math/// Namespace for mathematical concepts/** In contrast to the ones in algebra the concepts can require basic implementation concepts like std::Assignable */  namespace math {#ifdef __GXX_CONCEPTS__// ==================================// Classification of Arithmetic Types// ==================================// In addtion to std::Integralconcept Float<typename T>   : std::DefaultConstructible<T>, std::CopyConstructible<T>,    std::LessThanComparable<T>, std::EqualityComparable<T>{  T operator+(T);  T operator+(T, T);  T& operator+=(T&, T);  T operator-(T, T);  T operator-(T);  T& operator-=(T&, T);  T operator*(T, T);  T& operator*=(T&, T);  T operator/(T, T);  T& operator/=(T&, T);  // TBD: Some day, these will come from LessThanComparable,  // EqualityComparable, etc.  bool operator>(T, T);  bool operator<=(T, T);  bool operator>=(T, T);  bool operator!=(T, T);  requires std::Assignable<T>, std::SameType<std::Assignable<T>::result_type, T&>;}concept_map Float<float> {}concept_map Float<double> {}concept_map Float<long double> {}// The difference to Float is the lack of LessThanComparableconcept Complex<typename T>   : std::DefaultConstructible<T>, std::CopyConstructible<T>,    std::EqualityComparable<T>{  T operator+(T);  T operator+(T, T);  T& operator+=(T&, T);  T operator-(T, T);  T operator-(T);  T& operator-=(T&, T);  T operator*(T, T);  T& operator*=(T&, T);  T operator/(T, T);  T& operator/=(T&, T);  // TBD: Some day, these will come from EqualityComparable  bool operator!=(T, T);  requires std::Assignable<T>, std::SameType<std::Assignable<T>::result_type, T&>;}template <typename T>  requires Float<T>concept_map Complex<std::complex<T> > {}// TBD: Concept Arithmetic is useless like this, it should have operations and then be the base for// Integral, Float and Complexconcept Arithmetic<typename T> {}template <typename T>  requires std::Integral<T>concept_map Arithmetic<T> {}template <typename T>  requires Float<T>concept_map Arithmetic<T> {}template <typename T>  requires Arithmetic<T>concept_map Arithmetic< std::complex<T> > {}// ================// Utility Concepts// ================// Concepts for functions mapping to same type or convertibleauto concept UnaryIsoFunction<typename Operation, typename Element>{    requires std::Callable1<Operation, Element>;    requires std::Convertible<std::Callable1<Operation, Element>::result_type, Element>;    typename result_type = std::Callable1<Operation, Element>::result_type;};auto concept BinaryIsoFunction<typename Operation, typename Element>{    requires std::Callable2<Operation, Element, Element>;    requires std::Convertible<std::Callable2<Operation, Element, Element>::result_type, Element>;    typename result_type = std::Callable2<Operation, Element, Element>::result_type;};#if 0auto concept CompatibleBinaryFunction<typename A1, typename A2, typename Result>{    typename result_type;    result_type F(A1, A2);    requires std::Convertible<result_type, Result>;}#endif// ==================// Algebraic Concepts// ==================auto concept Magma<typename Operation, typename Element>    : BinaryIsoFunction<Operation, Element>{    requires std::Assignable<Element>;    requires std::Assignable<Element, BinaryIsoFunction<Operation, Element>::result_type>;};// For algebraic structures that are commutative but not associative// As an example floating point numbers are commutative but not associative//   w.r.t. addition and multiplicationauto concept CommutativeMagma<typename Operation, typename Element>  : Magma<Operation, Element>,     algebra::Commutative<Operation, Element>{};// SemiGroup is a refinement which must be nominalauto concept SemiGroup<typename Operation, typename Element>  : Magma<Operation, Element>,     algebra::SemiGroup<Operation, Element>{};auto concept CommutativeSemiGroup<typename Operation, typename Element>  : SemiGroup<Operation, Element>,    CommutativeMagma<Operation, Element>{};// Adding identity// auto concept Monoid<typename Operation, typename Element>  : SemiGroup<Operation, Element>,     algebra::Monoid<Operation, Element> {    requires std::Convertible<identity_result_type, Element>;};auto concept CommutativeMonoid<typename Operation, typename Element>  : CommutativeSemiGroup<Operation, Element>,     Monoid<Operation, Element>{};concept PartiallyInvertibleMonoid<typename Operation, typename Element>  : Monoid<Operation, Element>,     algebra::Inversion<Operation, Element> {    typename is_invertible_result_type;    is_invertible_result_type is_invertible(Operation, Element);    requires std::Convertible<is_invertible_result_type, bool>;    requires std::Convertible<inverse_result_type, Element>;    // Does it overwrites the axiom from algebra::Inversion    axiom Inversion(Operation op, Element x)    {	// Only for invertible elements:	if (is_invertible(op, x))	    op( x, inverse(op, x) ) == identity(op, x); 	if ( is_invertible(op, x) )	    op( inverse(op, x), x ) == identity(op, x);     }};auto concept PartiallyInvertibleCommutativeMonoid<typename Operation, typename Element>  : PartiallyInvertibleMonoid<Operation, Element>,     CommutativeMonoid<Operation, Element>   {};concept Group<typename Operation, typename Element>  : PartiallyInvertibleMonoid<Operation, Element>,    algebra::Group<Operation, Element>{    axiom AlwaysInvertible(Operation op, Element x)    {	is_invertible(op, x);    }    axiom Inversion(Operation op, Element x)    {	// In fact this is implied by AlwaysInvertible and inherited Inversion axiom	// However, we don't rely on the compiler to deduce this	op( x, inverse(op, x) ) == identity(op, x);	op( inverse(op, x), x ) == identity(op, x);    }};auto concept AbelianGroup<typename Operation, typename Element>  : Group<Operation, Element>,     PartiallyInvertibleCommutativeMonoid<Operation, Element>,    algebra::AbelianGroup<Operation, Element>{};// ========================// Additive scalar concepts// ========================concept AdditiveMagma<typename Element>  : Magma< math::add<Element>, Element >{    typename plus_assign_result_type;      plus_assign_result_type operator+=(Element& x, Element y);    // requires std::Convertible<plus_assign_result_type, Element>;    // Operator + is by default defined with +=    typename addition_result_type;      addition_result_type operator+(Element x, Element y);#if 0    {	Element tmp(x);	return tmp += y;                      defaults NYS    }#endif     requires std::Convertible<addition_result_type, Element>;    // Type consistency with Magma    requires std::Convertible< addition_result_type,                               Magma< math::add<Element>, Element >::result_type>;    // SameType requires more rigorous specializations on pure algebraic functors    // requires std::SameType< addition_result_type,     // 	                    Magma< math::add<Element>, Element >::result_type>;    axiom Consistency(math::add<Element> op, Element x, Element y)    {	op(x, y) == x + y;                    	// Consistency definition between + and += might change later        x + y == x += y;	// Element tmp = x; tmp+= y; tmp == x + y; not proposal-compliant    }   }auto concept AdditiveCommutativeMagma<typename Element>  : AdditiveMagma<Element>,    CommutativeMagma< math::add<Element>, Element >{};auto concept AdditiveSemiGroup<typename Element>  : AdditiveMagma<Element>,     SemiGroup< math::add<Element>, Element >{};// We really need only one of the additive concepts for the requirements, // the requirements of the other would be implied.// Vice versa, to derive concept maps of nested concepts from// concept maps of refined concepts, they are needed all.auto concept AdditiveCommutativeSemiGroup<typename Element>  : AdditiveSemiGroup<Element>,    AdditiveCommutativeMagma<Element>,    CommutativeSemiGroup< math::add<Element>, Element >{};concept AdditiveMonoid<typename Element>  : AdditiveSemiGroup<Element>,    Monoid< math::add<Element>, Element >{    Element zero(Element v);    axiom Consistency (math::add<Element> op, Element x)    {	zero(x) == identity(op, x);    }};// We really need only one of the additive concepts for the requirements, // the requirements of the other would be implied.// Vice versa, to derive concept maps of nested concepts from// concept maps of refined concepts, they are needed all.auto concept AdditiveCommutativeMonoid<typename Element>  : AdditiveMonoid<Element>,    AdditiveCommutativeSemiGroup<Element>,    CommutativeMonoid< math::add<Element>, Element >{};concept AdditivePartiallyInvertibleMonoid<typename Element>  : AdditiveMonoid<Element>,    PartiallyInvertibleMonoid< math::add<Element>, Element >{    typename minus_assign_result_type;      minus_assign_result_type operator-=(Element& x, Element y);    // requires std::Convertible<minus_assign_result_type, Element>;         // Operator - by default defined with -=    typename subtraction_result_type;      subtraction_result_type operator-(Element& x, Element y);#if 0    {	Element tmp(x);	return tmp -= y;                      defaults NYS    }#endif     requires std::Convertible<subtraction_result_type, Element>;    typename unary_result_type;      unary_result_type operator-(Element x);#if 0    {	return zero(x) - x;      defaults NYS    }#endif     requires std::Convertible<unary_result_type, Element>;        axiom Consistency(math::add<Element> op, Element x, Element y)    {	// consistency between additive and pure algebraic concept	if ( is_invertible(op, y) )	    op(x, inverse(op, y)) == x - y;            	if ( is_invertible(op, y) )	    inverse(op, y) == -y;                      	// consistency between unary and binary -	if ( is_invertible(op, x) )	    identity(op, x) - x == -x;                 	// Might change later	if ( is_invertible(op, y) )	    x - y == x -= y;                                                       	// Element tmp = x; tmp-= y; tmp == x - y; not proposal-compliant    }  };auto concept AdditivePartiallyInvertibleCommutativeMonoid<typename Element>  : AdditivePartiallyInvertibleMonoid<Element>,    AdditiveCommutativeMonoid<Element>,     PartiallyInvertibleCommutativeMonoid< math::add<Element>, Element >{};auto concept AdditiveGroup<typename Element>  : AdditivePartiallyInvertibleMonoid<Element>,    Group< math::add<Element>, Element >{};auto concept AdditiveAbelianGroup<typename Element>  : AdditiveGroup<Element>,    AdditiveCommutativeMonoid<Element>,    AbelianGroup< math::add<Element>, Element >{};// ============================// Multiplitive scalar concepts// ============================concept MultiplicativeMagma<typename Element>  : Magma< math::mult<Element>, Element >{    typename mult_assign_result_type;      mult_assign_result_type operator*=(Element& x, Element y);    // requires std::Convertible<mult_assign_result_type, Element>;    // Operator * is by default defined with *=    typename mult_result_type;      mult_result_type operator*(Element x, Element y);#if 0    {	Element tmp(x);	return tmp *= y;                      defaults NYS    }#endif     requires std::Convertible<mult_result_type, Element>;        // Type consistency with Magma    requires std::Convertible< mult_result_type,                               Magma< math::mult<Element>, Element >::result_type>;    // SameType requires more rigorous specializations on pure algebraic functors    // requires std::SameType< mult_result_type,     // 	                    Magma< math::mult<Element>, Element >::result_type>;    axiom Consistency(math::mult<Element> op, Element x, Element y)    {	op(x, y) == x * y;                    	// Consistency definition between * and *= might change later        x * y == x *= y;	// Element tmp = x; tmp*= y; tmp == x * y; not proposal-compliant    }  }auto concept MultiplicativeSemiGroup<typename Element>  : MultiplicativeMagma<Element>,    SemiGroup< math::mult<Element>, Element >{};auto concept MultiplicativeCommutativeSemiGroup<typename Element>  : MultiplicativeSemiGroup<Element>,    CommutativeSemiGroup< math::mult<Element>, Element >{};concept MultiplicativeMonoid<typename Element>  : MultiplicativeSemiGroup<Element>,    Monoid< math::mult<Element>, Element >{    Element one(Element v);    axiom Consistency (math::mult<Element> op, Element x)    {	one(x) == identity(op, x);    }};auto concept MultiplicativeCommutativeMonoid<typename Element>  : MultiplicativeMonoid<Element>,    MultiplicativeCommutativeSemiGroup<Element>,    CommutativeMonoid< math::mult<Element>, Element >{};concept MultiplicativePartiallyInvertibleMonoid<typename Element>  : MultiplicativeMonoid<Element>,    PartiallyInvertibleMonoid< math::mult<Element>, Element >{    typename divide_assign_result_type;      divide_assign_result_type operator/=(Element& x, Element y);    // requires std::Convertible<divide_assign_result_type, Element>;         // Operator / by default defined with /=    typename division_result_type = Element;      division_result_type operator/(Element x, Element y);#if 0    {	Element tmp(x);	return tmp /= y;                      defaults NYS    }#endif     requires std::Convertible<division_result_type, Element>;        axiom Consistency(math::mult<Element> op, Element x, Element y)    {	// consistency between multiplicative and pure algebraic concept	if ( is_invertible(op, y) )	    op(x, inverse(op, y)) == x / y;            	// Consistency between / and /=, might change later	if ( is_invertible(op, y) )	    x / y == x /= y;              	// Element tmp = x; tmp/= y; tmp == x / y; not proposal-compliant     }  }; auto concept MultiplicativePartiallyInvertibleCommutativeMonoid<typename Element>  : MultiplicativePartiallyInvertibleMonoid<Element>,    MultiplicativeCommutativeMonoid<Element>,    PartiallyInvertibleCommutativeMonoid< math::mult<Element>, Element >{}; auto concept MultiplicativeGroup<typename Element>  : MultiplicativeMonoid<Element>,    Group< math::mult<Element>, Element >{};auto concept MultiplicativeAbelianGroup<typename Element>  : MultiplicativeGroup<Element>,    MultiplicativeCommutativeMonoid<Element>,    AbelianGroup< math::mult<Element>, Element >{};