corsika/detail/framework/process/InteractionProcess.hpp

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+   *  traits test for InteractionProcess::doInteraction method.
+   */
+
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  struct has_method_doInteract : public detail::has_method_signature<TReturn, TArgs...> {
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! signature of templated method
+    template <class T>
+    static decltype(testSignature(&T::template doInteraction<TTemplate>)) test(
+        std::nullptr_t);
+
+    //! signature of non-templated method
+    template <class T>
+    static decltype(testSignature(&T::doInteraction)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  bool constexpr has_method_doInteract_v =
+      has_method_doInteract<TProcess, TReturn, TTemplate, TArgs...>::value;
+
+  /**
+   *  traits test for TEMPLATED InteractionProcess::getCrossSection method (PROPOSAL).
+   */
+
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  struct has_method_getCrossSectionTemplate
+      : public detail::has_method_signature<TReturn, TArgs...> {
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template getCrossSection<TTemplate>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::getCrossSection)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type shortcut
+  template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+  bool constexpr has_method_getCrossSectionTemplate_v =
+      has_method_getCrossSectionTemplate<TProcess, TReturn, TTemplate, TArgs...>::value;
+
+  /**
+   *  traits test for InteractionProcess::getCrossSection method.
+   */
+
+  template <class TProcess, typename TReturn, typename... TArgs>
+  struct has_method_getCrossSection
+      : public detail::has_method_signature<TReturn, TArgs...> {
+
+    ///! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template getCrossSection<TArgs...>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::getCrossSection)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     * @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! value traits type shortcut
+  template <class TProcess, typename TReturn, typename... TArgs>
+  bool constexpr has_method_getCrossSection_v =
+      has_method_getCrossSection<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/ProcessSequence.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/process/BaseProcess.hpp>
+#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
+#include <corsika/framework/process/DecayProcess.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/process/SecondariesProcess.hpp>
+#include <corsika/framework/process/StackProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+#include <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  inline ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                         IndexProcess2>::ProcessSequence(TProcess1 in_A, TProcess2 in_B)
+      : A_(in_A)
+      , B_(in_B) {
+
+    // make sure only BaseProcess types TProcess1/2 are passed
+
+    static_assert(is_process_v<TProcess1>,
+                  "can only use process derived from BaseProcess in "
+                  "ProcessSequence, for Process 1");
+    static_assert(is_process_v<TProcess2>,
+                  "can only use process derived from BaseProcess in "
+                  "ProcessSequence, for Process 2");
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn ProcessSequence<
+      TProcess1, TProcess2, IndexStart, IndexProcess1,
+      IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+                                         typename TParticle::node_type const& from,
+                                         typename TParticle::node_type const& to) {
+
+    ProcessReturn ret = ProcessReturn::Ok;
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_boundary_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TProcess1
+        static_assert(
+            has_method_doBoundaryCrossing_v<TProcess1, ProcessReturn, TParticle>,
+            "TDerived has no method with correct signature \"ProcessReturn "
+            "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+            "required for "
+            "BoundaryCrossingProcess<TDerived>. ");
+
+        ret |= A_.doBoundaryCrossing(particle, from, to);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_boundary_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on TProcess2
+        static_assert(
+            has_method_doBoundaryCrossing_v<TProcess2, ProcessReturn, TParticle>,
+            "TDerived has no method with correct signature \"ProcessReturn "
+            "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+            "required for "
+            "BoundaryCrossingProcess<TDerived>. ");
+
+        ret |= B_.doBoundaryCrossing(particle, from, to);
+      }
+    }
+
+    return ret;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+      doContinuous(Step<TParticle>& step,
+                   [[maybe_unused]] ContinuousProcessIndex const limitId) {
+    ProcessReturn ret = ProcessReturn::Ok;
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence) {
+
+        ret |= A_.doContinuous(step, limitId);
+      } else if constexpr (is_continuous_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        //~ static_assert(
+        //~ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle&, TTrack&> ||
+        //~ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle&,
+        //~ TTrack const&> ||
+        //~ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle const&,
+        //~ TTrack const&>,
+        //~ "TDerived has no method with correct signature \"ProcessReturn "
+        //~ "doContinuous(TParticle [const]&,TTrack [const]&,bool)\" required for "
+        //~ "ContinuousProcess<TDerived>. ");
+
+        ret |= A_.doContinuous(
+            step, limitId == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(A_))));
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (process2_type::is_process_sequence) {
+        ret |= B_.doContinuous(step, limitId);
+      } else if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        //~ static_assert(
+        //~ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle&, TTrack&> ||
+        //~ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle&,
+        //~ TTrack const&> ||
+        //~ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle const&,
+        //~ TTrack const&>,
+        //~ "TDerived has no method with correct signature \"ProcessReturn "
+        //~ "doContinuous(TParticle [const]&,TTrack [const]&,bool)\" required for "
+        //~ "ContinuousProcess<TDerived>. ");
+
+        ret |= B_.doContinuous(
+            step, limitId == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(B_))));
+      }
+    }
+
+    return ret;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TSecondaries>
+  inline void
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                  IndexProcess2>::doSecondaries([[maybe_unused]] TSecondaries& vS) {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_secondaries_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TProcess1
+        static_assert(
+            has_method_doSecondaries_v<TProcess1, void, TSecondaries&> ||
+                has_method_doSecondaries_v<TProcess1, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        A_.doSecondaries(vS);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_secondaries_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on TProcess2
+        static_assert(
+            has_method_doSecondaries_v<TProcess2, void, TSecondaries&> ||
+                has_method_doSecondaries_v<TProcess2, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        B_.doSecondaries(vS);
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  inline bool ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                              IndexProcess2>::checkStep() {
+    bool ret = false;
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_stack_process_v<process1_type> ||
+                    (process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence)) {
+        ret |= A_.checkStep();
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_stack_process_v<process2_type> ||
+                    (process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence)) {
+        ret |= B_.checkStep();
+      }
+    }
+    return ret;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TStack>
+  inline void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                              IndexProcess2>::doStack(TStack& stack) {
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_stack_process_v<process1_type> ||
+                    (process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence)) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doStack_v<TProcess1, void, TStack&> ||
+                          has_method_doStack_v<TProcess1, void, TStack const&>,
+                      "TDerived has no method with correct signature \"void "
+                      "doStack(TStack [const]&)\" required for "
+                      "StackProcess<TDerived>. ");
+
+        if (A_.checkStep()) { A_.doStack(stack); }
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_stack_process_v<process2_type> ||
+                    (process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence)) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doStack_v<TProcess2, void, TStack&> ||
+                          has_method_doStack_v<TProcess2, void, TStack const&>,
+                      "TDerived has no method with correct signature \"void "
+                      "doStack(TStack [const]&)\" required for "
+                      "StackProcess<TDerived>. ");
+
+        if (B_.checkStep()) { B_.doStack(stack); }
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle, typename TTrack>
+  inline ContinuousProcessStepLength
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                  IndexProcess2>::getMaxStepLength(TParticle&& particle,
+                                                   TTrack&& vTrack) {
+    // if no other process in the sequence implements it
+    ContinuousProcessStepLength max_length(std::numeric_limits<double>::infinity() *
+                                           meter);
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence) {
+        ContinuousProcessStepLength const step = A_.getMaxStepLength(particle, vTrack);
+        max_length = std::min(max_length, step);
+      } else if constexpr (is_continuous_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_getMaxStepLength_v<TProcess1, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        ContinuousProcessStepLength const step(
+            A_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(A_))));
+        max_length = std::min(max_length, step);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (process2_type::is_process_sequence) {
+        ContinuousProcessStepLength const step = B_.getMaxStepLength(particle, vTrack);
+        max_length = std::min(max_length, step);
+      } else if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        static_assert(has_method_getMaxStepLength_v<TProcess2, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        ContinuousProcessStepLength const step(
+            B_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(B_))));
+        max_length = std::min(max_length, step);
+      }
+    }
+    return max_length;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline CrossSectionType
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+      getCrossSection([[maybe_unused]] TParticle const& projectile,
+                      [[maybe_unused]] Code const targetId,
+                      [[maybe_unused]] FourMomentum const& targetP4) const {
+
+    CrossSectionType tot = CrossSectionType::zero();
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_interaction_process_v<process1_type>) {
+
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess1,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+
+        if constexpr (has_signature_cx1) {
+          tot += A_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else { // for PROPOSAL
+          tot += A_.getCrossSection(projectile, projectile.getPID(),
+                                    {projectile.getEnergy(), projectile.getMomentum()});
+        }
+
+      } else if constexpr (process1_type::is_process_sequence) {
+        tot += A_.getCrossSection(projectile, targetId, targetP4);
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_interaction_process_v<process2_type>) {
+
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess2,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+
+        if constexpr (has_signature_cx1) {
+          tot += B_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else { // for PROPOSAL
+          tot += B_.getCrossSection(projectile, projectile.getPID(),
+                                    {projectile.getEnergy(), projectile.getMomentum()});
+        }
+
+      } else if constexpr (process2_type::is_process_sequence) {
+        tot += B_.getCrossSection(projectile, targetId, targetP4);
+      }
+    }
+    return tot;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TStack>
+  void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                       IndexProcess2>::doCascadeEquations(TStack& stack) {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence &&
+                    !process1_type::is_switch_process_sequence) {
+
+        A_.doCascadeEquations(stack);
+
+      } else if constexpr (is_cascade_equations_process_v<process1_type>) {
+
+        // interface checking on TProcess1
+        static_assert(has_method_doCascadeEquations_v<TProcess1,
+                                                      void,     // return type
+                                                      TStack&>, // parameter
+                      "TDerived has no method with correct signature \"void "
+                      "doCascadeEquations(TStack&)\" required for "
+                      "CascadeEquationsProcess<TDerived>. ");
+
+        A_.doCascadeEquations(stack);
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (process2_type::is_process_sequence &&
+                    !process2_type::is_switch_process_sequence) {
+
+        B_.doCascadeEquations(stack);
+
+      } else if constexpr (is_cascade_equations_process_v<process2_type>) {
+
+        // interface checking on TProcess2
+        static_assert(has_method_doCascadeEquations_v<TProcess2,
+                                                      void,     // return type
+                                                      TStack&>, // parameter
+                      "TDerived has no method with correct signature \"void "
+                      "doCascadeEquations(TStack&)\" required for "
+                      "CascadeEquationsProcess<TDerived>. ");
+
+        B_.doCascadeEquations(stack);
+      }
+    }
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                       IndexProcess2>::initCascadeEquations() {
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr ((process1_type::is_process_sequence &&
+                     !process1_type::is_switch_process_sequence) ||
+                    is_cascade_equations_process_v<process1_type>) {
+        A_.initCascadeEquations();
+      }
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr ((process2_type::is_process_sequence &&
+                     !process2_type::is_switch_process_sequence) ||
+                    is_cascade_equations_process_v<process2_type>) {
+        B_.initCascadeEquations();
+      }
+    }
+  } // namespace corsika
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TSecondaryView, typename TRNG>
+  inline ProcessReturn
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+      selectInteraction(TSecondaryView&& view, FourMomentum const& projectileP4,
+                        [[maybe_unused]] NuclearComposition const& composition,
+                        [[maybe_unused]] TRNG&& rng,
+                        [[maybe_unused]] CrossSectionType const cx_select,
+                        [[maybe_unused]] CrossSectionType cx_sum) {
+
+    // TODO: add check for cx_select > cx_tot
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence) {
+        // if A is a process sequence --> check inside
+        ProcessReturn const ret =
+            A_.selectInteraction(view, projectileP4, composition, rng, cx_select, cx_sum);
+        // if A_ did succeed, stop routine. Not checking other static branch B_.
+        if (ret != ProcessReturn::Ok) { return ret; }
+      } else if constexpr (is_interaction_process_v<process1_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components
+        // for selected process A
+
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess1,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                TProcess1,                   // process object
+                CrossSectionType,            // return type
+                decltype(projectile) const&, // template argument
+                decltype(projectile) const&, // parameters
+                Code, FourMomentum const&>;
+
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "TProcess1 has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<TProcess1>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /*std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.cbegin(),
+                              weightedCrossSections.cend(), CrossSectionType::zero());
+
+        } else { // this is for PROPOSAL
+          cx_sum += A_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        // check if we should execute THIS process and then EXIT
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess1
+            static_assert(
+                has_method_doInteract_v<TProcess1,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "TProcess1 has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<TProcess1>. ");
+
+            A_.template doInteraction(view, projectileId, targetId, projectileP4,
+                                      targetP4);
+
+          } else { // this is for PROPOSAL
+            A_.template doInteraction(view, projectileId, projectileP4);
+          }
+
+          return ProcessReturn::Interacted;
+        }
+      }
+    } // end branch A
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+
+      if constexpr (process2_type::is_process_sequence) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectInteraction(view, projectileP4, composition, rng, cx_select,
+                                    cx_sum);
+      } else if constexpr (is_interaction_process_v<process2_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components, for selected process B
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TProcess2,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                TProcess2,                    // process object
+                CrossSectionType,             // return type
+                decltype(*projectile) const&, // template argument
+                decltype(*projectile) const&, // parameters
+                Code,                         // parameters
+                FourMomentum const&>;
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "TProcess2 has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<TProcess1>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /* std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.begin(), weightedCrossSections.end(),
+                              CrossSectionType::zero());
+        } else { // this is for PROPOSAL
+          cx_sum += B_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        // check if we should execute THIS process and then EXIT
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess2
+            static_assert(
+                has_method_doInteract_v<TProcess2,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "TProcess1 has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<TProcess2>. ");
+
+            B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+          } else { // this is for PROPOSAL
+            B_.doInteraction(view, projectileId, projectileP4);
+          }
+          return ProcessReturn::Interacted;
+        }
+      }
+    } // end branch B_
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  template <typename TParticle>
+  inline InverseTimeType
+  ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+                  IndexProcess2>::getInverseLifetime(TParticle&& particle) {
+
+    InverseTimeType tot = 0 / second; // default value
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (is_decay_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+        tot += A_.getInverseLifetime(particle);
+      }
+    }
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (is_decay_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+        tot += B_.getInverseLifetime(particle);
+      }
+    }
+    return tot;
+  }
+
+  template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+            int IndexProcess2>
+  // select decay process
+  template <typename TSecondaryView>
+  inline ProcessReturn ProcessSequence<
+      TProcess1, TProcess2, IndexStart, IndexProcess1,
+      IndexProcess2>::selectDecay(TSecondaryView&& view,
+                                  [[maybe_unused]] InverseTimeType decay_inv_select,
+                                  [[maybe_unused]] InverseTimeType decay_inv_sum) {
+
+    // TODO: add check for decay_inv_select>decay_inv_tot
+
+    if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+                                                 // errors if process1_type is invalid
+      if constexpr (process1_type::is_process_sequence) {
+        // if A_ is a process sequence --> check inside
+        ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum);
+        // if A_ did succeed, stop routine here (not checking other static branch B_)
+        if (ret != ProcessReturn::Ok) { return ret; }
+      } else if constexpr (is_decay_process_v<process1_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += A_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
+                                                // decay_inv_sum / decay_inv_tot
+          // interface checking on TProcess1
+          static_assert(has_method_doDecay_v<TProcess1, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          A_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch A_
+    }
+
+    if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+                                                 // errors if process2_type is invalid
+      if constexpr (process2_type::is_process_sequence) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectDecay(view, decay_inv_select, decay_inv_sum);
+      } else if constexpr (is_decay_process_v<process2_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += B_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) {
+
+          // interface checking on TProcess1
+          static_assert(has_method_doDecay_v<TProcess2, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          B_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch B_
+    }
+
+    return ProcessReturn::Ok;
+  }
+
+  /**
+   * traits marker to identify objects containing any StackProcesses.
+   */
+  namespace detail {
+    // need helper alias to achieve this:
+    template <
+        typename TProcess1, typename TProcess2,
+        typename = typename std::enable_if_t<
+            contains_stack_process_v<TProcess1> || is_stack_process_v<TProcess1> ||
+                contains_stack_process_v<TProcess2> || is_stack_process_v<TProcess2>,
+            int>>
+    using enable_if_stack = ProcessSequence<TProcess1, TProcess2>;
+  } // namespace detail
+
+  template <typename TProcess1, typename TProcess2>
+  struct contains_stack_process<detail::enable_if_stack<TProcess1, TProcess2>>
+      : std::true_type {};
+
+} // namespace corsika

