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: Brice Boyer (briceboyer) <boyer.brice@gmail.com>
8//          A. Breust (taken from FFLAS-FFPACK)
9// ==========================================================================
10
11#ifndef __GIVARO_modular_float_INL
12#define __GIVARO_modular_float_INL
13
14namespace Givaro {
15
16#define MOD Modular<Storage_t, Compute_t, typename std::enable_if<std::is_floating_point<Storage_t>::value>::type>
17#define TMPL template<typename Storage_t, typename Compute_t>
18#define COND_TMPL(T, ...) \
19    template<typename T, \
20    typename std::enable_if<(__VA_ARGS__), int>::type*>
21
22    // --------------------
23    // ----- Initialisation
24
25    TMPL
26    inline typename MOD::Element&
27    MOD::init(Element& a) const
28    {
29        return a = this->zero;
30    }
31
32    TMPL
33    COND_TMPL(Source,
34              std::is_same<Source, double>::value && std::is_same<Storage_t, float>::value)
35    inline typename MOD::Element&
36    MOD::init(Element& r, const Source a) const
37    {
38        r = Caster<Element>(std::fmod(a, _pc));
39        if (r < 0.f) r += _pc;
40        return r;
41    }
42
43    TMPL
44    COND_TMPL(Source,
45              std::is_integral<Source>::value && std::is_signed<Source>::value && sizeof(Source) >= sizeof(Storage_t))
46    inline typename MOD::Element&
47    MOD::init(Element& r, const Source a) const
48    {
49        r = Caster<Element>(std::abs(a) % Caster<Source>(_p));
50        if (a < 0) negin(r);
51        return r;
52    }
53
54    TMPL
55    COND_TMPL(Source,
56              std::is_integral<Source>::value && std::is_unsigned<Source>::value && sizeof(Source) >= sizeof(Storage_t))
57    inline typename MOD::Element&
58    MOD::init(Element& r, const Source a) const
59    {
60        return r = Caster<Element>(a % Caster<Source>(_p));
61    }
62
63    TMPL
64    inline typename MOD::Element&
65    MOD::init(Element& r, const Integer& a) const
66    {
67        r = Caster<Element>(a % _p);
68        if (r < 0) r += _pc;
69        return r;
70    }
71
72    // ------------
73    // ----- Reduce
74
75    TMPL
76    inline typename MOD::Element&  MOD::reduce (Element& x) const
77    {
78        x = std::fmod(x, Caster<Element>(_pc));
79        if (x < 0) x += Caster<Element>(_pc);
80        return x;
81    }
82
83    TMPL
84    inline typename MOD::Element&  MOD::reduce (Element& x, const Element& y) const
85    {
86        x = std::fmod(y, Caster<Element>(_pc));
87        if (x < 0.f) x += Caster<Element>(_pc);
88        return x;
89    }
90
91    // ------------------------
92    // ----- Classic arithmetic
93
94    TMPL
95    inline typename MOD::Element & MOD::add
96    (Element &x, const Element &y, const Element &z) const
97    {
98        Compute_t tmp = Caster<Compute_t>(y) + Caster<Compute_t>(z);
99        return x = Caster<Element>(tmp<_pc?tmp:tmp-_pc);
100    }
101
102    TMPL
103    inline typename MOD::Element & MOD::sub
104    (Element &x, const Element &y, const Element &z) const
105    {
106        return x = (y>=z)? y-z:(Caster<Element>(_pc)-z)+y;
107    }
108
109    TMPL
110    inline typename MOD::Element & MOD::mul
111    (Element &x, const Element &y, const Element &z) const
112    {
113        return x = Caster<Element>(std::fmod(Caster<Compute_t>(y)*Caster<Compute_t>(z), _pc));
114    }
115
116    TMPL
117    inline typename MOD::Element & MOD::div
118    (Element &x, const Element &y, const Element &z) const
119    {
120        return mulin(inv(x, z), y);
121    }
122
123    TMPL
124    inline typename MOD::Element & MOD::neg
125    (Element &x, const Element &y) const
126    {
127        return x = (y==Caster<Element>(0)?Caster<Element>(0):Caster<Element>(_pc)-y);
128    }
129
130    TMPL
131    inline typename MOD::Element & MOD::inv
132    (Element &x, const Element &y) const
133    {
134        invext(x,y,Caster<Element>(_pc));
135        return (x<Caster<Element>(0) ? x+=Caster<Element>(_pc) : x);
136    }
137
138    TMPL
139    inline typename MOD::Element & MOD::addin
140    (Element &x, const Element &y) const
141    {
142        Compute_t tmp = Caster<Compute_t>(x) + Caster<Compute_t>(y);
143        return x = Caster<Element>(tmp < _pc? tmp : tmp - _pc);
144    }
145
146    TMPL
147    inline typename MOD::Element & MOD::subin
148    (Element &x, const Element &y) const
149    {
150        if (x<y) x += (Caster<Element>(_pc)-y);
151        else     x -= y;
152        return x;
153    }
154
155    TMPL
156    inline typename MOD::Element & MOD::mulin
157    (Element &x, const Element &y) const
158    {
159        return x = Caster<Element>(std::fmod(Caster<Compute_t>(x)*Caster<Compute_t>(y), _pc));
160    }
161
162    TMPL
163    inline typename MOD::Element & MOD::divin
164    (Element &x, const Element &y) const
165    {
166        Element iy;
167        return mulin(x, inv(iy, y));
168    }
169
170    TMPL
171    inline typename MOD::Element & MOD::negin
172    (Element &x) const
173    {
174        return x = (x==Caster<Element>(0)?Caster<Element>(0):Caster<Element>(_pc)-x);
175    }
176
177    TMPL
178    inline typename MOD::Element & MOD::invin
179    (Element &x) const
180    {
181        return inv(x, x);
182    }
183
184    // -- axpy: r <- a * x + y
185    TMPL
186    inline typename MOD::Element & MOD::axpy
187    (Element &r, const Element &a, const Element &x, const Element &y) const
188    {
189        return r = Caster<Element>(std::fmod(Caster<Compute_t>(a)*Caster<Compute_t>(x)+Caster<Compute_t>(y), _pc));
190    }
191
192    TMPL
193    inline typename MOD::Element & MOD::axpyin
194    (Element &r, const Element &a, const Element &x) const
195    {
196        return r = Caster<Element>(std::fmod(Caster<Compute_t>(a)*Caster<Compute_t>(x)+Caster<Compute_t>(r), _pc));
197    }
198
199    // -- axmy: r <- a * x - y
200    TMPL
201    inline typename MOD::Element & MOD::axmy
202    (Element& r, const Element &a, const Element &x, const Element &y) const
203    {
204        return r = Caster<Element>(std::fmod(Caster<Compute_t>(a)*Caster<Compute_t>(x) + (_pc-Caster<Compute_t>(y)), _pc));
205    }
206
207    TMPL
208    inline typename MOD::Element & MOD::axmyin
209    (Element& r, const Element &a, const Element &x) const
210    {
211        maxpyin(r,a,x);
212        return negin(r);
213    }
214
215    // -- maxpy:   r <- y - a * x
216    TMPL
217    inline typename MOD::Element&  MOD::maxpy
218    (Element& r, const Element& a, const Element& x, const Element& y) const
219    {
220        r = Caster<Element>(std::fmod(Caster<Compute_t>(a) * Caster<Compute_t>(x)
221                                      + (_pc - Caster<Compute_t>(y)), _pc));
222        return negin(r);
223    }
224
225    TMPL
226    inline typename MOD::Element&  MOD::maxpyin
227    (Element& r, const Element& a, const Element& x) const
228    {
229        Compute_t tmp = Caster<Compute_t>(a) * Caster<Compute_t>(x) + (_pc - Caster<Compute_t>(r));
230        r = (tmp < _pc) ? Caster<Element>(tmp) : Caster<Element>(std::fmod(tmp, _pc));
231        return negin(r);
232    }
233
234#undef MOD
235#undef TMPL
236#undef COND_TMPL
237
238} // namespace Givaro
239
240#endif // __GIVARO_modular_float_INL
241
242/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
243// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
244