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