Show
Ignore:
Timestamp:
06/06/08 19:09:43 (6 months ago)
Author:
pernet
Message:

* Upgrade to fflas-ffpack v 1.3.3
* rename all the modular-balance into modular-balanced
* fix the few remaining compilation warning with gcc-4.3

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • trunk/linbox/linbox/fflas/fflas_bounds.inl

    r2954 r2962  
    22 
    33/* fflas/fflas_bounds.inl 
    4  * Copyright (C) 2007 Clement Pernet 
     4 * Copyright (C) 2008 Clement Pernet 
    55 * 
    66 * Written by Clement Pernet <Clement.Pernet@imag.fr> 
    … …  
    1515#endif 
    1616 
    17 /* From M Abshoff */ 
    18 #if defined(__sun)  
    19 #define lround(x) my_lround(x) 
    20 static long my_lround(double x) 
    21 { 
    22        return (long) ((x) >= 0 ? (x) + 0.5 : (x) - 0.5); 
    23 } 
    24 #endif 
    25  
    26 //--------------------------------------------------------------------- 
    27 // DotProdBound : 
    28 // Computes the maximal dimension k so that the product A*B+beta.C over Z, 
    29 // where A is m*k and B is k*n can be performed correctly with w Winograd 
    30 // recursion levels on the number of bits of the floating point mantissa 
    31 //--------------------------------------------------------------------- 
    32 template  <class Field>  
    33 inline size_t FFLAS::DotProdBoundCompute (const Field& F, const size_t w,  
    34                                           const typename Field::Element& beta, 
    35                                           const FFLAS_BASE base){ 
     17/** 
     18 * MatMulParameters 
     19 * 
     20 * \brief Computes the threshold parameters for the cascade 
     21 *        Matmul algorithm 
     22 * 
     23 *  
     24 * \param F Finite Field/Ring of the computation. 
     25 * \param k Common dimension of A and B, in the product A x B 
     26 * \param bet Computing AB + beta C 
     27 * \param delayedDim Returns the size of blocks that can be multiplied 
     28 *                   over Z with no overflow 
     29 * \param base Returns the type of BLAS representation to use 
     30 * \param winoRecLevel Returns the number of recursion levels of 
     31 *                     Strassen-Winograd's algorithm to perform 
     32 * \param winoLevelProvided tells whether the user forced the number of 
     33 *                          recursive level of Winograd's algorithm 
     34 */ 
     35template <class Field> 
     36inline void FFLAS::MatMulParameters (const Field& F, 
     37                                     const size_t k, 
     38                                     const typename Field::Element& beta, 
     39                                     size_t& delayedDim, 
     40                                     FFLAS_BASE& base, 
     41                                     size_t& winoRecLevel, 
     42                                     bool winoLevelProvided) { 
     43 
     44        // Strategy : determine Winograd's recursion first, then choose appropriate 
     45        // floating point representation, and finally the blocking dimension. 
     46        // Can be improved for some cases. 
     47 
     48        if (!winoLevelProvided) 
     49                winoRecLevel = WinoSteps (k); 
     50        base = BaseCompute (F, winoRecLevel); 
     51        delayedDim = DotProdBound (F, winoRecLevel, beta, base); 
     52 
     53        size_t n = k; 
     54        size_t winoDel = winoRecLevel; 
     55 
     56        // Computes the delayedDim, only depending on the recursive levels 
     57        // that must be performed over Z 
     58        while (winoDel > 0 && delayedDim < n) { 
     59                winoDel--; 
     60                delayedDim = DotProdBound (F, winoDel, beta, base); 
     61                n >>= 1; 
     62        } 
     63        delayedDim = MIN (n, delayedDim); 
     64} 
     65 
     66/** 
     67 * DotProdBound 
     68 * 
     69 * \brief  computes the maximal size for delaying the modular reduction 
     70 *         in a dotproduct 
     71 * 
     72 * This is the default version assuming a conversion to a positive modular representation 
     73 *  
     74 * \param F Finite Field/Ring of the computation 
     75 * \param winoRecLevel Number of recusrive Strassen-Winograd levels (if any, 0 otherwise) 
     76 * \param beta Computing AB + beta C 
     77 * \param base Type of floating point representation for delayed modular computations 
     78 *  
     79 */ 
     80template <class Field> 
     81inline size_t FFLAS::DotProdBound (const Field& F, 
     82                                   const size_t w,  
     83                                   const typename Field::Element& beta, 
     84                                   const FFLAS_BASE base) { 
    3685         
     86        FFLAS_INT_TYPE p; 
     87        F.characteristic(p); 
    3788        typename Field::Element mone; 
    38         static FFLAS_INT_TYPE p; 
     89        F.init (mone, -1.0); 
     90 
     91        unsigned long mantissa = 
     92                (base == FflasDouble) ? DOUBLE_MANTISSA : FLOAT_MANTISSA; 
     93 
     94        if (p == 0) 
     95                return 1; 
     96         
     97        double kmax; 
     98        if (w > 0) { 
     99                double c = computeFactor (F,w); 
     100                double d = (double (1ULL << mantissa) /(c*c) + 1); 
     101                if (d < 2) 
     102                        return 1; 
     103                kmax = floor (d * (1ULL << w)); 
     104        } else { 
     105                ////// A fixer: (p-1)/2 si balanced 
     106 
     107                double c = p-1; 
     108                double cplt=0; 
     109                if (!F.isZero (beta)){ 
     110                        if (F.isOne (beta) || F.areEqual (beta, mone)) cplt = c; 
     111                        else cplt = c*c; 
     112                } 
     113                kmax = floor ( (double ((1ULL << mantissa) - cplt)) / (c*c)); 
     114                if (kmax  <= 1) 
     115                        return 1; 
     116                } 
     117        //kmax--; // we computed a strict upper bound 
     118        return  (size_t) MIN (kmax, 1ULL << 31); 
     119} 
     120 
     121/** 
     122 * Internal function for the bound computation 
     123 * Generic implementation for positive representations 
     124 */ 
     125template <class Field> 
     126inline double FFLAS::computeFactor (const Field& F, const size_t w){ 
     127        FFLAS_INT_TYPE p; 
    39128        F.characteristic(p); 
    40         F.init (mone, -1.0); 
    41         size_t kmax; 
    42  
    43         unsigned long mantissa = (base == FflasDouble)? DOUBLE_MANTISSA : FLOAT_MANTISSA; 
     129        size_t ex=1; 
     130        for (size_t i=0; i < w; ++i)    ex *= 3; 
     131        return double(p - 1) * (1 + ex) / 2; 
     132} 
     133 
     134/** 
     135 * WinoSteps 
     136 * 
     137 * \brief Computes the number of recursive levels to perform 
     138 * 
     139 * \param m the common dimension in the product AxB 
     140 *  
     141 */ 
     142inline size_t FFLAS::WinoSteps (const size_t m) { 
     143        size_t w = 0; 
     144        size_t mt = m; 
     145        while (mt >= WINOTHRESHOLD) {w++; mt >>= 1;} 
     146        return w; 
     147} 
     148 
     149/** 
     150 * BaseCompute 
     151 * 
     152 * \brief Determines the type of floating point representation to convert to, 
     153 *        for BLAS computations 
     154 * \param F Finite Field/Ring of the computation 
     155 * \param w Number of recursive levels in Winograd's algorithm 
     156 *  
     157 */ 
     158template <class Field> 
     159inline FFLAS::FFLAS_BASE FFLAS::BaseCompute (const Field& F, const size_t w){ 
    44160         
    45         if (p == 0) 
    46                 kmax = 2; 
    47         else 
    48                 if (w > 0) { 
    49                         size_t ex=1; 
    50                         for (size_t i=0; i < w; ++i)    ex *= 3; 
    51                         //FFLAS_INT_TYPE c = (p-1)*(ex)/2; //bound for a centered representation 
    52                         FFLAS_INT_TYPE c; 
    53 #ifndef _LINBOX_LINBOX_CONFIG_H 
    54                         if (F.balanced) 
    55                                 c = (p-1)*(ex)/2; // balanced representation 
    56                         else 
    57 #endif 
    58                                 c = (p-1)*(1+ex)/2; 
    59                         kmax =  lround(( double(1ULL << mantissa) /double(c*c) + 1)*(1ULL << w)); 
    60                         if (kmax ==  ( 1ULL << w)) 
    61                                 kmax = 2; 
    62                 } 
    63                 else{ 
    64                         FFLAS_INT_TYPE c = p-1; 
    65                         FFLAS_INT_TYPE cplt=0; 
    66                         if (!F.isZero (beta)) { 
    67                                 if (F.isOne (beta) || F.areEqual (beta, mone)) 
    68                                         cplt = c; 
    69                                 else cplt = c*c; 
    70                         } 
    71                         kmax =  lround(( double((1ULL << mantissa) - cplt)) /double(c*c)); 
    72                         if (kmax  < 2) 
    73                                 kmax = 2; 
    74                 } 
    75         //cerr<<"kmax = "<<kmax<<endl; 
    76         return  (size_t) MIN(kmax,1ULL<<31); 
    77 } 
    78  
    79  
    80  
    81 template  < class Field >  
    82 inline size_t FFLAS::DotProdBound (const Field& F, const size_t w,  
    83                                    const typename Field::Element& beta, 
    84                                    const FFLAS_BASE base) 
    85 { 
    86         static Field G = F; 
    87         static FFLAS_INT_TYPE pig; 
    88         G.characteristic(pig); 
    89         FFLAS_INT_TYPE pif; 
    90         F.characteristic(pif); 
    91         static typename Field::Element b = beta; 
    92         static size_t w2 = w; 
    93         static FFLAS_BASE base2 = base; 
    94         static size_t kmax = DotProdBoundCompute (F, w, beta, base); 
    95         if ((b != beta) || (pif != pig) ||  (w2 != w) || (base2!=base)) { 
    96                 G = F; 
    97                 w2 = w; 
    98                 b = beta; 
    99                 base2 = base; 
    100                 kmax =  DotProdBoundCompute (F, w, beta, base); 
    101         }        
    102         return kmax; 
    103 } 
    104  
    105 //--------------------------------------------------------------------- 
    106 // TRSMBound 
    107 // Computes nmax s.t. (p-1)/2*(p^{nmax-1} + (p-2)^{nmax-1}) < 2^53 
    108 //--------------------------------------------------------------------- 
    109 inline size_t bound_compute_double(const FFLAS_INT_TYPE pi) { 
    110          
    111         FFLAS_INT_TYPE p=pi,p1=1,p2=1; 
    112         size_t nmax = 1; 
    113         FFLAS_INT_TYPE max = ( FFLAS_INT_TYPE(  1ULL<<(DOUBLE_MANTISSA+1) )/(p-1)); 
     161        FFLAS_INT_TYPE pi; 
     162        F.characteristic(pi); 
     163        FFLAS_BASE base; 
     164        switch (w) { 
     165        case 0: base = (pi < FLOAT_DOUBLE_THRESHOLD_0)? FflasFloat : FflasDouble; 
     166                break; 
     167        case 1:  base = (pi < FLOAT_DOUBLE_THRESHOLD_1)? FflasFloat : FflasDouble; 
     168                break; 
     169        case 2:  base = (pi < FLOAT_DOUBLE_THRESHOLD_2)? FflasFloat : FflasDouble; 
     170                break; 
     171        default: base = FflasDouble; 
     172                break; 
     173        } 
     174        return base; 
     175} 
     176 
     177 
     178/************************************************************************************* 
     179 * Specializations for ModularPositive and ModularBalanced over double and float 
     180 *************************************************************************************/ 
     181 
     182template <class Element> 
     183inline double computeFactor (const ModularBalanced<Element>& F, const size_t w){ 
     184        FFLAS_INT_TYPE p; 
     185        F.characteristic(p); 
     186        size_t ex=1; 
     187        for (size_t i=0; i < w; ++i)    ex *= 3; 
     188        return  (p - 1) * ex / 2;  
     189} 
     190 
     191template <> 
     192inline FFLAS::FFLAS_BASE FFLAS::BaseCompute (const Modular<double>& F, 
     193                                             const size_t w){ 
     194        return FflasDouble; 
     195} 
     196 
     197template <> 
     198inline FFLAS::FFLAS_BASE FFLAS::BaseCompute (const Modular<float>& F, 
     199                                             const size_t w){ 
     200        return FflasFloat; 
     201} 
     202 
     203template <> 
     204inline FFLAS::FFLAS_BASE FFLAS::BaseCompute (const ModularBalanced<double>& F, 
     205                                             const size_t w){ 
     206        return FflasDouble; 
     207} 
     208 
     209template <> 
     210inline FFLAS::FFLAS_BASE FFLAS::BaseCompute (const ModularBalanced<float>& F, 
     211                                             const size_t w){ 
     212        return FflasFloat; 
     213} 
     214 
     215/** 
     216 * TRSMBound 
     217 * 
     218 * \brief  computes the maximal size for delaying the modular reduction 
     219 *         in a triangular system resolution 
     220 * 
     221 * This is the default version over an arbitrary field. 
     222 * It is currently never used (the recursive algorithm is run until n=1 in this case) 
     223 *  
     224 * \param F Finite Field/Ring of the computation 
     225 *  
     226 */ 
     227template <class Field> 
     228inline size_t FFLAS::TRSMBound (const Field& F) { 
     229        return 1;        
     230} 
     231 
     232/** 
     233 * Specialization for positive modular representation over double 
     234 * Computes nmax s.t. (p-1)/2*(p^{nmax-1} + (p-2)^{nmax-1}) < 2^53 
     235 * See [Dumas Giorgi Pernet 06, arXiv:cs/0601133] 
     236 */ 
     237template<> 
     238inline size_t FFLAS::TRSMBound (const Modular<double>& F){ 
     239 
     240        FFLAS_INT_TYPE pi; 
     241        F.characteristic(pi); 
     242        unsigned long long p = pi, p1 = 1, p2 = 1; 
     243        size_t nmax = 0; 
     244        unsigned long long max = ( (1ULL << (DOUBLE_MANTISSA + 1) ) / (p - 1)); 
    114245        while ( (p1 + p2) < max ){ 
    115246                p1*=p; 
    … …  
    117248                nmax++; 
    118249        } 
    119         nmax--; 
    120250        return nmax; 
    121251} 
    122 inline size_t bound_compute_double_balanced(const FFLAS_INT_TYPE pi) { 
    123          
    124         FFLAS_INT_TYPE p=(pi+1)/2,p1=1; 
    125         size_t nmax=0; 
    126         FFLAS_INT_TYPE max = ( FFLAS_INT_TYPE(  1ULL<<(DOUBLE_MANTISSA))/(p-1)); 
    127         while ( (p1) < max ){ 
    128                 p1*=p; 
    129                 nmax++; 
    130         } 
    131         return nmax; 
    132 } 
    133 inline size_t bound_compute_float(const FFLAS_INT_TYPE pi) { 
    134          
    135         FFLAS_INT_TYPE p=pi,p1=1,p2=1; 
    136         size_t nmax=1; 
    137         FFLAS_INT_TYPE max = (FFLAS_INT_TYPE(   1ULL<<(FLOAT_MANTISSA+1) )/(p-1)); 
     252 
     253 
     254/** 
     255 * Specialization for positive modular representation over float 
     256 * Computes nmax s.t. (p-1)/2*(p^{nmax-1} + (p-2)^{nmax-1}) < 2^24 
     257 * See [Dumas Giorgi Pernet 06, arXiv:cs/0601133] 
     258 */ 
     259template<> 
     260inline size_t FFLAS::TRSMBound (const Modular<float>& F){ 
     261 
     262        FFLAS_INT_TYPE pi; 
     263        F.characteristic(pi); 
     264        unsigned long long p = pi, p1 = 1, p2 = 1; 
     265        size_t nmax = 0; 
     266        unsigned long long max = ( (1ULL << (FLOAT_MANTISSA + 1) ) / (p - 1)); 
    138267        while ( (p1 + p2) < max ){ 
    139268                p1*=p; 
    … …  
    141270                nmax++; 
    142271        } 
    143         nmax--; 
    144272        return nmax; 
    145273} 
    146 inline size_t bound_compute_float_balanced(const FFLAS_INT_TYPE pi) { 
    147          
    148         FFLAS_INT_TYPE p=(pi+1)/2,p1=1; 
    149         size_t nmax=0; 
    150         FFLAS_INT_TYPE max = ( FFLAS_INT_TYPE(  1ULL<<(FLOAT_MANTISSA))/(p-1)); 
    151         while ( (p1) < max ){ 
    152                 p1*=p; 
     274 
     275/** 
     276 * Specialization for balanced modular representation over double 
     277 * Computes nmax s.t. (p-1)/2*(((p+1)/2)^{nmax-1}) < 2^53 
     278 * See [Dumas Giorgi Pernet 06, arXiv:cs/0601133] 
     279 */ 
     280template<> 
     281inline size_t FFLAS::TRSMBound (const ModularBalanced<double>& F){ 
     282 
     283        FFLAS_INT_TYPE pi; 
     284        F.characteristic (pi); 
     285        unsigned long long p = (pi + 1) / 2, p1 = 1; 
     286        size_t nmax = 0; 
     287        unsigned long long max = ((1ULL << (DOUBLE_MANTISSA + 1)) / ((unsigned long long)(p - 1))); 
     288        while (p1 < max){ 
     289                p1 *= p; 
    153290                nmax++; 
    154291        } 
    … …  
    156293} 
    157294 
    158 template <class Field> 
    159 inline size_t 
    160 FFLAS::TRSMBound (const Field& F) { 
    161         return callTRSMBound<typename Field::Element> () (F); 
    162 } 
    163  
    164 template<class Element> 
    165 class FFLAS::callTRSMBound { 
    166 public: 
    167         template <class Field> 
    168         size_t operator () (const Field& F) { 
    169                 FFLAS_INT_TYPE pi; 
    170                 F.characteristic(pi); 
    171                 static long unsigned int p=pi; 
    172                 static size_t nmax=bound_compute_double(p); 
    173                 if (p == pi)  
    174                         return nmax; 
    175                 else  
    176                         return nmax=bound_compute_double(p=pi); 
    177         } 
    178 }; 
    179  
     295/** 
     296 * Specialization for balanced modular representation over float 
     297 * Computes nmax s.t. (p-1)/2*(((p+1)/2)^{nmax-1}) < 2^24 
     298 * See [Dumas Giorgi Pernet 06, arXiv:cs/0601133] 
     299 */ 
    180300template<> 
    181 class FFLAS::callTRSMBound<double> { 
    182 public: 
    183         template <class Field> 
    184         size_t operator () (const Field& F) { 
    185                 FFLAS_INT_TYPE pi; 
    186                 F.characteristic(pi); 
    187                 static FFLAS_INT_TYPE p=pi; 
    188 #ifdef _LINBOX_LINBOX_CONFIG_H 
    189                 static size_t nmax = bound_compute_double(pi); 
    190 #else 
    191                 static size_t nmax = (F.balanced) ? bound_compute_double_balanced(pi) : bound_compute_double(pi); 
    192 #endif                         
    193                 if (p == pi)  
    194                         return nmax; 
    195                 else 
    196 #ifdef _LINBOX_LINBOX_CONFIG_H 
    197                         return nmax= bound_compute_double (p=pi); //(F.balanced) ? bound_compute_balanced(p=pi) : bound_compute(p=pi); 
    198 #else 
    199                 return nmax = (F.balanced) ? bound_compute_double_balanced(p=pi) : bound_compute_double(p=pi); 
    200 #endif 
    201         } 
    202 }; 
    203  
    204 template<> 
    205 class FFLAS::callTRSMBound<float> { 
    206 public: 
    207         template <class Field> 
    208         size_t operator () (const Field& F) { 
    209                 FFLAS_INT_TYPE pi; 
    210                 F.characteristic(pi); 
    211                 static FFLAS_INT_TYPE p=pi; 
    212 #ifdef _LINBOX_LINBOX_CONFIG_H 
    213                 static size_t nmax = bound_compute_float(pi); 
    214 #else 
    215                 static size_t nmax = (F.balanced) ? bound_compute_float_balanced(pi) : bound_compute_float(pi); 
    216 #endif                         
    217                 if (p == pi)  
    218                         return nmax; 
    219                 else 
    220 #ifdef _LINBOX_LINBOX_CONFIG_H 
    221                         return nmax= bound_compute_float (p=pi); //(F.balanced) ? bound_compute_balanced(p=pi) : bound_compute(p=pi); 
    222 #else 
    223                 return nmax = (F.balanced) ? bound_compute_float_balanced(p=pi) : bound_compute_float(p=pi); 
    224 #endif 
    225         } 
    226 }; 
    227  
    228  
    229  
    230 /** Automatic computation of the parameters for Matrix Multiplication 
    231  * 
    232  */ 
    233  
    234 // To be tuned up: use statics to perform the computations only once for a given set of parameters 
    235 template <class Element> 
    236 class FFLAS::setMatMulParam { 
    237  public: 
    238         template <class Field> 
    239                 void operator() (const Field& F, const size_t m, 
    240                                  const typename Field::Element& beta, 
    241                                  size_t& w, FFLAS_BASE& base, size_t& kmax){ 
    242  
    243                 // Strategy : determine Winograd's recursion first, then choose appropriate 
    244                 // floating point representation, and finally the blocking dimension. 
    245                 // Can be improved for some cases. 
    246                 static size_t ms = m; 
    247                 static size_t ks = WinoSteps (m); 
    248  
    249                 if (m != ms) 
    250                         ks = WinoSteps (ms = m); 
    251  
    252                 w = ks; 
    253                  
    254                 base = BaseCompute<typename Field::Element> ()(F, w); 
    255                 kmax = DotProdBound (F, w, beta, base); 
    256         } 
    257 }; 
    258  
    259 template<> 
    260 class FFLAS::setMatMulParam<double> { 
    261 public: 
    262         template <class Field> 
    263         void operator () (const Field& F, const size_t m, 
    264                           const typename Field::Element& beta, 
    265                           size_t& w, FFLAS_BASE& base, size_t& kmax){ 
    266  
    267                 static size_t ms = m; 
    268                 static size_t ks = WinoSteps (m); 
    269  
    270                 if (m != ms) 
    271                         ks = WinoSteps (ms = m); 
    272  
    273                 w = ks; 
    274                 base = FflasDouble; 
    275                 kmax = DotProdBound (F, w, beta, base); 
    276         } 
    277 }; 
    278  
    279 template<> 
    280 class FFLAS::setMatMulParam<float> { 
    281 public: 
    282         template <class Field> 
    283         void operator () (const Field& F, const size_t m, 
    284                           const typename Field::Element& beta, 
    285                           size_t& w, FFLAS_BASE& base, size_t& kmax){ 
    286  
    287                 static size_t ms = m; 
    288                 static size_t ks = WinoSteps (m); 
    289  
    290                 if (m != ms) 
    291                         ks = WinoSteps (ms = m); 
    292  
    293                 w = ks; 
    294                 base = FflasFloat; 
    295                 kmax = DotProdBound (F, w, beta, base); 
    296         } 
    297 }; 
    298  
    299 inline size_t FFLAS::WinoSteps (const size_t m) { 
    300  
    301         size_t w=0; 
    302         size_t mt = m; 
    303         while (mt >= WINOTHRESHOLD){ 
    304                 w++; 
    305                 mt/=2; 
    306         } 
    307         //cerr<<"Winostep = "<<w<<endl; 
    308         return w; 
    309 } 
    310  
    311  
    312 template <class Element> 
    313 class FFLAS::BaseCompute { 
    314 public: 
    315         template <class Field> 
    316         FFLAS::FFLAS_BASE operator() (const Field& F, const size_t w){ 
    317                  
    318                 FFLAS_INT_TYPE pi; 
    319                 F.characteristic(pi); 
    320                 FFLAS_BASE base; 
    321                 switch (w) { 
    322                 case 0: base = (pi < FLOAT_DOUBLE_THRESHOLD_0)? FflasFloat : FflasDouble; 
    323                         break; 
    324                 case 1:  base = (pi < FLOAT_DOUBLE_THRESHOLD_1)? FflasFloat : FflasDouble; 
    325                         break; 
    326                 case 2:  base = (pi < FLOAT_DOUBLE_THRESHOLD_2)? FflasFloat : FflasDouble; 
    327                         break; 
    328                 default: base = FflasDouble; 
    329                         break; 
    330                 } 
    331                 //cerr<<"base = "<<base<<endl; 
    332                 return base; 
    333                 //      case 0: return (pi < FLOAT_DOUBLE_THRESHOLD_0)? FflasFloat : FflasDouble; 
    334                 //      case 1: return (pi < FLOAT_DOUBLE_THRESHOLD_1)? FflasFloat : FflasDouble; 
    335                 //      case 2: return (pi < FLOAT_DOUBLE_THRESHOLD_2)? FflasFloat : FflasDouble; 
    336                 //      default: return FflasDouble; 
    337         } 
    338 }; 
    339  
    340 template<> 
    341 class FFLAS::BaseCompute<double> { 
    342 public: 
    343         template <class Field> 
    344         FFLAS::FFLAS_BASE operator() (const Field& F, const size_t w){ 
    345                 return FflasDouble; 
    346         } 
    347 }; 
    348  
    349 template<> 
    350 class FFLAS::BaseCompute<float> { 
    351 public: 
    352         template <class Field> 
    353         FFLAS::FFLAS_BASE operator() (const Field& F, const size_t w){ 
    354                 return FflasFloat; 
    355         } 
    356 }; 
     301inline size_t FFLAS::TRSMBound (const ModularBalanced<float>& F){ 
     302 
     303        FFLAS_INT_TYPE pi; 
     304        F.characteristic (pi); 
     305        unsigned long long p = (pi + 1) / 2, p1 = 1; 
     306        size_t nmax = 0; 
     307        unsigned long long max = ((1ULL << (FLOAT_MANTISSA + 1)) / ((unsigned long long) (pi - 1))); 
     308        while (p1 < max){ 
     309                p1 *= p; 
     310                nmax++; 
     311        } 
     312        return nmax; 
     313 
     314} 
     315