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