1. C++ / Говнокод #14536

    +55

    1. 01
    2. 02
    3. 03
    4. 04
    5. 05
    6. 06
    7. 07
    8. 08
    9. 09
    10. 10
    11. 11
    12. 12
    13. 13
    14. 14
    15. 15
    16. 16
    17. 17
    18. 18
    19. 19
    20. 20
    21. 21
    22. 22
    23. 23
    24. 24
    25. 25
    26. 26
    27. 27
    28. 28
    29. 29
    30. 30
    31. 31
    32. 32
    33. 33
    34. 34
    35. 35
    36. 36
    37. 37
    38. 38
    39. 39
    40. 40
    41. 41
    42. 42
    43. 43
    44. 44
    45. 45
    46. 46
    47. 47
    48. 48
    49. 49
    50. 50
    51. 51
    52. 52
    53. 53
    #include <iostream>
    #include <thread>
    #include <cmath>
    #include <cstdlib>
    #include <algorithm>
    #include <iterator>
    #include <future>
    
    constexpr std::size_t array_size = 1e7;
    constexpr std::size_t chunks_number = 1e3;
    
    double local_pi(const std::size_t array_length)
    {
        std::vector<std::pair<double, double>> array_xy;       
        array_xy.reserve(array_length);      
        std::generate_n(std::back_inserter(array_xy), array_length, [](){
          return std::make_pair(rand()/static_cast<double>(RAND_MAX), 
    			    rand()/static_cast<double>(RAND_MAX));
          });
        std::vector<bool> bool_z;
        bool_z.reserve(array_length);
        auto bool_z_insterter = std::back_inserter(bool_z);
        std::for_each(array_xy.cbegin(), array_xy.cend(), [&](const std::pair<double, double>& ref){
          const auto func = [](const std::pair<double, double>& ref){
    	return sqrt(pow(ref.first, 2) + pow(ref.second, 2)) <= 1.0;
          };
          bool_z_insterter = func(ref);
        });
        const auto sum = 4.0 * (std::accumulate(bool_z.cbegin(), bool_z.cend(), 0.0)/array_length);
        return sum;
    }
    int main(){
      const auto concurrency_level = std::thread::hardware_concurrency();
      std::vector< double > chunks_vector;
      chunks_vector.reserve(chunks_number);
      for (std::size_t j = 0; j < chunks_number; ++j) {
        std::vector<std::future<double>> futures;
        futures.reserve(concurrency_level-1);
        for (std::size_t i = 0; i < concurrency_level-1; ++i){
          auto f = std::async(std::launch::async, local_pi, array_size/concurrency_level);
          futures.push_back(std::move(f));
        }
        auto pi = local_pi(array_size/concurrency_level);
        for (auto&& x : futures ){
          pi += x.get();
        }
        pi /= concurrency_level;
        chunks_vector.push_back(pi);
      }
      const auto total_pi = std::accumulate(chunks_vector.cbegin(), chunks_vector.cend(), 0.0)/chunks_vector.size();
      std::cout << "Pi equals = " << std::scientific << total_pi << std::endl;
      return 0;
    }

    Вычисляем pi методом Монте-Карло

    Запостил: Abbath, 10 Февраля 2014

    Комментарии (63) RSS

    • Неплохой очерк "месть преподу". Был бы отличный, если б не строки 17-18 с неканоничным рандомом :)
      Ответить
      • > неканоничным рандомом
        Кстати потестил сейчас каноничный рандом (std::mt19937_64 + std::uniform_real_distribution<double>) - и у меня челюсть отпала: он в 100(!) раз быстрее rand(). (Тестил на упрощенном коде, с выпиленными массивами).

        UPD: А если убрать всю эту шелуху с многопоточностью (concurrency_level = 1) - то вдвое.

        Значит или в rand() упрятан мутекс, или это thread-local переменная для seed'а в нем так ужасно сказывается на перфомансе.

        Мораль: rand() - сраное говно, как и любая функция, которая пользуется глобальным состоянием. А генераторы с локальным состоянием рулят ;)
        Ответить
        • P.S. Массивы замедляют код всего вдвое. Все-таки большую часть тормозов создает именно rand().
          Ответить
        • А еще rand не uniform, потому что использует конгруэнтный генератор. И гарантированно возвращает значение только из диапазона 0..32767.
          То есть rand() - вдвойне УГ
          Ответить
          • uniform? А по русски математически это как?
            Ответить
          • паскальный ранд круче, потому что там внутри не деление, а длинное умножение с возвратом result.hi32
            Ответить
            • > паскальный ранд круче
              Среди прочего линейно-конгруэнтного говна с неявным глобальным состоянием (скажи "нет" многопоточности) - он крутой ;)

              Но согласись, кресто11блядский <random> (он же boost::random) няшней на порядки :3
              Ответить
              • в бусте даже табличка была с описанием какой алгоритм насколько быстрее и няшнее, их там десятки же
                Ответить
                • Ну и всегда можно замутить свой tblib::nostalgia::pascal_gen с длинным умножением и возвратом старшей половины.
                  Ответить
                  • первое, что сделал
                    надеюсь, что это поможет точно портировать некоторые вещи, тонко подогнанные под конкретный генератор
                    Ответить
                    • boost::random::linear_congruential_engine<>
                      подставь константы из дельфячего рандома да перепроверь

                      насколько я понимаю, это
                      boost::random::linear_congruental_engine<uint32_t, 134775813, 1, 4294967296>
                      Ответить
                      • либо я слегка напиздел с 2^32 в uint32_t :)
                        Ответить
                        • Ответить
                        • If the template parameter m is 0, the modulus m used is std::numeric_limits<_UIntType>::max() plus 1.

                          Короче 0 туда надо передавать ;)

                          4294967296 по сути тоже 0, но выдаст ворнинг. Некрасиво.
                          Ответить
                          • тут надо посмотреть, этот генератор же выдаст младшую часть от результата, а в борланде берут старшую
                            надо где-то ещё чего-нить подмешать
                            Ответить
                            • Там тоже младшая, см. исследование ниже по треду. Там где старшую часть берут - это уже не сам генератор, а преобразование интервала в меньший, аля distribution в крестах-бустах.

                              В целом сам генератор подходит, а вот uniform_int_distribution - не совпадает, и нужно мутить tblib::almost_uniform_int_distribution через умножение и старшую часть ;)
                              Ответить
                    • > первое, что сделал
                      О, так ты сдался и все-таки поставил буст?
                      Ответить
                      • Нет, я просто навелосипедил генератор
                        только в класс обернул

                        Кстати, вопрос: какого хуя гцц ругается на такаую перегрузку:

                        int random(int n);
                        float random();
                        Ответить
                        • Не ругается. У тебя поди n с дефолтом там?

                          P.S. А, оно у тебя с random() из <stdlib.h> конфликтануло. Лишний повод юзать namespace :)
                          Ответить
                          • Так оно же в классе обёрнуто.
                            А теперь вызови float f = pasrnd.random();
                            int n = pasrnd.random(5);
                            Ответить
                            • http://ideone.com/XSU33C

                              Все нормально, как и ожидалось.

                              P.S. 146%, что у тебя n с дефолтом:
                              int random(int n = 42);
                              Тогда действительно конфликтанут.
                              Ответить
                              • нахуй дефолт
                                я вчера весь день угорал на крестофоруме, мне объясняли, что дефолтные аргументы в конструкторах Point3D - это нормально и писать
                                Point3D(1) намного лучше, чем Point3D(1.0f, 0.0f, 0.0f)
                                Ответить
                                • В GLSL правила немного другие. vec3(1.0) == vec3(1.0, 1.0, 1.0). mat4(1.0) - единичная матрица.
                                  https://www.opengl.org/wiki/Data_Type_%28GLSL%29#Constructors
                                  С векторами, соглашусь, спорно, а вот с матрицами весьма удобно.
                                  Ответить
                                  • > vec3(1.0) == vec3(1.0, 1.0, 1.0)
                                    А не vec3(1.0, 0.0, 0.0) разве?
                                    Ответить
                                    • Там по моей ссылке есть примеры. Вот ещё цитата из спецификации языка:
                                      >If there is a single scalar parameter to a vector constructor, it is used to initialize all components of the
                                      constructed vector to that scalar’s value. If there is a single scalar parameter to a matrix constructor, it is
                                      used to initialize all the components on the matrix’s diagonal, with the remaining components initialized
                                      to 0.0.
                                      Ответить
                                    • Мне объясняли, что это же пиздец как очевидно должно быть, что
                                      vec3(1) - это может быть только vec3(1.0f, 0.0f, 0.0f) и ничем другим это быть не может и что я тормоз, что не понимаю этой очевидности.
                                      Ответить
                                  • с матрицами - ну может быть, в том плане, что умножение на такую матрицу этот как умножение на число, хотя тоже лучше сделать не конструктор, а что-то типа SetDiagonal(float v)

                                    с векторами не то что спорно, это просто крестоблядское мудоёбство экономить пару символов и ради чего
                                    Ответить
                        • Кстати, кинь, пожалуйста, реализацию делфёвого RNG. Не могу нигде ее нагуглить, а делфи под рукой нет.
                          Ответить
                          • #pragma once
                            #include "time.h"
                            
                            struct Random {
                            	unsigned int randseed;
                            
                            	void NextSeed () 
                            	{
                            		randseed = randseed * 0x08088405 + 1;
                            	}
                            
                            	Random () : randseed ((unsigned int)(time(NULL))) {}
                            
                            	Random (int seed) : randseed(seed) {}
                            
                            	inline int random (int range)
                            	{
                            		NextSeed ();
                            		return int(((long long)(randseed)*(long long)(range)) >> 32);
                            	}
                            
                            	inline double frandom ()
                            	{
                            		const double rev_2e32 = 1.0/(65536.0*65536.0);
                            		NextSeed ();
                            		return randseed * rev_2e32;
                            	}
                            };

                            inline не потому что я царь, а потому что линковщик матерится
                            Ответить
                            • оригинал такой:
                              procedure       _RandInt;
                              asm
                              {     ->EAX     Range   }
                              {     <-EAX     Result  }
                                      PUSH    EBX
                              {$IFDEF PIC}
                                      PUSH    EAX
                                      CALL    GetGOT
                                      MOV     EBX,EAX
                                      POP     EAX
                                      MOV     ECX,[EBX].OFFSET RandSeed
                                      IMUL    EDX,[ECX],08088405H
                                      INC     EDX
                                      MOV     [ECX],EDX
                              {$ELSE}
                                      XOR     EBX, EBX
                                      IMUL    EDX,[EBX].RandSeed,08088405H
                                      INC     EDX
                                      MOV     [EBX].RandSeed,EDX
                              {$ENDIF}
                                      MUL     EDX
                                      MOV     EAX,EDX
                                      POP     EBX
                              end;
                              
                              procedure       _RandExt;
                              const two2neg32: double = ((1.0/$10000) / $10000);  // 2^-32
                              asm
                              {       FUNCTION _RandExt: Extended;    }
                              
                                      PUSH    EBX
                              {$IFDEF PIC}
                                      CALL    GetGOT
                                      MOV     EBX,EAX
                                      MOV     ECX,[EBX].OFFSET RandSeed
                                      IMUL    EDX,[ECX],08088405H
                                      INC     EDX
                                      MOV     [ECX],EDX
                              {$ELSE}
                                      XOR     EBX, EBX
                                      IMUL    EDX,[EBX].RandSeed,08088405H
                                      INC     EDX
                                      MOV     [EBX].RandSeed,EDX
                              {$ENDIF}
                              
                                      FLD     [EBX].two2neg32
                                      PUSH    0
                                      PUSH    EDX
                                      FILD    qword ptr [ESP]
                                      ADD     ESP,8
                                      FMULP   ST(1), ST(0)
                                      POP     EBX
                              end;
                              Ответить
                            • Ну в общем понятно. Первая половина генератора представляет собой классический конгруэнтник, такой же, как и написанный выше defecate++:
                              std::linear_congruential_engine<uint32_t, 134775813, 1, 0> gen(seed);
                              А вторая половина (где mul edx; mov eax, edx) как раз то самое
                              uniform_int_distribution<int> dist(0, 99);
                              int rnd = dist(gen);
                              Вот только на тестах странно получилось - 8 чисел из 10кк отличаются от твоих на единичку. Надо поизучать, почему.
                              Ответить
                              • Первый такой сид - 3598178048.

                                После NextSeed() получаем 3006477056, которое по показаниям питона и хаскеля должно дать (3006477056 * 100) `div` (2^32) = 69. Код Тараса возвращает нам 69. Крестокод какого-то хуя - 70.
                                Ответить
                                • хихи, проверяю на своём калькуляторе (внутри 10-байтный float)
                                  > 3006477056*100 / (2**32)
                                  > 69.99999880790710449
                                  походу говнобуст зачем-то воткнул плавучку
                                  Ответить
                                  • > походу говнобуст зачем-то воткнул плавучку
                                    Вот такой код там для целых типов: http://pastebin.com/CuchWzhY

                                    Будем разбираться ;)
                                    Ответить
                                    • > Будем разбираться ;)
                                      Да, в этом надо долго разбираться...
                                      Ответить
                                      • > Да, в этом надо долго разбираться...
                                        Все, разобрался. В общем тут юзается точный алгоритм с отбрасыванием (например для генерации чисел от 0 до 99, последние 97 чисел от 4294967200 до 4294967296 тупо выкидываются и у генератора запрашивается новое, остальные же нарезаются на равные интервалы по 100).
                                        // rngrange - диапазон RNG
                                        // range - диапазон, который надо получить
                                        scaling = rngrange / range;
                                        past = scaling * range;
                                        do {
                                            ret = rng() - rngmin;
                                        } while (ret >= past);
                                        return ret / scaling;
                                        Вот из-за этого / scaling и получается несовпадение на границах интервалов (здесь все интервалы равные, в делфи - 3 штуки на 1 число меньше). А если бы попалось одно из тех 96 чисел - произошел бы рассинхрон с твоим алгоритмом.

                                        Так что если нужен в точности дельфиний генератор - заюзать uniform_int_distribution<int> dist(0, 99) нельзя.
                                        Ответить
                                        • > Все, разобрался.
                                          errata:
                                          последние 97 чисел от 4294967200 до 4294967296 тупо выкидываются читаем как последние 96 чисел от 4294967200 до 4294967295 тупо выкидываются

                                          в делфи - 3 штуки на 1 число меньше читаем как в делфи - 4 штуки на 1 число меньше
                                          Ответить
                                    • > Вот такой код там для целых типов: http://pastebin.com/CuchWzhY
                                      Что-то херня какая-то. Судя по этому коду он должен отбрасывать неудачные рандомы (собственно один из кошерных и точных алгоритмов по преобразованию диапазонов рандома). Но он же не отбрасывает (т.к. рассинхрона этого и твоего рандома после таких ситуаций не происходит), а работает как на флоатах ;(

                                      Ничего не понимаю. Ну не может же оно попадать во вторую специализацию, где true_type и работа с флоатами.

                                      UPD: Не тот класс смотрел, сейчас раскурю тот, который надо.
                                      Ответить
                        • Разобрался с перегрузкой, или ошибка еще осталась?
                          Ответить
                          • Откуда я знаю, у меня гцц только для ведроида, а эклипс я запускаю лишь по крайней нужде.
                            Очевидно, что дефолтного значения там быть не могло, я не настолько ебанулся.
                            Ответить
                      • Я как-то жабий Random портанул на плюсцы, чтобы буст не тянуть.
                        Ответить
      • Код был написан для:
        1)Смотри как я могу
        2)STL дрочерство
        3)саморазвитие
        4)тренировка парного программирования

        А преподу сдается матлабовский код.
        function p = f(n)
          summ = 0;
          total_length = n/1e6;
          for i=1:total_length
            x = rand(n/total_length,1);
            y = rand(n/total_length,1);
            z = sqrt(x.^2 + y.^2);
            summ += sum(z < 1);
          end
          p = (4*summ)/n;
        end
        Ответить
        • А-а-а, от оно чё...
          Я просто сам над оуеими халявщиками издевался: препод ожидает простую сишную программу, демонстрирующую алгоритм, а ему преподносят приплюсплюснутое нечто, полное СТЛа :)
          Ответить
        • А что должно произойти при итерации от целого к вещественному?
          Я бы на месте преподавателя за такие названия и форматирование заставлял пистать гусиным пером в прописи.
          Ответить
          • Будет считать до наибольшего целого меньше чем total_length же. Препод - доктор физико математических наук, он сам так пишет))
            Ответить
        • function m_pi = montecarlo_pi(num_runs)
            approximation = 0;
            range = 10;
            radius_squared = 1;
            for i = 1:num_runs
              x = rand(range, 1);
              y = rand(range, 1);
              approximation += sum(abs(floor(x .^ 2 + y .^ 2) .- radius_squared)) ./ range;
            end
            m_pi = (4 * approximation) / num_runs;
          end


          Без извлечения корней в каждой итерации.
          Ответить
          • function m_pi = montecarlo_pi(num_runs)
              range = 10;
              approximation = zeros(range, 1);
              radius_squared = 1;
              for i = 1:num_runs
                x = rand(range, 1);
                y = rand(range, 1);
                approximation = approximation + abs(floor(x .^ 2 + y .^ 2) .- radius_squared);
              end
              m_pi = (4 * sum(approximation)) / (num_runs * range);
            end


            Лол, не, тоже говно, вот так надо распараллеливать!

            PS. Я видел и не таких докторов, которые программируют на уровне бухгалтерши заведующей 1С-предприятием. Степень совсем не показатель.
            Ответить
            • А доктор наук разве обязан программировать лучше бухгалтерши-1сницы? Это же не его основная специализация. Может что-то полезное наваять для своей работы - и ладно. Ну кроме случаев, когда он доктор CS :)
              Ответить
              • Ну как бы а зачем такой доктор берется объяснять студентам, как код писать? Мне чет тяжело представить в какой отрасли кроме программирования нужно знать как реализовать (в отличие от использования) метода Монте Карло.
                Ответить
                • Любая естественная наука тащемта, физика например, или биология, или даже экономика, рассчитать что-то сложное в лоб. Этот метод еще при создании атомной бомбы применяли (собственно тогда и создали).

                  Кто же для тебя скрипт модели напишет если не сам/а.
                  Ответить
            • > распараллеливать!
              упоротый? в параллельном варианте там вообще цикла нет. Просто двумерный массив вместо одномерного.
              Ответить
              • И матрица, миллион на миллион, а чо?
                Ответить
                • >И матрица, миллион на 10
                  fixed. Да, будет быстрее чем цикл до миллиона.
                  Ответить
                  • Лень проверять (на работе Октавы нет), но предполагаю, что десять миллионов рандомов дадут точность ну так... до третьего знака, а как известно, на каждый следующий разряд нужно на порядок больше рандомов.
                    Хотя, конечно, это лаба, и точность тут не уперлась.
                    Ответить
    • > total_pi
      Это точно. Полный пи.
      Ответить
      • Юзать 2 вектора, лямбду и три stl алгоритма там, где хватило бы одного простейшего цикла на 4 строчки... stl головного мозга, однако ;)
        Ответить
    • > bool_z_insterter
      !!!
      Ответить

    Добавить комментарий