1 /* field/modular-balanced-float.h 2 * Copyright (C) 2003 Pascal Giorgi 3 * 2005,2008 Clement Pernet 4 * Written by Clement Pernet <clement.pernet@gmail.com> 5 * Pascal Giorgi <pascal.giorgi@ens-lyon.fr> 6 * Modified Brice Boyer (briceboyer) <boyer.brice@gmail.com> 7 * ------------------------------------ 8 * 9 * 10 * ========LICENCE======== 11 * This file is part of the library LinBox. 12 * 13 * LinBox is free software: you can redistribute it and/or modify 14 * it under the terms of the GNU Lesser General Public 15 * License as published by the Free Software Foundation; either 16 * version 2.1 of the License, or (at your option) any later version. 17 * 18 * This library is distributed in the hope that it will be useful, 19 * but WITHOUT ANY WARRANTY; without even the implied warranty of 20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 21 * Lesser General Public License for more details. 22 * 23 * You should have received a copy of the GNU Lesser General Public 24 * License along with this library; if not, write to the Free Software 25 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 26 * ========LICENCE======== 27 *. 28 */ 29 30 /*! @file field/modular/modular-balanced-float.h 31 * @ingroup field 32 * @brief Balanced representation of <code>Z/mZ</code> over \c float . 33 */ 34 35 #ifndef __LINBOX_modular_balanced_float_H 36 #define __LINBOX_modular_balanced_float_H 37 38 #ifdef __INTEL_COMPILER 39 #define FmodF fmodf 40 #else 41 #define FmodF fmod 42 #endif 43 44 #include "linbox/linbox-config.h" 45 #include "linbox/integer.h" 46 #include "linbox/ring/modular.h" 47 #include "linbox/vector/vector-domain.h" 48 #include "linbox/field/field-traits.h" 49 #include "linbox/util/field-axpy.h" 50 #include "linbox/util/debug.h" 51 #include <cmath> 52 #include "linbox/field/field-traits.h" 53 #include "linbox/randiter/modular-balanced.h" 54 55 #include <givaro/modular-balanced-float.h> 56 57 58 // Namespace in which all LinBox code resides 59 namespace LinBox 60 { 61 62 class MultiModFloat; 63 64 template <> 65 class FieldAXPY<Givaro::ModularBalanced<float> > { 66 public: 67 typedef float Element; 68 typedef float Abnormal; 69 typedef Givaro::ModularBalanced<Element> Field; 70 FieldAXPY(const Field & F)71 FieldAXPY (const Field &F) : 72 _field (&F), 73 _y(0.) , _bound( (Element) (((1UL << 24) - (unsigned long) (field().characteristic()*field().characteristic())))) 74 {} 75 FieldAXPY(const FieldAXPY & faxpy)76 FieldAXPY (const FieldAXPY &faxpy) : 77 _field (faxpy._field), 78 _y(faxpy._y), _bound(faxpy._bound) 79 {} 80 81 FieldAXPY<Givaro::ModularBalanced<Element> > &operator = (const FieldAXPY &faxpy) { 82 _field = faxpy._field; 83 _y= faxpy._y; 84 _bound= faxpy._bound; 85 return *this; 86 } 87 mulacc(const Element & a,const Element & x)88 inline Element& mulacc (const Element &a, const Element &x) { 89 // Element tmp= a*x; 90 // return accumulate(tmp); 91 return accumulate(a*x); 92 } 93 accumulate(const Element & tmp)94 inline Element& accumulate (const Element &tmp) { 95 _y += tmp; 96 if (_y > _bound) 97 return _y = fmodf (_y, field().characteristic()); 98 else 99 return _y; 100 } subumulate(const Element & tmp)101 inline Element& subumulate (const Element &tmp) { 102 _y -= tmp; 103 if (_y < 0) 104 return _y += field().characteristic(); 105 else 106 return _y; 107 } 108 get(Element & y)109 inline Element& get (Element &y) { 110 _y = fmodf (_y, field().characteristic()); 111 return y=_y ; 112 } 113 assign(const Element y)114 inline FieldAXPY &assign (const Element y) { 115 _y = y; 116 return *this; 117 } 118 reset()119 inline void reset() { 120 _y = 0.; 121 } 122 set(const Element & tmp)123 inline Element& set (const Element &tmp) { 124 _y = tmp; 125 if (_y > _bound) 126 return _y = fmodf (_y, field().characteristic()); 127 else 128 return _y; 129 } 130 field()131 inline const Field & field() const { return *_field; } 132 133 private: 134 const Field *_field; 135 Element _y; 136 Element _bound; 137 }; 138 139 140 template <> 141 class DotProductDomain<Givaro::ModularBalanced<float> > : public VectorDomainBase<Givaro::ModularBalanced<float> > { 142 public: 143 typedef float Element; DotProductDomain()144 DotProductDomain(){} DotProductDomain(const Givaro::ModularBalanced<Element> & F)145 DotProductDomain (const Givaro::ModularBalanced<Element> &F) : 146 VectorDomainBase<Givaro::ModularBalanced<Element> > (F), _bound( (Element) ( (1UL<<24) - (unsigned long) (field().characteristic()*field().characteristic()))) 147 { 148 _nmax= (size_t)floor((Element(1<<11)* Element(1<<11)*2.)/ (field().characteristic() * field().characteristic())); 149 } 150 151 using VectorDomainBase<Givaro::ModularBalanced<Element> >::field; 152 protected: 153 template <class Vector1, class Vector2> dotSpecializedDD(Element & res,const Vector1 & v1,const Vector2 & v2)154 inline Element &dotSpecializedDD (Element &res, const Vector1 &v1, const Vector2 &v2) const 155 { 156 157 Element y = 0.; 158 if (v1.size() < _nmax) { 159 for (size_t i = 0; i< v1.size();++i) 160 y += v1[i] * v2[i] ; 161 y = fmodf(y, field().characteristic()); 162 } 163 else{ 164 Element t = 0.; 165 size_t i=0; 166 for (;i< v1.size()- _nmax ;i=i+_nmax){ 167 for (size_t j=i;j<i+_nmax;++j) 168 y += v1[j] * v2[j]; 169 t+= fmodf(y, field().characteristic()); 170 y=0.; 171 } 172 for (;i < v1.size();++i) 173 y += v1[i] * v2[i]; 174 t+= fmodf(y, field().characteristic()); 175 y = fmodf(t, field().characteristic()); 176 } 177 //!@bug should not be neccessary (use assign) 178 return field().init(res, y); 179 } 180 181 template <class Vector1, class Vector2> dotSpecializedDSP(Element & res,const Vector1 & v1,const Vector2 & v2)182 inline Element &dotSpecializedDSP (Element &res, const Vector1 &v1, const Vector2 &v2) const 183 { 184 185 Element y = 0.; 186 187 188 if (v1.first.size() < _nmax) { 189 for (size_t i=0;i<v1.first.size();++i) 190 y+= v1.second[i] * v2[v1.first[i]]; 191 y = fmodf(y, field().characteristic()); 192 } 193 else { 194 Element t =0.; 195 size_t i=0; 196 for (;i< v1.first.size()- _nmax ;i=i+_nmax){ 197 for (size_t j=i;j<i+_nmax;++j) 198 y += v1.second[j] * v2[v1.first[j]]; 199 t+=fmodf(y, field().characteristic()); 200 y=0.; 201 } 202 for (;i < v1.first.size();++i) 203 y += v1.second[i] * v2[v1.first[i]]; 204 t+= fmodf(y, field().characteristic()); 205 y = fmodf(t, field().characteristic()); 206 } 207 //!@bug should not be neccessary (use assign) 208 return field().init(res, y); 209 } 210 private: 211 Element _bound; 212 size_t _nmax; 213 214 }; 215 } // Namespace LinBox 216 217 #include "linbox/randiter/modular-balanced.h" 218 219 #undef FmodF 220 221 #endif //__LINBOX_modular_balanced_float_H 222 223 // Local Variables: 224 // mode: C++ 225 // tab-width: 4 226 // indent-tabs-mode: nil 227 // c-basic-offset: 4 228 // End: 229 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 230