Show
Ignore:
Timestamp:
08/05/09 14:40:03 (15 years ago)
Author:
mido
Message:

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • library/bdm/itpp_ext.cpp

    r404 r477  
    2020extern "C" {  /* QR factorization of a general matrix A  */ 
    2121        void dgeqrf_ ( int *m, int *n, double *a, int *lda, double *tau, double *work, 
    22                        int *lwork, int *info ); 
     22        int *lwork, int *info ); 
    2323}; 
    2424 
    2525namespace itpp { 
    26         Array<int> to_Arr ( const ivec &indices ) { 
    27                 Array<int> a ( indices.size() ); 
    28                 for ( int i = 0; i < a.size(); i++ ) { 
    29                         a ( i ) = indices ( i ); 
    30                 } 
    31                 return a; 
    32         } 
    33  
    34         ivec linspace ( int from, int to ) { 
    35                 int n=to-from+1; 
    36                 int i; 
    37                 it_assert_debug ( n>0,"wrong linspace" ); 
    38                 ivec iv ( n ); for ( i=0;i<n;i++ ) iv ( i ) =from+i; 
    39                 return iv; 
    40         }; 
    41  
    42         void set_subvector ( vec &ov, const ivec &iv, const vec &v ) 
    43         { 
    44                 it_assert_debug ( ( iv.length() <=v.length() ), 
    45                                                         "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 
    46                                                                         "of range of v" ); 
    47                 for ( int i = 0; i < iv.length(); i++ ) { 
    48                         it_assert_debug ( iv ( i ) <ov.length(), 
    49                                                           "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 
    50                                                                           "of range of v" ); 
    51                         ov ( iv ( i ) ) = v ( i ); 
    52                 } 
    53         } 
    54  
    55         vec get_vec(const vec &v, const ivec &indexlist){ 
    56                 int size = indexlist.size(); 
    57                 vec temp(size); 
    58                 for (int i = 0; i < size; ++i) { 
    59                         temp(i) = v._data()[indexlist(i)]; 
    60                 } 
    61                 return temp; 
    62         } 
     26Array<int> to_Arr ( const ivec &indices ) { 
     27        Array<int> a ( indices.size() ); 
     28        for ( int i = 0; i < a.size(); i++ ) { 
     29                a ( i ) = indices ( i ); 
     30        } 
     31        return a; 
     32} 
     33 
     34ivec linspace ( int from, int to ) { 
     35        int n = to - from + 1; 
     36        int i; 
     37        it_assert_debug ( n > 0, "wrong linspace" ); 
     38        ivec iv ( n ); 
     39        for ( i = 0; i < n; i++ ) iv ( i ) = from + i; 
     40        return iv; 
     41}; 
     42 
     43void set_subvector ( vec &ov, const ivec &iv, const vec &v ) { 
     44        it_assert_debug ( ( iv.length() <= v.length() ), 
     45                          "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 
     46                          "of range of v" ); 
     47        for ( int i = 0; i < iv.length(); i++ ) { 
     48                it_assert_debug ( iv ( i ) < ov.length(), 
     49                                  "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 
     50                                  "of range of v" ); 
     51                ov ( iv ( i ) ) = v ( i ); 
     52        } 
     53} 
     54 
     55vec get_vec ( const vec &v, const ivec &indexlist ) { 
     56        int size = indexlist.size(); 
     57        vec temp ( size ); 
     58        for ( int i = 0; i < size; ++i ) { 
     59                temp ( i ) = v._data() [indexlist ( i ) ]; 
     60        } 
     61        return temp; 
     62} 
    6363 
    6464// Gamma 
    … …  
    6868#define R_FINITE std::isfinite 
    6969 
    70         bvec operator>(const vec &t1, const vec &t2) { 
    71                 it_assert_debug(t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 
    72                 bvec temp(t1.length()); 
    73                 for (int i = 0; i < t1.length(); i++) 
    74                         temp(i) = (t1[i] > t2[i]); 
    75                 return temp; 
    76         } 
    77          
    78         bvec operator<(const vec &t1, const vec &t2) { 
    79                 it_assert_debug(t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 
    80                 bvec temp(t1.length()); 
    81                 for (int i = 0; i < t1.length(); i++) 
    82                         temp(i) = (t1[i] < t2[i]); 
    83                 return temp; 
    84         } 
    85  
    86  
    87         bvec operator& ( const bvec &a, const bvec &b ) { 
    88                 it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" ); 
    89  
    90                 bvec temp ( a.size() ); 
    91                 for ( int i = 0;i < a.size();i++ ) { 
    92                         temp ( i ) = a ( i ) & b ( i ); 
    93                 } 
    94                 return temp; 
    95         } 
    96  
    97         bvec operator| ( const bvec &a, const bvec &b ) { 
    98                 it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" ); 
    99  
    100                 bvec temp ( a.size() ); 
    101                 for ( int i = 0;i < a.size();i++ ) { 
    102                         temp ( i ) = a ( i ) | b ( i ); 
    103                 } 
    104                 return temp; 
    105         } 
     70bvec operator> ( const vec &t1, const vec &t2 ) { 
     71        it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" ); 
     72        bvec temp ( t1.length() ); 
     73        for ( int i = 0; i < t1.length(); i++ ) 
     74                temp ( i ) = ( t1[i] > t2[i] ); 
     75        return temp; 
     76} 
     77 
     78bvec operator< ( const vec &t1, const vec &t2 ) { 
     79        it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" ); 
     80        bvec temp ( t1.length() ); 
     81        for ( int i = 0; i < t1.length(); i++ ) 
     82                temp ( i ) = ( t1[i] < t2[i] ); 
     83        return temp; 
     84} 
     85 
     86 
     87bvec operator& ( const bvec &a, const bvec &b ) { 
     88        it_assert_debug ( b.size() == a.size(), "operator&(): Vectors of different lengths" ); 
     89 
     90        bvec temp ( a.size() ); 
     91        for ( int i = 0; i < a.size(); i++ ) { 
     92                temp ( i ) = a ( i ) & b ( i ); 
     93        } 
     94        return temp; 
     95} 
     96 
     97bvec operator| ( const bvec &a, const bvec &b ) { 
     98        it_assert_debug ( b.size() != a.size(), "operator&(): Vectors of different lengths" ); 
     99 
     100        bvec temp ( a.size() ); 
     101        for ( int i = 0; i < a.size(); i++ ) { 
     102                temp ( i ) = a ( i ) | b ( i ); 
     103        } 
     104        return temp; 
     105} 
    106106 
    107107//#if 0 
    108         Gamma_RNG::Gamma_RNG ( double a, double b ) { 
    109         setup ( a,b ); 
    110 } 
    111         double  Gamma_RNG::sample() { 
    112                 //A copy of rgamma code from the R package!! 
    113                 // 
    114  
    115                 /* Constants : */ 
    116                 const static double sqrt32 = 5.656854; 
    117                 const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ 
    118  
    119                 /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) 
    120                  * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) 
    121                  * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) 
    122                  */ 
    123                 const static double q1 = 0.04166669; 
    124                 const static double q2 = 0.02083148; 
    125                 const static double q3 = 0.00801191; 
    126                 const static double q4 = 0.00144121; 
    127                 const static double q5 = -7.388e-5; 
    128                 const static double q6 = 2.4511e-4; 
    129                 const static double q7 = 2.424e-4; 
    130  
    131                 const static double a1 = 0.3333333; 
    132                 const static double a2 = -0.250003; 
    133                 const static double a3 = 0.2000062; 
    134                 const static double a4 = -0.1662921; 
    135                 const static double a5 = 0.1423657; 
    136                 const static double a6 = -0.1367177; 
    137                 const static double a7 = 0.1233795; 
    138  
    139                 /* State variables [FIXME for threading!] :*/ 
    140                 static double aa = 0.; 
    141                 static double aaa = 0.; 
    142                 static double s, s2, d;    /* no. 1 (step 1) */ 
    143                 static double q0, b, si, c;/* no. 2 (step 4) */ 
    144  
    145                 double e, p, q, r, t, u, v, w, x, ret_val; 
    146                 double a=alpha; 
    147                 double scale=1.0/beta; 
    148  
    149                 if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 ) 
    150                         {it_error ( "Gamma_RNG wrong parameters" );} 
    151  
    152                 if ( a < 1. ) { /* GS algorithm for parameters a < 1 */ 
    153                         if ( a == 0 ) 
    154                                 return 0.; 
    155                         e = 1.0 + exp_m1 * a; 
    156                         for ( ;; ) {  //VS repeat 
    157                                 p = e * unif_rand(); 
    158                                 if ( p >= 1.0 ) { 
    159                                         x = -log ( ( e - p ) / a ); 
    160                                         if ( exp_rand() >= ( 1.0 - a ) * log ( x ) ) 
    161                                                 break; 
    162                                 } 
    163                                 else { 
    164                                         x = exp ( log ( p ) / a ); 
    165                                         if ( exp_rand() >= x ) 
    166                                                 break; 
    167                                 } 
     108Gamma_RNG::Gamma_RNG ( double a, double b ) { 
     109        setup ( a, b ); 
     110} 
     111double  Gamma_RNG::sample() { 
     112        //A copy of rgamma code from the R package!! 
     113        // 
     114 
     115        /* Constants : */ 
     116        const static double sqrt32 = 5.656854; 
     117        const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ 
     118 
     119        /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) 
     120         * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) 
     121         * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) 
     122         */ 
     123        const static double q1 = 0.04166669; 
     124        const static double q2 = 0.02083148; 
     125        const static double q3 = 0.00801191; 
     126        const static double q4 = 0.00144121; 
     127        const static double q5 = -7.388e-5; 
     128        const static double q6 = 2.4511e-4; 
     129        const static double q7 = 2.424e-4; 
     130 
     131        const static double a1 = 0.3333333; 
     132        const static double a2 = -0.250003; 
     133        const static double a3 = 0.2000062; 
     134        const static double a4 = -0.1662921; 
     135        const static double a5 = 0.1423657; 
     136        const static double a6 = -0.1367177; 
     137        const static double a7 = 0.1233795; 
     138 
     139        /* State variables [FIXME for threading!] :*/ 
     140        static double aa = 0.; 
     141        static double aaa = 0.; 
     142        static double s, s2, d;    /* no. 1 (step 1) */ 
     143        static double q0, b, si, c;/* no. 2 (step 4) */ 
     144 
     145        double e, p, q, r, t, u, v, w, x, ret_val; 
     146        double a = alpha; 
     147        double scale = 1.0 / beta; 
     148 
     149        if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 ) { 
     150                it_error ( "Gamma_RNG wrong parameters" ); 
     151        } 
     152 
     153        if ( a < 1. ) { /* GS algorithm for parameters a < 1 */ 
     154                if ( a == 0 ) 
     155                        return 0.; 
     156                e = 1.0 + exp_m1 * a; 
     157                for ( ;; ) {  //VS repeat 
     158                        p = e * unif_rand(); 
     159                        if ( p >= 1.0 ) { 
     160                                x = -log ( ( e - p ) / a ); 
     161                                if ( exp_rand() >= ( 1.0 - a ) * log ( x ) ) 
     162                                        break; 
     163                        } else { 
     164                                x = exp ( log ( p ) / a ); 
     165                                if ( exp_rand() >= x ) 
     166                                        break; 
    168167                        } 
    169                         return scale * x; 
    170                 } 
    171  
    172                 /* --- a >= 1 : GD algorithm --- */ 
    173  
    174                 /* Step 1: Recalculations of s2, s, d if a has changed */ 
    175                 if ( a != aa ) { 
    176                         aa = a; 
    177                         s2 = a - 0.5; 
    178                         s = sqrt ( s2 ); 
    179                         d = sqrt32 - s * 12.0; 
    180                 } 
    181                 /* Step 2: t = standard normal deviate, 
    182                            x = (s,1/2) -normal deviate. */ 
    183  
    184                 /* immediate acceptance (i) */ 
    185                 t = norm_rand(); 
    186                 x = s + 0.5 * t; 
    187                 ret_val = x * x; 
    188                 if ( t >= 0.0 ) 
     168                } 
     169                return scale * x; 
     170        } 
     171 
     172        /* --- a >= 1 : GD algorithm --- */ 
     173 
     174        /* Step 1: Recalculations of s2, s, d if a has changed */ 
     175        if ( a != aa ) { 
     176                aa = a; 
     177                s2 = a - 0.5; 
     178                s = sqrt ( s2 ); 
     179                d = sqrt32 - s * 12.0; 
     180        } 
     181        /* Step 2: t = standard normal deviate, 
     182                   x = (s,1/2) -normal deviate. */ 
     183 
     184        /* immediate acceptance (i) */ 
     185        t = norm_rand(); 
     186        x = s + 0.5 * t; 
     187        ret_val = x * x; 
     188        if ( t >= 0.0 ) 
     189                return scale * ret_val; 
     190 
     191        /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 
     192        u = unif_rand(); 
     193        if ( ( d * u ) <= ( t * t * t ) ) 
     194                return scale * ret_val; 
     195 
     196        /* Step 4: recalculations of q0, b, si, c if necessary */ 
     197 
     198        if ( a != aaa ) { 
     199                aaa = a; 
     200                r = 1.0 / a; 
     201                q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r 
     202                         + q2 ) * r + q1 ) * r; 
     203 
     204                /* Approximation depending on size of parameter a */ 
     205                /* The constants in the expressions for b, si and c */ 
     206                /* were established by numerical experiments */ 
     207 
     208                if ( a <= 3.686 ) { 
     209                        b = 0.463 + s + 0.178 * s2; 
     210                        si = 1.235; 
     211                        c = 0.195 / s - 0.079 + 0.16 * s; 
     212                } else if ( a <= 13.022 ) { 
     213                        b = 1.654 + 0.0076 * s2; 
     214                        si = 1.68 / s + 0.275; 
     215                        c = 0.062 / s + 0.024; 
     216                } else { 
     217                        b = 1.77; 
     218                        si = 0.75; 
     219                        c = 0.1515 / s; 
     220                } 
     221        } 
     222        /* Step 5: no quotient test if x not positive */ 
     223 
     224        if ( x > 0.0 ) { 
     225                /* Step 6: calculation of v and quotient q */ 
     226                v = t / ( s + s ); 
     227                if ( fabs ( v ) <= 0.25 ) 
     228                        q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v 
     229                                                     + a3 ) * v + a2 ) * v + a1 ) * v; 
     230                else 
     231                        q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
     232 
     233 
     234                /* Step 7: quotient acceptance (q) */ 
     235                if ( log ( 1.0 - u ) <= q ) 
    189236                        return scale * ret_val; 
    190  
    191                 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 
     237        } 
     238 
     239        for ( ;; ) { //VS repeat 
     240                /* Step 8: e = standard exponential deviate 
     241                 *      u =  0,1 -uniform deviate 
     242                 *      t = (b,si)-double exponential (laplace) sample */ 
     243                e = exp_rand(); 
    192244                u = unif_rand(); 
    193                 if ( ( d * u ) <= ( t * t * t ) ) 
    194                         return scale * ret_val; 
    195  
    196                 /* Step 4: recalculations of q0, b, si, c if necessary */ 
    197  
    198                 if ( a != aaa ) { 
    199                         aaa = a; 
    200                         r = 1.0 / a; 
    201                         q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r 
    202                                  + q2 ) * r + q1 ) * r; 
    203  
    204                         /* Approximation depending on size of parameter a */ 
    205                         /* The constants in the expressions for b, si and c */ 
    206                         /* were established by numerical experiments */ 
    207  
    208                         if ( a <= 3.686 ) { 
    209                                 b = 0.463 + s + 0.178 * s2; 
    210                                 si = 1.235; 
    211                                 c = 0.195 / s - 0.079 + 0.16 * s; 
    212                         } 
    213                         else if ( a <= 13.022 ) { 
    214                                 b = 1.654 + 0.0076 * s2; 
    215                                 si = 1.68 / s + 0.275; 
    216                                 c = 0.062 / s + 0.024; 
    217                         } 
    218                         else { 
    219                                 b = 1.77; 
    220                                 si = 0.75; 
    221                                 c = 0.1515 / s; 
    222                         } 
    223                 } 
    224                 /* Step 5: no quotient test if x not positive */ 
    225  
    226                 if ( x > 0.0 ) { 
    227                         /* Step 6: calculation of v and quotient q */ 
     245                u = u + u - 1.0; 
     246                if ( u < 0.0 ) 
     247                        t = b - si * e; 
     248                else 
     249                        t = b + si * e; 
     250                /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */ 
     251                if ( t >= -0.71874483771719 ) { 
     252                        /* Step 10:      calculation of v and quotient q */ 
    228253                        v = t / ( s + s ); 
    229254                        if ( fabs ( v ) <= 0.25 ) 
    230                                 q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v 
    231                                                              + a3 ) * v + a2 ) * v + a1 ) * v; 
     255                                q = q0 + 0.5 * t * t * 
     256                                    ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v 
     257                                        + a2 ) * v + a1 ) * v; 
    232258                        else 
    233259                                q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
    234  
    235  
    236                         /* Step 7: quotient acceptance (q) */ 
    237                         if ( log ( 1.0 - u ) <= q ) 
    238                                 return scale * ret_val; 
    239                 } 
    240  
    241                 for ( ;; ) { //VS repeat 
    242                         /* Step 8: e = standard exponential deviate 
    243                          *      u =  0,1 -uniform deviate 
    244                          *      t = (b,si)-double exponential (laplace) sample */ 
    245                         e = exp_rand(); 
    246                         u = unif_rand(); 
    247                         u = u + u - 1.0; 
    248                         if ( u < 0.0 ) 
    249                                 t = b - si * e; 
    250                         else 
    251                                 t = b + si * e; 
    252                         /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */ 
    253                         if ( t >= -0.71874483771719 ) { 
    254                                 /* Step 10:      calculation of v and quotient q */ 
    255                                 v = t / ( s + s ); 
    256                                 if ( fabs ( v ) <= 0.25 ) 
    257                                         q = q0 + 0.5 * t * t * 
    258                                             ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v 
    259                                                 + a2 ) * v + a1 ) * v; 
    260                                 else 
    261                                         q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
    262                                 /* Step 11:      hat acceptance (h) */ 
    263                                 /* (if q not positive go to step 8) */ 
    264                                 if ( q > 0.0 ) { 
    265                                         // TODO: w = expm1(q); 
    266                                         w = exp ( q )-1; 
    267                                         /*  ^^^^^ original code had approximation with rel.err < 2e-7 */ 
    268                                         /* if t is rejected sample again at step 8 */ 
    269                                         if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) ) 
    270                                                 break; 
    271                                 } 
     260                        /* Step 11:      hat acceptance (h) */ 
     261                        /* (if q not positive go to step 8) */ 
     262                        if ( q > 0.0 ) { 
     263                                // TODO: w = expm1(q); 
     264                                w = exp ( q ) - 1; 
     265                                /*  ^^^^^ original code had approximation with rel.err < 2e-7 */ 
     266                                /* if t is rejected sample again at step 8 */ 
     267                                if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) ) 
     268                                        break; 
    272269                        } 
    273                 } /* repeat .. until  `t' is accepted */ 
    274                 x = s + 0.5 * t; 
    275                 return scale * x * x; 
    276         } 
    277  
    278  
    279         bool qr ( const mat &A, mat &R ) { 
    280                 int info; 
    281                 int m = A.rows(); 
    282                 int n = A.cols(); 
    283                 int lwork = n; 
    284                 int k = std::min ( m, n ); 
    285                 vec tau ( k ); 
    286                 vec work ( lwork ); 
    287  
    288                 R = A; 
    289  
    290                 // perform workspace query for optimum lwork value 
    291                 int lwork_tmp = -1; 
    292                 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 
    293                           &info ); 
    294                 if ( info == 0 ) { 
    295                         lwork = static_cast<int> ( work ( 0 ) ); 
    296                         work.set_size ( lwork, false ); 
    297                 } 
    298                 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info ); 
    299  
    300                 // construct R 
    301                 for ( int i=0; i<m; i++ ) 
    302                         for ( int j=0; j<std::min ( i,n ); j++ ) 
    303                                 R ( i,j ) = 0; 
    304  
    305                 return ( info==0 ); 
    306         } 
     270                } 
     271        } /* repeat .. until  `t' is accepted */ 
     272        x = s + 0.5 * t; 
     273        return scale * x * x; 
     274} 
     275 
     276 
     277bool qr ( const mat &A, mat &R ) { 
     278        int info; 
     279        int m = A.rows(); 
     280        int n = A.cols(); 
     281        int lwork = n; 
     282        int k = std::min ( m, n ); 
     283        vec tau ( k ); 
     284        vec work ( lwork ); 
     285 
     286        R = A; 
     287 
     288        // perform workspace query for optimum lwork value 
     289        int lwork_tmp = -1; 
     290        dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 
     291                  &info ); 
     292        if ( info == 0 ) { 
     293                lwork = static_cast<int> ( work ( 0 ) ); 
     294                work.set_size ( lwork, false ); 
     295        } 
     296        dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info ); 
     297 
     298        // construct R 
     299        for ( int i = 0; i < m; i++ ) 
     300                for ( int j = 0; j < std::min ( i, n ); j++ ) 
     301                        R ( i, j ) = 0; 
     302 
     303        return ( info == 0 ); 
     304} 
    307305 
    308306//#endif 
    309         std::string num2str(double d){ 
    310                 char tmp[20];//that should do 
    311                 sprintf(tmp,"%f",d); 
    312                 return std::string(tmp); 
    313         }; 
    314         std::string num2str(int i){ 
    315                 char tmp[10];//that should do 
    316                 sprintf(tmp,"%d",i); 
    317                 return std::string(tmp); 
    318         }; 
     307std::string num2str ( double d ) { 
     308        char tmp[20];//that should do 
     309        sprintf ( tmp, "%f", d ); 
     310        return std::string ( tmp ); 
     311}; 
     312std::string num2str ( int i ) { 
     313        char tmp[10];//that should do 
     314        sprintf ( tmp, "%d", i ); 
     315        return std::string ( tmp ); 
     316}; 
    319317 
    320318// digamma 
    … …  
    323321#define el 0.5772156649015329 
    324322 
    325 double psi(double x) { 
    326     double s,ps,xa,x2; 
    327     int n,k; 
    328     static double a[] = { 
    329         -0.8333333333333e-01, 
    330          0.83333333333333333e-02, 
    331         -0.39682539682539683e-02, 
    332          0.41666666666666667e-02, 
    333         -0.75757575757575758e-02, 
    334          0.21092796092796093e-01, 
    335         -0.83333333333333333e-01, 
    336          0.4432598039215686}; 
    337  
    338     xa = fabs(x); 
    339     s = 0.0; 
    340     if ((x == (int)x) && (x <= 0.0)) { 
    341         ps = 1e308; 
    342         return ps; 
    343     } 
    344     if (xa == (int)xa) { 
    345         n = xa; 
    346         for (k=1;k<n;k++) { 
    347             s += 1.0/k; 
    348         } 
    349         ps =  s-el; 
    350     } 
    351     else if ((xa+0.5) == ((int)(xa+0.5))) { 
    352         n = xa-0.5; 
    353         for (k=1;k<=n;k++) { 
    354             s += 1.0/(2.0*k-1.0); 
    355         } 
    356         ps = 2.0*s-el-1.386294361119891; 
    357     } 
    358     else { 
    359         if (xa < 10.0) { 
    360             n = 10-(int)xa; 
    361             for (k=0;k<n;k++) { 
    362                 s += 1.0/(xa+k); 
    363             } 
    364             xa += n; 
    365         } 
    366         x2 = 1.0/(xa*xa); 
    367         ps = log(xa)-0.5/xa+x2*(((((((a[7]*x2+a[6])*x2+a[5])*x2+ 
    368             a[4])*x2+a[3])*x2+a[2])*x2+a[1])*x2+a[0]); 
    369         ps -= s; 
    370     } 
    371     if (x < 0.0) 
    372         ps = ps - M_PI*std::cos(M_PI*x)/std::sin(M_PI*x)-1.0/x; 
    373         return ps; 
    374 } 
    375  
    376 } 
     323double psi ( double x ) { 
     324        double s, ps, xa, x2; 
     325        int n, k; 
     326        static double a[] = { 
     327                -0.8333333333333e-01, 
     328                0.83333333333333333e-02, 
     329                -0.39682539682539683e-02, 
     330                0.41666666666666667e-02, 
     331                -0.75757575757575758e-02, 
     332                0.21092796092796093e-01, 
     333                -0.83333333333333333e-01, 
     334                0.4432598039215686 
     335        }; 
     336 
     337        xa = fabs ( x ); 
     338        s = 0.0; 
     339        if ( ( x == ( int ) x ) && ( x <= 0.0 ) ) { 
     340                ps = 1e308; 
     341                return ps; 
     342        } 
     343        if ( xa == ( int ) xa ) { 
     344                n = xa; 
     345                for ( k = 1; k < n; k++ ) { 
     346                        s += 1.0 / k; 
     347                } 
     348                ps =  s - el; 
     349        } else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) ) ) { 
     350                n = xa - 0.5; 
     351                for ( k = 1; k <= n; k++ ) { 
     352                        s += 1.0 / ( 2.0 * k - 1.0 ); 
     353                } 
     354                ps = 2.0 * s - el - 1.386294361119891; 
     355        } else { 
     356                if ( xa < 10.0 ) { 
     357                        n = 10 - ( int ) xa; 
     358                        for ( k = 0; k < n; k++ ) { 
     359                                s += 1.0 / ( xa + k ); 
     360                        } 
     361                        xa += n; 
     362                } 
     363                x2 = 1.0 / ( xa * xa ); 
     364                ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5] ) * x2 + 
     365                                                        a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0] ); 
     366                ps -= s; 
     367        } 
     368        if ( x < 0.0 ) 
     369                ps = ps - M_PI * std::cos ( M_PI * x ) / std::sin ( M_PI * x ) - 1.0 / x; 
     370        return ps; 
     371} 
     372 
     373}