Functional Programming in C++14

lost_templar

这一阵子折腾 C++14，练手的时候把自己之前写过的一个 生成非均匀随机变量的库重构了一下，还没有搞完，记录了重构的大概过程，先发过来大家瞧瞧。

https://github.com/fengwang/vg

标准组件 (Stardard Components)

Seed

写一个函数，这个函数返回一个随机整数。

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. }
    

复制代码

引擎 -- Engine

写一个函数，这个函数返回另外一个函数，返回的函数被调用后返回一个在 [0,1] 区间服从均匀分布的随机数。

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. }
    

复制代码

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
正态分布 -- Normal

写一个函数，这个函数返回一个函数 g, g 返回一个函数 h，h 接收一个函数 e 作为参数，h 返回一个服从 N(0,1) 分布的随机数，而作为函数 h 参数的函数 e 被调用后会返回一个在区间 [0,1] 上服从均匀分布的随机数。

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. }
    

复制代码

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

仿照正态分布的生成函数，写一个在区间 [lower, upper] 上服从均匀分布的函数。

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. }
    

复制代码

Generator

写一个函数，这个函数接收两个参数，第一个参数是一个生成具体随机变量分布 F 的算法的函数 f（分布），第二个参数是一个生成在 [0,1] 区间服从均匀分布随机数的函数 g（引擎）， 返回一个函数 h，h 接收任意数量的参数并返回一个函数 i，i 被调用后组合 f 与 g 生成服从 F 分布的随机变量。

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. }   
    

复制代码

如果 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. }
   

复制代码

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. }
    

复制代码

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  **
    

复制代码

其他引擎 -- 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. }
    

复制代码

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. }
    

复制代码

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. }
    

复制代码

当均值为 0， 方差为 1 的时候，高斯分布退化为正态分布。

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. }
    

复制代码

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  *
    

复制代码

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. }
    

复制代码

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. }
    

复制代码

Output:

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

复制代码

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. }
    

复制代码

如果 G_a 和 G_b 不相关，并且分别服从 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

珍爱生命，远离C++

lost_templar

cokeboL 发表于 2015-06-15 21:21
珍爱生命，远离C++

其实， C++11/14 已经重新发明了 C++；
函数式编程，放弃了 class/struct，远离了虚拟，继承，多态，深拷贝，浅拷贝等阴暗角落，尤其是对 GP 依赖的减少，使得智商因素对代码质量的影响没有那么大了。
建议重新审视一番。

Life is too short to learn C++ 已经不再总是成立了。

fender0107401

我只能说：C++真强大、lz真会玩。

C++是个筐，啥都能往里装。

cokeboL

回复 3# lost_templar


远离的原因不是难学，而是为了学语言花费的时间太久，审视成本太大
而且，c++这些东西，不用c++都能搞定而且花费的成本更少

windoze

回复 5# cokeboL

汽车能跑飞机能飞船能下水，但问题是有时候你就需要在一个项目里上天入地下海挨个干一遍，此时除了搞一辆会飞的水陆两用车你还能怎么办？

stevenyi

functional programming体现在什么地方？

cokeboL

回复 6# windoze


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

windoze

回复 8# cokeboL

C++最适合的场合就是JVM之类的项目，既有机器级底层操作，又有上层的高度抽象，外加对性能还有变态的需求，代码规模还特别大。

你说这一类的项目，除了C++你还有什么可以用？Rust？

lxyscls

lost_templar 发表于 2015-06-15 21:55
其实， C++11/14 已经重新发明了 C++；
函数式编程，放弃了 class/struct，远离了虚拟，继承，多态，深 ...

我靠，这么说正在学的过时了