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