Functional Programming in C++14

lost_templar

這一陣子折騰 C++14，練手的時(shí)候把自己之前寫過(guò)的一個(gè) 生成非均勻隨機(jī)變量的庫(kù)重構(gòu)了一下，還沒有搞完，記錄了重構(gòu)的大概過(guò)程，先發(fā)過(guò)來(lái)大家瞧瞧。

https://github.com/fengwang/vg

標(biāo)準(zhǔn)組件 (Stardard Components)

Seed

寫一個(gè)函數(shù)，這個(gè)函數(shù)返回一個(gè)隨機(jī)整數(shù)。

Write a function that returns a random integer.

1. 
   
2. inline auto make_seed() noexcept
   
3. {
   
4.     typedef std::uint_fast64_t        seed_type;
   
5.     auto const now                  = std::chrono::system_clock::now();
   
6.     auto const duration             = now.time_since_epoch();
   
7.     seed_type  const nano_seed      = std::chrono::duration_cast<std::chrono::nanoseconds>(duration).count();
   
8.     seed_type const random_address  = reinterpret_cast<seed_type>( &nano_seed );
   
9.     seed_type ans                   = nano_seed + ( random_address | (random_address << 32) );
   
10.     ans                             = ( ans & 0x5555555555555555ULL ) <<  1 | ( ans & 0xAAAAAAAAAAAAAAAAULL ) >>  1;
    
11.     ans                             = ( ans & 0x3333333333333333ULL ) <<  2 | ( ans & 0xCCCCCCCCCCCCCCCCULL ) >>  2;
    
12.     ans                             = ( ans & 0x0F0F0F0F0F0F0F0FULL ) <<  4 | ( ans & 0xF0F0F0F0F0F0F0F0ULL ) >>  4;
    
13.     ans                             = ( ans & 0x00FF00FF00FF00FFULL ) <<  8 | ( ans & 0xFF00FF00FF00FF00ULL ) >>  8;
    
14.     ans                             = ( ans & 0x0000FFFF0000FFFFULL ) << 16 | ( ans & 0xFFFF0000FFFF0000ULL ) >> 16;
    
15.     ans                             = ( ans & 0x00000000FFFFFFFFULL ) << 32 | ( ans & 0xFFFFFFFF00000000ULL ) >> 32;
    
16.     ans                            ^= nano_seed;
    
17.     return                            ans;
    
18. }
    

復(fù)制代碼

引擎 -- Engine

寫一個(gè)函數(shù)，這個(gè)函數(shù)返回另外一個(gè)函數(shù)，返回的函數(shù)被調(diào)用后返回一個(gè)在 [0,1] 區(qū)間服從均勻分布的隨機(jī)數(shù)。

Write a function that return a function that returns a random variate that uniformly distributed in range [0,1].

1. 
   
2. inline auto linear_congruential() noexcept
   
3. {
   
4.     static const std::uint_fast64_t a_ = 6364136223846793005ULL;
   
5.     static const std::uint_fast64_t c_ = 1442695040888963407ULL;
   
6.     std::uint_fast64_t              x_ = make_seed();
   
7. 
   
8.     return [=]() mutable noexcept
   
9.     {
   
10.         x_ *= a_;
    
11.         x_ += c_;
    
12.         return static_cast<long double>( x_ ) / static_cast<long double>( std::numeric_limits<std::uint_fast64_t>::max() );
    
13.     };
    
14. }
    

復(fù)制代碼

Reference

    D. E. Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorithms, Third Edition. Addison-Wesley, 1997. ISBN 0-201-89684-2. Section 3.2.1: The Linear Congruential Method, pp. 10–26.

分布 -- Distribution
正態(tài)分布 -- Normal

寫一個(gè)函數(shù)，這個(gè)函數(shù)返回一個(gè)函數(shù) g, g 返回一個(gè)函數(shù) h，h 接收一個(gè)函數(shù) e 作為參數(shù)，h 返回一個(gè)服從 N(0,1) 分布的隨機(jī)數(shù)，而作為函數(shù) h 參數(shù)的函數(shù) e 被調(diào)用后會(huì)返回一個(gè)在區(qū)間 [0,1] 上服從均勻分布的隨機(jī)數(shù)。

Write a function that returns a function that takes a function, which returns a random variate uniform in range [0,1], and returns a variate that is subject to N(0,1).

