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: Clement Pernet <clement.pernet@gmail.com>
8//          A. Breust (taken from FFLAS-FFPACK)
9// ==========================================================================
10
11#ifndef __GIVARO_modular_balanced_float_INL
12#define __GIVARO_modular_balanced_float_INL
13
14#include <cmath> // fmod
15
16#define NORMALISE(x)				\
17{						\
18    if (x < _mhalfp) x += _p;		\
19    else if (x > _halfp) x -= _p;		\
20}
21
22#define NORMALISE_HI(x)				\
23{						\
24    if (x > _halfp) x -= _p;		\
25}
26
27namespace Givaro
28{
29
30    //----- Classic arithmetic
31
32    inline ModularBalanced<float>::Element&
33    ModularBalanced<float>::mul(Element& r, const Element& a, const Element& b) const
34    {
35        return reduce(r = a * b);
36    }
37
38    inline ModularBalanced<float>::Element&
39    ModularBalanced<float>::div(Element& r, const Element& a, const Element& b) const
40    {
41        Element tmp;
42        return mul (r, a, inv(tmp, b));
43    }
44
45    inline ModularBalanced<float>::Element&
46    ModularBalanced<float>::add(Element& r, const Element& a, const Element& b) const
47    {
48        r = a + b;
49        NORMALISE(r);
50        return r;
51    }
52
53    inline ModularBalanced<float>::Element&
54    ModularBalanced<float>::sub(Element& r, const Element& a, const Element& b) const
55    {
56        r = a - b;
57        NORMALISE(r);
58        return r;
59    }
60
61    inline ModularBalanced<float>::Element&
62    ModularBalanced<float>::neg(Element& r, const Element& a) const
63    {
64        return r = -a;
65    }
66
67    inline ModularBalanced<float>::Element&
68    ModularBalanced<float>::inv(Element& r, const Element& a) const
69    {
70        r = invext(a, _p);
71        NORMALISE(r);
72        return r;
73    }
74
75    inline bool ModularBalanced<float>::isUnit(const Element& a) const
76    {
77        Element u,d;
78        extended_euclid(u,d,a,_p);
79        return isOne(d) || isMOne(d);
80    }
81
82    inline ModularBalanced<float>::Element&
83    ModularBalanced<float>::mulin(Element& r, const Element& a) const
84    {
85        return mul(r, r, a);
86    }
87
88    inline ModularBalanced<float>::Element&
89    ModularBalanced<float>::divin(Element& r, const Element& a) const
90    {
91        return div(r, r, a);
92    }
93
94    inline ModularBalanced<float>::Element&
95    ModularBalanced<float>::addin(Element& r, const Element& a) const
96    {
97        return add(r, r, a);
98    }
99
100    inline ModularBalanced<float>::Element&
101    ModularBalanced<float>::subin(Element& r, const Element& a) const
102    {
103        return sub(r, r, a);
104    }
105
106    inline ModularBalanced<float>::Element&
107    ModularBalanced<float>::negin(Element& r) const
108    {
109        return neg(r, r);
110    }
111
112    inline ModularBalanced<float>::Element&
113    ModularBalanced<float>::invin(Element& r) const
114    {
115        return inv(r, r);
116    }
117
118    //----- Special arithmetic
119
120    inline ModularBalanced<float>::Element&
121    ModularBalanced<float>::axpy(Element& r, const Element& a, const Element& x, const Element& y) const
122    {
123        return reduce(r = a * x + y);
124    }
125
126    inline ModularBalanced<float>::Element&
127    ModularBalanced<float>::axpyin(Element& r, const Element& a, const Element& x) const
128    {
129        return reduce(r += a * x);
130    }
131
132    inline ModularBalanced<float>::Element&
133    ModularBalanced<float>::axmy(Element& r, const Element& a, const Element& x, const Element& y) const
134    {
135        return reduce(r = a * x - y);
136    }
137
138    inline ModularBalanced<float>::Element&
139    ModularBalanced<float>::axmyin(Element& r, const Element& a, const Element& x) const
140    {
141        return reduce(r = a * x - r);
142    }
143
144    inline ModularBalanced<float>::Element&
145    ModularBalanced<float>::maxpy(Element& r, const Element& a, const Element& x, const Element& y) const
146    {
147        return reduce(r = y - a * x);
148    }
149
150    inline ModularBalanced<float>::Element&
151    ModularBalanced<float>:: maxpyin(Element& r, const Element& a, const Element& x) const
152    {
153        return reduce(r -= a * x);
154    }
155
156    //----- Init
157
158    inline ModularBalanced<float>::Element&
159    ModularBalanced<float>::init(Element& x) const
160    {
161        return x = 0;
162    }
163
164    inline ModularBalanced<float>::Element&
165    ModularBalanced<float>::init(Element& x, const float y) const
166    {
167        x = static_cast<Element>(fmodf(y, _p));
168        NORMALISE(x);
169        return x;
170    }
171
172    inline ModularBalanced<float>::Element&
173    ModularBalanced<float>::init(Element& x, const double y) const
174    {
175        x = static_cast<Element>(fmod(y, double(_p)));
176        NORMALISE(x);
177        return x;
178    }
179
180    inline ModularBalanced<float>::Element&
181    ModularBalanced<float>::init(Element& x, const int32_t y) const
182    {
183        x = static_cast<Element>(y % static_cast<int32_t>(_up));
184        NORMALISE(x);
185        return x;
186    }
187
188    inline ModularBalanced<float>::Element&
189    ModularBalanced<float>::init(Element& x, const uint32_t y) const
190    {
191        x = static_cast<Element>(y %  static_cast<uint32_t>(_up));
192        NORMALISE_HI(x);
193        return x;
194    }
195
196    inline ModularBalanced<float>::Element&
197    ModularBalanced<float>::init(Element& x, const int64_t y) const
198    {
199        x = static_cast<Element>(y % static_cast<int64_t>(_up));
200        NORMALISE(x);
201        return x;
202    }
203
204    inline ModularBalanced<float>::Element&
205    ModularBalanced<float>::init(Element& x, const uint64_t y) const
206    {
207        x = static_cast<Element>(y % static_cast<uint64_t>(_up));
208        NORMALISE_HI(x);
209        return x;
210    }
211
212    inline ModularBalanced<float>::Element&
213    ModularBalanced<float>::init(Element& x, const Integer& y) const
214    {
215        x = static_cast<Element>(y % _p);
216        NORMALISE_HI(x);
217        return x;
218    }
219
220    inline ModularBalanced<float>::Element&
221    ModularBalanced<float>::assign(Element& x, const Element& y) const
222    {
223        return x = y;
224    }
225
226    //----- Reduce
227
228    inline ModularBalanced<float>::Element&
229    ModularBalanced<float>::reduce(Element& x, const Element& y) const
230    {
231        x = fmodf(y, _p);
232        NORMALISE(x);
233        return x;
234    }
235
236    inline ModularBalanced<float>::Element&
237    ModularBalanced<float>::reduce(Element& x) const
238    {
239        x = fmodf(x, _p);
240        NORMALISE(x);
241        return x;
242    }
243
244    //----- IO
245
246    inline std::ostream&
247    ModularBalanced<float>::write(std::ostream& os) const
248    {
249        return os << "ModularBalanced<float> mod " << _p;
250    }
251
252    inline std::ostream&
253    ModularBalanced<float>::write(std::ostream& os, const Element& x) const
254    {
255        return os << x;
256    }
257
258    inline std::istream&
259    ModularBalanced<float>::read(std::istream& is, Element& x) const
260    {
261        Element tmp;
262        is >> tmp;
263        init(x, tmp);
264        return is;
265    }
266
267}
268
269#undef NORMALISE
270#undef NORMALISE_HI
271
272#endif
273/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
274// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
275