corsika/detail/framework/process/SecondariesProcess.hpp

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+   * traits test for SecondariesProcess::doSecondaries method.
+   */
+  template <class TProcess, typename TReturn, typename... TArg>
+  struct has_method_doSecondaries
+      : public detail::has_method_signature<TReturn, TArg...> {
+
+    //! method signature
+    using detail::has_method_signature<TReturn, TArg...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template doSecondaries<TArg...>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::doSecondaries)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  /**
+   * @file SecondariesProcess.hpp
+   *
+   * @brief value traits type.
+   */
+  template <class TProcess, typename TReturn, typename... TArg>
+  bool constexpr has_method_doSecondaries_v =
+      has_method_doSecondaries<TProcess, TReturn, TArg...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/StackProcess.hpp

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+     traits test for StackProcess::doStack method
+  */
+
+  template <class TProcess, typename TReturn, typename... TArgs>
+  struct has_method_doStack : public detail::has_method_signature<TReturn, TArgs...> {
+
+    //! method signature
+    using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template doStack<TArgs...>)) test(std::nullptr_t);
+
+    //! templated-template parameter option
+    template <template <typename> typename T>
+    static decltype(testSignature(&T<TArgs...>::template doStack)) test(std::nullptr_t);
+
+    //! non-templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::doStack)) test(std::nullptr_t);
+
+  public:
+    /**
+        @name traits results
+        @{
+    */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! @file StackProcess.hpp
+  //! value traits type
+  template <class TProcess, typename TReturn, typename... TArgs>
+  bool constexpr has_method_doStack_v =
+      has_method_doStack<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika

corsika/detail/framework/process/SwitchProcessSequence.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/BaseProcess.hpp>
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
+#include <corsika/framework/process/DecayProcess.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/process/SecondariesProcess.hpp>
+#include <corsika/framework/process/StackProcess.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn SwitchProcessSequence<
+      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+      IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+                                         typename TParticle::node_type const& from,
+                                         typename TParticle::node_type const& to) {
+    if (select_(particle)) {
+      if constexpr (is_boundary_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TSequence
+        if constexpr (is_boundary_process_v<process1_type>) {
+
+          static_assert(
+              has_method_doBoundaryCrossing_v<TSequence, ProcessReturn, TParticle&>,
+              "TDerived has no method with correct signature \"ProcessReturn "
+              "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+              "required for "
+              "BoundaryCrossingProcess<TDerived>. ");
+        }
+
+        return A_.doBoundaryCrossing(particle, from, to);
+      }
+    } else {
+
+      if constexpr (is_boundary_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on USequence
+        if constexpr (is_boundary_process_v<process2_type>) {
+
+          static_assert(
+              has_method_doBoundaryCrossing_v<USequence, ProcessReturn, TParticle>,
+              "TDerived has no method with correct signature \"ProcessReturn "
+              "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+              "required for "
+              "BoundaryCrossingProcess<TDerived>. ");
+        }
+
+        return B_.doBoundaryCrossing(particle, from, to);
+      }
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  inline ProcessReturn SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+                                             IndexProcess1, IndexProcess2>::
+      doContinuous(Step<TParticle>& step,
+                   [[maybe_unused]] ContinuousProcessIndex const idLimit) {
+    if (select_(step.getParticlePre())) {
+      if constexpr (process1_type::is_process_sequence) {
+        return A_.doContinuous(step, idLimit);
+      }
+      if constexpr (is_continuous_process_v<process1_type>) {
+
+        //        static_assert(
+        //            has_method_doContinuous_v<TSequence, ProcessReturn, TParticle&,
+        //            TTrack&> ||
+        //                has_method_doContinuous_v<TSequence, ProcessReturn, TParticle&,
+        //                                          TTrack const&> ||
+        //                has_method_doContinuous_v<TSequence, ProcessReturn, TParticle
+        //                const&,
+        //                                          TTrack const&>,
+        //            "TDerived has no method with correct signature \"ProcessReturn "
+        //            "doContinuous(TParticle[const]&,TTrack[const]&,bool)\" required for
+        //            " "ContinuousProcess<TDerived>. ");
+
+        return A_.doContinuous(
+            step, idLimit == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(A_))));
+      }
+    } else {
+      if constexpr (process2_type::is_process_sequence) {
+        return B_.doContinuous(step, idLimit);
+      }
+      if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on USequence
+        //        static_assert(
+        //            has_method_doContinuous_v<USequence, ProcessReturn, TParticle&,
+        //            TTrack&> ||
+        //                has_method_doContinuous_v<USequence, ProcessReturn, TParticle&,
+        //                                          TTrack const&> ||
+        //                has_method_doContinuous_v<USequence, ProcessReturn, TParticle
+        //                const&,
+        //                                          TTrack const&>,
+        //            "TDerived has no method with correct signature \"ProcessReturn "
+        //            "doContinuous(TParticle [const]&,TTrack[const]&,bool)\" required for
+        //            " "ContinuousProcess<TDerived>. ");
+
+        return B_.doContinuous(
+            step, idLimit == ContinuousProcessIndex(
+                                 static_cast<void const*>(std::addressof(B_))));
+      }
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TSecondaries>
+  inline void
+  SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+                        IndexProcess2>::doSecondaries(TSecondaries& vS) {
+    const auto& particle = vS.parent();
+    if (select_(particle)) {
+      if constexpr (is_secondaries_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+
+        // interface checking on TSequence
+        static_assert(
+            has_method_doSecondaries_v<TSequence, void, TSecondaries&> ||
+                has_method_doSecondaries_v<TSequence, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        A_.doSecondaries(vS);
+      }
+    } else {
+      if constexpr (is_secondaries_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+
+        // interface checking on USequence
+        static_assert(
+            has_method_doSecondaries_v<USequence, void, TSecondaries&> ||
+                has_method_doSecondaries_v<USequence, void, TSecondaries const&>,
+            "TDerived has no method with correct signature \"void "
+            "doSecondaries(TStackView [const]&)\" required for "
+            "SecondariesProcessProcess<TDerived>. ");
+
+        B_.doSecondaries(vS);
+      }
+    }
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle, typename TTrack>
+  inline ContinuousProcessStepLength
+  SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+                        IndexProcess2>::getMaxStepLength(TParticle& particle,
+                                                         TTrack& vTrack) {
+    if (select_(particle)) {
+      if constexpr (process1_type::is_process_sequence) {
+        return A_.getMaxStepLength(particle, vTrack);
+      }
+      if constexpr (is_continuous_process_v<process1_type>) {
+
+        // interface checking on TSequence
+        static_assert(has_method_getMaxStepLength_v<TSequence, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        return ContinuousProcessStepLength(
+            A_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(A_))));
+      }
+    } else {
+      if constexpr (process2_type::is_process_sequence) {
+        return B_.getMaxStepLength(particle, vTrack);
+      }
+      if constexpr (is_continuous_process_v<process2_type>) {
+
+        // interface checking on USequence
+        static_assert(has_method_getMaxStepLength_v<USequence, LengthType,
+                                                    TParticle const&, TTrack const&>,
+                      "TDerived has no method with correct signature \"LengthType "
+                      "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+                      "ContinuousProcess<TDerived>. ");
+
+        return ContinuousProcessStepLength(
+            B_.getMaxStepLength(particle, vTrack),
+            ContinuousProcessIndex(static_cast<void const*>(std::addressof(B_))));
+      }
+    }
+
+    // if no other process in the sequence implements it
+    return ContinuousProcessStepLength(std::numeric_limits<double>::infinity() * meter);
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  CrossSectionType SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+                                         IndexProcess1, IndexProcess2>::
+      getCrossSection(TParticle const& projectile, [[maybe_unused]] Code const targetId,
+                      [[maybe_unused]] FourMomentum const& targetP4) const {
+
+    if (select_(projectile)) {
+      if constexpr (is_interaction_process_v<process1_type>) {
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TSequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        if constexpr (has_signature_cx1) {
+
+          return A_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else {
+          return A_.getCrossSection(projectile, projectile.getPID(),
+                                    {projectile.getEnergy(), projectile.getMomentum()});
+        }
+      } else if constexpr (process1_type::is_process_sequence) {
+        return A_.getCrossSection(projectile, targetId, targetP4);
+      }
+
+    } else {
+      if constexpr (is_interaction_process_v<process2_type>) {
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<USequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        if constexpr (has_signature_cx1) {
+
+          return B_.getCrossSection(projectile.getPID(), targetId,
+                                    {projectile.getEnergy(), projectile.getMomentum()},
+                                    targetP4);
+        } else {
+          return B_.getCrossSection(projectile, targetId, targetP4);
+        }
+      } else if constexpr (process2_type::is_process_sequence) {
+        return B_.getCrossSection(projectile, targetId, targetP4);
+      }
+    }
+    return CrossSectionType::zero(); // default value
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TSecondaryView, typename TRNG>
+  inline ProcessReturn SwitchProcessSequence<
+      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+      IndexProcess2>::selectInteraction(TSecondaryView& view,
+                                        FourMomentum const& projectileP4,
+                                        NuclearComposition const& composition, TRNG& rng,
+                                        [[maybe_unused]] CrossSectionType const cx_select,
+                                        [[maybe_unused]] CrossSectionType cx_sum) {
+
+    if (select_(view.parent())) {
+      if constexpr (process1_type::is_process_sequence) {
+        // if A_ is a process sequence --> check inside
+        return A_.selectInteraction(view, projectileP4, composition, rng, cx_select,
+                                    cx_sum);
+      } else if constexpr (is_interaction_process_v<process1_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components
+        // for selected process A
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<TSequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                TSequence,                   // process object
+                CrossSectionType,            // return type
+                decltype(projectile) const&, // template argument
+                decltype(projectile) const&, // parameters
+                Code, FourMomentum const&>;
+
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "TSequence has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<TSequence>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /*std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.cbegin(),
+                              weightedCrossSections.cend(), CrossSectionType::zero());
+        } else { // this is for PROPOSAL
+          cx_sum += A_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess1
+            static_assert(
+                has_method_doInteract_v<TSequence,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "TSequence has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<TSequence>. ");
+
+            A_.template doInteraction(view, projectileId, targetId, projectileP4,
+                                      targetP4);
+          } else { // this is for PROPOSAL
+            A_.template doInteraction(view, projectileId, projectileP4);
+          }
+
+          return ProcessReturn::Interacted;
+        } // end collision branch A
+      }
+
+    } else { // selection: end branch A, start branch B
+
+      if constexpr (process2_type::is_process_sequence) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectInteraction(view, projectileP4, composition, rng, cx_select,
+                                    cx_sum);
+      } else if constexpr (is_interaction_process_v<process2_type>) {
+
+        auto const& projectile = view.parent();
+        Code const projectileId = projectile.getPID();
+
+        // get cross section vector for all material components, for selected process B
+        bool constexpr has_signature_cx1 =
+            has_method_getCrossSection_v<USequence,        // process object
+                                         CrossSectionType, // return type
+                                         Code, Code,       // parameters
+                                         FourMomentum const&, FourMomentum const&>;
+        bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+            has_method_getCrossSectionTemplate_v<
+                USequence,                   // process object
+                CrossSectionType,            // return type
+                decltype(projectile) const&, // template argument
+                decltype(projectile) const&, // parameters
+                Code, FourMomentum const&>;
+
+        static_assert((has_signature_cx1 || has_signature_cx2),
+                      "USequence has no method with correct signature \"CrossSectionType "
+                      "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+                      "const&)\" required by "
+                      "InteractionProcess<USequence>. ");
+
+        std::vector<CrossSectionType> weightedCrossSections;
+        if constexpr (has_signature_cx1) {
+          /* std::vector<CrossSectionType> const*/ weightedCrossSections =
+              composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+                FourMomentum const targetP4(
+                    get_mass(targetId),
+                    MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                                   {0_GeV, 0_GeV, 0_GeV}));
+                return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+              });
+
+          cx_sum +=
+              std::accumulate(weightedCrossSections.begin(), weightedCrossSections.end(),
+                              CrossSectionType::zero());
+        } else { // this is for PROPOSAL
+          cx_sum += B_.template getCrossSection(projectile, projectileId, projectileP4);
+        }
+
+        // check if we should execute THIS process and then EXIT
+        if (cx_select < cx_sum) {
+
+          if constexpr (has_signature_cx1) {
+
+            // now also sample targetId from weighted cross sections
+            Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+            FourMomentum const targetP4(
+                get_mass(targetId),
+                MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+                               {0_GeV, 0_GeV, 0_GeV}));
+
+            // interface checking on TProcess2
+            static_assert(
+                has_method_doInteract_v<USequence,       // process object
+                                        void,            // return type
+                                        TSecondaryView,  // template argument
+                                        TSecondaryView&, // method parameters
+                                        Code, Code, FourMomentum const&,
+                                        FourMomentum const&>,
+                "USequence has no method with correct signature \"void "
+                "doInteraction<TSecondaryView>(TSecondaryView&, "
+                "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+                "InteractionProcess<USequence>. ");
+
+            B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+          } else { // this is for PROPOSAL
+            B_.doInteraction(view, projectileId, projectileP4);
+          }
+
+          return ProcessReturn::Interacted;
+        } // end collision in branch B
+      }
+    } // end branch B_
+
+    return ProcessReturn::Ok;
+  } // namespace corsika
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  template <typename TParticle>
+  inline InverseTimeType
+  SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+                        IndexProcess2>::getInverseLifetime(TParticle&& particle) {
+
+    if (select_(particle)) {
+      if constexpr (is_decay_process_v<process1_type> ||
+                    process1_type::is_process_sequence) {
+        return A_.getInverseLifetime(particle);
+      }
+
+    } else {
+
+      if constexpr (is_decay_process_v<process2_type> ||
+                    process2_type::is_process_sequence) {
+        return B_.getInverseLifetime(particle);
+      }
+    }
+    return 0 / second; // default value
+  }
+
+  template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
+            int IndexProcess1, int IndexProcess2>
+  // select decay process
+  template <typename TSecondaryView>
+  inline ProcessReturn SwitchProcessSequence<
+      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+      IndexProcess2>::selectDecay(TSecondaryView& view,
+                                  [[maybe_unused]] InverseTimeType decay_inv_select,
+                                  [[maybe_unused]] InverseTimeType decay_inv_sum) {
+    if (select_(view.parent())) {
+      if constexpr (process1_type::is_process_sequence) {
+        // if A_ is a process sequence --> check inside
+        ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum);
+        // if A_ did succeed, stop routine here (not checking other static branch B_)
+        if (ret != ProcessReturn::Ok) { return ret; }
+      } else if constexpr (is_decay_process_v<process1_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += A_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) {
+          // more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
+
+          // interface checking on TSequence
+          static_assert(has_method_doDecay_v<TSequence, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          A_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch A_
+
+    } else {
+
+      if constexpr (process2_type::is_process_sequence) {
+        // if B_ is a process sequence --> check inside
+        return B_.selectDecay(view, decay_inv_select, decay_inv_sum);
+      } else if constexpr (is_decay_process_v<process2_type>) {
+        // if this is not a ContinuousProcess --> evaluate probability
+        decay_inv_sum += B_.getInverseLifetime(view.parent());
+        // check if we should execute THIS process and then EXIT
+        if (decay_inv_select < decay_inv_sum) {
+
+          // interface checking on TSequence
+          static_assert(has_method_doDecay_v<USequence, void, TSecondaryView&>,
+                        "TDerived has no method with correct signature \"void "
+                        "doDecay(TSecondaryView&)\" required for "
+                        "DecayProcess<TDerived>. ");
+
+          B_.doDecay(view);
+          return ProcessReturn::Decayed;
+        }
+      } // end branch B_
+    }
+    return ProcessReturn::Ok;
+  }
+
+} // namespace corsika