1. 
   
2. inline auto normal() noexcept
   
3. {
   
4.     return []() noexcept
   
5.     {
   
6.         return []( auto& engine ) noexcept
   
7.         {
   
8.             long double const one = 1;
   
9. 
   
10.             for (;;)
    
11.             {
    
12.                 auto const u1 = engine();
    
13.                 auto const u2 = engine();
    
14.                 auto const v1 = u1 + u1 - one;
    
15.                 auto const v2 = u2 + u2 - one;
    
16.                 auto const s = v1 * v1 + v2 * v2;
    
17. 
    
18.                 if ( s < one )
    
19.                 {
    
20.                     auto const tmp = -std::log( s );
    
21.                     return v1 * std::sqrt(( tmp + tmp ) / s );
    
22.                 }
    
23.             }
    
24.         };
    
25.     };
    
26. }
    

復(fù)制代碼

Reference

    G. E. P. Box and Mervin E. Muller, A Note on the Generation of Random Normal Deviates, The Annals of Mathematical Statistics (195, Vol. 29, No. 2 pp. 610-611

均勻分布 -- Uniform

仿照正態(tài)分布的生成函數(shù)，寫一個(gè)在區(qū)間 [lower, upper] 上服從均勻分布的函數(shù)。

Write a uniform ditribution function defined in range [lower, upper].

1. 
   
2. auto inline uniform() noexcept
   
3. {
   
4.     return []( auto lower, auto upper ) noexcept
   
5.     {
   
6.         return [=]( auto& engine ) noexcept
   
7.         {
   
8.             return engine() * ( upper - lower ) + lower;
   
9.         };
   
10.     };
    
11. }
    

復(fù)制代碼

Generator

寫一個(gè)函數(shù)，這個(gè)函數(shù)接收兩個(gè)參數(shù)，第一個(gè)參數(shù)是一個(gè)生成具體隨機(jī)變量分布 F 的算法的函數(shù) f（分布），第二個(gè)參數(shù)是一個(gè)生成在 [0,1] 區(qū)間服從均勻分布隨機(jī)數(shù)的函數(shù) g（引擎）， 返回一個(gè)函數(shù) h，h 接收任意數(shù)量的參數(shù)并返回一個(gè)函數(shù) i，i 被調(diào)用后組合 f 與 g 生成服從 F 分布的隨機(jī)變量。

Write a function that takes two functions as parameters, the first one for a specified distribution F， the second one is an engine generating uniform random variates in range [0,1], returns a function that takes arbitrary parameters and returns a function that returns a random variate subject to distribution F.

1. 
   
2. template< typename Distribution, typename Engine >
   
3. auto make_generator( Distribution distribution, Engine engine ) noexcept
   
4. {
   
5.     auto engine_ = engine();
   
6.     return [=]( auto ... args ) mutable noexcept
   
7.     {
   
8.         return [=]() mutable noexcept
   
9.         {
   
10.             return distribution()( args... )( engine_ );
    
11.         };
    
12.     };
    
13. }   
    

復(fù)制代碼

如果 Generator 的引擎沒有指定，缺省使用 linear_congruential 引擎。

If the engine parameter for the make_generator is not specified, use linear_congruential.

1. 
   
2. template< typename Distribution >
   
3. auto make_generator( Distribution distribution ) noexcept
   
4. {
   
5.     return make_generator( distribution, mitchell_moore );
   
6. }
   

復(fù)制代碼

Test:

1. 
   
2. //normal_test.cc
   
3. auto v = make_generator(normal)();
   
4. 
   
5. map<int, int> sample;
   
6. 
   
7. for ( unsigned long int i = 0; i != 500; ++i )
   
8.     sample[static_cast<int>(round(4.0*v()))]++;
   
9. 
   
10. for ( auto element : sample )
    
11. {
    
12.     cout << "\n" << element.first << "\t";
    
13.     for ( auto j = 0; j < element.second; ++j ) cout << "*";
    
14. }
    

復(fù)制代碼

Output:

1. 
   
2. -13 *
   
3. -12 *
   
4. -10 **
   
5. -9  *******
   
6. -8  ******
   
7. -7  ***************
   
8. -6  ********************
   
9. -5  ******************
   
10. -4  *******************************
    
11. -3  ************************************
    
12. -2  *****************************************
    
13. -1  *********************************************
    
14. 0   *******************************************************
    
15. 1   *************************************************
    
16. 2   ******************************************
    
17. 3   *******************************************
    
18. 4   *****************************
    
19. 5   ***************************
    
20. 6   ********
    
21. 7   ****
    
22. 8   ***********
    
23. 9   ****
    
24. 10  **
    

復(fù)制代碼

其他引擎 -- More Engines
mt19937

1. 
   
2. inline auto mt19937() noexcept
   
3. {
   
4.     std::uint_fast64_t                      mt[312];
   
5.     std::uint_fast64_t                      mag01[2];
   
6.     std::size_t                             mti;
   
7. 
   
8.     mt[0]    = make_seed();
   
9.     mag01[0] = 0ULL;
   
10.     mag01[1] = 0xB5026F5AA96619E9ULL;
    
11. 
    
12.     for ( mti = 1; mti < 312; ++mti )
    
13.     {
    
14.         mt[mti] =  6364136223846793005ULL * ( mt[mti-1] ^( mt[mti-1] >> 62 ) );
    
15.         mt[mti] += mti;
    
16.     }
    
17. 
    
18.     return [=]() mutable noexcept
    
19.     {
    
20.         std::uint_fast64_t x;
    
21. 
    
22.         if ( mti > 311 )
    
23.         {
    
24.             for ( std::size_t i = 0; i < 156; ++i )
    
25.             {
    
26.                 x = ( mt[i] & 0xFFFFFFFF80000000ULL ) | ( mt[i+1] & 0x7FFFFFFFULL );
    
27.                 mt[i] = mt[i+156] ^( x >> 1 ) ^ mag01[x&1];
    
28.             }
    
29. 
    
30.             for ( std::size_t i = 156; i < 311; ++i )
    
31.             {
    
32.                 x = ( mt[i] & 0xFFFFFFFF80000000ULL ) | ( mt[i+1] & 0x7FFFFFFFULL );
    
33.                 mt[i] = mt[i-156] ^( x >> 1 ) ^ mag01[x&1];
    
34.             }
    
35. 
    
36.             x = ( mt[311] & 0xFFFFFFFF80000000ULL ) | ( mt[0] & 0x7FFFFFFFULL );
    
37.             mt[311] = mt[155] ^( x >> 1 ) ^ mag01[x&1];
    
38.             mti = 0;
    
39.         }
    
40. 
    
41.         x  = mt[mti++];
    
42.         x ^= ( x >> 29 ) & 0x5555555555555555ULL;
    
43.         x ^= ( x << 17 ) & 0x71D67FFFEDA60000ULL;
    
44.         x ^= ( x << 37 ) & 0xFFF7EEE000000000ULL;
    
45.         x ^= ( x >> 43 );
    
46. 
    
47.         return static_cast<long double>( x ) / static_cast<long double>( std::numeric_limits<std::uint_fast64_t>::max() );
    
48.     };
    
49. }
    

復(fù)制代碼

Reference

    Makoto Matsumoto, Takuji Nishimura: Mersenne twister. A 623-dimensionally equidistributed uniform pseudorandom number generator. In: ACM Transactions on Modeling and Computer Simulation. 8, 1998, ISSN 1049-3301, S. 3–30.

lagged_fibonacci

1. 
   
2. inline auto lagged_fibonacci() noexcept
   
3. {
   
4.     std::uint_fast64_t data[44497];
   
5.     std::size_t mti = make_seed() % 21034;
   
6. 
   
7.     data[0] = 1;
   
8.     data[1] = 1;
   
9. 
   
10.     for ( std::size_t index = 2; index != 44497; ++index )
    
11.         data[index] = data[index-2] + data[index-1];
    
12. 
    
13.     return[=]() mutable noexcept
    
14.     {
    
15.         if ( mti > 21034 )
    
16.         {
    
17.             mti = 0;
    
18.             std::copy( data+21034, data+44497, data );
    
19.             for ( std::size_t index = 23463; index != 44497; ++index ) data[index] = data[index-2] + data[index-1];
    
20.         }
    
21. 
    
22.         auto const mtj  = mti + 23463;
    
23.         auto const x    = data[mti++] + data[mtj];
    
24.         auto const ans  = static_cast<long double>( x ) / static_cast<long double>( std::numeric_limits<std::uint_fast64_t>::max() );
    
25.         return ans;
    
26.     };
    
27. }
    

復(fù)制代碼

Reference

    A New Class of Random Number Generators, George Marsaglia and Arif Zaman, The Annals of Applied Probability, Vol. 1, No. 3, 1991

其他非均勻分布 -- More Non-Uniform Distributions

If not specified, the reference for the distribution generator is

    Devroye, Luc 1986, Non-Uniform Random Variate Generation. New York: Springer-Verlag.

gaussian

1. 
   
2. inline auto gaussian() noexcept
   
3. {
   
4.     return []( auto mean, auto variance ) noexcept
   
5.     {
   
6.         return [=]( auto& engine ) noexcept
   
7.         {
   
8.             long double const normal_ = normal()()( engine );
   
9.             return normal_ * variance + mean;
   
10.         };
    
11.     };
    
12. }
    

復(fù)制代碼

當(dāng)均值為 0， 方差為 1 的時(shí)候，高斯分布退化為正態(tài)分布。

Normal distribution is a special case of Gaussian distribution when the average is 0 and the variance is 1.

Test:

1. 
   
2. auto v = vg::make_generator(vg::gaussian, vg::mt19937)( 0.0, 4.0 );
   
3. 
   
4. map< int, int > sample;
   
5. 
   
6. for ( unsigned long int i = 0; i != 500; ++i )
   
7.     sample[static_cast<int>(round(v()))]++;
   
8. 
   
9. for ( auto i = sample.begin(); i != sample.end(); ++i )
   
10. {
    
11.     cout << "\n" << (*i).first << "\t";
    
12.     for ( auto j = 0; j < (*i).second; ++j ) cout << "*";
    
13. }
    

復(fù)制代碼

Output:

1. 
   
2. -12 *
   
3. -11 *
   
4. -10 **
   
5. -9  *********
   
6. -8  *******
   
7. -7  *********
   
8. -6  *******************
   
9. -5  **************************
   
10. -4  **************************
    
11. -3  *********************************************
    
12. -2  *****************************************
    
13. -1  **********************************************
    
14. 0   **********************************************************
    
15. 1   ****************************************
    
16. 2   ******************************************
    
17. 3   *************************************
    
18. 4   *********************************
    
19. 5   *****************************
    
20. 6   *********
    
21. 7   *****
    
22. 8   *****
    
23. 9   *****
    
24. 10  ***
    
25. 11  *
    

復(fù)制代碼

gamma

1. 
   
2. inline auto gamma() noexcept
   
3. {
   
4.     return []( auto alpha, auto beta ) noexcept
   
5.     {
   
6.         assert( alpha > 0.0 );
   
7.         assert( beta > 0.0 );
   
8. 
   
9.         return [=]( auto& engine ) noexcept
   
10.         {
    
11.             long double const zero  = 0.0;
    
12.             long double const one   = 1.0;
    
13.             long double const three = 3.0;
    
14.             long double const a     = alpha < one ? alpha + one : alpha;
    
15.             long double const d     = a - one / three;
    
16.             long double const c     = one / (three * std::sqrt( d ));
    
17.             long double u           = zero;
    
18.             long double v           = zero;
    
19.             long double x           = zero;
    
20. 
    
21.             for ( ;; )
    
22.             {
    
23.                 for ( ;; )
    
24.                 {
    
25.                     x = normal()()( engine );
    
26.                     v = one + c * x;
    
27. 
    
28.                     if ( v > zero ) break;
    
29.                 }
    
30. 
    
31.                 const long double xx = x * x;
    
32. 
    
33.                 v *= v * v;
    
34. 
    
35.                 u = engine();
    
36. 
    
37.                 if ( u < one - 0.0331 * xx * xx )
    
38.                     break;
    
39. 
    
40.                 if ( std::log( u ) < 0.5 * xx + d * ( one - v + std::log( v ) ) )
    
41.                     break;
    
42.             }
    
43. 
    
44.             long double tmp = d * v;
    
45. 
    
46.             if ( alpha < one )
    
47.                 tmp *= std::pow( engine(), one / alpha );
    
48. 
    
49.             return tmp * beta;
    
50.         };
    
51.     };
    
52. }
    

復(fù)制代碼

Reference

    Marsaglia and Tsang, "A Simple Method for generating gamma variables", ACM Transactions on Mathematical Software, Vol 26, No 3 (2000), p363-372.

Test:

1. 
   
2. auto v = vg::make_generator(vg::gamma)( 1.0, 2.0 );
   
3. 
   
4. map<int, int> sample;
   
5. 
   
6. for ( unsigned long int i = 0; i != 100; ++i )
   
7.     sample[static_cast<int>(std::round(v()))]++;
   
8. 
   
9. for ( auto element : sample )
   
10. {
    
11.     cout << "\n" << element.first << "\t";
    
12.     for ( auto j = 0; j < element.second; ++j ) cout << "*";
    
13. }
    

復(fù)制代碼

Output:

1. 
   
2. 0   ***********************************************
   
3. 1   **************************
   
4. 2   **********
   
5. 3   ***
   
6. 4   ***
   
7. 5   ***
   
8. 6   *****
   
9. 8   *
   
10. 9   *
    
11. 11  *      
    

復(fù)制代碼

beta

1. 
   
2. inline auto beta() noexcept
   
3. {
   
4.     return []( auto alpha, auto beta ) noexcept
   
5.     {
   
6.         return [=]( auto& engine ) noexcept
   
7.         {
   
8.             long double one = 1.0;
   
9.             auto const a = gamma()( alpha, one )( engine );
   
10.             auto const b = gamma()( beta, one )( engine );
    
11.             return a / ( a + b );
    
12.         };
    
13.     };
    
14. }
    

復(fù)制代碼

如果 G_a 和 G_b 不相關(guān)，并且分別服從 gamma(a,1), gamma(b,1) 分布，那么 G_a/(G_a+G_b) 服從 beta(a,b) 分布。

If G_a and G_b are indenpendent gamma(a,1), gamma(b,1) random variables, then G_a/(G_a+G_b) is beta(a,b) distributed.


更多在 https://github.com/fengwang/vg

cokeboL

珍愛生命，遠(yuǎn)離C++

lost_templar

cokeboL 發(fā)表于 2015-06-15 21:21
珍愛生命，遠(yuǎn)離C++

其實(shí)， C++11/14 已經(jīng)重新發(fā)明了 C++；
函數(shù)式編程，放棄了 class/struct，遠(yuǎn)離了虛擬，繼承，多態(tài)，深拷貝，淺拷貝等陰暗角落，尤其是對(duì) GP 依賴的減少，使得智商因素對(duì)代碼質(zhì)量的影響沒有那么大了。
建議重新審視一番。

Life is too short to learn C++ 已經(jīng)不再總是成立了。

fender0107401

我只能說(shuō)：C++真強(qiáng)大、lz真會(huì)玩。

C++是個(gè)筐，啥都能往里裝。

cokeboL

回復(fù) 3# lost_templar


遠(yuǎn)離的原因不是難學(xué)，而是為了學(xué)語(yǔ)言花費(fèi)的時(shí)間太久，審視成本太大
而且，c++這些東西，不用c++都能搞定而且花費(fèi)的成本更少

windoze

回復(fù) 5# cokeboL

汽車能跑飛機(jī)能飛船能下水，但問(wèn)題是有時(shí)候你就需要在一個(gè)項(xiàng)目里上天入地下海挨個(gè)干一遍，此時(shí)除了搞一輛會(huì)飛的水陸兩用車你還能怎么辦？

stevenyi

functional programming體現(xiàn)在什么地方？

cokeboL

回復(fù) 6# windoze


    如果一定需要啥，就搞啥。。

windoze

回復(fù) 8# cokeboL

C++最適合的場(chǎng)合就是JVM之類的項(xiàng)目，既有機(jī)器級(jí)底層操作，又有上層的高度抽象，外加對(duì)性能還有變態(tài)的需求，代碼規(guī)模還特別大。

你說(shuō)這一類的項(xiàng)目，除了C++你還有什么可以用？Rust？

lxyscls

lost_templar 發(fā)表于 2015-06-15 21:55
其實(shí)， C++11/14 已經(jīng)重新發(fā)明了 C++；
函數(shù)式編程，放棄了 class/struct，遠(yuǎn)離了虛擬，繼承，多態(tài)，深 ...

我靠，這么說(shuō)正在學(xué)的過(guò)時(shí)了