1 // ========================================================================== 2 // Copyright(c)'1994-2015 by The Givaro group 3 // This file is part of Givaro. 4 // Givaro is governed by the CeCILL-B license under French law 5 // and abiding by the rules of distribution of free software. 6 // see the COPYRIGHT file for more details. 7 // Authors: Pascal Giorgi <pascal.giorgi@ens-lyon.fr> 8 // Clement Pernet <clement.pernet@gmail.com> 9 // Brice Boyer (briceboyer) <boyer.brice@gmail.com> (modified) 10 // A. Breust (taken from FFLAS-FFPACK) 11 // ========================================================================== 12 13 #ifndef __GIVARO_modular_balanced_double_H 14 #define __GIVARO_modular_balanced_double_H 15 16 #include "givaro/givinteger.h" 17 #include "givaro/givcaster.h" 18 #include "givaro/givranditer.h" 19 #include "givaro/ring-interface.h" 20 #include "givaro/modular-general.h" 21 22 #include <iostream> 23 24 namespace Givaro 25 { 26 template<class TAG> class ModularBalanced; 27 28 template <> 29 class ModularBalanced<double> 30 { 31 public: 32 33 // ----- Exported types 34 using Self_t = ModularBalanced<double>; 35 using Element = double; 36 using Element_ptr = Element*; 37 using ConstElement = const Element; 38 using ConstElement_ptr = const Element*; 39 using Residu_t = double; 40 enum { size_rep = sizeof(Element) }; 41 42 // ----- Constantes 43 const Element zero = 0.0; 44 const Element one = 1.0; 45 const Element mOne = -1.0; 46 47 // ----- Constructors ModularBalanced()48 ModularBalanced() 49 : _p(0), _halfp(0), _mhalfp(0), _up(0.0) 50 {} 51 ModularBalanced(double p)52 ModularBalanced(double p) 53 : _p(p) 54 , _halfp((_p - 1.f) / 2.f) 55 , _mhalfp(_halfp - _p + 1.f) 56 , _up(static_cast<uint32_t>(_p)) 57 { 58 assert(_p >= minCardinality()); 59 assert(_p <= maxCardinality()); 60 } 61 ModularBalanced(const Self_t & F)62 ModularBalanced(const Self_t& F) 63 : _p(F._p), _halfp(F._halfp), _mhalfp(F._mhalfp), _up(F._up) 64 {} 65 66 // ----- Accessors minElement()67 inline Element minElement() const { return _mhalfp; } maxElement()68 inline Element maxElement() const { return _halfp; } 69 70 // ----- Access to the modulus residu()71 inline Residu_t residu() const { return _p; } size()72 inline Residu_t size() const { return _p; } characteristic()73 inline Residu_t characteristic() const { return _p; } cardinality()74 inline Residu_t cardinality() const { return _p; } characteristic(T & p)75 template<class T> inline T& characteristic(T& p) const { return Caster(p,_p); } cardinality(T & p)76 template<class T> inline T& cardinality(T& p) const { return Caster(p,_p); } 77 maxCardinality()78 static inline Residu_t maxCardinality() { return 134217727; } // 2^12.5 minCardinality()79 static inline Residu_t minCardinality() { return 3.f; } 80 81 // ----- Checkers isZero(const Element & a)82 inline bool isZero(const Element& a) const { return a == zero; } isOne(const Element & a)83 inline bool isOne (const Element& a) const { return a == one; } isMOne(const Element & a)84 inline bool isMOne(const Element& a) const { return a == mOne; } 85 inline bool isUnit(const Element& a) const; areEqual(const Element & a,const Element & b)86 inline bool areEqual(const Element& a, const Element& b) const { return a == b; } length(const Element a)87 inline size_t length(const Element a) const { return size_rep; } 88 89 // ----- Ring-wise operators 90 inline bool operator==(const Self_t& F) const { return _p == F._p; } 91 inline bool operator!=(const Self_t& F) const { return _p != F._p; } 92 inline Self_t& operator=(const Self_t& F) 93 { 94 F.assign(const_cast<Element&>(one), F.one); 95 F.assign(const_cast<Element&>(zero), F.zero); 96 F.assign(const_cast<Element&>(mOne), F.mOne); 97 _p = F._p; 98 _halfp = F._halfp; 99 _mhalfp = F._mhalfp; 100 _up = F._up; 101 return *this; 102 } 103 104 // ----- Initialisation 105 Element& init(Element& a) const; 106 Element& init(Element& r, const float a) const; 107 Element& init(Element& r, const double a) const; 108 Element& init(Element& r, const int64_t a) const; 109 Element& init(Element& r, const uint64_t a) const; 110 Element& init(Element& r, const Integer& a) const; init(Element & r,const T & a)111 template<typename T> Element& init(Element& r, const T& a) const 112 { r = Caster<Element>(a); return reduce(r); } 113 114 Element& assign(Element& r, const Element& a) const; 115 116 // ----- Convert convert(T & r,const Element & a)117 template<typename T> T& convert(T& r, const Element& a) const 118 { return r = static_cast<T>(a); } 119 120 Element& reduce(Element& r, const Element& a) const; 121 Element& reduce(Element& r) const; 122 123 // ----- Classic arithmetic 124 Element& mul(Element& r, const Element& a, const Element& b) const; 125 Element& div(Element& r, const Element& a, const Element& b) const; 126 Element& add(Element& r, const Element& a, const Element& b) const; 127 Element& sub(Element& r, const Element& a, const Element& b) const; 128 Element& neg(Element& r, const Element& a) const; 129 Element& inv(Element& r, const Element& a) const; 130 131 Element& mulin(Element& r, const Element& a) const; 132 Element& divin(Element& r, const Element& a) const; 133 Element& addin(Element& r, const Element& a) const; 134 Element& subin(Element& r, const Element& a) const; 135 Element& negin(Element& r) const; 136 Element& invin(Element& r) const; 137 138 // -- axpy: r <- a * x + y 139 // -- axpyin: r <- a * x + r 140 Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const; 141 Element& axpyin(Element& r, const Element& a, const Element& x) const; 142 143 // -- axmy: r <- a * x - y 144 // -- axmyin: r <- a * x - r 145 Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const; 146 Element& axmyin(Element& r, const Element& a, const Element& x) const; 147 148 // -- maxpy: r <- y - a * x 149 // -- maxpyin: r <- r - a * x 150 Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const; 151 Element& maxpyin(Element& r, const Element& a, const Element& x) const; 152 153 // ----- Random generators 154 typedef ModularRandIter<Self_t> RandIter; 155 typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter; random(Random & g,Element & r)156 template< class Random > Element& random(Random& g, Element& r) const 157 { return init(r, g()); } nonzerorandom(Random & g,Element & a)158 template< class Random > Element& nonzerorandom(Random& g, Element& a) const 159 { while (isZero(init(a, g()))) 160 ; 161 return a; } 162 163 // --- IO methods 164 std::ostream& write(std::ostream& s) const; 165 std::istream& read (std::istream& s, Element& a) const; 166 std::ostream& write(std::ostream& s, const Element& a) const; 167 168 protected: 169 170 Element _p; 171 Element _halfp; 172 Element _mhalfp; 173 uint64_t _up; 174 }; 175 } 176 177 #include "givaro/modular-balanced-double.inl" 178 179 #endif 180 181 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 182 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 183