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_int64_INL 12#define __GIVARO_modular_balanced_int64_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<int64_t>::Element& 33 ModularBalanced<int64_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<int64_t>::Element& 42 ModularBalanced<int64_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<int64_t>::Element& 49 ModularBalanced<int64_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<int64_t>::Element& 57 ModularBalanced<int64_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<int64_t>::Element& 65 ModularBalanced<int64_t>::neg(Element& r, const Element& a) const 66 { 67 return r = -a; 68 } 69 70 inline ModularBalanced<int64_t>::Element& 71 ModularBalanced<int64_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<int64_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<int64_t>::Element& 86 ModularBalanced<int64_t>::mulin(Element& r, const Element& a) const 87 { 88 return mul(r, r, a); 89 } 90 91 inline ModularBalanced<int64_t>::Element& 92 ModularBalanced<int64_t>::divin(Element& r, const Element& a) const 93 { 94 return div(r, r, a); 95 } 96 97 inline ModularBalanced<int64_t>::Element& 98 ModularBalanced<int64_t>::addin(Element& r, const Element& a) const 99 { 100 return add(r, r, a); 101 } 102 103 inline ModularBalanced<int64_t>::Element& 104 ModularBalanced<int64_t>::subin(Element& r, const Element& a) const 105 { 106 return sub(r, r, a); 107 } 108 109 inline ModularBalanced<int64_t>::Element& 110 ModularBalanced<int64_t>::negin(Element& r) const 111 { 112 return neg(r, r); 113 } 114 115 inline ModularBalanced<int64_t>::Element& 116 ModularBalanced<int64_t>::invin(Element& r) const 117 { 118 return inv(r, r); 119 } 120 121 //----- Special arithmetic 122 123 inline ModularBalanced<int64_t>::Element& 124 ModularBalanced<int64_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<int64_t>::Element& 134 ModularBalanced<int64_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<int64_t>::Element& 144 ModularBalanced<int64_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<int64_t>::Element& 154 ModularBalanced<int64_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<int64_t>::Element& 164 ModularBalanced<int64_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<int64_t>::Element& 170 ModularBalanced<int64_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<int64_t>::Element& 178 ModularBalanced<int64_t>::init(Element& x) const 179 { 180 return x = 0; 181 } 182 183 inline ModularBalanced<int64_t>::Element& 184 ModularBalanced<int64_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<int64_t>::Element& 192 ModularBalanced<int64_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<int64_t>::Element& 200 ModularBalanced<int64_t>::init(Element& x, const Integer& y) const 201 { 202 x = static_cast<Element>(y % _p); 203 NORMALISE_HI(x); 204 return x; 205 } 206 207 inline ModularBalanced<int64_t>::Element& 208 ModularBalanced<int64_t>::assign(Element& x, const Element& y) const 209 { 210 return x = y; 211 } 212 213 //----- Reduce 214 215 inline ModularBalanced<int64_t>::Element& 216 ModularBalanced<int64_t>::reduce(Element& x, const Element& y) const 217 { 218 x = y % _p; 219 NORMALISE(x); 220 return x; 221 } 222 223 inline ModularBalanced<int64_t>::Element& 224 ModularBalanced<int64_t>::reduce(Element& x) const 225 { 226 x %= _p; 227 NORMALISE(x); 228 return x; 229 } 230 231 //----- IO 232 233 inline std::ostream& 234 ModularBalanced<int64_t>::write(std::ostream& os) const 235 { 236 return os << "ModularBalanced<int64_t> mod " << _p; 237 } 238 239 inline std::ostream& 240 ModularBalanced<int64_t>::write(std::ostream& os, const Element& x) const 241 { 242 return os << x; 243 } 244 245 inline std::istream& 246 ModularBalanced<int64_t>::read(std::istream& is, Element& x) const 247 { 248 Element tmp; 249 is >> tmp; 250 init(x, tmp); 251 return is; 252 } 253 254} 255 256#undef NORMALISE 257#undef NORMALISE_HI 258 259#endif 260/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 261// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 262