corsika/detail/framework/random/RNGManager.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <iterator>
+#include <random>
+#include <sstream>
+#include <tuple>
+#include <utility>
+
+namespace corsika {
+
+  template <typename CBPRNG>
+  inline void RNGManager<CBPRNG>::registerRandomStream(string_type const& pStreamName) {
+
+    auto const& it = rngs_.find(pStreamName);
+
+    if (it == rngs_.end()) // key not in container, so create one and initialize the value
+      rngs_.emplace(std::piecewise_construct, std::forward_as_tuple(pStreamName.c_str()),
+                    std::forward_as_tuple(seed_, uint32_t(rngs_.size())));
+  }
+
+  template <typename CBPRNG>
+  inline typename RNGManager<CBPRNG>::prng_type& RNGManager<CBPRNG>::getRandomStream(
+      string_type const& pStreamName) {
+
+    if (isRegistered(pStreamName)) {
+      return rngs_.at(pStreamName);
+    } else { // this stream name is not in the map
+      throw std::runtime_error("'" + pStreamName + "' is not a registered stream.");
+    }
+  }
+
+  template <typename CBPRNG>
+  inline bool RNGManager<CBPRNG>::isRegistered(string_type const& pStreamName) const {
+    return rngs_.count(pStreamName) > 0;
+  }
+
+  template <typename CBPRNG>
+  inline std::stringstream RNGManager<CBPRNG>::dumpState() const {
+    std::stringstream buffer;
+    for (auto const& [streamName, rng] : rngs_) {
+      buffer << '"' << streamName << "\" = \"" << rng << '"' << std::endl;
+    }
+
+    return buffer;
+  }
+
+} // namespace corsika

corsika/detail/framework/random/random_iterator/RandomIterator.h

+/*----------------------------------------------------------------------------
+*
+* Copyright (C) 2021 - 2024 Antonio Augusto Alves Junior
+*
+* This file is part of RandomIterator.
+
+* Redistribution and use in source and binary forms, with or without modification,
+* are permitted provided that the following conditions are met:
+*
+* 1. Redistributions of source code must retain the above copyright notice,
+* this list of conditions and the following disclaimer.
+*
+* 2. Redistributions in binary form must reproduce the above copyright notice,
+* this list of conditions and the following disclaimer in the documentation and/or
+* other materials provided with the distribution.
+*
+* 3. Neither the name of the copyright holder nor the names of its contributors
+* may be used to endorse or promote products derived from this software without
+* specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND
+* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+* INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+* STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+* OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+* ----------------------------------------------------------------------------*/
+
+/*
+ * RandomIterator.h
+ *
+ *  Created on: 21/02/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+//  This is the only RandomIterator header that is guaranteed to
+//  change with every RandomIterator release.
+//
+//  RandomIterator_VERSION % 100 is the sub-minor version
+//  RandomIterator_VERSION / 100 % 1000 is the minor version
+//  RandomIterator_VERSION / 100000 is the major version
+//
+//  Because this header does not #include <RandomIterator/detail/Config.h>,
+//  it is the only RandomIterator header that does not cause
+//  RandomIterator_HOST_SYSTEM and RandomIterator_DEVICE_SYSTEM to be defined.
+
+/*! \def RandomIterator_VERSION
+ *  \brief The preprocessor macro \p RandomIterator_VERSION encodes the version
+ *         number of the RandomIterator.
+ *
+ *         <tt>RandomIterator_VERSION % 100</tt> is the patch version.
+ *         <tt>RandomIterator_VERSION / 100 % 1000</tt> is the feature version.
+ *         <tt>RandomIterator_VERSION / 100000</tt> is the major version.
+ */
+#define RandomIterator_VERSION 100001
+
+
+/*! \def RandomIterator_MAJOR_VERSION
+ *  \brief The preprocessor macro \p RandomIterator_MAJOR_VERSION encodes the
+ *         major version number of RandomIterator.
+ */
+#define RandomIterator_MAJOR_VERSION     (RandomIterator_VERSION / 100000)
+
+/*! \def RandomIterator_MINOR_VERSION
+ *  \brief The preprocessor macro \p RandomIterator_MINOR_VERSION encodes the
+ *         minor version number of RandomIterator.
+ */
+#define RandomIterator_MINOR_VERSION     (RandomIterator_VERSION / 100 % 1000)
+
+/*! \def RandomIterator_PATCH_NUMBER
+ *  \brief The preprocessor macro \p RandomIterator_PATCH_NUMBER encodes the
+ *         patch number of the RandomIterator library.
+ */
+#define RandomIterator_PATCH_NUMBER  0
+
+
+
+// Declare these namespaces here for the purpose of Doxygenating them
+
+/*!
+ *  \brief \p random_iterator is the top-level namespace which contains all RandomIterator
+ *         functions and types.
+ */
+namespace random_iterator{ }
+
+

corsika/detail/framework/random/random_iterator/Stream.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/*
+ * Stream.hpp
+ *
+ *  Created on: 23/02/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <stdint.h>
+
+#include "detail/tbb/iterators.h"
+#include "detail/Engine.hpp"
+#include "detail/functors/EngineCaller.hpp"
+
+namespace random_iterator {
+
+  /// philox engine
+  typedef detail::philox philox;
+  /// threefry  engine
+  typedef detail::threefry threefry;
+  /// squares3_128  engine
+  typedef detail::squares3_128 squares3_128;
+  /// squares4_128  engine
+  typedef detail::squares4_128 squares4_128;
+  /// squares3_64  engine
+  typedef detail::squares3_64 squares3_64;
+  /// squares4_64  engine
+  typedef detail::squares4_64 squares4_64;
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+  /// ars  engine
+  typedef detail::ars ars;
+#endif
+
+  /**
+   * \class Stream
+   *
+   * \brief representation for streams of pseudorandom numbers
+   *
+   * Objects of this class represent collections of up to \f$ {2}^{32} \f$  streams of
+   * random numbers for each seed. Each of such streams has a length of \f$ {2}^{64} \f$.
+   * If the chosen distribution generates data with 64bits, like ``double``, then each
+   * stream can handle up to 128 ExaByte of data. Stream objects are thread-safe and
+   * lazy-evaluated. Output is produced only when requested. It is user responsibility to
+   * retrieve and store, if needed, the generated output.
+   *
+   * \tparam DistibutionType STL compliant distribution type. For example:
+   * ```std::uniform_real_distribution```, ```std::exponential_distribution```.
+   *
+   * \tparam EngineType RNG type.
+   *
+   * There are seven options:
+   *
+   *  * random_iterator::philox
+   *  * random_iterator::threefry
+   *  * random_iterator::ars
+   *  * random_iterator::squares3_128
+   *  * random_iterator::squares4_128
+   *  * random_iterator::squares3_64
+   *  * random_iterator::squares4_64
+   *
+   *  Usage:
+   *
+   * 	1) Direct class instantiation:
+   *
+   *  \code
+   *  \\this code will print the 10 doubles randomly distributed in the range [0.0, 1.0].
+   *
+   *  \\instantiate real uniform distribution from C++ standard libray
+   *  std::uniform_real_distribution<double> udist(0.0, 1.0);
+   *
+   *  \\
+   *  Stream<std::uniform_real_distribution<double>, random_iterator::philox>
+   * uniform_stream(udist, 0x1a2b3c4d, 1 );
+   *
+   *  for( size_t i=0; i<10; ++i )
+   *   std::cout << uniform_stream[i] << std::endl;
+   *  \endcode
+   *
+   * 	2) Using ```random_iterator::make_stream(...)```:
+   *
+   * 	\code
+   *  \\this code will print the 10 doubles randomly distributed in the range [0.0, 1.0].
+   *
+   *  \\instantiate real uniform distribution from C++ standard libray
+   *  std::uniform_real_distribution<double> udist(0.0, 1.0);
+   *
+   *  \\instantiate the generator
+   *  random_iterator::philox   rng(0x1a2b3c4d);
+   *
+   *  \\create stream
+   *  auto uniform_stream = random_iterator::make_stream(udist, rng, 1 );
+   *
+   *  for( size_t i=0; i<10; ++i )
+   *   std::cout << uniform_stream[i] << std::endl;
+   * 	\endcode
+   */
+  template <typename DistibutionType, typename EngineType>
+  class Stream {
+    typedef detail::EngineCaller<DistibutionType, EngineType> caller_type; //!
+    typedef random_iterator_tbb::counting_iterator<typename EngineType::advance_type>
+        counting_iterator; //!
+
+  public:
+    typedef EngineType engine_type;                          //! Engine type
+    typedef typename engine_type::state_type state_type;     //! Type of RNG state
+    typedef typename engine_type::seed_type seed_type;       //! Type of RNG seed
+    typedef typename engine_type::advance_type advance_type; //! Type of RNG displacement
+    typedef typename engine_type::init_type init_type;       //! Type of RNG initializer
+
+    typedef DistibutionType distribution_type; //! type of the STL compliant distribution
+    typedef typename distribution_type::result_type result_type; //! type of the result
+
+    typedef random_iterator_tbb::transform_iterator<caller_type, counting_iterator>
+        iterator; //! type of stream iterator
+
+    Stream() = delete;
+
+    /**
+     * Constructor.
+     *
+     * @param distribution Distribution with stl compliant interface.
+     * @param seed Seed for pseudorandom number generation.
+     * @param stream Pseudorandom number stream.
+     */
+    Stream(distribution_type const& distribution, seed_type seed, uint32_t stream = 0)
+        : distribution_(distribution)
+        , engine_(seed, stream)
+        , seed_(seed)
+        , stream_(stream) {}
+
+    /**
+     * Copy constructor
+     *
+     * @param other
+     */
+    Stream(Stream<DistibutionType, EngineType> const& other)
+        : distribution_(other.getDistribution())
+        , engine_(other.getEngine())
+        , seed_(other.getSeed())
+        , stream_(other.getStream()) {}
+
+    /**
+     * Returns an iterator to the first element of the stream.
+     * @return  iterator to the first element of the stream.
+     */
+    inline iterator begin() const {
+
+      counting_iterator first(0);
+      caller_type caller(distribution_, seed_, stream_);
+
+      return iterator(first, caller);
+    }
+
+    /**
+     * Returns an iterator to the past last element of the stream.
+     * @return iterator to the past last element of the stream.
+     */
+    inline iterator end() const {
+
+      counting_iterator last(std::numeric_limits<advance_type>::max());
+      caller_type caller(distribution_, seed_, stream_);
+
+      return iterator(last, caller);
+    }
+
+    /**
+     * Returns the element at specified location n. No bounds checking is performed.
+     *
+     * @param n position in the sequence.
+     * @return element at location n
+     */
+    inline result_type operator[](advance_type n) const {
+
+      caller_type caller(distribution_, seed_, stream_);
+
+      return caller(n);
+    }
+
+    /**
+     *  At each call, this function will produce a pseudorandom number and advance the
+     * engine state.
+     *
+     * @return pseudorandom number distributed according with `DistributionType`
+     */
+    inline result_type operator()(void) { return distribution_(engine_); }
+
+    /**
+     * Get the distribution
+     *
+     * @return Copy of the object's distribution.
+     */
+    inline distribution_type getDistribution() const { return distribution_; }
+
+    /**
+     * Set the distribution
+     * @param dist
+     *
+     */
+    inline void setDistribution(distribution_type dist) { distribution_ = dist; }
+
+    /**
+     * Get the associated engine
+     *
+     * @return Copy of the object's engine
+     */
+    inline engine_type getEngine() const { return engine_; }
+
+    /**
+     * Set the object's engine
+     *
+     * @param engine
+     */
+    inline void setEngine(engine_type engine) { engine_ = engine; }
+
+    /**
+     * Get the seed
+     *
+     * @return Copy of the object's seed
+     */
+    inline const seed_type& getSeed() const { return seed_; }
+
+    /**
+     * Set the seed
+     *
+     * @param seed
+     */
+    inline void setSeed(const seed_type& seed) { seed_ = seed; }
+
+    /**
+     * Get the object's stream number
+     *
+     * @return stream number
+     */
+    inline uint32_t getStream() const { return stream_; }
+
+    /**
+     * Set the object's stream number
+     *
+     * @param stream
+     */
+    inline void setStream(uint32_t stream) { stream_ = stream; }
+
+    friend inline std::ostream& operator<<(
+        std::ostream& os, const Stream<DistibutionType, EngineType>& be) {
+      return os << "engine: " << be.getEngine() << " stream: " << be.getStream()
+                << " distribution: " << be.getDistribution();
+    }
+
+  private:
+    distribution_type distribution_;
+    engine_type engine_;
+    seed_type seed_;
+    uint32_t stream_;
+  };
+
+  /**
+   * Stream
+   *
+   * @param dist STL compliant distribution instance.
+   * @param eng Engine instance
+   * @param stream Stream number, in the range [0, 2^{32} - 1]
+   * @return Stream object.
+   */
+  template <typename DistibutionType, typename EngineType>
+  Stream<DistibutionType, EngineType> make_stream(DistibutionType const& dist,
+                                                  EngineType eng, uint32_t stream) {
+
+    return Stream<DistibutionType, EngineType>(dist, eng.getSeed(), stream);
+  }
+
+} // namespace random_iterator

