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