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