corsika/detail/framework/random/random_iterator/detail/Engine.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/*
+ * Engine.hpp
+ *
+ *  Created on: 22 de fev. de 2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <limits>
+
+#include "EngineTraits.hpp"
+#include "SplitMix.hpp"
+#include "Squares3_128.hpp"
+#include "Squares4_128.hpp"
+#include "Squares3_64.hpp"
+#include "Squares4_64.hpp"
+
+namespace random_iterator {
+
+  namespace detail {
+
+    template <typename Engine123>
+    class Engine {
+
+      typedef unsigned trigger_type;
+
+    public:
+      typedef Engine123 engine_type;
+      typedef typename detail::random_traits<engine_type>::state_type state_type;
+      typedef typename detail::random_traits<engine_type>::seed_type seed_type;
+      typedef typename detail::random_traits<engine_type>::advance_type advance_type;
+      typedef typename detail::random_traits<engine_type>::init_type init_type;
+      typedef typename detail::random_traits<engine_type>::result_type result_type;
+
+      static const unsigned arity = detail::random_traits<engine_type>::arity;
+
+      Engine() = delete;
+
+      Engine(result_type seed)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed_type{})
+          , trigger_(arity) {
+        init_type temp{{}};
+
+        for (unsigned i = 0; i < temp.size(); ++i) {
+          temp[i] = detail::splitmix<result_type>(seed);
+        }
+
+        seed_ = temp;
+      }
+
+      Engine(result_type seed, uint32_t stream)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed_type{})
+          , trigger_(arity) {
+        init_type temp{{}};
+
+        for (unsigned i = 0; i < temp.size(); ++i) {
+          temp[i] = detail::splitmix<result_type>(seed);
+        }
+
+        state_[arity - 1] = stream;
+        state_[arity - 1] = state_[arity - 1] << 32;
+
+        seed_ = temp;
+      }
+
+      Engine(init_type seed)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed)
+          , trigger_(arity) {}
+
+      Engine(init_type seed, uint32_t stream)
+          : engine_(engine_type{})
+          , cache_(state_type{})
+          , state_(state_type{})
+          , seed_(seed)
+          , trigger_(arity) {
+        state_[arity - 1] = stream;
+        state_[arity - 1] << 32;
+      }
+
+      Engine(Engine<Engine123> const& other)
+          : engine_(engine_type{})
+          , cache_(other.getCache())
+          , state_(other.getState())
+          , seed_(other.getSeed())
+          , trigger_(other.getTrigger()) {}
+
+      inline Engine<Engine123>& operator=(Engine<Engine123> const& other) {
+
+        if (this == &other) return *this;
+
+        engine_ = engine_type{};
+        cache_ = other.getCache();
+        state_ = other.getState();
+        seed_ = other.getSeed();
+        trigger_ = other.getTrigger();
+        return *this;
+      }
+
+      inline result_type operator()(void) {
+        result_type result = 0;
+
+        if (trigger_ == arity) {
+
+          trigger_ = 0;
+          cache_ = engine_(state_.incr(), seed_);
+          result = cache_[trigger_];
+          ++trigger_;
+        } else {
+          result = cache_[trigger_];
+          ++trigger_;
+        }
+
+        return result;
+      }
+
+      inline void discard(advance_type n) {
+
+        state_.incr(n);
+        trigger_ = arity;
+      }
+
+      inline void reset(void) {
+        trigger_ = arity;
+        state_ = state_type{};
+        cache_ = state_type{};
+      }
+
+      inline const seed_type& getSeed() const { return seed_; }
+
+      inline void setSeed(seed_type seed) { seed_ = seed; }
+
+      inline void setSeed(result_type seed) {
+
+        init_type temp{{}};
+        for (unsigned i = 0; i < temp.size(); ++i) {
+          temp[i] = detail::splitmix<result_type>(seed);
+        }
+        seed_ = temp;
+      }
+
+      inline const state_type& getState() const { return state_; }
+
+      inline void setState(const state_type& state) { state_ = state; }
+
+      inline trigger_type getTrigger() const { return trigger_; }
+
+      inline void setTrigger(trigger_type trigger) { trigger_ = trigger; }
+      static constexpr result_type min() { return 0; }
+
+      static constexpr result_type max() {
+        return std::numeric_limits<result_type>::max();
+      }
+
+      friend inline std::ostream& operator<<(std::ostream& os,
+                                             const Engine<Engine123>& be) {
+        return os << "state: " << be.getState() << " seed: " << be.getSeed()
+                  << " trigger: " << be.getTrigger();
+      }
+
+    private:
+      inline state_type getCache() const { return cache_; }
+
+      engine_type engine_;
+      state_type cache_;
+      state_type state_;
+      seed_type seed_;
+      trigger_type trigger_;
+    };
+
+    typedef Engine<random_iterator_r123::Philox2x64> philox;
+    typedef Engine<random_iterator_r123::Threefry2x64> threefry;
+    typedef detail::Squares3_64 squares3_64;
+    typedef detail::Squares4_64 squares4_64;
+    typedef detail::Squares3_128 squares3_128;
+    typedef detail::Squares4_128 squares4_128;
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+    typedef Engine<random_iterator_r123::ARS2x64> ars;
+#endif
+
+  } // namespace detail
+
+} // namespace random_iterator

corsika/detail/framework/random/random_iterator/detail/EngineTraits.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/*
+ * EngineTraits.hpp
+ *
+ *  Created on: 23 de fev. de 2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <stdint.h>
+#include "Random123/array.h"
+#include "Random123/philox.h"
+#include "Random123/threefry.h"
+#include "Random123/ars.h"
+#include "Random123/ReinterpretCtr.hpp"
+
+namespace random_iterator {
+
+  namespace detail {
+
+    template <typename Engine>
+    struct random_traits;
+    /*
+     * random_traits<T>::state_type { counter, state}
+     * random_traits<T>::advance_type;
+     * random_traits<T>::init_type;
+     * random_traits<T>::result_type;
+     */
+
+    // philox
+    template <>
+    struct random_traits<random_iterator_r123::Philox2x64> {
+      typedef typename random_iterator_r123::Philox2x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Philox2x64::key_type seed_type;
+      typedef typename random_iterator_r123::Philox2x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 2 };
+    };
+
+    template <>
+    struct random_traits<random_iterator_r123::Philox4x64> {
+      typedef typename random_iterator_r123::Philox4x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Philox4x64::key_type seed_type;
+      typedef typename random_iterator_r123::Philox4x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 4 };
+    };
+    //
+    template <>
+    struct random_traits<random_iterator_r123::Threefry2x64> {
+      typedef typename random_iterator_r123::Threefry2x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Threefry2x64::key_type seed_type;
+      typedef typename random_iterator_r123::Threefry2x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 2 };
+    };
+
+    //
+    template <>
+    struct random_traits<random_iterator_r123::Threefry4x64> {
+      typedef typename random_iterator_r123::Threefry4x64::ctr_type state_type;
+      typedef typename random_iterator_r123::Threefry4x64::key_type seed_type;
+      typedef typename random_iterator_r123::Threefry4x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 4 };
+    };
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+    template <>
+    struct random_traits<random_iterator_r123::ARS4x32> {
+      typedef typename random_iterator_r123::ARS4x32::ctr_type state_type;
+      typedef typename random_iterator_r123::ARS4x32::key_type seed_type;
+      typedef typename random_iterator_r123::ARS4x32::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 4 };
+    };
+
+    template <>
+    struct random_traits<random_iterator_r123::ARS2x64> {
+      typedef typename random_iterator_r123::ARS2x64::ctr_type state_type;
+      typedef typename random_iterator_r123::ARS2x64::key_type seed_type;
+      typedef typename random_iterator_r123::ARS2x64::ukey_type init_type;
+      typedef uint64_t advance_type;
+      typedef state_type::value_type result_type;
+
+      enum { arity = 2 };
+    };
+
+#endif
+
+  } // namespace detail
+
+} // namespace random_iterator

corsika/detail/framework/random/random_iterator/detail/Macros.h

+/*----------------------------------------------------------------------------
+*
+* Copyright (C) 2021 - 2024 Antonio Augusto Alves Junior
+*
+* This file is part of RandomIterator.
+
+* Redistribution and use in source and binary forms, with or without modification,
+* are permitted provided that the following conditions are met:
+*
+* 1. Redistributions of source code must retain the above copyright notice,
+* this list of conditions and the following disclaimer.
+*
+* 2. Redistributions in binary form must reproduce the above copyright notice,
+* this list of conditions and the following disclaimer in the documentation and/or
+* other materials provided with the distribution.
+*
+* 3. Neither the name of the copyright holder nor the names of its contributors
+* may be used to endorse or promote products derived from this software without
+* specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND
+* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+* INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+* STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+* OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+* ----------------------------------------------------------------------------*/
+/*
+ * Macros.h
+ *
+ *  Created on: 04/03/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#if defined(__has_builtin)
+    //__has_builtin(__builtin_expect)
+	#if __has_builtin(__builtin_expect)
+		 #define   _LIKELY_(x) __builtin_expect(x, 1)
+		 #define _UNLIKELY_(x) __builtin_expect(x, 0)
+	#else
+		 #define   _LIKELY_(x) x
+		 #define _UNLIKELY_(x) x
+	#endif //__has_builtin(__builtin_expect)
+#else
+	#define   _LIKELY_(x) x
+	#define _UNLIKELY_(x) x
+#endif //defined(__has_builtin)
+
+
+//detect endianess
+#if defined(__BYTE_ORDER) && __BYTE_ORDER == __BIG_ENDIAN || \
+    defined(__BIG_ENDIAN__) ||                               \
+    defined(__ARMEB__) ||                                    \
+    defined(__THUMBEB__) ||                                  \
+    defined(__AARCH64EB__) ||                                \
+    defined(_MIBSEB) || defined(__MIBSEB) || defined(__MIBSEB__)
+#ifndef __BIG_ENDIAN__
+#define __BIG_ENDIAN__
+#endif
+#elif defined(__BYTE_ORDER) && __BYTE_ORDER == __LITTLE_ENDIAN || \
+    defined(__LITTLE_ENDIAN__) ||                                 \
+    defined(__ARMEL__) ||                                         \
+    defined(__THUMBEL__) ||                                       \
+    defined(__AARCH64EL__) ||                                     \
+    defined(_MIPSEL) || defined(__MIPSEL) || defined(__MIPSEL__) || \
+    defined(_WIN32) || defined(__i386__) || defined(__x86_64__) || \
+    defined(_X86_) || defined(_IA64_)
+#ifndef __LITTLE_ENDIAN__
+#define __LITTLE_ENDIAN__
+#endif
+#else
+#error "I don't know what architecture this is!"
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/MicroURNG.hpp

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __MicroURNG_dot_hpp__
+#define __MicroURNG_dot_hpp__
+
+#include <stdexcept>
+#include <limits>
+
+namespace random_iterator_r123 {
+  /**
+      Given a CBRNG whose ctr_type has an unsigned integral value_type,
+      MicroURNG<CBRNG>(c, k) is a type that satisfies the
+      requirements of a C++11 Uniform Random Number Generator.
+
+      The intended purpose is for a MicroURNG to be passed
+      as an argument to a C++11 Distribution, e.g.,
+      std::normal_distribution.  See examples/MicroURNG.cpp.
+
+      The MicroURNG functor has a period of "only"
+
+         ctr_type.size()*2^32,
+
+      after which it will silently repeat.
+
+      The high 32 bits of the highest word in the counter c, passed to
+      the constructor must be zero.  MicroURNG uses these bits to
+      "count".
+
+      Older versions of the library permitted a second template
+      parameter by which the caller could control the number of
+      bits devoted to the URNG's internal counter.  This flexibility
+      has been disabled because URNGs created with different
+      numbers of counter bits could, conceivably "collide".
+
+  \code
+         typedef ?someCBRNG? RNG;
+         RNG::ctr_type c = ...; // under application control
+         RNG::key_type k = ...; //
+         std::normal_distribution<float> nd;
+         MicroURNG<RNG> urng(c, k);
+         for(???){
+           ...
+           nd(urng);  // may be called several hundred times with BITS=10
+           ...
+         }
+  \endcode
+  */
+
+  template <typename CBRNG>
+  class MicroURNG {
+    // According to C++11, a URNG requires only a result_type,
+    // operator()(), min() and max() methods.  Everything else
+    // (ctr_type, key_type, reset() method, etc.) is "value added"
+    // for the benefit of users that "know" that they're dealing with
+    // a MicroURNG.
+  public:
+    typedef CBRNG cbrng_type;
+    static const int BITS = 32;
+    typedef typename cbrng_type::ctr_type ctr_type;
+    typedef typename cbrng_type::key_type key_type;
+    typedef typename cbrng_type::ukey_type ukey_type;
+    typedef typename ctr_type::value_type result_type;
+
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(std::numeric_limits<result_type>::digits >= BITS,
+                                       "The result_type must have at least 32 bits");
+
+    result_type operator()() {
+      if (last_elem == 0) {
+        // jam n into the high bits of c
+        const size_t W = std::numeric_limits<result_type>::digits;
+        ctr_type c = c0;
+        c[c0.size() - 1] |= n << (W - BITS);
+        rdata = b(c, k);
+        n++;
+        last_elem = rdata.size();
+      }
+      return rdata[--last_elem];
+    }
+    MicroURNG(cbrng_type _b, ctr_type _c0, ukey_type _uk)
+        : b(_b)
+        , c0(_c0)
+        , k(_uk)
+        , n(0)
+        , last_elem(0) {
+      chkhighbits();
+    }
+    MicroURNG(ctr_type _c0, ukey_type _uk)
+        : b()
+        , c0(_c0)
+        , k(_uk)
+        , n(0)
+        , last_elem(0) {
+      chkhighbits();
+    }
+
+    // _Min and _Max work around a bug in the library shipped with MacOS Xcode 4.5.2.
+    // See the commment in conventional/Engine.hpp.
+    const static result_type _Min = 0;
+    const static result_type _Max = ~((result_type)0);
+
+    static RANDOM_ITERATOR_R123_CONSTEXPR result_type min
+    RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
+      return _Min;
+    }
+    static RANDOM_ITERATOR_R123_CONSTEXPR result_type max
+    RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
+      return _Max;
+    }
+    // extra methods:
+    const ctr_type& counter() const { return c0; }
+    void reset(ctr_type _c0, ukey_type _uk) {
+      c0 = _c0;
+      chkhighbits();
+      k = _uk;
+      n = 0;
+      last_elem = 0;
+    }
+
+  private:
+    cbrng_type b;
+    ctr_type c0;
+    key_type k;
+    RANDOM_ITERATOR_R123_ULONG_LONG n;
+    size_t last_elem;
+    ctr_type rdata;
+    void chkhighbits() {
+      result_type r = c0[c0.size() - 1];
+      result_type mask = ((uint64_t)std::numeric_limits<result_type>::max
+                              RANDOM_ITERATOR_R123_NO_MACRO_SUBST()) >>
+                         BITS;
+      if ((r & mask) != r)
+        throw std::runtime_error("MicroURNG: c0, does not have high bits clear");
+    }
+  };
+} // namespace random_iterator_r123
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/ReinterpretCtr.hpp

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __ReinterpretCtr_dot_hpp__
+#define __ReinterpretCtr_dot_hpp__
+
+#include "features/compilerfeatures.h"
+#include <cstring>
+
+namespace random_iterator_r123 {
+  /*!
+    ReinterpretCtr uses memcpy to map back and forth
+    between a CBRNG's ctr_type and the specified ToType.  For example,
+    after:
+
+      typedef ReinterpretCtr<r123array4x32, Philox2x64> G;
+
+    G is a bona fide CBRNG with ctr_type r123array4x32.
+
+    WARNING:  ReinterpretCtr is endian dependent.  The
+    values returned by G, declared as above,
+    will depend on the endianness of the machine on which it runs.
+   */
+
+  template <typename ToType, typename CBRNG>
+  struct ReinterpretCtr {
+    typedef ToType ctr_type;
+    typedef typename CBRNG::key_type key_type;
+    typedef typename CBRNG::ctr_type bctype;
+    typedef typename CBRNG::ukey_type ukey_type;
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(
+        sizeof(ToType) == sizeof(bctype) && sizeof(typename bctype::value_type) != 16,
+        "ReinterpretCtr:  sizeof(ToType) is not the same as sizeof(CBRNG::ctr_type) or "
+        "CBRNG::ctr_type::value_type looks like it might be __m128i");
+    // It's amazingly difficult to safely do conversions with __m128i.
+    // If we use the operator() implementation below with a CBRNG
+    // whose ctr_type is r123array1xm128i, gcc4.6 optimizes away the
+    // memcpys, inlines the operator()(c,k), and produces assembly
+    // language that ends with an aesenclast instruction with a
+    // destination operand pointing to an unaligned memory address ...
+    // Segfault!  See:  http://gcc.gnu.org/bugzilla/show_bug.cgi?id=50444
+    // MSVC also produces code that crashes.  We suspect a
+    // similar mechanism but haven't done the debugging necessary to
+    // be sure.  We were able to 'fix' gcc4.6 by making bc a mutable
+    // data member rather than declaring it in the scope of
+    // operator().  That didn't fix the MSVC problems, though.
+    //
+    // Conclusion - don't touch __m128i, at least for now.  The
+    // easiest (but highly imprecise) way to do that is the static
+    // assertion above that rejects bctype::value_types of size 16. -
+    // Sep 2011.
+    ctr_type operator()(ctr_type c, key_type k) {
+      bctype bc;
+      std::memcpy(&bc, &c, sizeof(c));
+      CBRNG b;
+      bc = b(bc, k);
+      std::memcpy(&c, &bc, sizeof(bc));
+      return c;
+    }
+  };
+} // namespace random_iterator_r123
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/aes.h

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __Random123_aes_dot_hpp__
+#define __Random123_aes_dot_hpp__
+
+#include "features/compilerfeatures.h"
+#include "array.h"
+
+/* Implement a bona fide AES block cipher.  It's minimally
+// checked against the test vector in FIPS-197 in ut_aes.cpp. */
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+
+/** @ingroup AESNI */
+typedef struct r123array1xm128i aesni1xm128i_ctr_t;
+/** @ingroup AESNI */
+typedef struct r123array1xm128i aesni1xm128i_ukey_t;
+/** @ingroup AESNI */
+typedef struct r123array4x32 aesni4x32_ukey_t;
+/** @ingroup AESNI */
+enum r123_enum_aesni1xm128i { aesni1xm128i_rounds = 10 };
+
+/** \cond HIDDEN_FROM_DOXYGEN */
+RANDOM_ITERATOR_R123_STATIC_INLINE __m128i AES_128_ASSIST (__m128i temp1, __m128i temp2) { 
+    __m128i temp3; 
+    temp2 = _mm_shuffle_epi32 (temp2 ,0xff); 
+    temp3 = _mm_slli_si128 (temp1, 0x4);
+    temp1 = _mm_xor_si128 (temp1, temp3);
+    temp3 = _mm_slli_si128 (temp3, 0x4);
+    temp1 = _mm_xor_si128 (temp1, temp3);
+    temp3 = _mm_slli_si128 (temp3, 0x4);
+    temp1 = _mm_xor_si128 (temp1, temp3);
+    temp1 = _mm_xor_si128 (temp1, temp2); 
+    return temp1; 
+}
+
+RANDOM_ITERATOR_R123_STATIC_INLINE void aesni1xm128iexpand(aesni1xm128i_ukey_t uk, __m128i ret[11])
+{
+    __m128i rkey = uk.v[0].m;
+    __m128i tmp2;
+
+    ret[0] = rkey;
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x1);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[1] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x2);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[2] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x4);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[3] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x8);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[4] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x10);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[5] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x20);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[6] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x40);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[7] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x80);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[8] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x1b);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[9] = rkey;
+
+    tmp2 = _mm_aeskeygenassist_si128(rkey, 0x36);
+    rkey = AES_128_ASSIST(rkey, tmp2);
+    ret[10] = rkey;
+}
+/** \endcond */
+    
+#ifdef __cplusplus
+/** @ingroup AESNI */
+struct aesni1xm128i_key_t{ 
+    __m128i k[11]; 
+    aesni1xm128i_key_t(){
+        aesni1xm128i_ukey_t uk;
+        uk.v[0].m = _mm_setzero_si128();
+        aesni1xm128iexpand(uk, k);
+    }
+    aesni1xm128i_key_t(const aesni1xm128i_ukey_t& uk){
+        aesni1xm128iexpand(uk, k);
+    }
+    aesni1xm128i_key_t(const aesni4x32_ukey_t& uk){
+        aesni1xm128i_ukey_t uk128;
+        uk128.v[0].m = _mm_set_epi32(uk.v[3], uk.v[2], uk.v[1], uk.v[0]);
+        aesni1xm128iexpand(uk128, k);
+    }
+    aesni1xm128i_key_t& operator=(const aesni1xm128i_ukey_t& uk){
+        aesni1xm128iexpand(uk, k);
+        return *this;
+    }
+    aesni1xm128i_key_t& operator=(const aesni4x32_ukey_t& uk){
+        aesni1xm128i_ukey_t uk128;
+        uk128.v[0].m = _mm_set_epi32(uk.v[3], uk.v[2], uk.v[1], uk.v[0]);
+        aesni1xm128iexpand(uk128, k);
+        return *this;
+    }
+    bool operator==(const aesni1xm128i_key_t& rhs) const{
+        for(int i=0; i<11; ++i){
+            // Sigh... No r123m128i(__m128i) constructor!
+            r123m128i li; li.m = k[i];
+            r123m128i ri; ri.m = rhs.k[i];
+            if( li != ri ) return false;
+        }
+        return true;
+    }
+    bool operator!=(const aesni1xm128i_key_t& rhs) const{
+        return !(*this == rhs);
+    }
+    friend std::ostream& operator<<(std::ostream& os, const aesni1xm128i_key_t& v){
+        r123m128i ki;
+        for(int i=0; i<10; ++i){
+            ki.m = v.k[i];
+            os << ki << " ";
+        }
+        ki.m = v.k[10];
+        return os << ki;
+    }
+    friend std::istream& operator>>(std::istream& is, aesni1xm128i_key_t& v){
+        r123m128i ki;
+        for(int i=0; i<11; ++i){
+            is >> ki;
+            v.k[i] = ki;
+        }
+        return is;
+    }
+};
+#else
+typedef struct { 
+    __m128i k[11]; 
+}aesni1xm128i_key_t;
+
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni1xm128i_key_t aesni1xm128ikeyinit(aesni1xm128i_ukey_t uk){
+    aesni1xm128i_key_t ret;
+    aesni1xm128iexpand(uk, ret.k);
+    return ret;
+}
+#endif
+
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni1xm128i_ctr_t aesni1xm128i(aesni1xm128i_ctr_t in, aesni1xm128i_key_t k) {
+    __m128i x = _mm_xor_si128(k.k[0], in.v[0].m);
+    x = _mm_aesenc_si128(x, k.k[1]);
+    x = _mm_aesenc_si128(x, k.k[2]);
+    x = _mm_aesenc_si128(x, k.k[3]);
+    x = _mm_aesenc_si128(x, k.k[4]);
+    x = _mm_aesenc_si128(x, k.k[5]);
+    x = _mm_aesenc_si128(x, k.k[6]);
+    x = _mm_aesenc_si128(x, k.k[7]);
+    x = _mm_aesenc_si128(x, k.k[8]);
+    x = _mm_aesenc_si128(x, k.k[9]);
+    x = _mm_aesenclast_si128(x, k.k[10]);
+    {
+      aesni1xm128i_ctr_t ret;
+      ret.v[0].m = x;
+      return ret;
+    }
+}
+
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni1xm128i_ctr_t aesni1xm128i_R(unsigned R, aesni1xm128i_ctr_t in, aesni1xm128i_key_t k){
+    RANDOM_ITERATOR_R123_ASSERT(R==10);
+    return aesni1xm128i(in, k);
+}
+
+
+/** @ingroup AESNI */
+typedef struct r123array4x32 aesni4x32_ctr_t;
+/** @ingroup AESNI */
+typedef aesni1xm128i_key_t aesni4x32_key_t;
+/** @ingroup AESNI */
+enum r123_enum_aesni4x32 { aesni4x32_rounds = 10 };
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni4x32_key_t aesni4x32keyinit(aesni4x32_ukey_t uk){
+    aesni1xm128i_ukey_t uk128;
+    aesni4x32_key_t ret;
+    uk128.v[0].m = _mm_set_epi32(uk.v[3], uk.v[2], uk.v[1], uk.v[0]);
+    aesni1xm128iexpand(uk128, ret.k);
+    return ret;
+}
+
+/** @ingroup AESNI */
+/** The aesni4x32_R function provides a C API to the @ref AESNI "AESNI" CBRNG, allowing the number of rounds to be specified explicitly **/
+RANDOM_ITERATOR_R123_STATIC_INLINE aesni4x32_ctr_t aesni4x32_R(unsigned int Nrounds, aesni4x32_ctr_t c, aesni4x32_key_t k){
+    aesni1xm128i_ctr_t c128;
+    c128.v[0].m = _mm_set_epi32(c.v[3], c.v[2], c.v[1], c.v[0]);
+    c128 = aesni1xm128i_R(Nrounds, c128, k);
+    _mm_storeu_si128((__m128i*)&c.v[0], c128.v[0].m);
+    return c;
+}
+
+#define aesni4x32_rounds aesni1xm128i_rounds
+
+/** The aesni4x32 macro provides a C API to the @ref AESNI "AESNI" CBRNG, uses the default number of rounds i.e. \c aesni4x32_rounds **/
+/** @ingroup AESNI */
+#define aesni4x32(c,k) aesni4x32_R(aesni4x32_rounds, c, k)
+
+#ifdef __cplusplus
+namespace random_iterator_r123{
+/** 
+@defgroup AESNI ARS and AESNI Classes and Typedefs
+
+The ARS4x32, ARS1xm128i, AESNI4x32 and AESNI1xm128i classes export the member functions, typedefs and
+operator overloads required by a @ref CBRNG "CBRNG" class.
+
+ARS1xm128i and AESNI1xm128i are based on the AES block cipher and rely on the AES-NI hardware instructions
+available on some some new (2011) CPUs.
+
+The ARS1xm128i CBRNG and the use of AES for random number generation are described in 
+<a href="http://dl.acm.org/citation.cfm?doid=2063405"><i>Parallel Random Numbers:  As Easy as 1, 2, 3</i> </a>.
+Although it uses some cryptographic primitives, ARS1xm128i uses a cryptographically weak key schedule and is \b not suitable for cryptographic use.
+
+@class AESNI1xm128i
+@ingroup AESNI
+AESNI exports the member functions, typedefs and operator overloads required by a @ref CBRNG class.
+
+AESNI1xm128i uses the crypotgraphic AES round function, including the cryptographic key schedule.
+
+In contrast to the other CBRNGs in the Random123 library, the AESNI1xm128i_R::key_type is opaque
+and is \b not identical to the AESNI1xm128i_R::ukey_type.  Creating a key_type, using either the constructor
+or assignment operator, is significantly more time-consuming than running the bijection (hundreds
+of clock cycles vs. tens of clock cycles).
+
+AESNI1xm128i is only available when the feature-test macro RANDOM_ITERATOR_R123_USE_AES_NI is true, which
+should occur only when the compiler is configured to generate AES-NI instructions (or
+when defaults are overridden by compile-time, compiler-command-line options).
+
+As of September 2011, the authors know of no statistical flaws with AESNI1xm128i.  It
+would be an event of major cryptographic note if any such flaws were ever found.
+*/
+struct AESNI1xm128i{
+    typedef aesni1xm128i_ctr_t ctr_type;
+    typedef aesni1xm128i_ukey_t ukey_type;
+    typedef aesni1xm128i_key_t key_type;
+    static const unsigned int rounds=10;
+    ctr_type operator()(ctr_type ctr, key_type key) const{
+        return aesni1xm128i(ctr, key);
+    }
+};
+
+/* @class AESNI4x32 */
+struct AESNI4x32{
+    typedef aesni4x32_ctr_t ctr_type;
+    typedef aesni4x32_ukey_t ukey_type;
+    typedef aesni4x32_key_t key_type;
+    static const unsigned int rounds=10;
+    ctr_type operator()(ctr_type ctr, key_type key) const{
+        return aesni4x32(ctr, key);
+    }
+};
+
+/** @ingroup AESNI
+    @class AESNI1xm128i_R
+
+AESNI1xm128i_R is provided for completeness, but is only instantiable with ROUNDS=10, in
+which case it is identical to AESNI1xm128i */
+template <unsigned ROUNDS=10> 
+struct AESNI1xm128i_R : public AESNI1xm128i{
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(ROUNDS==10, "AESNI1xm128i_R<R> is only valid with R=10");
+};
+
+/** @class AESNI4x32_R **/
+template <unsigned ROUNDS=10> 
+struct AESNI4x32_R : public AESNI4x32{
+    RANDOM_ITERATOR_R123_STATIC_ASSERT(ROUNDS==10, "AESNI4x32_R<R> is only valid with R=10");
+};
+} // namespace random_iterator_r123
+#endif /* __cplusplus */
+
+#endif /* RANDOM_ITERATOR_R123_USE_AES_NI */
+
+#if RANDOM_ITERATOR_R123_USE_AES_OPENSSL
+#include "string.h"
+#include <openssl/aes.h>
+typedef struct r123array16x8 aesopenssl16x8_ctr_t;
+typedef struct r123array16x8 aesopenssl16x8_ukey_t;
+#ifdef __cplusplus
+struct aesopenssl16x8_key_t{
+    AES_KEY k;
+    aesopenssl16x8_key_t(){
+        aesopenssl16x8_ukey_t ukey={{}};
+        AES_set_encrypt_key((const unsigned char *)&ukey.v[0], 128, &k);
+    }
+    aesopenssl16x8_key_t(const aesopenssl16x8_ukey_t& ukey){
+        AES_set_encrypt_key((const unsigned char *)&ukey.v[0], 128, &k);
+    }
+    aesopenssl16x8_key_t& operator=(const aesopenssl16x8_ukey_t& ukey){
+        AES_set_encrypt_key((const unsigned char *)&ukey.v[0], 128, &k);
+        return *this;
+    }
+    bool operator==(const aesopenssl16x8_key_t& rhs) const{
+        return (k.rounds == rhs.k.rounds) && 0==::memcmp(&k.rd_key[0], &rhs.k.rd_key[0], (k.rounds+1) * 4 * sizeof(uint32_t));
+    }
+    bool operator!=(const aesopenssl16x8_key_t& rhs) const{
+        return !(*this == rhs);
+    }
+    friend std::ostream& operator<<(std::ostream& os, const aesopenssl16x8_key_t& v){
+        os << v.k.rounds;
+        const unsigned int *p = &v.k.rd_key[0];
+        for(int i=0; i<(v.k.rounds+1); ++i){
+            os << " " << p[0] << " " << p[1] << " " << p[2] << " " << p[3];
+            p += 4;
+        }
+        return os;
+    }
+    friend std::istream& operator>>(std::istream& is, aesopenssl16x8_key_t& v){
+        is >> v.k.rounds;
+        unsigned int *p = &v.k.rd_key[0];
+        for(int i=0; i<(v.k.rounds+1); ++i){
+            is >> p[0] >> p[1] >> p[2] >> p[3];
+            p += 4;
+        }
+        return is;
+    }
+};
+#else
+typedef struct aesopenssl16x8_key_t{
+    AES_KEY k;
+}aesopenssl16x8_key_t;
+RANDOM_ITERATOR_R123_STATIC_INLINE struct aesopenssl16x8_key_t aesopenssl16x8keyinit(aesopenssl16x8_ukey_t uk){
+    aesopenssl16x8_key_t ret;
+    AES_set_encrypt_key((const unsigned char *)&uk.v[0], 128, &ret.k);
+    return ret;
+}
+#endif
+
+RANDOM_ITERATOR_R123_STATIC_INLINE RANDOM_ITERATOR_R123_FORCE_INLINE(aesopenssl16x8_ctr_t aesopenssl16x8_R(aesopenssl16x8_ctr_t ctr, aesopenssl16x8_key_t key));
+RANDOM_ITERATOR_R123_STATIC_INLINE
+aesopenssl16x8_ctr_t aesopenssl16x8_R(aesopenssl16x8_ctr_t ctr, aesopenssl16x8_key_t key){
+    aesopenssl16x8_ctr_t ret;
+    AES_encrypt((const unsigned char*)&ctr.v[0], (unsigned char *)&ret.v[0], &key.k);
+    return ret;
+}
+
+#define aesopenssl16x8_rounds aesni4x32_rounds
+#define aesopenssl16x8(c,k) aesopenssl16x8_R(aesopenssl16x8_rounds)
+
+#ifdef __cplusplus
+namespace random_iterator_r123{
+struct AESOpenSSL16x8{
+    typedef aesopenssl16x8_ctr_t ctr_type;
+    typedef aesopenssl16x8_key_t key_type;
+    typedef aesopenssl16x8_ukey_t ukey_type;
+    static const unsigned int rounds=10;
+    ctr_type operator()(const ctr_type& in, const key_type& k){
+        ctr_type out;
+        AES_encrypt((const unsigned char *)&in[0], (unsigned char *)&out[0], &k.k);
+        return out;
+    }
+};
+} // namespace random_iterator_r123
+#endif /* __cplusplus */
+#endif /* RANDOM_ITERATOR_R123_USE_AES_OPENSSL */
+
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/array.h

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef _r123array_dot_h__
+#define _r123array_dot_h__
+#include "features/compilerfeatures.h"
+#include "features/sse.h"
+
+#if !defined(__cplusplus) || defined(__METAL_MACOS__)
+#define CXXMETHODS(_N, W, T)
+#define CXXOVERLOADS(_N, W, T)
+#define CXXMETHODS_REQUIRING_STL
+#else
+
+#include <stddef.h>
+#include <algorithm>
+#include <stdexcept>
+#include <iterator>
+#include <limits>
+#include <iostream>
+
+/** @defgroup arrayNxW The r123arrayNxW classes 
+
+    Each of the r123arrayNxW is a fixed size array of N W-bit unsigned integers.
+    It is functionally equivalent to the C++11 std::array<N, uintW_t>,
+    but does not require C++11 features or libraries.
+
+    In addition to meeting most of the requirements of a Container,
+    it also has a member function, incr(), which increments the zero-th
+    element and carrys overflows into higher indexed elements.  Thus,
+    by using incr(), sequences of up to 2^(N*W) distinct values
+    can be produced. 
+
+    If SSE is supported by the compiler, then the class
+    r123array1xm128i is also defined, in which the data member is an
+    array of one r123m128i object.
+
+    When compiling with __CUDA_ARCH__ defined, the reverse iterator
+    methods (rbegin, rend, crbegin, crend) are not defined because
+    CUDA does not support std::reverse_iterator.
+
+*/
+
+/** @cond HIDDEN_FROM_DOXYGEN */
+
+template <typename value_type>
+inline RANDOM_ITERATOR_R123_CUDA_DEVICE value_type assemble_from_u32(uint32_t *p32){
+    value_type v=0;
+    for(size_t i=0; i<(3+sizeof(value_type))/4; ++i)
+        v |= ((value_type)(*p32++)) << (32*i);
+    return v;
+}
+
+/** @endcond */
+
+#ifdef __CUDA_ARCH__
+/* CUDA can't handle std::reverse_iterator.  We *could* implement it
+   ourselves, but let's not bother until somebody really feels a need
+   to reverse-iterate through an r123array */
+#define CXXMETHODS_REQUIRING_STL
+#else
+#define  CXXMETHODS_REQUIRING_STL \
+    public: \
+    typedef std::reverse_iterator<iterator> reverse_iterator;           \
+    typedef std::reverse_iterator<const_iterator> const_reverse_iterator; \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE reverse_iterator rbegin(){ return reverse_iterator(end()); }                         \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reverse_iterator rbegin() const{ return const_reverse_iterator(end()); } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE reverse_iterator rend(){ return reverse_iterator(begin()); }        \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reverse_iterator rend() const{ return const_reverse_iterator(begin()); } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reverse_iterator crbegin() const{ return const_reverse_iterator(cend()); } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reverse_iterator crend() const{ return const_reverse_iterator(cbegin()); } 
+#endif
+
+// Work-alike methods and typedefs modeled on std::array:
+#define CXXMETHODS(_N, W, T)                                            \
+    typedef T value_type;                                               \
+    typedef T* iterator;                                                \
+    typedef const T* const_iterator;                                    \
+    typedef value_type& reference;                                      \
+    typedef const value_type& const_reference;                          \
+    typedef size_t size_type;                                           \
+    typedef ptrdiff_t difference_type;                                  \
+    typedef T* pointer;                                                 \
+    typedef const T* const_pointer;                                     \
+    /* Boost.array has static_size.  C++11 specializes tuple_size */    \
+    enum {static_size = _N};                                            \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE reference operator[](size_type i){return v[i];}                     \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reference operator[](size_type i) const {return v[i];}        \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE reference at(size_type i){ if(i >=  _N) RANDOM_ITERATOR_R123_THROW(std::out_of_range("array index out of range")); return (*this)[i]; } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reference at(size_type i) const { if(i >=  _N) RANDOM_ITERATOR_R123_THROW(std::out_of_range("array index out of range")); return (*this)[i]; } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE size_type size() const { return  _N; }                              \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE size_type max_size() const { return _N; }                           \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE bool empty() const { return _N==0; };                               \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE iterator begin() { return &v[0]; }                                  \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE iterator end() { return &v[_N]; }                                   \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_iterator begin() const { return &v[0]; }                      \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_iterator end() const { return &v[_N]; }                       \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_iterator cbegin() const { return &v[0]; }                     \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_iterator cend() const { return &v[_N]; }                      \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE pointer data(){ return &v[0]; }                                     \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_pointer data() const{ return &v[0]; }                         \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE reference front(){ return v[0]; }                                   \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reference front() const{ return v[0]; }                       \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE reference back(){ return v[_N-1]; }                                 \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE const_reference back() const{ return v[_N-1]; }                     \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE bool operator==(const r123array##_N##x##W& rhs) const{ \
+	/* CUDA3 does not have std::equal */ \
+	for (size_t i = 0; i < _N; ++i) \
+	    if (v[i] != rhs.v[i]) return false; \
+	return true; \
+    } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE bool operator!=(const r123array##_N##x##W& rhs) const{ return !(*this == rhs); } \
+    /* CUDA3 does not have std::fill_n */ \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE void fill(const value_type& val){ for (size_t i = 0; i < _N; ++i) v[i] = val; } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE void swap(r123array##_N##x##W& rhs){ \
+	/* CUDA3 does not have std::swap_ranges */ \
+	for (size_t i = 0; i < _N; ++i) { \
+	    T tmp = v[i]; \
+	    v[i] = rhs.v[i]; \
+	    rhs.v[i] = tmp; \
+	} \
+    } \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE r123array##_N##x##W& incr(RANDOM_ITERATOR_R123_ULONG_LONG n=1){                         \
+        /* This test is tricky because we're trying to avoid spurious   \
+           complaints about illegal shifts, yet still be compile-time   \
+           evaulated. */                                                \
+        if(sizeof(T)<sizeof(n) && n>>((sizeof(T)<sizeof(n))?8*sizeof(T):0) ) \
+            return incr_carefully(n);                                   \
+        if(n==1){                                                       \
+            ++v[0];                                                     \
+            if(_N==1 || RANDOM_ITERATOR_R123_BUILTIN_EXPECT(!!v[0], 1)) return *this;   \
+        }else{                                                          \
+            v[0] += n;                                                  \
+            if(_N==1 || RANDOM_ITERATOR_R123_BUILTIN_EXPECT(n<=v[0], 1)) return *this;  \
+        }                                                               \
+        /* We expect that the N==?? tests will be                       \
+           constant-folded/optimized away by the compiler, so only the  \
+           overflow tests (!!v[i]) remain to be done at runtime.  For  \
+           small values of N, it would be better to do this as an       \
+           uncondtional sequence of adc.  An experiment/optimization    \
+           for another day...                                           \
+           N.B.  The weird subscripting: v[_N>3?3:0] is to silence      \
+           a spurious error from icpc                                   \
+           */                                                           \
+        ++v[_N>1?1:0];                                                  \
+        if(_N==2 || RANDOM_ITERATOR_R123_BUILTIN_EXPECT(!!v[_N>1?1:0], 1)) return *this; \
+        ++v[_N>2?2:0];                                                  \
+        if(_N==3 || RANDOM_ITERATOR_R123_BUILTIN_EXPECT(!!v[_N>2?2:0], 1)) return *this;  \
+        ++v[_N>3?3:0];                                                  \
+        for(size_t i=4; i<_N; ++i){                                     \
+            if( RANDOM_ITERATOR_R123_BUILTIN_EXPECT(!!v[i-1], 1) ) return *this;        \
+            ++v[i];                                                     \
+        }                                                               \
+        return *this;                                                   \
+    }                                                                   \
+    /* seed(SeedSeq) would be a constructor if having a constructor */  \
+    /* didn't cause headaches with defaults */                          \
+    template <typename SeedSeq>                                         \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE static r123array##_N##x##W seed(SeedSeq &ss){      \
+        r123array##_N##x##W ret;                                        \
+        const size_t Ngen = _N*((3+sizeof(value_type))/4);              \
+        uint32_t u32[Ngen];                                             \
+        uint32_t *p32 = &u32[0];                                        \
+        ss.generate(&u32[0], &u32[Ngen]);                               \
+        for(size_t i=0; i<_N; ++i){                                     \
+            ret.v[i] = assemble_from_u32<value_type>(p32);              \
+            p32 += (3+sizeof(value_type))/4;                            \
+        }                                                               \
+        return ret;                                                     \
+    }                                                                   \
+protected:                                                              \
+    RANDOM_ITERATOR_R123_CUDA_DEVICE r123array##_N##x##W& incr_carefully(RANDOM_ITERATOR_R123_ULONG_LONG n){ \
+        /* n may be greater than the maximum value of a single value_type */ \
+        value_type vtn;                                                 \
+        vtn = n;                                                        \
+        v[0] += n;                                                      \
+        const unsigned rshift = 8* ((sizeof(n)>sizeof(value_type))? sizeof(value_type) : 0); \
+        for(size_t i=1; i<_N; ++i){                                     \
+            if(rshift){                                                 \
+                n >>= rshift;                                           \
+            }else{                                                      \
+                n=0;                                                    \
+            }                                                           \
+            if( v[i-1] < vtn )                                          \
+                ++n;                                                    \
+            if( n==0 ) break;                                           \
+            vtn = n;                                                    \
+            v[i] += n;                                                  \
+        }                                                               \
+        return *this;                                                   \
+    }                                                                   \
+
+/** @cond HIDDEN_FROM_DOXYGEN */
+
+// There are several tricky considerations for the insertion and extraction
+// operators:
+// - we would like to be able to print r123array16x8 as a sequence of 16 integers,
+//   not as 16 bytes.
+// - we would like to be able to print r123array1xm128i.
+// - we do not want an int conversion operator in r123m128i because it causes
+//   lots of ambiguity problems with automatic promotions.
+// Solution: r123arrayinsertable and r123arrayextractable
+
+template<typename T>
+struct r123arrayinsertable{
+    const T& v;
+    r123arrayinsertable(const T& t_) : v(t_) {} 
+    friend std::ostream& operator<<(std::ostream& os, const r123arrayinsertable<T>& t){
+        return os << t.v;
+    }
+};
+
+template<>
+struct r123arrayinsertable<uint8_t>{
+    const uint8_t& v;
+    r123arrayinsertable(const uint8_t& t_) : v(t_) {} 
+    friend std::ostream& operator<<(std::ostream& os, const r123arrayinsertable<uint8_t>& t){
+        return os << (int)t.v;
+    }
+};
+
+template<typename T>
+struct r123arrayextractable{
+    T& v;
+    r123arrayextractable(T& t_) : v(t_) {}
+    friend std::istream& operator>>(std::istream& is, r123arrayextractable<T>& t){
+        return is >> t.v;
+    }
+};
+
+template<>
+struct r123arrayextractable<uint8_t>{
+    uint8_t& v;
+    r123arrayextractable(uint8_t& t_) : v(t_) {} 
+    friend std::istream& operator>>(std::istream& is, r123arrayextractable<uint8_t>& t){
+        int i;
+        is >>  i;
+        t.v = i;
+        return is;
+    }
+};
+/** @endcond */
+
+#define CXXOVERLOADS(_N, W, T)                                          \
+                                                                        \
+inline std::ostream& operator<<(std::ostream& os, const r123array##_N##x##W& a){   \
+    os << r123arrayinsertable<T>(a.v[0]);                                  \
+    for(size_t i=1; i<_N; ++i)                                          \
+        os << " " << r123arrayinsertable<T>(a.v[i]);                       \
+    return os;                                                          \
+}                                                                       \
+                                                                        \
+inline std::istream& operator>>(std::istream& is, r123array##_N##x##W& a){         \
+    for(size_t i=0; i<_N; ++i){                                         \
+        r123arrayextractable<T> x(a.v[i]);                                 \
+        is >> x;                                                        \
+    }                                                                   \
+    return is;                                                          \
+}                                                                       \
+                                                                        \
+namespace random_iterator_r123{                                                        \
+ typedef r123array##_N##x##W Array##_N##x##W;                          \
+}
+                                                                        
+#endif /* __cplusplus */
+
+/* _r123array_tpl expands to a declaration of struct r123arrayNxW.  
+
+   In C, it's nothing more than a struct containing an array of N
+   objects of type T.
+
+   In C++ it's the same, but endowed with an assortment of member
+   functions, typedefs and friends.  In C++, r123arrayNxW looks a lot
+   like std::array<T,N>, has most of the capabilities of a container,
+   and satisfies the requirements outlined in compat/Engine.hpp for
+   counter and key types.  ArrayNxW, in the r123 namespace is
+   a typedef equivalent to r123arrayNxW.
+*/
+
+#define _r123array_tpl(_N, W, T)                    \
+    /** @ingroup arrayNxW */                        \
+    /** @see arrayNxW */                            \
+struct r123array##_N##x##W{                         \
+ T v[_N];                                       \
+ CXXMETHODS(_N, W, T)                           \
+ CXXMETHODS_REQUIRING_STL                       \
+};                                              \
+                                                \
+CXXOVERLOADS(_N, W, T)
+
+
+_r123array_tpl(1, 32, uint32_t)  /* r123array1x32 */
+_r123array_tpl(2, 32, uint32_t)  /* r123array2x32 */
+_r123array_tpl(4, 32, uint32_t)  /* r123array4x32 */
+_r123array_tpl(8, 32, uint32_t)  /* r123array8x32 */
+
+#if RANDOM_ITERATOR_R123_USE_64BIT
+_r123array_tpl(1, 64, uint64_t)  /* r123array1x64 */
+_r123array_tpl(2, 64, uint64_t)  /* r123array2x64 */
+_r123array_tpl(4, 64, uint64_t)  /* r123array4x64 */
+#endif
+
+_r123array_tpl(16, 8, uint8_t)  /* r123array16x8 for ARSsw, AESsw */
+
+#if RANDOM_ITERATOR_R123_USE_SSE
+_r123array_tpl(1, m128i, r123m128i) /* r123array1x128i for ARSni, AESni */
+#endif
+
+/* In C++, it's natural to use sizeof(a::value_type), but in C it's
+   pretty convoluted to figure out the width of the value_type of an
+   r123arrayNxW:
+*/
+#define RANDOM_ITERATOR_R123_W(a)   (8*sizeof(((a *)0)->v[0]))
+
+/** @namespace random_iterator_r123
+  Most of the Random123 C++ API is contained in the r123 namespace. 
+*/
+
+#endif
+

corsika/detail/framework/random/random_iterator/detail/Random123/ars.h

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __Random123_ars_dot_hpp__
+#define __Random123_ars_dot_hpp__
+
+#include "features/compilerfeatures.h"
+#include "array.h"
+
+#if RANDOM_ITERATOR_R123_USE_AES_NI
+
+#ifndef ARS1xm128i_DEFAULT_ROUNDS
+#define ARS1xm128i_DEFAULT_ROUNDS 7
+#endif
+
+/** @ingroup AESNI */
+enum r123_enum_ars1xm128i {ars1xm128i_rounds = ARS1xm128i_DEFAULT_ROUNDS};
+
+/* ARS1xm128i with Weyl keys.  Fast, and Crush-resistant, but NOT CRYPTO. */
+/** @ingroup AESNI */
+typedef struct r123array1xm128i ars1xm128i_ctr_t;
+/** @ingroup AESNI */
+typedef struct r123array1xm128i ars1xm128i_key_t;
+/** @ingroup AESNI */
+typedef struct r123array1xm128i ars1xm128i_ukey_t;
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE ars1xm128i_key_t ars1xm128ikeyinit(ars1xm128i_ukey_t uk) { return uk; }
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE ars1xm128i_ctr_t ars1xm128i_R(unsigned int Nrounds, ars1xm128i_ctr_t in, ars1xm128i_key_t k){
+    __m128i kweyl = _mm_set_epi64x(RANDOM_ITERATOR_R123_64BIT(0xBB67AE8584CAA73B), /* sqrt(3) - 1.0 */
+                                   RANDOM_ITERATOR_R123_64BIT(0x9E3779B97F4A7C15)); /* golden ratio */
+    /* N.B.  the aesenc instructions do the xor *after*
+    // so if we want to follow the AES pattern, we
+    // have to do the initial xor explicitly */
+    __m128i kk = k.v[0].m;
+    __m128i v = _mm_xor_si128(in.v[0].m, kk);
+    ars1xm128i_ctr_t ret;
+    RANDOM_ITERATOR_R123_ASSERT(Nrounds<=10);
+    if( Nrounds>1 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>2 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>3 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>4 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>5 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>6 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>7 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>8 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    if( Nrounds>9 ){
+        kk = _mm_add_epi64(kk, kweyl);
+        v = _mm_aesenc_si128(v, kk);
+    }
+    kk = _mm_add_epi64(kk, kweyl);
+    v = _mm_aesenclast_si128(v, kk);
+    ret.v[0].m = v;
+    return ret;
+}
+
+/** @def ars1xm128i
+@ingroup AESNI
+The ars1mx128i macro provides a C API interface to the @ref AESNI "ARS" CBRNG with the default number of rounds i.e. \c ars1xm128i_rounds **/
+#define ars1xm128i(c,k) ars1xm128i_R(ars1xm128i_rounds, c, k)
+
+/** @ingroup AESNI */
+typedef struct r123array4x32 ars4x32_ctr_t;
+/** @ingroup AESNI */
+typedef struct r123array4x32 ars4x32_key_t;
+/** @ingroup AESNI */
+typedef struct r123array4x32 ars4x32_ukey_t;
+
+typedef struct r123array2x64 ars2x64_ctr_t;
+/** @ingroup AESNI */
+typedef struct r123array2x64 ars2x64_key_t;
+/** @ingroup AESNI */
+typedef struct r123array2x64 ars2x64_ukey_t;
+
+/** @ingroup AESNI */
+enum r123_enum_ars4x32 {ars4x32_rounds = ARS1xm128i_DEFAULT_ROUNDS};
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE ars4x32_key_t ars4x32keyinit(ars4x32_ukey_t uk) { return uk; }
+/** @ingroup AESNI */
+RANDOM_ITERATOR_R123_STATIC_INLINE ars4x32_ctr_t ars4x32_R(unsigned int Nrounds, ars4x32_ctr_t c, ars4x32_key_t k){
+    ars1xm128i_ctr_t c128;
+    ars1xm128i_key_t k128;
+    c128.v[0].m = _mm_set_epi32(c.v[3], c.v[2], c.v[1], c.v[0]);
+    k128.v[0].m = _mm_set_epi32(k.v[3], k.v[2], k.v[1], k.v[0]);
+    c128 = ars1xm128i_R(Nrounds, c128, k128);
+    _mm_storeu_si128((__m128i*)&c.v[0], c128.v[0].m);
+    return c;
+}
+
+/** @def ars4x32
+@ingroup AESNI
+The ars4x32 macro provides a C API interface to the @ref AESNI "ARS" CBRNG with the default number of rounds i.e. \c ars4x32_rounds **/
+#define ars4x32(c,k) ars4x32_R(ars4x32_rounds, c, k)
+
+#ifdef __cplusplus
+namespace random_iterator_r123{
+/** 
+@ingroup AESNI
+
+ARS1xm128i_R exports the member functions, typedefs and operator overloads required by a @ref CBRNG class.
+
+ARS1xm128i uses the crypotgraphic AES round function, but a @b non-cryptographc key schedule
+to save time and space.
+
+ARS1xm128i is only available when the feature-test macro RANDOM_ITERATOR_R123_USE_AES_NI is true, which
+should occur only when the compiler is configured to generate AES-NI instructions (or
+when defaults are overridden by compile-time, compiler-command-line options).
+
+The template argument, ROUNDS, is the number of times the ARS round
+functions will be applied.
+
+As of September 2011, the authors know of no statistical flaws with
+ROUNDS=5 or more.
+
+@class ARS1xm128i_R
+
+*/
+template<unsigned int ROUNDS>
+struct ARS1xm128i_R{
+    typedef ars1xm128i_ctr_t ctr_type;
+    typedef ars1xm128i_key_t key_type;
+    typedef ars1xm128i_key_t ukey_type;
+    static const unsigned int rounds=ROUNDS;
+    RANDOM_ITERATOR_R123_FORCE_INLINE(ctr_type operator()(ctr_type ctr, key_type key) const){
+        return ars1xm128i_R(ROUNDS, ctr, key);
+    }
+};
+
+/** @class ARS4x32_R
+    @ingroup AESNI
+*/
+
+template<unsigned int ROUNDS>
+struct ARS4x32_R{
+    typedef ars4x32_ctr_t ctr_type;
+    typedef ars4x32_key_t key_type;
+    typedef ars4x32_key_t ukey_type;
+    static const unsigned int rounds=ROUNDS;
+    RANDOM_ITERATOR_R123_FORCE_INLINE(ctr_type operator()(ctr_type ctr, key_type key) const){
+        return ars4x32_R(ROUNDS, ctr, key);
+    }
+
+
+};
+
+template<unsigned int ROUNDS>
+struct ARS2x64_R{
+
+    typedef ars2x64_ctr_t ctr_type;
+    typedef ars2x64_key_t key_type;
+    typedef ars2x64_key_t ukey_type;
+
+    static const unsigned int rounds=ROUNDS;
+
+    RANDOM_ITERATOR_R123_FORCE_INLINE(ctr_type operator()(ctr_type ctr, key_type key) const){
+
+    	ars4x32_ctr_t ctr_{ctr.v[0], ctr.v[0]>>32 , ctr.v[1], ctr.v[1]>>32};
+    	ars4x32_key_t key_{key.v[0], key.v[0]>>32 , key.v[1], key.v[1]>>32};
+    	ars4x32_ctr_t res_ = ars4x32_R(ROUNDS, ctr_, key_);
+    	ctr_type res{{}};
+    	res.v[0] = (res.v[0] | res_[1])<<32;
+    	res.v[0] = (res.v[0] | res_[0]);
+    	res.v[1] = (res.v[1] | res_[3])<<32;
+    	res.v[1] = (res.v[1] | res_[2]);
+
+        return res;
+    }
+};
+/**
+@ingroup AESNI
+
+@class ARS1xm128i_R
+  ARS1xm128i is equivalent to ARS1xm128i_R<7>.    With 7 rounds,
+  the ARS1xm128i CBRNG  has a considerable safety margin over the minimum number
+  of rounds with no known statistical flaws, but still has excellent
+  performance. */
+typedef ARS1xm128i_R<ars1xm128i_rounds> ARS1xm128i;
+typedef ARS4x32_R<ars4x32_rounds> ARS4x32;
+typedef ARS2x64_R<ars4x32_rounds> ARS2x64;
+
+} // namespace random_iterator_r123
+
+#endif /* __cplusplus */
+
+#endif /* RANDOM_ITERATOR_R123_USE_AES_NI */
+
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/boxmuller.hpp

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+// This file implements the Box-Muller method for generating gaussian
+// random variables (GRVs).  Box-Muller has the advantage of
+// deterministically requiring exactly two uniform random variables as
+// input and producing exactly two GRVs as output, which makes it
+// especially well-suited to the counter-based generators in
+// Random123.  Other methods (e.g., Ziggurat, polar) require an
+// indeterminate number of inputs for each output and so require a
+// 'MicroURNG' to be used with Random123.  The down side of Box-Muller
+// is that it calls sincos, log and sqrt, which may be slow.  However,
+// on GPUs, these functions are remarkably fast, which makes
+// Box-Muller the fastest GRV generator we know of on GPUs.
+//
+// This file exports two structs and one overloaded function,
+// all in the r123 namespace:
+//   struct r123::float2{ float x,y; }
+//   struct r123::double2{ double x,y; }
+//
+//   r123::float2  r123::boxmuller(uint32_t u0, uint32_t u1);
+//   r123::double2 r123::boxmuller(uint64_t u0, uint64_t u1);
+//
+// float2 and double2 are identical to their synonymous global-
+// namespace structures in CUDA.
+//
+// This file may not be as portable, and has not been tested as
+// rigorously as other files in the library, e.g., the generators.
+// Nevertheless, we hope it is useful and we encourage developers to
+// copy it and modify it for their own use.  We invite comments and
+// improvements.
+
+#ifndef _r123_BOXMULLER_HPP__
+#define _r123_BOXMULLER_HPP__
+
+#include <Random123/features/compilerfeatures.h>
+#include <Random123/uniform.hpp>
+#include <math.h>
+
+namespace random_iterator_r123 {
+
+#if !defined(__CUDACC__)
+  typedef struct {
+    float x, y;
+  } float2;
+  typedef struct {
+    double x, y;
+  } double2;
+#else
+  typedef ::float2 float2;
+  typedef ::double2 double2;
+#endif
+
+#if !defined(RANDOM_ITERATOR_R123_NO_SINCOS) && defined(__APPLE__)
+/* MacOS X 10.10.5 (2015) doesn't have sincosf */
+#define RANDOM_ITERATOR_R123_NO_SINCOS 1
+#endif
+
+#if RANDOM_ITERATOR_R123_NO_SINCOS /* enable this if sincos and sincosf are not in the \
+                                      math library */
+  RANDOM_ITERATOR_R123_CUDA_DEVICE RANDOM_ITERATOR_R123_STATIC_INLINE void sincosf(
+      float x, float* s, float* c) {
+    *s = sinf(x);
+    *c = cosf(x);
+  }
+
+  RANDOM_ITERATOR_R123_CUDA_DEVICE RANDOM_ITERATOR_R123_STATIC_INLINE void sincos(
+      double x, double* s, double* c) {
+    *s = sin(x);
+    *c = cos(x);
+  }
+#endif /* sincos is not in the math library */
+
+#if !defined(CUDART_VERSION) || \
+    CUDART_VERSION < 5000 /* enabled if sincospi and sincospif are not in math lib */
+
+  RANDOM_ITERATOR_R123_CUDA_DEVICE RANDOM_ITERATOR_R123_STATIC_INLINE void sincospif(
+      float x, float* s, float* c) {
+    const float PIf = 3.1415926535897932f;
+    sincosf(PIf * x, s, c);
+  }
+
+  RANDOM_ITERATOR_R123_CUDA_DEVICE RANDOM_ITERATOR_R123_STATIC_INLINE void sincospi(
+      double x, double* s, double* c) {
+    const double PI = 3.1415926535897932;
+    sincos(PI * x, s, c);
+  }
+#endif /* sincospi is not in math lib */
+
+  /*
+   * take two 32bit unsigned random values and return a float2 with
+   * two random floats in a normal distribution via a Box-Muller transform
+   */
+  RANDOM_ITERATOR_R123_CUDA_DEVICE RANDOM_ITERATOR_R123_STATIC_INLINE float2
+  boxmuller(uint32_t u0, uint32_t u1) {
+    float r;
+    float2 f;
+    sincospif(uneg11<float>(u0), &f.x, &f.y);
+    r = sqrtf(-2.f * logf(u01<float>(u1))); // u01 is guaranteed to avoid 0.
+    f.x *= r;
+    f.y *= r;
+    return f;
+  }
+
+  /*
+   * take two 64bit unsigned random values and return a double2 with
+   * two random doubles in a normal distribution via a Box-Muller transform
+   */
+  RANDOM_ITERATOR_R123_CUDA_DEVICE RANDOM_ITERATOR_R123_STATIC_INLINE double2
+  boxmuller(uint64_t u0, uint64_t u1) {
+    double r;
+    double2 f;
+
+    sincospi(uneg11<double>(u0), &f.x, &f.y);
+    r = sqrt(-2. * log(u01<double>(u1))); // u01 is guaranteed to avoid 0.
+    f.x *= r;
+    f.y *= r;
+    return f;
+  }
+} // namespace random_iterator_r123
+
+#endif /* BOXMULLER_H__ */

corsika/detail/framework/random/random_iterator/detail/Random123/conventional/Engine.hpp

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __Engine_dot_hpp_
+#define __Engine_dot_hpp_
+
+#include "../features/compilerfeatures.h"
+#include "../array.h"
+#include <limits>
+#include <stdexcept>
+#include <sstream>
+#include <algorithm>
+#include <vector>
+#if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
+#include <type_traits>
+#endif
+
+namespace random_iterator_r123 {
+  /**
+    If G satisfies the requirements of a CBRNG, and has a ctr_type whose
+    value_type is an unsigned integral type, then Engine<G> satisfies
+    the requirements of a C++11 "Uniform Random Number Engine" and can
+    be used in any context where such an object is expected.
+
+    Note that wrapping a counter based RNG with a traditional API in
+    this way obscures much of the power of counter based PRNGs.
+    Nevertheless, it may be of value in applications that are already
+    coded to work with the C++11 random number engines.
+
+    The MicroURNG template in MicroURNG.hpp
+    provides the more limited functionality of a C++11 "Uniform
+    Random Number Generator", but leaves the application in control
+    of counters and keys and hence may be preferable to the Engine template.
+    For example, a MicroURNG allows one to use C++11 "Random Number
+    Distributions"  without giving up control over the counters
+    and keys.
+  */
+
+  template <typename CBRNG>
+  struct Engine {
+    typedef CBRNG cbrng_type;
+    typedef typename CBRNG::ctr_type ctr_type;
+    typedef typename CBRNG::key_type key_type;
+    typedef typename CBRNG::ukey_type ukey_type;
+    typedef typename ctr_type::value_type result_type;
+
+  protected:
+    cbrng_type b;
+    key_type key;
+    ctr_type c;
+    ctr_type v;
+
+    void fix_invariant() {
+      if (v.back() != 0) {
+        result_type vv = v.back();
+        v = b(c, key);
+        v.back() = vv;
+      }
+    }
+
+  public:
+    explicit Engine()
+        : b()
+        , c() {
+      ukey_type x = {{}};
+      v.back() = 0;
+      key = x;
+    }
+    explicit Engine(result_type r)
+        : b()
+        , c() {
+      ukey_type x = {{typename ukey_type::value_type(r)}};
+      v.back() = 0;
+      key = x;
+    }
+    // 26.5.3 says that the SeedSeq templates shouldn't particpate in
+    // overload resolution unless the type qualifies as a SeedSeq.
+    // How that is determined is unspecified, except that "as a
+    // minimum a type shall not qualify as a SeedSeq if it is
+    // implicitly convertible to a result_type."
+    //
+    // First, we make sure that even the non-const copy constructor
+    // works as expected.  In addition, if we've got C++11
+    // type_traits, we use enable_if and is_convertible to implement
+    // the convertible-to-result_type restriction.  Otherwise, the
+    // template is unconditional and will match in some surpirsing
+    // and undesirable situations.
+    Engine(Engine& e)
+        : b(e.b)
+        , key(e.key)
+        , c(e.c) {
+      v.back() = e.v.back();
+      fix_invariant();
+    }
+    Engine(const Engine& e)
+        : b(e.b)
+        , key(e.key)
+        , c(e.c) {
+      v.back() = e.v.back();
+      fix_invariant();
+    }
+
+    template <typename SeedSeq>
+    explicit Engine(SeedSeq& s
+#if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
+                    ,
+                    typename std::enable_if<
+                        !std::is_convertible<SeedSeq, result_type>::value>::type* = 0
+#endif
+                    )
+        : b()
+        , c() {
+      ukey_type ukey = ukey_type::seed(s);
+      key = ukey;
+      v.back() = 0;
+    }
+    void seed(result_type r) { *this = Engine(r); }
+    template <typename SeedSeq>
+    void seed(SeedSeq& s
+#if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
+              ,
+              typename std::enable_if<
+                  !std::is_convertible<SeedSeq, result_type>::value>::type* = 0
+#endif
+    ) {
+      *this = Engine(s);
+    }
+    void seed() { *this = Engine(); }
+    friend bool operator==(const Engine& lhs, const Engine& rhs) {
+      return lhs.c == rhs.c && lhs.v.back() == rhs.v.back() && lhs.key == rhs.key;
+    }
+    friend bool operator!=(const Engine& lhs, const Engine& rhs) {
+      return lhs.c != rhs.c || lhs.v.back() != rhs.v.back() || lhs.key != rhs.key;
+    }
+
+    friend std::ostream& operator<<(std::ostream& os, const Engine& be) {
+      return os << be.c << " " << be.key << " " << be.v.back();
+    }
+
+    friend std::istream& operator>>(std::istream& is, Engine& be) {
+      is >> be.c >> be.key >> be.v.back();
+      be.fix_invariant();
+      return is;
+    }
+
+    // The <random> shipped with MacOS Xcode 4.5.2 imposes a
+    // non-standard requirement that URNGs also have static data
+    // members: _Min and _Max.  Later versions of libc++ impose the
+    // requirement only when constexpr isn't supported.  Although the
+    // Xcode 4.5.2 requirement is clearly non-standard, it is unlikely
+    // to be fixed and it is very easy work around.  We certainly
+    // don't want to go to great lengths to accommodate every buggy
+    // library we come across, but in this particular case, the effort
+    // is low and the benefit is high, so it's worth doing.  Thanks to
+    // Yan Zhou for pointing this out to us.  See similar code in
+    // ../MicroURNG.hpp
+    const static result_type _Min = 0;
+    const static result_type _Max = ~((result_type)0);
+
+    static RANDOM_ITERATOR_R123_CONSTEXPR result_type min
+    RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
+      return _Min;
+    }
+    static RANDOM_ITERATOR_R123_CONSTEXPR result_type max
+    RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
+      return _Max;
+    }
+
+    result_type operator()() {
+      if (c.size() == 1) // short-circuit the scalar case.  Compilers aren't mind-readers.
+        return b(c.incr(), key)[0];
+      result_type& elem = v.back();
+      if (elem == 0) {
+        v = b(c.incr(), key);
+        result_type ret = v.back();
+        elem = c.size() - 1;
+        return ret;
+      }
+      return v[--elem];
+    }
+
+    void discard(RANDOM_ITERATOR_R123_ULONG_LONG skip) {
+      // don't forget:  elem counts down
+      size_t nelem = c.size();
+      size_t sub = skip % nelem;
+      result_type& elem = v.back();
+      skip /= nelem;
+      if (elem < sub) {
+        elem += nelem;
+        skip++;
+      }
+      elem -= sub;
+      c.incr(skip);
+      fix_invariant();
+    }
+
+    //--------------------------
+    // Some bonus methods, not required for a Random Number
+    // Engine
+
+    // Constructors and seed() method for ukey_type seem useful
+    // We need const and non-const to supersede the SeedSeq template.
+    explicit Engine(const ukey_type& uk)
+        : key(uk)
+        , c() {
+      v.back() = 0;
+    }
+    explicit Engine(ukey_type& uk)
+        : key(uk)
+        , c() {
+      v.back() = 0;
+    }
+    void seed(const ukey_type& uk) { *this = Engine(uk); }
+    void seed(ukey_type& uk) { *this = Engine(uk); }
+
+#if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
+    template <typename DUMMY = void>
+    explicit Engine(const key_type& k,
+                    typename std::enable_if<!std::is_same<ukey_type, key_type>::value,
+                                            DUMMY>::type* = 0)
+        : key(k)
+        , c() {
+      v.back() = 0;
+    }
+
+    template <typename DUMMY = void>
+    void seed(const key_type& k,
+              typename std::enable_if<!std::is_same<ukey_type, key_type>::value,
+                                      DUMMY>::type* = 0) {
+      *this = Engine(k);
+    }
+#endif
+
+    // Forward the e(counter) to the CBRNG we are templated
+    // on, using the current value of the key.
+    ctr_type operator()(const ctr_type& c) const { return b(c, key); }
+
+    key_type getkey() const { return key; }
+
+    // N.B.  setkey(k) is different from seed(k) because seed(k) zeros
+    // the counter (per the C++11 requirements for an Engine), whereas
+    // setkey does not.
+    void setkey(const key_type& k) {
+      key = k;
+      fix_invariant();
+    }
+
+    // Maybe the caller want's to know the details of
+    // the internal state, e.g., so it can call a different
+    // bijection with the same counter.
+    std::pair<ctr_type, result_type> getcounter() const {
+      return std::make_pair(c, v.back());
+    }
+
+    // And the inverse.
+    void setcounter(const ctr_type& _c, result_type _elem) {
+      static const size_t nelem = c.size();
+      if (_elem >= nelem)
+        throw std::range_error("Engine::setcounter called  with elem out of range");
+      c = _c;
+      v.back() = _elem;
+      fix_invariant();
+    }
+
+    void setcounter(const std::pair<ctr_type, result_type>& ce) {
+      setcounter(ce.first, ce.second);
+    }
+  };
+} // namespace random_iterator_r123
+
+#endif

corsika/detail/framework/random/random_iterator/detail/Random123/conventional/gsl_cbrng.h

+/*
+Copyright 2010-2011, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __r123_compat_gslrng_dot_h__
+#define __r123_compat_gslrng_dot_h__
+
+#include <gsl/gsl_rng.h>
+#include <string.h>
+
+/**
+   The macro:  GSL_CBRNG(NAME, CBRNGNAME)
+   declares the necessary structs and  constants that define a
+   gsl_rng_NAME type based on the counter-based RNG CBRNGNAME.  For example:
+
+   Usage:
+
+   @code
+   #include <Random123/threefry.h>
+   #include <Random123/conventional/gsl_cbrng.h>  // this file
+   GSL_CBRNG(cbrng, threefry4x32); // creates gsl_rng_cbrng
+
+   int main(int argc, char **argv){
+       gsl_rng *r = gsl_rng_alloc(gsl_rng_cbrng);
+       ... use r as you would use any other gsl_rng ...   
+    }
+    @endcode
+
+    It requires that NAME be the name of a CBRNG that follows the
+    naming and stylistic conventions of the Random123 library.
+
+    Note that wrapping a \ref CBRNG "counter-based PRNG" with a traditional API in
+    this way obscures much of the power of the CBRNG API.
+    Nevertheless, it may be of value to applications that are already
+    coded to work with GSL random number generators, and that wish
+    to use the RNGs in the Random123 library.
+
+ */ 
+
+#define GSL_CBRNG(NAME, CBRNGNAME)                                      \
+const gsl_rng_type *gsl_rng_##NAME;                                     \
+                                                                        \
+typedef struct{                                                         \
+    CBRNGNAME##_ctr_t ctr;                                                   \
+    CBRNGNAME##_ctr_t r;                                                     \
+    CBRNGNAME##_key_t key;                                                   \
+    int elem;                                                           \
+} NAME##_state;                                                         \
+                                                                        \
+static unsigned long int NAME##_get(void *vstate){                      \
+    NAME##_state *st = (NAME##_state *)vstate;                          \
+    const int N=sizeof(st->ctr.v)/sizeof(st->ctr.v[0]);                 \
+    if( st->elem == 0 ){                                                \
+        ++st->ctr.v[0];                                                 \
+        if( N>1 && st->ctr.v[0] == 0 ) ++st->ctr.v[1];                  \
+        if( N>2 && st->ctr.v[1] == 0 ) ++st->ctr.v[2];                  \
+        if( N>3 && st->ctr.v[2] == 0 ) ++st->ctr.v[3];                  \
+        st->r = CBRNGNAME(st->ctr, st->key);                                 \
+        st->elem = N;                                                   \
+    }                                                                   \
+    return 0xffffffffUL & st->r.v[--st->elem];                          \
+}                                                                       \
+                                                                        \
+static double                                                           \
+NAME##_get_double (void * vstate)                                       \
+{                                                                       \
+    return NAME##_get (vstate)/4294967296.0;                            \
+}                                                                       \
+                                                                        \
+static void NAME##_set(void *vstate, unsigned long int s){              \
+    NAME##_state *st = (NAME##_state *)vstate;                          \
+    st->elem = 0;                                                       \
+    /* Assume that key and ctr have an array member, v,                 \
+       as if they are r123arrayNxW.  If not, this will fail             \
+       to compile.  In particular, this macro fails to compile          \
+       when the underlying CBRNG requires use of keyinit */             \
+    memset(&st->ctr.v[0], 0, sizeof(st->ctr.v));                        \
+    memset(&st->key.v[0], 0, sizeof(st->key.v));                        \
+    /* GSL 1.15 documentation says this about gsl_rng_set:              \
+         Note that the most generators only accept 32-bit seeds, with higher \
+         values being reduced modulo 2^32.  For generators with smaller \
+         ranges the maximum seed value will typically be lower.         \
+     so we won't jump through any hoops here to deal with               \
+     high bits if sizeof(unsigned long) > sizeof(uint32_t). */          \
+    st->key.v[0] = s;                                                   \
+}                                                                       \
+                                                                        \
+static const gsl_rng_type NAME##_type = {                               \
+    #NAME,                                                              \
+    0xffffffffUL,                                                       \
+    0,                                                                  \
+    sizeof(NAME##_state),                                               \
+    &NAME##_set,                                                        \
+    &NAME##_get,                                                        \
+    &NAME##_get_double                                                  \
+};                                                                      \
+                                                                        \
+const gsl_rng_type *gsl_rng_##NAME = &NAME##_type
+
+#endif
+

corsika/detail/framework/random/random_iterator/detail/Random123/features/clangfeatures.h

+/*
+Copyright 2010-2016, D. E. Shaw Research.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of D. E. Shaw Research nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef __clangfeatures_dot_hpp
+#define __clangfeatures_dot_hpp
+
+#ifndef RANDOM_ITERATOR_R123_USE_X86INTRIN_H
+#if (defined(__x86_64__)||defined(__i386__))
+#define RANDOM_ITERATOR_R123_USE_X86INTRIN_H 1
+#else
+#define RANDOM_ITERATOR_R123_USE_X86INTRIN_H 0
+#endif
+#endif
+
+#ifndef RANDOM_ITERATOR_R123_USE_CXX11_UNRESTRICTED_UNIONS
+#define RANDOM_ITERATOR_R123_USE_CXX11_UNRESTRICTED_UNIONS __has_feature(cxx_unrestricted_unions)
+#endif
+
+#ifndef RANDOM_ITERATOR_R123_USE_CXX11_STATIC_ASSERT
+#define RANDOM_ITERATOR_R123_USE_CXX11_STATIC_ASSERT __has_feature(cxx_static_assert)
+#endif
+
+// With clang-3.6, -Wall warns about unused-local-typedefs.
+// The "obvious" thing to do is to ignore -Wunused-local-typedefs,
+// but that doesn't work because earlier versions of clang blow
+// up on an 'unknown warning group'.  So we briefly ignore -Wall...
+// It's tempting to just give up on static assertions in pre-c++11 code.
+#if !RANDOM_ITERATOR_R123_USE_CXX11_STATIC_ASSERT && !defined(RANDOM_ITERATOR_R123_STATIC_ASSERT)
+#define RANDOM_ITERATOR_R123_STATIC_ASSERT(expr, msg) \
+_Pragma("clang diagnostic push")                      \
+_Pragma("clang diagnostic ignored \"-Wall\"")     \
+typedef char static_assertion[(!!(expr))*2-1] \
+_Pragma("clang diagnostic pop")
+#endif
+
+#ifndef RANDOM_ITERATOR_R123_USE_CXX11_CONSTEXPR
+#define RANDOM_ITERATOR_R123_USE_CXX11_CONSTEXPR __has_feature(cxx_constexpr)
+#endif
+
+#ifndef RANDOM_ITERATOR_R123_USE_CXX11_EXPLICIT_CONVERSIONS
+#define RANDOM_ITERATOR_R123_USE_CXX11_EXPLICIT_CONVERSIONS __has_feature(cxx_explicit_conversions)
+#endif
+
+// With clang-3.0, the apparently simpler:
+//  #define RANDOM_ITERATOR_R123_USE_CXX11_RANDOM __has_include(<random>)
+// dumps core.
+#ifndef RANDOM_ITERATOR_R123_USE_CXX11_RANDOM
+#if __cplusplus>=201103L && __has_include(<random>)
+#define RANDOM_ITERATOR_R123_USE_CXX11_RANDOM 1
+#else
+#define RANDOM_ITERATOR_R123_USE_CXX11_RANDOM 0
+#endif
+#endif
+
+#ifndef RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
+#if __cplusplus>=201103L && __has_include(<type_traits>)
+#define RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS 1
+#else
+#define RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS 0
+#endif
+#endif
+
+#include "gccfeatures.h"
+
+#endif