corsika/detail/framework/random/random_iterator/detail/tbb/tbb_stddef.h

+/*
+    Copyright (c) 2005-2020 Intel Corporation
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+        http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+*/
+
+#ifndef __RANDOM_ITERATOR_TBB_tbb_stddef_H
+#define __RANDOM_ITERATOR_TBB_tbb_stddef_H
+
+// Marketing-driven product version
+#define RANDOM_ITERATOR_TBB_VERSION_MAJOR 2020
+#define RANDOM_ITERATOR_TBB_VERSION_MINOR 3
+
+// Engineering-focused interface version
+#define RANDOM_ITERATOR_TBB_INTERFACE_VERSION 11103
+#define RANDOM_ITERATOR_TBB_INTERFACE_VERSION_MAJOR TBB_INTERFACE_VERSION/1000
+
+// The oldest major interface version still supported
+// To be used in SONAME, manifests, etc.
+#define RANDOM_ITERATOR_TBB_COMPATIBLE_INTERFACE_VERSION 2
+
+#define __RANDOM_ITERATOR_TBB_STRING_AUX(x) #x
+#define __RANDOM_ITERATOR_TBB_STRING(x) __RANDOM_ITERATOR_TBB_STRING_AUX(x)
+
+// We do not need defines below for resource processing on windows
+#if !defined RC_INVOKED
+
+// Define groups for Doxygen documentation
+/**
+ * @defgroup algorithms         Algorithms
+ * @defgroup containers         Containers
+ * @defgroup memory_allocation  Memory Allocation
+ * @defgroup synchronization    Synchronization
+ * @defgroup timing             Timing
+ * @defgroup task_scheduling    Task Scheduling
+ */
+
+// Simple text that is displayed on the main page of Doxygen documentation.
+/**
+ * \mainpage Main Page
+ *
+ * Click the tabs above for information about the
+ * - <a href="./modules.html">Modules</a> (groups of functionality) implemented by the library
+ * - <a href="./annotated.html">Classes</a> provided by the library
+ * - <a href="./files.html">Files</a> constituting the library.
+ * .
+ * Please note that significant part of TBB functionality is implemented in the form of
+ * template functions, descriptions of which are not accessible on the <a href="./annotated.html">Classes</a>
+ * tab. Use <a href="./modules.html">Modules</a> or <a href="./namespacemembers.html">Namespace/Namespace Members</a>
+ * tabs to find them.
+ *
+ * Additional pieces of information can be found here
+ * - \subpage concepts
+ * .
+ */
+
+/** \page concepts TBB concepts
+
+    A concept is a set of requirements to a type, which are necessary and sufficient
+    for the type to model a particular behavior or a set of behaviors. Some concepts
+    are specific to a particular algorithm (e.g. algorithm body), while other ones
+    are common to several algorithms (e.g. range concept).
+
+    All TBB algorithms make use of different classes implementing various concepts.
+    Implementation classes are supplied by the user as type arguments of template
+    parameters and/or as objects passed as function call arguments. The library
+    provides predefined  implementations of some concepts (e.g. several kinds of
+    \ref range_req "ranges"), while other ones must always be implemented by the user.
+
+    TBB defines a set of minimal requirements each concept must conform to. Here is
+    the list of different concepts hyperlinked to the corresponding requirements specifications:
+    - \subpage range_req
+    - \subpage parallel_do_body_req
+    - \subpage parallel_for_body_req
+    - \subpage parallel_reduce_body_req
+    - \subpage parallel_scan_body_req
+    - \subpage parallel_sort_iter_req
+**/
+
+// tbb_config.h should be included the first since it contains macro definitions used in other headers
+#include "tbb_config.h"
+
+#if _MSC_VER >=1400
+    #define __RANDOM_ITERATOR_TBB_EXPORTED_FUNC   __cdecl
+    #define __RANDOM_ITERATOR_TBB_EXPORTED_METHOD __thiscall
+#else
+    #define __RANDOM_ITERATOR_TBB_EXPORTED_FUNC
+    #define __RANDOM_ITERATOR_TBB_EXPORTED_METHOD
+#endif
+
+#if __INTEL_COMPILER || _MSC_VER
+#define __RANDOM_ITERATOR_TBB_NOINLINE(decl) __declspec(noinline) decl
+#elif __GNUC__
+#define __RANDOM_ITERATOR_TBB_NOINLINE(decl) decl __attribute__ ((noinline))
+#else
+#define __RANDOM_ITERATOR_TBB_NOINLINE(decl) decl
+#endif
+
+#if __RANDOM_ITERATOR_TBB_NOEXCEPT_PRESENT
+#define __RANDOM_ITERATOR_TBB_NOEXCEPT(expression) noexcept(expression)
+#else
+#define __RANDOM_ITERATOR_TBB_NOEXCEPT(expression)
+#endif
+
+#include <cstddef>      /* Need size_t and ptrdiff_t */
+
+#if _MSC_VER
+    #define __RANDOM_ITERATOR_TBB_tbb_windef_H
+    #include "internal/_tbb_windef.h"
+    #undef __RANDOM_ITERATOR_TBB_tbb_windef_H
+#endif
+#if !defined(_MSC_VER) || _MSC_VER>=1600
+    #include <stdint.h>
+#endif
+
+//! Type for an assertion handler
+typedef void(*assertion_handler_type)( const char* filename, int line, const char* expression, const char * comment );
+
+#if __RANDOM_ITERATOR_TBBMALLOC_BUILD
+namespace rml { namespace internal {
+ #define __RANDOM_ITERATOR_TBB_ASSERT_RELEASE(predicate,message) ((predicate)?((void)0) : rml::internal::assertion_failure(__FILE__,__LINE__,#predicate,message))
+#else
+namespace random_iterator_tbb {
+ #define __RANDOM_ITERATOR_TBB_ASSERT_RELEASE(predicate,message) ((predicate)?((void)0) : tbb::assertion_failure(__FILE__,__LINE__,#predicate,message))
+#endif
+
+    //! Set assertion handler and return previous value of it.
+    assertion_handler_type __RANDOM_ITERATOR_TBB_EXPORTED_FUNC set_assertion_handler( assertion_handler_type new_handler );
+
+    //! Process an assertion failure.
+    /** Normally called from __RANDOM_ITERATOR_TBB_ASSERT macro.
+        If assertion handler is null, print message for assertion failure and abort.
+        Otherwise call the assertion handler. */
+    void __RANDOM_ITERATOR_TBB_EXPORTED_FUNC assertion_failure( const char* filename, int line, const char* expression, const char* comment );
+
+#if __RANDOM_ITERATOR_TBBMALLOC_BUILD
+}}  // namespace rml::internal
+#else
+} // namespace random_iterator_tbb
+#endif
+
+#if RANDOM_ITERATOR_TBB_USE_ASSERT
+
+    //! Assert that predicate is true.
+    /** If predicate is false, print assertion failure message.
+        If the comment argument is not NULL, it is printed as part of the failure message.
+        The comment argument has no other effect. */
+    #define __RANDOM_ITERATOR_TBB_ASSERT(predicate,message) __RANDOM_ITERATOR_TBB_ASSERT_RELEASE(predicate,message)
+
+    #define __RANDOM_ITERATOR_TBB_ASSERT_EX __RANDOM_ITERATOR_TBB_ASSERT
+
+#else /* !RANDOM_ITERATOR_TBB_USE_ASSERT */
+
+    //! No-op version of __RANDOM_ITERATOR_TBB_ASSERT.
+    #define __RANDOM_ITERATOR_TBB_ASSERT(predicate,comment) ((void)0)
+    //! "Extended" version is useful to suppress warnings if a variable is only used with an assert
+    #define __RANDOM_ITERATOR_TBB_ASSERT_EX(predicate,comment) ((void)(1 && (predicate)))
+
+#endif /* !RANDOM_ITERATOR_TBB_USE_ASSERT */
+
+//! The namespace random_iterator_tbb contains all components of the library.
+namespace random_iterator_tbb {
+
+    namespace internal {
+#if _MSC_VER && _MSC_VER<1600
+        typedef __int8 int8_t;
+        typedef __int16 int16_t;
+        typedef __int32 int32_t;
+        typedef __int64 int64_t;
+        typedef unsigned __int8 uint8_t;
+        typedef unsigned __int16 uint16_t;
+        typedef unsigned __int32 uint32_t;
+        typedef unsigned __int64 uint64_t;
+#else /* Posix */
+        using ::int8_t;
+        using ::int16_t;
+        using ::int32_t;
+        using ::int64_t;
+        using ::uint8_t;
+        using ::uint16_t;
+        using ::uint32_t;
+        using ::uint64_t;
+#endif /* Posix */
+    } // namespace internal
+
+    using std::size_t;
+    using std::ptrdiff_t;
+
+//! The function returns the interface version of the TBB shared library being used.
+/**
+ * The version it returns is determined at runtime, not at compile/link time.
+ * So it can be different than the value of TBB_INTERFACE_VERSION obtained at compile time.
+ */
+extern "C" int __RANDOM_ITERATOR_TBB_EXPORTED_FUNC TBB_runtime_interface_version();
+
+/**
+ * @cond INTERNAL
+ * @brief Identifiers declared inside namespace internal should never be used directly by client code.
+ */
+namespace internal {
+
+//! Compile-time constant that is upper bound on cache line/sector size.
+/** It should be used only in situations where having a compile-time upper
+    bound is more useful than a run-time exact answer.
+    @ingroup memory_allocation */
+const size_t NFS_MaxLineSize = 128;
+
+/** Label for data that may be accessed from different threads, and that may eventually become wrapped
+    in a formal atomic type.
+
+    Note that no problems have yet been observed relating to the definition currently being empty,
+    even if at least "volatile" would seem to be in order to avoid data sometimes temporarily hiding
+    in a register (although "volatile" as a "poor man's atomic" lacks several other features of a proper
+    atomic, some of which are now provided instead through specialized functions).
+
+    Note that usage is intentionally compatible with a definition as qualifier "volatile",
+    both as a way to have the compiler help enforce use of the label and to quickly rule out
+    one potential issue.
+
+    Note however that, with some architecture/compiler combinations, e.g. on IA-64 architecture, "volatile"
+    also has non-portable memory semantics that are needlessly expensive for "relaxed" operations.
+
+    Note that this must only be applied to data that will not change bit patterns when cast to/from
+    an integral type of the same length; tbb::atomic must be used instead for, e.g., floating-point types.
+
+    TODO: apply wherever relevant **/
+#define __RANDOM_ITERATOR_TBB_atomic // intentionally empty, see above
+
+#if __RANDOM_ITERATOR_TBB_OVERRIDE_PRESENT
+#define __RANDOM_ITERATOR_TBB_override override
+#else
+#define __RANDOM_ITERATOR_TBB_override // formal comment only
+#endif
+
+#if __RANDOM_ITERATOR_TBB_CPP17_FALLTHROUGH_PRESENT
+#define __RANDOM_ITERATOR_TBB_fallthrough [[fallthrough]]
+#elif __RANDOM_ITERATOR_TBB_FALLTHROUGH_PRESENT
+#define __RANDOM_ITERATOR_TBB_fallthrough __attribute__ ((fallthrough))
+#else
+#define __RANDOM_ITERATOR_TBB_fallthrough
+#endif
+
+template<class T, size_t S, size_t R>
+struct padded_base : T {
+    char pad[S - R];
+};
+template<class T, size_t S> struct padded_base<T, S, 0> : T {};
+
+//! Pads type T to fill out to a multiple of cache line size.
+template<class T, size_t S = NFS_MaxLineSize>
+struct padded : padded_base<T, S, sizeof(T) % S> {};
+
+//! Extended variant of the standard offsetof macro
+/** The standard offsetof macro is not sufficient for TBB as it can be used for
+    POD-types only. The constant 0x1000 (not NULL) is necessary to appease GCC. **/
+#define __RANDOM_ITERATOR_TBB_offsetof(class_name, member_name) \
+    ((ptrdiff_t)&(reinterpret_cast<class_name*>(0x1000)->member_name) - 0x1000)
+
+//! Returns address of the object containing a member with the given name and address
+#define __RANDOM_ITERATOR_TBB_get_object_ref(class_name, member_name, member_addr) \
+    (*reinterpret_cast<class_name*>((char*)member_addr - __RANDOM_ITERATOR_TBB_offsetof(class_name, member_name)))
+
+//! Throws std::runtime_error with what() returning error_code description prefixed with aux_info
+void __RANDOM_ITERATOR_TBB_EXPORTED_FUNC handle_perror( int error_code, const char* aux_info );
+
+#if RANDOM_ITERATOR_TBB_USE_EXCEPTIONS
+    #define __RANDOM_ITERATOR_TBB_TRY try
+    #define __RANDOM_ITERATOR_TBB_CATCH(e) catch(e)
+    #define __RANDOM_ITERATOR_TBB_THROW(e) throw e
+    #define __RANDOM_ITERATOR_TBB_RETHROW() throw
+#else /* !RANDOM_ITERATOR_TBB_USE_EXCEPTIONS */
+    inline bool __RANDOM_ITERATOR_TBB_false() { return false; }
+    #define __RANDOM_ITERATOR_TBB_TRY
+    #define __RANDOM_ITERATOR_TBB_CATCH(e) if ( tbb::internal::__RANDOM_ITERATOR_TBB_false() )
+    #define __RANDOM_ITERATOR_TBB_THROW(e) tbb::internal::suppress_unused_warning(e)
+    #define __RANDOM_ITERATOR_TBB_RETHROW() ((void)0)
+#endif /* !RANDOM_ITERATOR_TBB_USE_EXCEPTIONS */
+
+//! Report a runtime warning.
+void __RANDOM_ITERATOR_TBB_EXPORTED_FUNC runtime_warning( const char* format, ... );
+
+#if RANDOM_ITERATOR_TBB_USE_ASSERT
+static void* const poisoned_ptr = reinterpret_cast<void*>(-1);
+
+//! Set p to invalid pointer value.
+//  Also works for regular (non-__RANDOM_ITERATOR_TBB_atomic) pointers.
+template<typename T>
+inline void poison_pointer( T* __RANDOM_ITERATOR_TBB_atomic & p ) { p = reinterpret_cast<T*>(poisoned_ptr); }
+
+/** Expected to be used in assertions only, thus no empty form is defined. **/
+template<typename T>
+inline bool is_poisoned( T* p ) { return p == reinterpret_cast<T*>(poisoned_ptr); }
+#else
+template<typename T>
+inline void poison_pointer( T* __RANDOM_ITERATOR_TBB_atomic & ) {/*do nothing*/}
+#endif /* !RANDOM_ITERATOR_TBB_USE_ASSERT */
+
+//! Cast between unrelated pointer types.
+/** This method should be used sparingly as a last resort for dealing with
+    situations that inherently break strict ISO C++ aliasing rules. */
+// T is a pointer type because it will be explicitly provided by the programmer as a template argument;
+// U is a referent type to enable the compiler to check that "ptr" is a pointer, deducing U in the process.
+template<typename T, typename U>
+inline T punned_cast( U* ptr ) {
+    uintptr_t x = reinterpret_cast<uintptr_t>(ptr);
+    return reinterpret_cast<T>(x);
+}
+
+#if __RANDOM_ITERATOR_TBB_DEFAULTED_AND_DELETED_FUNC_PRESENT
+
+//! Base class for types that should not be assigned.
+class no_assign {
+public:
+    void operator=( const no_assign& ) = delete;
+    no_assign( const no_assign& ) = default;
+    no_assign() = default;
+};
+
+//! Base class for types that should not be copied or assigned.
+class no_copy: no_assign {
+public:
+    no_copy( const no_copy& ) = delete;
+    no_copy() = default;
+};
+
+#else /*__RANDOM_ITERATOR_TBB_DEFAULTED_AND_DELETED_FUNC_PRESENT*/
+
+//! Base class for types that should not be assigned.
+class no_assign {
+    // Deny assignment
+    void operator=( const no_assign& );
+public:
+#if __GNUC__
+    //! Explicitly define default construction, because otherwise gcc issues gratuitous warning.
+    no_assign() {}
+#endif /* __GNUC__ */
+};
+
+//! Base class for types that should not be copied or assigned.
+class no_copy: no_assign {
+    //! Deny copy construction
+    no_copy( const no_copy& );
+public:
+    //! Allow default construction
+    no_copy() {}
+};
+
+#endif /*__RANDOM_ITERATOR_TBB_DEFAULTED_AND_DELETED_FUNC_PRESENT*/
+
+#if TBB_DEPRECATED_MUTEX_COPYING
+class mutex_copy_deprecated_and_disabled {};
+#else
+// By default various implementations of mutexes are not copy constructible
+// and not copy assignable.
+class mutex_copy_deprecated_and_disabled : no_copy {};
+#endif
+
+//! A function to check if passed in pointer is aligned on a specific border
+template<typename T>
+inline bool is_aligned(T* pointer, uintptr_t alignment) {
+    return 0==((uintptr_t)pointer & (alignment-1));
+}
+
+//! A function to check if passed integer is a power of 2
+template<typename integer_type>
+inline bool is_power_of_two(integer_type arg) {
+    return arg && (0 == (arg & (arg - 1)));
+}
+
+//! A function to compute arg modulo divisor where divisor is a power of 2.
+template<typename argument_integer_type, typename divisor_integer_type>
+inline argument_integer_type modulo_power_of_two(argument_integer_type arg, divisor_integer_type divisor) {
+    __RANDOM_ITERATOR_TBB_ASSERT( is_power_of_two(divisor), "Divisor should be a power of two" );
+    return (arg & (divisor - 1));
+}
+
+
+//! A function to determine if arg is a power of 2 at least as big as another power of 2.
+// i.e. for strictly positive i and j, with j being a power of 2,
+// determines whether i==j<<k for some nonnegative k (so i==j yields true).
+template<typename argument_integer_type, typename power2_integer_type>
+inline bool is_power_of_two_at_least(argument_integer_type arg, power2_integer_type power2) {
+    __RANDOM_ITERATOR_TBB_ASSERT( is_power_of_two(power2), "Divisor should be a power of two" );
+    return 0 == (arg & (arg - power2));
+}
+
+//! Utility template function to prevent "unused" warnings by various compilers.
+template<typename T1> void suppress_unused_warning( const T1& ) {}
+template<typename T1, typename T2> void suppress_unused_warning( const T1&, const T2& ) {}
+template<typename T1, typename T2, typename T3> void suppress_unused_warning( const T1&, const T2&, const T3& ) {}
+
+// Struct to be used as a version tag for inline functions.
+/** Version tag can be necessary to prevent loader on Linux from using the wrong
+    symbol in debug builds (when inline functions are compiled as out-of-line). **/
+struct version_tag_v3 {};
+
+typedef version_tag_v3 version_tag;
+
+} // internal
+
+//! Dummy type that distinguishes splitting constructor from copy constructor.
+/**
+ * See description of parallel_for and parallel_reduce for example usages.
+ * @ingroup algorithms
+ */
+class split {
+};
+
+//! Type enables transmission of splitting proportion from partitioners to range objects
+/**
+ * In order to make use of such facility Range objects must implement
+ * splitting constructor with this type passed and initialize static
+ * constant boolean field 'is_splittable_in_proportion' with the value
+ * of 'true'
+ */
+class proportional_split: internal::no_assign {
+public:
+    proportional_split(size_t _left = 1, size_t _right = 1) : my_left(_left), my_right(_right) { }
+
+    size_t left() const { return my_left; }
+    size_t right() const { return my_right; }
+
+    // used when range does not support proportional split
+    operator split() const { return split(); }
+
+#if __RANDOM_ITERATOR_TBB_ENABLE_RANGE_FEEDBACK
+    void set_proportion(size_t _left, size_t _right) {
+        my_left = _left;
+        my_right = _right;
+    }
+#endif
+private:
+    size_t my_left, my_right;
+};
+
+} // tbb
+
+// Following is a set of classes and functions typically used in compile-time "metaprogramming".
+// TODO: move all that to a separate header
+
+#if __RANDOM_ITERATOR_TBB_CPP11_SMART_POINTERS_PRESENT
+#include <memory> // for unique_ptr
+#endif
+
+#if __RANDOM_ITERATOR_TBB_CPP11_RVALUE_REF_PRESENT || __RANDOM_ITERATOR_TBB_CPP11_DECLTYPE_PRESENT || _LIBCPP_VERSION
+#include <utility> // for std::move, std::forward, std::declval
+#endif
+
+namespace random_iterator_tbb {
+namespace internal {
+
+#if __RANDOM_ITERATOR_TBB_CPP11_SMART_POINTERS_PRESENT && __RANDOM_ITERATOR_TBB_CPP11_RVALUE_REF_PRESENT && __RANDOM_ITERATOR_TBB_CPP11_VARIADIC_TEMPLATES_PRESENT
+    template<typename T, typename... Args>
+    std::unique_ptr<T> make_unique(Args&&... args) {
+        return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
+    }
+#endif
+
+//! Class for determining type of std::allocator<T>::value_type.
+template<typename T>
+struct allocator_type {
+    typedef T value_type;
+};
+
+#if _MSC_VER
+//! Microsoft std::allocator has non-standard extension that strips const from a type.
+template<typename T>
+struct allocator_type<const T> {
+    typedef T value_type;
+};
+#endif
+
+// Ad-hoc implementation of true_type & false_type
+// Intended strictly for internal use! For public APIs (traits etc), use C++11 analogues.
+template <bool v>
+struct bool_constant {
+    static /*constexpr*/ const bool value = v;
+};
+typedef bool_constant<true> true_type;
+typedef bool_constant<false> false_type;
+
+//! A template to select either 32-bit or 64-bit constant as compile time, depending on machine word size.
+template <unsigned u, unsigned long long ull >
+struct select_size_t_constant {
+    //Explicit cast is needed to avoid compiler warnings about possible truncation.
+    //The value of the right size,   which is selected by ?:, is anyway not truncated or promoted.
+    static const size_t value = (size_t)((sizeof(size_t)==sizeof(u)) ? u : ull);
+};
+
+#if __RANDOM_ITERATOR_TBB_CPP11_RVALUE_REF_PRESENT
+using std::move;
+using std::forward;
+#elif defined(_LIBCPP_NAMESPACE)
+// libc++ defines "pre-C++11 move and forward" similarly to ours; use it to avoid name conflicts in some cases.
+using std::_LIBCPP_NAMESPACE::move;
+using std::_LIBCPP_NAMESPACE::forward;
+#else
+// It is assumed that cv qualifiers, if any, are part of the deduced type.
+template <typename T>
+T& move( T& x ) { return x; }
+template <typename T>
+T& forward( T& x ) { return x; }
+#endif /* __RANDOM_ITERATOR_TBB_CPP11_RVALUE_REF_PRESENT */
+
+// Helper macros to simplify writing templates working with both C++03 and C++11.
+#if __RANDOM_ITERATOR_TBB_CPP11_RVALUE_REF_PRESENT
+#define  __RANDOM_ITERATOR_TBB_FORWARDING_REF(A) A&&
+#else
+// It is assumed that cv qualifiers, if any, are part of a deduced type.
+// Thus this macro should not be used in public interfaces.
+#define  __RANDOM_ITERATOR_TBB_FORWARDING_REF(A) A&
+#endif
+#if __RANDOM_ITERATOR_TBB_CPP11_VARIADIC_TEMPLATES_PRESENT
+#define __RANDOM_ITERATOR_TBB_PARAMETER_PACK ...
+#define __RANDOM_ITERATOR_TBB_PACK_EXPANSION(A) A...
+#else
+#define __RANDOM_ITERATOR_TBB_PARAMETER_PACK
+#define __RANDOM_ITERATOR_TBB_PACK_EXPANSION(A) A
+#endif /* __RANDOM_ITERATOR_TBB_CPP11_VARIADIC_TEMPLATES_PRESENT */
+
+#if __RANDOM_ITERATOR_TBB_CPP11_DECLTYPE_PRESENT
+#if __RANDOM_ITERATOR_TBB_CPP11_DECLVAL_BROKEN
+// Ad-hoc implementation of std::declval
+template <class T> __RANDOM_ITERATOR_TBB_FORWARDING_REF(T) declval() /*noexcept*/;
+#else
+using std::declval;
+#endif
+#endif
+
+template <bool condition>
+struct STATIC_ASSERTION_FAILED;
+
+template <>
+struct STATIC_ASSERTION_FAILED<false> { enum {value=1};};
+
+template<>
+struct STATIC_ASSERTION_FAILED<true>; //intentionally left undefined to cause compile time error
+
+//! @endcond
+}} // namespace random_iterator_tbb::internal
+
+#if __RANDOM_ITERATOR_TBB_STATIC_ASSERT_PRESENT
+#define __RANDOM_ITERATOR_TBB_STATIC_ASSERT(condition,msg) static_assert(condition,msg)
+#else
+//please note condition is intentionally inverted to get a bit more understandable error msg
+#define __RANDOM_ITERATOR_TBB_STATIC_ASSERT_IMPL1(condition,msg,line)       \
+    enum {static_assert_on_line_##line = tbb::internal::STATIC_ASSERTION_FAILED<!(condition)>::value}
+
+#define __RANDOM_ITERATOR_TBB_STATIC_ASSERT_IMPL(condition,msg,line) __RANDOM_ITERATOR_TBB_STATIC_ASSERT_IMPL1(condition,msg,line)
+//! Verify condition, at compile time
+#define __RANDOM_ITERATOR_TBB_STATIC_ASSERT(condition,msg) __RANDOM_ITERATOR_TBB_STATIC_ASSERT_IMPL(condition,msg,__LINE__)
+#endif
+
+#endif /* RC_INVOKED */
+#endif /* __RANDOM_ITERATOR_TBB_tbb_stddef_H */

corsika/detail/framework/random/random_iterator/detail/uint128.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.
+ */
+
+/*
+ * uint128_t.hpp
+ *
+ *  Created on: 13/03/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <cstdint>
+#include <ostream>
+#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include "Macros.h"
+
+namespace random_iterator {
+
+  class uint128_t;
+}
+
+namespace std { // This is probably not a good idea
+  template <>
+  struct is_arithmetic<random_iterator::uint128_t> : std::true_type {};
+  template <>
+  struct is_integral<random_iterator::uint128_t> : std::true_type {};
+  template <>
+  struct is_unsigned<random_iterator::uint128_t> : std::true_type {};
+} // namespace std
+
+namespace random_iterator {
+
+  // Give uint128_t type traits
+
+  class uint128_t {
+
+  public:
+    // Constructors
+    uint128_t() = default;
+
+    uint128_t(uint128_t const& rhs) = default;
+
+    uint128_t(uint128_t&& rhs) = default;
+
+    uint128_t(std::string& s);
+
+    uint128_t(const char* s);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    uint128_t(const T& rhs)
+#ifdef __BIG_ENDIAN__
+        : UPPER(0)
+        , LOWER(rhs)
+#endif
+#ifdef __LITTLE_ENDIAN__
+        : LOWER(rhs)
+        , UPPER(0)
+#endif
+    {
+      if (std::is_signed<T>::value) {
+        if (rhs < T(0)) { UPPER = -1; }
+      }
+    }
+
+    template <typename S, typename T,
+              typename = typename std::enable_if<
+                  std::is_integral<S>::value && std::is_integral<T>::value, void>::type>
+    uint128_t(const S& upper_rhs, const T& lower_rhs)
+#ifdef __BIG_ENDIAN__
+        : UPPER(upper_rhs)
+        , LOWER(lower_rhs)
+#endif
+#ifdef __LITTLE_ENDIAN__
+        : LOWER(lower_rhs)
+        , UPPER(upper_rhs)
+#endif
+    {
+    }
+
+    //  RHS input args only
+
+    // Assignment Operator
+    inline uint128_t& operator=(const uint128_t& rhs) = default;
+
+    inline uint128_t& operator=(uint128_t&& rhs) = default;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator=(const T& rhs) {
+      UPPER = 0;
+
+      if (std::is_signed<T>::value) {
+        if (rhs < T(0)) { UPPER = -1; }
+      }
+
+      LOWER = rhs;
+      return *this;
+    }
+
+    // Typecast Operators
+    inline operator bool() const;
+    inline operator uint8_t() const;
+    inline operator uint16_t() const;
+    inline operator uint32_t() const;
+    inline operator uint64_t() const;
+
+    // Bitwise Operators
+    inline uint128_t operator&(const uint128_t& rhs) const;
+
+    // inline  void export_bits(std::vector<uint8_t> & ret) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator&(const T& rhs) const {
+      return uint128_t(0, LOWER & (uint64_t)rhs);
+    }
+
+    inline uint128_t& operator&=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator&=(const T& rhs) {
+      UPPER = 0;
+      LOWER &= rhs;
+      return *this;
+    }
+
+    inline uint128_t operator|(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator|(const T& rhs) const {
+      return uint128_t(UPPER, LOWER | (uint64_t)rhs);
+    }
+
+    inline uint128_t& operator|=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator|=(const T& rhs) {
+      LOWER |= (uint64_t)rhs;
+      return *this;
+    }
+
+    inline uint128_t operator^(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator^(const T& rhs) const {
+      return uint128_t(UPPER, LOWER ^ (uint64_t)rhs);
+    }
+
+    inline uint128_t& operator^=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator^=(const T& rhs) {
+      LOWER ^= (uint64_t)rhs;
+      return *this;
+    }
+
+    inline uint128_t operator~() const;
+
+    // Bit Shift Operators
+    inline uint128_t operator<<(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator<<(const T& rhs) const {
+      return *this << uint128_t(rhs);
+    }
+
+    inline uint128_t& operator<<=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator<<=(const T& rhs) {
+      *this = *this << uint128_t(rhs);
+      return *this;
+    }
+
+    inline uint128_t operator>>(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator>>(const T& rhs) const {
+      return *this >> uint128_t(rhs);
+    }
+
+    inline uint128_t& operator>>=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator>>=(const T& rhs) {
+      *this = *this >> uint128_t(rhs);
+      return *this;
+    }
+
+    // Logical Operators
+    inline bool operator!() const;
+    inline bool operator&&(const uint128_t& rhs) const;
+    inline bool operator||(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator&&(const T& rhs) const {
+      return static_cast<bool>(*this && rhs);
+    }
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator||(const T& rhs) const {
+      return static_cast<bool>(*this || rhs);
+    }
+
+    // Comparison Operators
+    inline bool operator==(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator==(const T& rhs) const {
+      return (!UPPER && (LOWER == (uint64_t)rhs));
+    }
+
+    inline bool operator!=(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator!=(const T& rhs) const {
+      return (UPPER | (LOWER != (uint64_t)rhs));
+    }
+
+    inline bool operator>(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator>(const T& rhs) const {
+      return (UPPER || (LOWER > (uint64_t)rhs));
+    }
+
+    inline bool operator<(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator<(const T& rhs) const {
+      return (!UPPER) ? (LOWER < (uint64_t)rhs) : false;
+    }
+
+    inline bool operator>=(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator>=(const T& rhs) const {
+      return ((*this > rhs) | (*this == rhs));
+    }
+
+    inline bool operator<=(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline bool operator<=(const T& rhs) const {
+      return ((*this < rhs) | (*this == rhs));
+    }
+
+    // Arithmetic Operators
+    inline uint128_t operator+(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator+(const T& rhs) const {
+      return uint128_t(UPPER + ((LOWER + (uint64_t)rhs) < LOWER), LOWER + (uint64_t)rhs);
+    }
+
+    inline uint128_t& operator+=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator+=(const T& rhs) {
+      return *this += uint128_t(rhs);
+    }
+
+    inline uint128_t operator-(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator-(const T& rhs) const {
+      return uint128_t((uint64_t)(UPPER - ((LOWER - rhs) > LOWER)),
+                       (uint64_t)(LOWER - rhs));
+    }
+
+    inline uint128_t& operator-=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator-=(const T& rhs) {
+      return *this = *this - uint128_t(rhs);
+    }
+
+    inline uint128_t operator*(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator*(const T& rhs) const {
+      return *this * uint128_t(rhs);
+    }
+
+    inline uint128_t& operator*=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator*=(const T& rhs) {
+      return *this = *this * uint128_t(rhs);
+    }
+
+    inline uint128_t operator/(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator/(const T& rhs) const {
+      return *this / uint128_t(rhs);
+    }
+
+    inline uint128_t& operator/=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator/=(const T& rhs) {
+      return *this = *this / uint128_t(rhs);
+    }
+
+    inline uint128_t operator%(const uint128_t& rhs) const;
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t operator%(const T& rhs) const {
+      return *this % uint128_t(rhs);
+    }
+
+    inline uint128_t& operator%=(const uint128_t& rhs);
+
+    template <typename T,
+              typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+    inline uint128_t& operator%=(const T& rhs) {
+      return *this = *this % uint128_t(rhs);
+    }
+
+    // Increment Operator
+    inline uint128_t& operator++();
+    inline uint128_t operator++(int);
+
+    // Decrement Operator
+    inline uint128_t& operator--();
+    inline uint128_t operator--(int);
+
+    // Nothing done since promotion doesn't work here
+    inline uint128_t operator+() const;
+
+    // two's complement
+    inline uint128_t operator-() const;
+
+    // Get private values
+    inline uint64_t upper() const;
+    inline uint64_t lower() const;
+
+    // Get bitsize of value
+    inline uint8_t bits() const;
+
+    // rotate right 64
+    inline uint128_t rotate_right() const;
+
+    // Get string representation of value
+    inline std::string str(uint8_t base = 10, const unsigned int& len = 0) const;
+
+  private:
+    std::pair<uint128_t, uint128_t> divmod(const uint128_t& lhs,
+                                           const uint128_t& rhs) const;
+    void init(const char* s);
+    void _init_hex(const char* s);
+    void _init_dec(const char* s);
+    void _init_oct(const char* s);
+
+#if defined(__powerpc64__) || defined(__x86_64__)
+    static inline uint128_t mul128(uint128_t const x, uint128_t const y) {
+      uint128_t z;
+#ifdef __powerpc64__
+      asm("mulhdu %3\n\t" : "=a"(z.LOWER), "=d"(z.UPPER) : "%0"(x.LOWER), "rm"(y.LOWER));
+#else
+      asm("mulq %3\n\t" : "=a"(z.LOWER), "=d"(z.UPPER) : "%0"(x.LOWER), "rm"(y.LOWER));
+#endif
+      z.UPPER += (x.UPPER * y.LOWER) + (x.LOWER * y.UPPER);
+      return z;
+    }
+#endif
+
+#ifdef __BIG_ENDIAN__
+    uint64_t UPPER, LOWER;
+#endif
+#ifdef __LITTLE_ENDIAN__
+    uint64_t LOWER, UPPER;
+#endif
+  };
+
+  // lhs type T as first arguemnt
+  // If the output is not a bool, casts to type T
+
+  // Bitwise Operators
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator&(const T& lhs, const uint128_t& rhs) {
+    return rhs & lhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator&=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(rhs & lhs);
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator|(const T& lhs, const uint128_t& rhs) {
+    return rhs | lhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator|=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(rhs | lhs);
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator^(const T& lhs, const uint128_t& rhs) {
+    return rhs ^ lhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator^=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(rhs ^ lhs);
+  }
+
+  // Bitshift operators
+  inline uint128_t operator<<(const bool& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const uint8_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const uint16_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const uint32_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const uint64_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const int8_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const int16_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const int32_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator<<(const int64_t& lhs, const uint128_t& rhs);
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator<<=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(uint128_t(lhs) << rhs);
+  }
+
+  inline uint128_t operator>>(const bool& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const uint8_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const uint16_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const uint32_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const uint64_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const int8_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const int16_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const int32_t& lhs, const uint128_t& rhs);
+  inline uint128_t operator>>(const int64_t& lhs, const uint128_t& rhs);
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator>>=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(uint128_t(lhs) >> rhs);
+  }
+
+  // Comparison Operators
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline bool operator==(const T& lhs, const uint128_t& rhs) {
+    return (!rhs.upper() && ((uint64_t)lhs == rhs.lower()));
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline bool operator!=(const T& lhs, const uint128_t& rhs) {
+    return (rhs.upper() | ((uint64_t)lhs != rhs.lower()));
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline bool operator>(const T& lhs, const uint128_t& rhs) {
+    return (!rhs.upper()) && ((uint64_t)lhs > rhs.lower());
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline bool operator<(const T& lhs, const uint128_t& rhs) {
+    if (rhs.upper()) { return true; }
+    return ((uint64_t)lhs < rhs.lower());
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline bool operator>=(const T& lhs, const uint128_t& rhs) {
+    if (rhs.upper()) { return false; }
+    return ((uint64_t)lhs >= rhs.lower());
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline bool operator<=(const T& lhs, const uint128_t& rhs) {
+    if (rhs.upper()) { return true; }
+    return ((uint64_t)lhs <= rhs.lower());
+  }
+
+  // Arithmetic Operators
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator+(const T& lhs, const uint128_t& rhs) {
+    return rhs + lhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator+=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(rhs + lhs);
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator-(const T& lhs, const uint128_t& rhs) {
+    return -(rhs - lhs);
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator-=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(-(rhs - lhs));
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator*(const T& lhs, const uint128_t& rhs) {
+    return rhs * lhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator*=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(rhs * lhs);
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator/(const T& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) / rhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator/=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(uint128_t(lhs) / rhs);
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline uint128_t operator%(const T& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) % rhs;
+  }
+
+  template <typename T,
+            typename = typename std::enable_if<std::is_integral<T>::value, T>::type>
+  inline T& operator%=(T& lhs, const uint128_t& rhs) {
+    return lhs = static_cast<T>(uint128_t(lhs) % rhs);
+  }
+
+  // IO Operator
+  inline std::ostream& operator<<(std::ostream& stream, const uint128_t& rhs);
+
+} // namespace random_iterator
+
+#include "uint128.inl"

corsika/detail/framework/random/random_iterator/detail/uint128.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.
+ */
+
+/*
+ * uint128_t.inl
+ *
+ *  Created on: 13/03/2021
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#pragma once
+
+#include <cstring>
+
+namespace random_iterator {
+
+  const uint128_t uint128_0(0);
+  const uint128_t uint128_1(1);
+
+  inline uint128_t::uint128_t(std::string& s) { init(s.c_str()); }
+
+  inline uint128_t::uint128_t(const char* s) { init(s); }
+
+  inline void uint128_t::init(const char* s) {
+    if (s == NULL || s[0] == 0) {
+      LOWER = UPPER = 0;
+      return;
+    }
+
+    while (*s == ' ') ++s;
+
+    if (std::memcmp(s, "0x", 2) == 0) {
+      _init_hex(&s[2]);
+    } else if (std::memcmp(s, "0o", 2) == 0) {
+      _init_oct(&s[2]);
+    } else {
+      _init_dec(s);
+    }
+  }
+
+  inline void uint128_t::_init_hex(const char* s) {
+    // 2**128 = 0x100000000000000000000000000000000.
+    LOWER = UPPER = 0;
+    int i;
+    for (i = 0; *s && i < 16; ++s, ++i) {
+      if ('0' <= *s && *s <= '9') {
+        LOWER *= 16;
+        LOWER += *s - '0';
+      } else if ('A' <= *s && *s <= 'F') {
+        LOWER *= 16;
+        LOWER += *s + (10 - 'A');
+      } else if ('a' <= *s && *s <= 'f') {
+        LOWER *= 16;
+        LOWER += *s + (10 - 'a');
+      } else {
+        return;
+      }
+    }
+    for (; *s && i < 32; ++s, ++i) {
+      if ('0' <= *s && *s <= '9') {
+        *this *= 16;
+        *this += *s - '0';
+      } else if ('A' <= *s && *s <= 'F') {
+        *this *= 16;
+        *this += *s + (10 - 'A');
+      } else if ('a' <= *s && *s <= 'f') {
+        *this *= 16;
+        *this += *s + (10 - 'a');
+      } else {
+        return;
+      }
+    }
+  }
+
+  inline void uint128_t::_init_dec(const char* s) {
+    // 2**128 = 340282366920938463463374607431768211456.
+    LOWER = UPPER = 0;
+    for (int i = 0; '0' <= *s && *s <= '9' && i < 39; ++s, ++i) {
+      *this *= 10;
+      *this += *s - '0';
+    }
+  }
+
+  inline void uint128_t::_init_oct(const char* s) {
+    // 2**128 = 0o4000000000000000000000000000000000000000000.
+    LOWER = UPPER = 0;
+    for (int i = 0; '0' <= *s && *s <= '7' && i < 43; ++s, ++i) {
+      *this *= 8;
+      *this += *s - '0';
+    }
+  }
+
+  inline uint128_t::operator bool() const { return (bool)(UPPER | LOWER); }
+
+  inline uint128_t::operator uint8_t() const { return (uint8_t)LOWER; }
+
+  inline uint128_t::operator uint16_t() const { return (uint16_t)LOWER; }
+
+  inline uint128_t::operator uint32_t() const { return (uint32_t)LOWER; }
+
+  inline uint128_t::operator uint64_t() const { return (uint64_t)LOWER; }
+
+  inline uint128_t uint128_t::operator&(const uint128_t& rhs) const {
+    return uint128_t(UPPER & rhs.UPPER, LOWER & rhs.LOWER);
+  }
+
+  inline uint128_t& uint128_t::operator&=(const uint128_t& rhs) {
+    UPPER &= rhs.UPPER;
+    LOWER &= rhs.LOWER;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator|(const uint128_t& rhs) const {
+    return uint128_t(UPPER | rhs.UPPER, LOWER | rhs.LOWER);
+  }
+
+  inline uint128_t& uint128_t::operator|=(const uint128_t& rhs) {
+    UPPER |= rhs.UPPER;
+    LOWER |= rhs.LOWER;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator^(const uint128_t& rhs) const {
+    return uint128_t(UPPER ^ rhs.UPPER, LOWER ^ rhs.LOWER);
+  }
+
+  inline uint128_t& uint128_t::operator^=(const uint128_t& rhs) {
+    UPPER ^= rhs.UPPER;
+    LOWER ^= rhs.LOWER;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator~() const { return uint128_t(~UPPER, ~LOWER); }
+
+  inline uint128_t uint128_t::operator<<(const uint128_t& rhs) const {
+    const uint64_t shift = rhs.LOWER;
+    if (((bool)rhs.UPPER) || (shift >= 128)) {
+      return uint128_0;
+    } else if (shift == 64) {
+      return uint128_t(LOWER, 0);
+    } else if (shift == 0) {
+      return *this;
+    } else if (shift < 64) {
+      return uint128_t((UPPER << shift) + (LOWER >> (64 - shift)), LOWER << shift);
+    } else if ((128 > shift) && (shift > 64)) {
+      return uint128_t(LOWER << (shift - 64), 0);
+    } else {
+      return uint128_0;
+    }
+  }
+
+  inline uint128_t& uint128_t::operator<<=(const uint128_t& rhs) {
+    *this = *this << rhs;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator>>(const uint128_t& rhs) const {
+
+    const uint64_t shift = rhs.LOWER;
+    if (((bool)rhs.UPPER) || (shift >= 128)) {
+      return uint128_0;
+    } else if (shift == 64) {
+      return uint128_t(0, UPPER);
+    } else if (shift == 0) {
+      return *this;
+    } else if (shift < 64) {
+      return uint128_t(UPPER >> shift, (UPPER << (64 - shift)) + (LOWER >> shift));
+    } else if ((128 > shift) && (shift > 64)) {
+      return uint128_t(0, (UPPER >> (shift - 64)));
+    } else {
+      return uint128_0;
+    }
+  }
+
+  inline uint128_t& uint128_t::operator>>=(const uint128_t& rhs) {
+    *this = *this >> rhs;
+    return *this;
+  }
+
+  inline bool uint128_t::operator!() const { return !(bool)(UPPER | LOWER); }
+
+  inline bool uint128_t::operator&&(const uint128_t& rhs) const {
+    return ((bool)*this && rhs);
+  }
+
+  inline bool uint128_t::operator||(const uint128_t& rhs) const {
+    return ((bool)*this || rhs);
+  }
+
+  inline bool uint128_t::operator==(const uint128_t& rhs) const {
+    return ((UPPER == rhs.UPPER) && (LOWER == rhs.LOWER));
+  }
+
+  inline bool uint128_t::operator!=(const uint128_t& rhs) const {
+    return ((UPPER != rhs.UPPER) | (LOWER != rhs.LOWER));
+  }
+
+  inline bool uint128_t::operator>(const uint128_t& rhs) const {
+    if (UPPER == rhs.UPPER) { return (LOWER > rhs.LOWER); }
+    return (UPPER > rhs.UPPER);
+  }
+
+  inline bool uint128_t::operator<(const uint128_t& rhs) const {
+    if (UPPER == rhs.UPPER) { return (LOWER < rhs.LOWER); }
+    return (UPPER < rhs.UPPER);
+  }
+
+  inline bool uint128_t::operator>=(const uint128_t& rhs) const {
+    return ((*this > rhs) || (*this == rhs));
+  }
+
+  inline bool uint128_t::operator<=(const uint128_t& rhs) const {
+    return ((*this < rhs) || (*this == rhs));
+  }
+
+  inline uint128_t uint128_t::operator+(const uint128_t& rhs) const {
+    return uint128_t(UPPER + rhs.UPPER + ((LOWER + rhs.LOWER) < LOWER),
+                     LOWER + rhs.LOWER);
+  }
+
+  inline uint128_t& uint128_t::operator+=(const uint128_t& rhs) {
+    UPPER += rhs.UPPER + ((LOWER + rhs.LOWER) < LOWER);
+    LOWER += rhs.LOWER;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator-(const uint128_t& rhs) const {
+    return uint128_t(UPPER - rhs.UPPER - ((LOWER - rhs.LOWER) > LOWER),
+                     LOWER - rhs.LOWER);
+  }
+
+  inline uint128_t& uint128_t::operator-=(const uint128_t& rhs) {
+    *this = *this - rhs;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator*(const uint128_t& rhs) const {
+
+#ifdef __powerpc64__
+    return mul128(*this, rhs);
+#elif __x86_64__
+    return mul128(*this, rhs);
+#else
+
+    // split values into 4 32-bit parts
+    uint64_t top[4] = {UPPER >> 32, UPPER & 0xffffffff, LOWER >> 32, LOWER & 0xffffffff};
+    uint64_t bottom[4] = {rhs.UPPER >> 32, rhs.UPPER & 0xffffffff, rhs.LOWER >> 32,
+                          rhs.LOWER & 0xffffffff};
+    uint64_t products[4][4];
+
+    // multiply each component of the values
+    for (int y = 3; y > -1; y--) {
+      for (int x = 3; x > -1; x--) { products[3 - x][y] = top[x] * bottom[y]; }
+    }
+
+    // first row
+    uint64_t fourth32 = (products[0][3] & 0xffffffff);
+    uint64_t third32 = (products[0][2] & 0xffffffff) + (products[0][3] >> 32);
+    uint64_t second32 = (products[0][1] & 0xffffffff) + (products[0][2] >> 32);
+    uint64_t first32 = (products[0][0] & 0xffffffff) + (products[0][1] >> 32);
+
+    // second row
+    third32 += (products[1][3] & 0xffffffff);
+    second32 += (products[1][2] & 0xffffffff) + (products[1][3] >> 32);
+    first32 += (products[1][1] & 0xffffffff) + (products[1][2] >> 32);
+
+    // third row
+    second32 += (products[2][3] & 0xffffffff);
+    first32 += (products[2][2] & 0xffffffff) + (products[2][3] >> 32);
+
+    // fourth row
+    first32 += (products[3][3] & 0xffffffff);
+
+    // move carry to next digit
+    third32 += fourth32 >> 32;
+    second32 += third32 >> 32;
+    first32 += second32 >> 32;
+
+    // remove carry from current digit
+    fourth32 &= 0xffffffff;
+    third32 &= 0xffffffff;
+    second32 &= 0xffffffff;
+    first32 &= 0xffffffff;
+
+    // combine components
+    return uint128_t((first32 << 32) | second32, (third32 << 32) | fourth32);
+#endif
+  }
+
+  inline uint128_t& uint128_t::operator*=(const uint128_t& rhs) {
+    *this = *this * rhs;
+    return *this;
+  }
+
+  inline std::pair<uint128_t, uint128_t> uint128_t::divmod(const uint128_t& lhs,
+                                                           const uint128_t& rhs) const {
+    // Save some calculations /////////////////////
+    if (rhs == uint128_0) {
+      throw std::domain_error("Error: division or modulus by 0");
+    } else if (rhs == uint128_1) {
+      return std::pair<uint128_t, uint128_t>(lhs, uint128_0);
+    } else if (lhs == rhs) {
+      return std::pair<uint128_t, uint128_t>(uint128_1, uint128_0);
+    } else if ((lhs == uint128_0) || (lhs < rhs)) {
+      return std::pair<uint128_t, uint128_t>(uint128_0, lhs);
+    }
+
+    std::pair<uint128_t, uint128_t> qr(uint128_0, uint128_0);
+    for (uint8_t x = lhs.bits(); x > 0; x--) {
+      qr.first <<= uint128_1;
+      qr.second <<= uint128_1;
+
+      if ((lhs >> (x - 1U)) & 1) { ++qr.second; }
+
+      if (qr.second >= rhs) {
+        qr.second -= rhs;
+        ++qr.first;
+      }
+    }
+    return qr;
+  }
+
+  inline uint128_t uint128_t::operator/(const uint128_t& rhs) const {
+    return divmod(*this, rhs).first;
+  }
+
+  inline uint128_t& uint128_t::operator/=(const uint128_t& rhs) {
+    *this = *this / rhs;
+    return *this;
+  }
+
+  inline uint128_t uint128_t::operator%(const uint128_t& rhs) const {
+    return divmod(*this, rhs).second;
+  }
+
+  inline uint128_t& uint128_t::operator%=(const uint128_t& rhs) {
+    *this = *this % rhs;
+    return *this;
+  }
+
+  inline uint128_t& uint128_t::operator++() { return *this += uint128_1; }
+
+  inline uint128_t uint128_t::operator++(int) {
+    uint128_t temp(*this);
+    ++*this;
+    return temp;
+  }
+
+  inline uint128_t& uint128_t::operator--() { return *this -= uint128_1; }
+
+  inline uint128_t uint128_t::operator--(int) {
+    uint128_t temp(*this);
+    --*this;
+    return temp;
+  }
+
+  inline uint128_t uint128_t::operator+() const { return *this; }
+
+  inline uint128_t uint128_t::operator-() const { return ~*this + uint128_1; }
+
+  inline uint64_t uint128_t::upper() const { return UPPER; }
+
+  inline uint64_t uint128_t::lower() const { return LOWER; }
+
+  inline uint8_t uint128_t::bits() const {
+    uint8_t out = 0;
+    if (UPPER) {
+      out = 64;
+      uint64_t up = UPPER;
+      while (up) {
+        up >>= 1;
+        out++;
+      }
+    } else {
+      uint64_t low = LOWER;
+      while (low) {
+        low >>= 1;
+        out++;
+      }
+    }
+    return out;
+  }
+
+  inline uint128_t uint128_t::rotate_right() const { return uint128_t(LOWER, UPPER); }
+
+  inline std::string uint128_t::str(uint8_t base, const unsigned int& len) const {
+    if ((base < 2) || (base > 16)) {
+      throw std::invalid_argument("Base must be in the range [2, 16]");
+    }
+    std::string out = "";
+    if (!(*this)) {
+      out = "0";
+    } else {
+      std::pair<uint128_t, uint128_t> qr(*this, uint128_0);
+      do {
+        qr = divmod(qr.first, base);
+        out = "0123456789abcdef"[(uint8_t)qr.second] + out;
+      } while (qr.first);
+    }
+    if (out.size() < len) { out = std::string(len - out.size(), '0') + out; }
+    return out;
+  }
+
+  inline uint128_t operator<<(const bool& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const uint8_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const uint16_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const uint32_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const uint64_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const int8_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const int16_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const int32_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator<<(const int64_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) << rhs;
+  }
+
+  inline uint128_t operator>>(const bool& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const uint8_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const uint16_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const uint32_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const uint64_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const int8_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const int16_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const int32_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline uint128_t operator>>(const int64_t& lhs, const uint128_t& rhs) {
+    return uint128_t(lhs) >> rhs;
+  }
+
+  inline std::ostream& operator<<(std::ostream& stream, const uint128_t& rhs) {
+    if (stream.flags() & stream.oct) {
+      stream << rhs.str(8);
+    } else if (stream.flags() & stream.dec) {
+      stream << rhs.str(10);
+    } else if (stream.flags() & stream.hex) {
+      stream << rhs.str(16);
+    }
+    return stream;
+  }
+
+} // namespace random_iterator

corsika/detail/framework/stack/CombinedStack.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/Logging.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  template <template <typename> class TParticleInterfaceA,
+            template <typename> class TParticleInterfaceB, typename TStackIterator>
+  template <typename... TArgs1>
+  inline void CombinedParticleInterface<
+      TParticleInterfaceA, TParticleInterfaceB,
+      TStackIterator>::setParticleData(std::tuple<TArgs1...> const vA) {
+    pi_a_type::setParticleData(vA);
+    pi_b_type::setParticleData();
+  }
+
+  template <template <typename> class TParticleInterfaceA,
+            template <typename> class TParticleInterfaceB, typename TStackIterator>
+  template <typename... TArgs1, typename... TArgs2>
+  inline void CombinedParticleInterface<
+      TParticleInterfaceA, TParticleInterfaceB,
+      TStackIterator>::setParticleData(std::tuple<TArgs1...> const vA,
+                                       std::tuple<TArgs2...> const vB) {
+    pi_a_type::setParticleData(vA);
+    pi_b_type::setParticleData(vB);
+  }
+
+  template <template <typename> class TParticleInterfaceA,
+            template <typename> class TParticleInterfaceB, typename TStackIterator>
+  template <typename... TArgs1>
+  inline void CombinedParticleInterface<
+      TParticleInterfaceA, TParticleInterfaceB,
+      TStackIterator>::setParticleData(pi_a_type& p, std::tuple<TArgs1...> const vA) {
+    // static_assert(MT<I>::has_not, "error");
+    pi_a_type::setParticleData(static_cast<pi_a_type&>(p), vA); // original stack
+    pi_b_type::setParticleData(static_cast<pi_b_type&>(p));     // addon stack
+  }
+
+  template <template <typename> class TParticleInterfaceA,
+            template <typename> class TParticleInterfaceB, typename TStackIterator>
+  template <typename... TArgs1, typename... TArgs2>
+  inline void CombinedParticleInterface<
+      TParticleInterfaceA, TParticleInterfaceB,
+      TStackIterator>::setParticleData(pi_c_type& p, std::tuple<TArgs1...> const vA,
+                                       std::tuple<TArgs2...> const vB) {
+
+    pi_a_type::setParticleData(static_cast<pi_a_type&>(p), vA);
+    pi_b_type::setParticleData(static_cast<pi_b_type&>(p), vB);
+  }
+
+  template <template <typename> class TParticleInterfaceA,
+            template <typename> class TParticleInterfaceB, typename TStackIterator>
+  inline std::string CombinedParticleInterface<TParticleInterfaceA, TParticleInterfaceB,
+                                               TStackIterator>::asString() const {
+    return fmt::format("[[{}][{}]]", pi_a_type::asString(), pi_b_type::asString());
+  }
+
+  template <typename Stack1Impl, typename Stack2Impl>
+  inline void CombinedStackImpl<Stack1Impl, Stack2Impl>::clear() {
+    Stack1Impl::clear();
+    Stack2Impl::clear();
+  }
+  template <typename Stack1Impl, typename Stack2Impl>
+  inline void CombinedStackImpl<Stack1Impl, Stack2Impl>::copy(const unsigned int i1,
+                                                              const unsigned int i2) {
+    if (i1 >= getSize() || i2 >= getSize()) {
+      std::ostringstream err;
+      err << "CombinedStack: trying to access data beyond size of stack!";
+      throw std::runtime_error(err.str());
+    }
+    Stack1Impl::copy(i1, i2);
+    Stack2Impl::copy(i1, i2);
+  }
+
+  template <typename Stack1Impl, typename Stack2Impl>
+  inline void CombinedStackImpl<Stack1Impl, Stack2Impl>::swap(const unsigned int i1,
+                                                              const unsigned int i2) {
+    if (i1 >= getSize() || i2 >= getSize()) {
+      std::ostringstream err;
+      err << "CombinedStack: trying to access data beyond size of stack!";
+      throw std::runtime_error(err.str());
+    }
+    Stack1Impl::swap(i1, i2);
+    Stack2Impl::swap(i1, i2);
+  }
+
+  template <typename Stack1Impl, typename Stack2Impl>
+  inline void CombinedStackImpl<Stack1Impl, Stack2Impl>::incrementSize() {
+    Stack1Impl::incrementSize();
+    Stack2Impl::incrementSize();
+  }
+
+  template <typename Stack1Impl, typename Stack2Impl>
+  inline void CombinedStackImpl<Stack1Impl, Stack2Impl>::decrementSize() {
+    Stack1Impl::decrementSize();
+    Stack2Impl::decrementSize();
+  }
+
+} // namespace corsika

corsika/detail/framework/stack/SecondaryView.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/stack/Stack.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <stdexcept>
+#include <vector>
+
+namespace corsika {
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... Args>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::addSecondary(
+      const Args... v) {
+    CORSIKA_LOG_TRACE("SecondaryView::addSecondary(Args&&)");
+    stack_view_iterator proj = getProjectile(); // make this const
+    return addSecondary(proj, v...);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::begin() {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+
+      if (!isErased(i)) break;
+    }
+    return stack_view_iterator(*this, i + 1);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::const_stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::begin() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!isErased(i)) break;
+    }
+
+    return const_stack_view_iterator(*this, i + 1);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::const_stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::cbegin() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!isErased(i)) break;
+    }
+
+    return const_stack_view_iterator(*this, i + 1);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::last() {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!isErased(getSize() - 1 - i)) break;
+    }
+    return stack_view_iterator(*this, getSize() - 1 - i + 1);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::const_stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::last() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!isErased(getSize() - 1 - i)) break;
+    }
+    return stack_view_iterator(*this, getSize() - 1 - i + 1);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::const_stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::clast() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+
+      if (!isErased(getSize() - 1 - i)) break;
+    }
+    return const_stack_view_iterator(*this, getSize() - 1 - i + 1);
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::swap(
+      stack_view_iterator a, stack_view_iterator b) {
+
+    CORSIKA_LOG_TRACE("View::swap");
+    inner_stack_.swap(getIndexFromIterator(a.getIndex()),
+                      getIndexFromIterator(b.getIndex()));
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::copy(
+      stack_view_iterator a, stack_view_iterator b) {
+
+    CORSIKA_LOG_TRACE("View::copy");
+    inner_stack_.copy(getIndexFromIterator(a.getIndex()),
+                      getIndexFromIterator(b.getIndex()));
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::copy(
+      const_stack_view_iterator a, stack_view_iterator b) {
+
+    CORSIKA_LOG_TRACE("View::copy");
+    inner_stack_.copy(getIndexFromIterator(a.getIndex()),
+                      getIndexFromIterator(b.getIndex()));
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void SecondaryView<TStackDataType, TParticleInterface,
+                            MSecondaryProducer>::erase(stack_view_iterator p) {
+
+    CORSIKA_LOG_TRACE("SecondaryView::Delete");
+    if (isEmpty()) { /*error*/
+      throw std::runtime_error("Stack, cannot delete entry since size is zero");
+    }
+    if (isErased(p.getIndex() - 1)) { /*error*/
+      throw std::runtime_error("Stack, cannot delete entry since already deleted");
+    }
+    inner_stack_.erase(getIndexFromIterator(p.getIndex()));
+    inner_stack_reference_type::nDeleted_++; // also count in SecondaryView
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool SecondaryView<TStackDataType, TParticleInterface,
+                            MSecondaryProducer>::purgeLastIfDeleted() {
+    CORSIKA_LOG_TRACE("SecondaryView::purgeLastIfDeleted");
+    if (!isErased(getSize() - 1))
+      return false; // the last particle is not marked for deletion. Do nothing.
+    inner_stack_.purge(getIndexFromIterator(getSize()));
+    inner_stack_reference_type::nDeleted_--;
+    indices_.pop_back();
+    return true;
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::purge() {
+    unsigned int iStack = 0;
+    unsigned int size = getSize();
+    while (iStack < size) {
+
+      if (isErased(iStack)) {
+        inner_stack_.purge(iStack);
+        indices_.erase(indices_.begin() + iStack);
+      }
+      size = getSize();
+      iStack++;
+    }
+    inner_stack_reference_type::nDeleted_ = 0;
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline std::string SecondaryView<TStackDataType, TParticleInterface,
+                                   MSecondaryProducer>::asString() const {
+    std::string str(fmt::format("size {}\n", getSize()));
+    // we make our own begin/end since we want ALL entries
+    std::string new_line = "     ";
+    for (unsigned int iPart = 0; iPart != getSize(); ++iPart) {
+      const_stack_view_iterator itPart(*this, iPart);
+      str += fmt::format(
+          "{}{}{}", new_line, itPart.asString(),
+          (inner_stack_.deleted_[getIndexFromIterator(itPart.getIndex())] ? " [deleted]"
+                                                                          : ""));
+      new_line = "\n     ";
+    }
+    return str;
+  }
+
+  template <typename TStackDataType, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... Args>
+  inline typename SecondaryView<TStackDataType, TParticleInterface,
+                                MSecondaryProducer>::stack_view_iterator
+  SecondaryView<TStackDataType, TParticleInterface, MSecondaryProducer>::addSecondary(
+      stack_view_iterator& proj, const Args... v) {
+    CORSIKA_LOG_TRACE("SecondaryView::addSecondary(stack_view_iterator&, Args&&)");
+    // make space on stack
+    inner_stack_reference_type::getStackData().incrementSize();
+    inner_stack_.deleted_.push_back(false);
+    // get current number of secondaries on stack
+    const unsigned int idSec = getSize();
+    // determine index on (inner) stack where new particle will be located
+    const unsigned int index = inner_stack_reference_type::getStackData().getSize() - 1;
+    indices_.push_back(index);
+    // NOTE: "+1" is since "0" is special marker here for PROJECTILE, see
+    // getIndexFromIterator
+    auto sec = stack_view_iterator(*this, idSec + 1, proj, v...);
+    MSecondaryProducer<TStackDataType, TParticleInterface>::new_secondary(sec);
+    return sec;
+  }
+
+} // namespace corsika

corsika/detail/framework/stack/Stack.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/core/Logging.hpp>
+#include <corsika/framework/stack/StackIteratorInterface.hpp>
+
+#include <stdexcept>
+#include <string>
+#include <vector>
+#include <utility>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... TArgs>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::clear(
+      TArgs... args) {
+    data_.clear(args...);
+    deleted_ = std::vector<bool>(data_.getSize(), false);
+    nDeleted_ = 0;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::begin() {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[i]) break;
+    }
+    return stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::end() {
+    return stack_iterator_type(*this, getSize());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::last() {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[getSize() - 1 - i]) break;
+    }
+    return stack_iterator_type(*this, getSize() - 1 - i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::begin() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[i]) break;
+    }
+    return const_stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::end() const {
+    return const_stack_iterator_type(*this, getSize());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface,
+                 MSecondaryProducer>::const_stack_iterator_type
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::last() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[getSize() - 1 - i]) break;
+    }
+    return const_stack_iterator_type(*this, getSize() - 1 - i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::cbegin() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[i]) break;
+    }
+    return const_stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::cend() const {
+    return const_stack_iterator_type(*this, getSize());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::clast() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[getSize() - 1 - i]) break;
+    }
+
+    return const_stack_iterator_type(*this, getSize() - 1 - i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::at(unsigned int i) {
+    return stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<
+          StackData, MParticleInterface, MSecondaryProducer>::at(unsigned int i) const {
+    return const_stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::first() {
+    return stack_iterator_type{*this, 0};
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::cfirst() const {
+    return const_stack_iterator_type{*this, 0};
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::getNextParticle() {
+    while (purgeLastIfDeleted()) {}
+    return last();
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... TArgs>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<
+          StackData, MParticleInterface,
+          MSecondaryProducer>::addParticle(const TArgs... v) {
+    data_.incrementSize();
+    deleted_.push_back(false);
+    return stack_iterator_type(*this, getSize() - 1, v...);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::swap(
+      stack_iterator_type a, stack_iterator_type b) {
+    swap(a.getIndex(), b.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::copy(
+      stack_iterator_type a, stack_iterator_type b) {
+    copy(a.getIndex(), b.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::copy(
+      const_stack_iterator_type a, stack_iterator_type b) {
+    data_.copy(a.getIndex(), b.getIndex());
+    if (deleted_[b.getIndex()] && !deleted_[a.getIndex()]) nDeleted_--;
+    if (!deleted_[b.getIndex()] && deleted_[a.getIndex()]) nDeleted_++;
+    deleted_[b.getIndex()] = deleted_[a.getIndex()];
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::erase(
+      stack_iterator_type p) {
+    if (this->isEmpty()) {
+      throw std::runtime_error("Stack, cannot delete entry since size is zero");
+    }
+    if (deleted_[p.getIndex()]) {
+      throw std::runtime_error("Stack, cannot delete entry since already deleted");
+    }
+    this->erase(p.getIndex());
+  }
+
+  /*
+   * delete this particle
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::erase(
+      particle_interface_type p) {
+    this->erase(p.getIterator());
+  }
+
+  /*
+   * check if there are no further non-deleted particles on stack
+   */
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isEmpty() {
+    return getEntries() == 0;
+  }
+
+  /*
+   * check if this particle was already deleted
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      const stack_iterator_type& p) const {
+    return isErased(p.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      const const_stack_iterator_type& p) const {
+    return isErased(p.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      const particle_interface_type& p) const {
+    return isErased(p.getIterator());
+  }
+
+  /*
+   * Function to ultimatively remove the last entry from the stack,
+   * if it was marked as deleted before. If this is not the case,
+   * the function will just return false and do nothing.
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::purgeLastIfDeleted() {
+    if (!deleted_.back())
+      return false; // the last particle is not marked for deletion. Do nothing.
+
+    CORSIKA_LOG_TRACE("stack: purgeLastIfDeleted: yes");
+    data_.decrementSize();
+    nDeleted_--;
+    deleted_.pop_back();
+    return true;
+  }
+
+  /*
+   * Function to ultimatively remove all entries from the stack
+   * marked as deleted.
+   *
+   * Careful: this will re-order the entries on the stack, since
+   * "gaps" in the stack are filled with entries from the back
+   * (copied).
+   */
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::purge() {
+    unsigned int iStackFront = 0;
+    unsigned int iStackBack = getSize() - 1;
+
+    for (unsigned int iDeleted = 0; iDeleted < getErased(); ++iDeleted) {
+      // search first delete entry on stack
+      while (!deleted_[iStackFront]) { iStackFront++; }
+      // search for last non-deleted particle on stack
+      while (deleted_[iStackBack]) { iStackBack--; }
+      // copy entry from iStackBack to iStackFront
+      data_.copy(iStackBack, iStackFront);
+      data_.decrementSize();
+    }
+    deleted_.clear();
+    nDeleted_ = 0;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline unsigned int Stack<StackData, MParticleInterface, MSecondaryProducer>::getSize()
+      const {
+    return data_.getSize();
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline std::string Stack<StackData, MParticleInterface, MSecondaryProducer>::asString()
+      const {
+    std::string str(fmt::format("size {}, entries {}, deleted {} \n", getSize(),
+                                getEntries(), getErased()));
+    // we make our own begin/end since we want ALL entries
+    std::string new_line = "     ";
+    for (unsigned int iPart = 0; iPart != getSize(); ++iPart) {
+      const_stack_iterator_type itPart(*this, iPart);
+      str += fmt::format("{}{}{}", new_line, itPart.asString(),
+                         (deleted_[itPart.getIndex()] ? " [deleted]" : ""));
+      new_line = "\n     ";
+    }
+    return str;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... TArgs>
+  inline typename Stack<StackData, MParticleInterface,
+                        MSecondaryProducer>::stack_iterator_type
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::addSecondary(
+      stack_iterator_type& parent, const TArgs... v) {
+    data_.incrementSize();
+    deleted_.push_back(false);
+    return stack_iterator_type(*this, getSize() - 1, parent, v...);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::swap(
+      unsigned int const a, unsigned int const b) {
+    data_.swap(a, b);
+    std::vector<bool>::swap(deleted_[a], deleted_[b]);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::copy(
+      unsigned int const a, unsigned int const b) {
+    data_.copy(a, b);
+    if (deleted_[b] && !deleted_[a]) nDeleted_--;
+    if (!deleted_[b] && deleted_[a]) nDeleted_++;
+    deleted_[b] = deleted_[a];
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      unsigned int const i) const {
+    if (i >= deleted_.size()) return false;
+    return deleted_.at(i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::erase(
+      unsigned int const i) {
+    deleted_[i] = true;
+    nDeleted_++;
+  }
+
+  /*
+   * will remove from storage the element i. This is a helper
+   * function for SecondaryView.
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::purge(
+      unsigned int i) {
+    unsigned int iStackBack = getSize() - 1;
+    // search for last non-deleted particle on stack
+    while (deleted_[iStackBack]) { iStackBack--; }
+    // copy entry from iStackBack to iStackFront
+    data_.copy(iStackBack, i);
+    if (deleted_[i]) nDeleted_--;
+    deleted_[i] = deleted_[iStackBack];
+    data_.decrementSize();
+    deleted_.pop_back();
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline unsigned int
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::getIndexFromIterator(
+      const unsigned int vI) const {
+    return vI;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename Stack<StackData, MParticleInterface, MSecondaryProducer>::value_type&
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::getStackData() {
+    return data_;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline const typename Stack<StackData, MParticleInterface,
+                              MSecondaryProducer>::value_type&
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::getStackData() const {
+    return data_;
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/COMBoost.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 <cmath>
+
+#include <Eigen/Dense>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  inline COMBoost::COMBoost(FourMomentum const& P4projectile,
+                            HEPMassType const massTarget)
+      : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
+      , rotatedCS_{make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents())} {
+    auto const& pProjectile = P4projectile.getSpaceLikeComponents();
+    auto const pProjNormSquared = pProjectile.getSquaredNorm();
+    auto const pProjNorm = sqrt(pProjNormSquared);
+
+    auto const eProjectile = P4projectile.getTimeLikeComponent();
+    auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
+    auto const s =
+        massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
+
+    auto const sqrtS = sqrt(s);
+    auto const sinhEta = -pProjNorm / sqrtS;
+    auto const coshEta = sqrt(1 + pProjNormSquared / s);
+
+    setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
+
+  inline COMBoost::COMBoost(FourMomentum const& P4projectile,
+                            FourMomentum const& P4target)
+      : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()} {
+
+    // this is the center-of-momentum CM frame
+    auto const pCM =
+        P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents();
+    auto const pCM2 = pCM.getSquaredNorm();
+    auto const pCMnorm = sqrt(pCM2);
+    if (pCMnorm == 0_eV) {
+      // CM is at reset
+      rotatedCS_ = originalCS_;
+    } else {
+      rotatedCS_ = make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents() +
+                                                     P4target.getSpaceLikeComponents());
+    }
+
+    auto const s = (P4projectile + P4target).getNormSqr();
+    auto const sqrtS = sqrt(s);
+    auto const sinhEta = -pCMnorm / sqrtS;
+    auto const coshEta = sqrt(1 + pCM2 / s);
+
+    setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
+
+  inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType const mass)
+      : originalCS_{momentum.getCoordinateSystem()}
+      , rotatedCS_{make_rotationToZ(momentum.getCoordinateSystem(), momentum)} {
+    auto const squaredNorm = momentum.getSquaredNorm();
+    auto const norm = sqrt(squaredNorm);
+    auto const sinhEta = -norm / mass;
+    auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
+    setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
+
+  template <typename FourVector>
+  inline FourVector COMBoost::toCoM(FourVector const& p4) const {
+    auto const pComponents = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+    Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
+    Eigen::Vector2d lab;
+
+    lab << (p4.getTimeLikeComponent() * (1 / 1_GeV)),
+        (eVecRotated(2) * (1 / 1_GeV).magnitude());
+
+    auto const boostedZ = boost_ * lab;
+    auto const E_CoM = boostedZ(0) * 1_GeV;
+
+    eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
+
+    CORSIKA_LOG_TRACE("E0={}, p={}, E0'={}, p'={}", p4.getTimeLikeComponent() / 1_GeV,
+                      eVecRotated(2) * (1 / 1_GeV).magnitude(), E_CoM / 1_GeV, boostedZ);
+
+    return FourVector(E_CoM, MomentumVector(rotatedCS_, eVecRotated));
+  }
+
+  template <typename FourVector>
+  inline FourVector COMBoost::fromCoM(FourVector const& p4) const {
+    auto pCM = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+    auto const Ecm = p4.getTimeLikeComponent();
+
+    Eigen::Vector2d com;
+    com << (Ecm * (1 / 1_GeV)), (pCM.getEigenVector()(2) * (1 / 1_GeV).magnitude());
+
+    CORSIKA_LOG_TRACE("Ecm={} GeV, pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
+                      Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV,
+                      p4.getNorm() / 1_GeV);
+
+    auto const boostedZ = inverseBoost_ * com;
+    auto const E_lab = boostedZ(0) * 1_GeV;
+
+    pCM.getEigenVector()[2] = boostedZ(1) * (1_GeV).magnitude();
+
+    Vector<typename decltype(pCM)::dimension_type> pLab{rotatedCS_, pCM};
+    pLab.rebase(originalCS_);
+
+    CORSIKA_LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
+                      " plab={} GeV (norm={} GeV) "
+                      " GeV), invariant mass = {}",
+                      E_lab / 1_GeV, FourVector{E_lab, pLab}.getNorm() / 1_GeV,
+                      pLab.getComponents(), pLab.getNorm() / 1_GeV);
+
+    return FourVector{E_lab, pLab};
+  }
+
+  inline void COMBoost::setBoost(double const coshEta, double const sinhEta) {
+    boost_ << coshEta, sinhEta, sinhEta, coshEta;
+    inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
+  }
+
+  inline CoordinateSystemPtr const& COMBoost::getRotatedCS() const { return rotatedCS_; }
+
+  inline CoordinateSystemPtr const& COMBoost::getOriginalCS() const {
+    return originalCS_;
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/CorsikaData.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/corsika.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <boost/filesystem/path.hpp>
+
+#include <cstdlib>
+#include <stdexcept>
+#include <string>
+
+namespace corsika {
+
+  inline boost::filesystem::path corsika_data(boost::filesystem::path const& key) {
+    boost::filesystem::path fname =
+        boost::filesystem::path(corsika::CORSIKA_DATA_DIR) / key;
+    // LCOV_EXCL_START, this cannot be easily tested system-independently
+    if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
+      fname = boost::filesystem::path(p) / key;
+    }
+    // LCOV_EXCL_STOP
+    CORSIKA_LOG_INFO("opening data file={}", fname);
+    return fname;
+  }
+} // namespace corsika

corsika/detail/framework/utility/CorsikaFenvFallback.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.
+ */
+
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <cfenv>
+
+extern "C" {
+#warning No enabling/disabling of floating point exceptions - platform needs better implementation
+
+inline int feenableexcept(int /*excepts*/) noexcept { return -1; }
+
+inline int fedisableexcept(int /*excepts*/) noexcept { return -1; }
+}

corsika/detail/framework/utility/CorsikaFenvOSX.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.
+ */
+
+/**
+ * Import public domain code
+ *
+ * Provide portable or fallback versions of feenableexcept() and fedisableexcept()
+ * Exist by default in glibc since version 2.2, but not in the standard
+ * fenv.h / cfenv headers for C 99 or C++ 11
+ *
+ * \author Lukas Nellen
+ * \date 14 Jan 2019
+ *
+ */
+
+#include <cfenv>
+
+// Implementation for OS X on intel X64_86
+// code from
+// https://stackoverflow.com/questions/37819235/how-do-you-enable-floating-point-exceptions-for-clang-in-os-x
+// based on
+// http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c
+
+extern "C" {
+
+inline int feenableexcept(int excepts) noexcept {
+  static fenv_t fenv;
+  int new_excepts = excepts & FE_ALL_EXCEPT;
+  // previous masks
+  int old_excepts;
+
+  if (fegetenv(&fenv)) { return -1; }
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
+
+  // unmask
+  fenv.__control &= ~new_excepts;
+  fenv.__mxcsr &= ~(new_excepts << 7);
+
+  return fesetenv(&fenv) ? -1 : old_excepts;
+}
+
+inline int fedisableexcept(int excepts) noexcept {
+  static fenv_t fenv;
+  int new_excepts = excepts & FE_ALL_EXCEPT;
+  // all previous masks
+  int old_excepts;
+
+  if (fegetenv(&fenv)) { return -1; }
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
+
+  // mask
+  fenv.__control |= new_excepts;
+  fenv.__mxcsr |= new_excepts << 7;
+
+  return fesetenv(&fenv) ? -1 : old_excepts;
+}
+}

corsika/detail/framework/utility/CubicSolver.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/utility/CubicSolver.hpp>
+#include <cmath>
+
+namespace corsika {
+
+  namespace andre {
+
+    //---------------------------------------------------------------------------
+    // solve cubic equation A x^3 + B*x^2 + C*x + D = 0
+    //                        x^3 + a*x^2 + b*x + c = 0
+    // mainly along WolframAlpha formulas
+    inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
+                                                         long double C, long double D,
+                                                         double const epsilon) {
+
+      if (std::abs(A) < epsilon) { return solve_quadratic_real(B, C, epsilon); }
+
+      long double a = B / A;
+      long double b = C / A;
+      long double c = D / A;
+
+      long double a2 = a * a;
+      long double q = (3 * b - a2) / 9;
+      long double r = (a * (9 * b - 2 * a2) - 27 * c) / 54;
+      long double q3 = q * q * q;
+
+      // disc = q**3 + r**2
+      // w:e = r*r exactly
+      long double w = r * r;
+      long double e = std::fma(r, r, -w);
+      // s:t =  q*q exactly
+      long double s = q * q;
+      long double t = std::fma(q, q, -s);
+      // s:t * q + w:e = s*q + w + t*q +e = s*q+w + u:v + e = f + u:v + e
+      long double f = std::fma(s, q, w);
+      long double u = t * q;
+      long double v = std::fma(t, q, -u);
+      // error-free sum f+u. See Knuth, TAOCP vol. 2
+      long double uf = u + f;
+      long double au = uf - u;
+      long double ab = uf - au;
+      au = f - au;
+      ab = u - ab;
+      // sum all terms into final result
+      long double const disc = (((e + uf) + au) + ab) + v;
+
+      CORSIKA_LOG_TRACE("disc={} {}", disc, q3 + r * r);
+
+      if (std::abs(disc) < epsilon) {
+
+        a /= 3;
+        long double const cbrtR = std::cbrt(r);
+        return {double(2 * cbrtR - a), double(-cbrtR - a)}; // 2nd solution is doublet
+
+      } else if (disc > 0) {
+
+        long double const S = std::cbrt(r + std::sqrt(disc));
+        long double const T = std::cbrt(r - std::sqrt(disc));
+        a /= 3;
+        return {double((S + T) - a)}; // plus two imaginary solution
+
+      } else { // disc < 0
+
+        long double t = r / std::sqrt(-q3);
+        if (t < -1) t = -1;
+        if (t > 1) t = 1;
+        t = std::acos(t);
+        a /= 3;
+        q = 2 * std::sqrt(-q);
+        return {double(q * std::cos(t / 3) - a),
+                double(q * std::cos((t + 2 * M_PI) / 3) - a),
+                double(q * std::cos((t + 4 * M_PI) / 3) - a)};
+      }
+    }
+  } // namespace andre
+
+  inline std::vector<double> solve_cubic_depressed_disciminant_real(
+      long double p, long double q, long double const disc, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("p={}, q={}, disc={}", p, q, disc);
+
+    if (std::abs(disc) < epsilon) { // disc==0   multiple roots !
+      if (std::abs(p) < epsilon) {  // tripple root
+        return {0};
+      }
+      // double root, single root
+      CORSIKA_LOG_TRACE("cubic double root");
+      return {double(-3 * q / (2 * p)), double(3 * q / p)};
+    }
+
+    if (std::abs(p) < epsilon) { // p==0  --> x^3 + q = 0
+      return {double(std::cbrt(-q))};
+    }
+
+    if (disc > 0) { // three real roots
+      CORSIKA_LOG_TRACE("cubic three roots");
+      long double const p_third = p / 3;
+      long double const sqrt_minus_p_third = std::sqrt(-p_third);
+
+      long double const arg = std::acos(q / (2 * p_third) / sqrt_minus_p_third) / 3;
+
+      long double constexpr pi = M_PI;
+      return {double(2 * sqrt_minus_p_third * std::cos(arg)),
+              double(2 * sqrt_minus_p_third * std::cos(arg - 2 * pi / 3)),
+              double(2 * sqrt_minus_p_third * std::cos(arg - 4 * pi / 3))};
+    }
+
+    // thus disc < 0  ->  one real root
+    if (p < 0) {
+      CORSIKA_LOG_TRACE("cubic cosh");
+      long double const abs_q = std::abs(q);
+      long double const p_third = p / 3;
+      long double const sqrt_minus_p_third = std::sqrt(-p_third);
+      CORSIKA_LOG_TRACE("sqrt_minus_p_third={}, arcosh({})={}", sqrt_minus_p_third,
+                        -abs_q / (2 * p_third) / sqrt_minus_p_third,
+                        std::acosh(-abs_q / (2 * p_third) / sqrt_minus_p_third));
+      CORSIKA_LOG_TRACE(
+          "complex: {}",
+          -2 * abs_q / q * sqrt_minus_p_third *
+              std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3));
+      double const z =
+          -2 * abs_q / q * sqrt_minus_p_third *
+          std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3);
+      return {z};
+    } else { // p > 0
+      CORSIKA_LOG_TRACE("cubic sinh");
+      long double const p_third = p / 3;
+      long double const sqrt_p_third = std::sqrt(p_third);
+      return {double(-2 * sqrt_p_third *
+                     std::sinh(std::asinh(q / (2 * p_third * sqrt_p_third)) / 3))};
+    }
+  }
+
+  inline std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+                                                        double const epsilon) {
+
+    // thanks!
+    // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+    // long double const disc = -(4 * p * p * p + 27 * q * q);
+    // disc = (p/3)**3 + (q/2)**2
+    long double p_third = p / 3;
+    long double q_half = q / 2;
+    // w:e = (q/2)*(q/2) exactly
+    long double w = q_half * q_half;
+    long double e = std::fma(q_half, q_half, -w);
+    // s:t =  (p/3)*(p/3) exactly
+    long double s = p_third * p_third;
+    long double t = std::fma(p_third, p_third, -s);
+    // s:t * p + w:e = s*p + w + t*p +e = s*p+w + u:v + e = f + u:v + e
+    long double f = std::fma(s, p_third, w);
+    long double u = t * p_third;
+    long double v = std::fma(t, p_third, -u);
+    // error-free sum f+u. See Knuth, TAOCP vol. 2
+    long double a = u + f;
+    long double b = a - u;
+    long double c = a - b;
+    b = f - b;
+    c = u - c;
+    // sum all terms into final result
+    long double const disc = -(((e + a) + b) + c) + v;
+    return solve_cubic_depressed_disciminant_real(p, q, disc, epsilon);
+  }
+
+  /**
+   * Analytical approach. Not very stable in all conditions.
+   */
+  inline std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+                                                       long double c, long double d,
+                                                       double const epsilon) {
+
+    CORSIKA_LOG_TRACE("cubic: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b, c,
+                      d, epsilon, (std::abs(a - 1) < epsilon), (std::abs(b) < epsilon));
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_quadratic_real(b, c, d, epsilon);
+    }
+
+    if ((std::abs(a - 1) < epsilon) &&
+        (std::abs(b) < epsilon)) { // this is a depressed cubic
+      return solve_cubic_depressed_real(c, d, epsilon);
+    }
+
+    // p = (3ac - b^2) / (3a^2) = 3 * ( c/(3*a) - b^2/(9*a^2) )
+    long double b_over_a = b / a;
+    long double const p_third = std::fma(-b_over_a, b_over_a / 9, c / (a * 3));
+
+    // p = std::fma(a * 3, c, -b2) / (3 * a2);
+    // q = (2*b^3 - 9*abc + 27*a^2*d) / (27a^3) = 2 * ( b^3/(27*a^3) - bc/(6*a^2) +
+    // d/(2*a) )
+    long double const q_half_term1 = std::fma(b_over_a, b_over_a / 27, -c / (a * 6));
+    long double const q_half = std::fma(b_over_a, q_half_term1, d / (a * 2));
+
+    std::vector<double> zs = solve_cubic_depressed_real(3 * p_third, 2 * q_half, epsilon);
+    CORSIKA_LOG_TRACE("cubic: solve_depressed={}, b/3a={}", fmt::join(zs, ", "),
+                      b / (a * 3));
+    for (auto& z : zs) {
+      z -= b / (a * 3);
+      double const b1 = z + b / a;
+      double const b0 = b1 * z + c / a;
+      std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+      CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+                        static_pow<3>(z) * a + static_pow<2>(z) * b + z * c + d);
+    }
+    CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(zs, ", "));
+    return zs;
+  }
+
+  template <typename T> // T must be floating point type
+  inline T cubic_function(T x, T a, T b, T c, T d) {
+    T x2 = x * x;
+    return x2 * x * a + x2 * b + x * c + d;
+  }
+
+  template <typename T> // T must be floating point type
+  inline T cubic_function_dfdx(T x, T a, T b, T c) {
+    T x2 = x * x;
+    return x2 * a * 3 + x * b * 2 + c;
+  }
+
+  template <typename T> // T must be floating point type
+  inline T cubic_function_d2fd2x(T x, T a, T b) {
+    return x * a * 6 + b * 2;
+  }
+
+  /**
+   * Iterative approach. https://en.wikipedia.org/wiki/Halley%27s_method
+   *  Halley's method
+   */
+
+  inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
+                                              long double d, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("cubic_iterative: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}",
+                      a, b, c, d, epsilon, (std::abs(a - 1) < epsilon),
+                      (std::abs(b) < epsilon));
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_quadratic_real(b, c, d, epsilon);
+    }
+
+    auto pre_opt = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+    long double x1 = 0; // start value
+
+    if (pre_opt.size()) {
+      x1 = pre_opt[0]; //*std::max_element(pre_opt.begin(), pre_opt.end());
+#ifdef _C8_DEBUG_
+      for (long double test_v : pre_opt) {
+        CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
+                          cubic_function(test_v, a, b, c, d));
+      }
+#endif
+    } else {
+      // this is only if the former solve_cubic_real_analytic would not result
+      // in any solution. We have no test case for this. This is excluded from tests:
+      // LCOV_EXCL_START
+      long double const dist = std::fma(b / a, b / a, -3 * c / a);
+      long double const xinfl = -b / (a * 3);
+
+      x1 = xinfl;
+      long double f_test = cubic_function(xinfl, a, b, c, d);
+
+      if (std::abs(f_test) > epsilon) {
+        if (std::abs(dist) < epsilon) {
+          x1 = xinfl - std::cbrt(f_test);
+        } else if (dist > 0) {
+          if (f_test > 0)
+            x1 = xinfl - 2 / 3 * std::sqrt(dist);
+          else
+            x1 = xinfl + 2 / 3 * std::sqrt(dist);
+        }
+      }
+      // LCOV_EXCL_STOP
+    }
+
+    long double f_x1 = cubic_function(x1, a, b, c, d);
+    long double dx1 = 0;
+
+    int niter = 0;
+    const int maxiter = 100;
+    do {
+      long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
+      long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
+      // if (potential) saddle point... avoid
+      if (std::abs(f_prime_x1) < epsilon) {
+        dx1 = std::cbrt(f_x1);
+      } else {
+        dx1 = f_x1 * f_prime_x1 * 2 / (f_prime_x1 * f_prime_x1 * 2 - f_x1 * f_prime2_x1);
+      }
+      x1 -= dx1;
+      f_x1 = cubic_function(x1, a, b, c, d);
+      CORSIKA_LOG_TRACE(
+          "niter={} x1={:.20f} f_x1={:.20f} f_prime={:.20f} f_prime2={:.20f} dx1={}",
+          niter, x1, f_x1, f_prime_x1, f_prime2_x1,
+          f_x1 * f_prime_x1 / (f_prime_x1 * f_prime_x1 - f_x1 * f_prime2_x1 * 0.5));
+    } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon * 1000) &&
+             (std::abs(dx1) > epsilon));
+
+    CORSIKA_LOG_TRACE("niter={}", niter);
+    if (niter >= maxiter) {
+      CORSIKA_LOG_DEBUG("niter reached max iterations {}", niter);
+      return andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+    }
+
+    CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1);
+
+    double const b1 = x1 + b / a;
+    double const b0 = b1 * x1 + c / a;
+    std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, 1e-3);
+    CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+                      cubic_function(x1, a, b, c, d));
+
+    quad_check.push_back(x1);
+    CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(quad_check, ", "));
+    return quad_check;
+  } // namespace corsika
+
+} // namespace corsika

corsika/detail/framework/utility/LinearSolver.inl

+/*
+ * (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 <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+  inline std::vector<double> solve_linear_real(double a, double b) {
+
+    if (a == 0) {
+      return {}; // no (b!=0), or infinite number (b==0) of solutions....
+    }
+
+    return {-b / a};
+  }
+
+  inline std::vector<std::complex<double>> solve_linear(double a, double b) {
+
+    if (std::abs(a) == 0) {
+      return {}; // no (b!=0), or infinite number (b==0) of solutions....
+    }
+
+    return {std::complex<double>(-b / a, 0)};
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/QuadraticSolver.inl

+/*
+ * (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.
+ */
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+                                                           long double c,
+                                                           double const epsilon) {
+    if (std::abs(a) < epsilon) { return solve_linear(b, c); }
+
+    if (std::abs(c) < epsilon) {
+      std::vector<std::complex<double>> lin_result = solve_linear(a, b);
+      lin_result.push_back({0.});
+      return lin_result;
+    }
+
+    long double const radicant = static_pow<2>(b) - a * c * 4;
+
+    if (radicant < -epsilon) { // two complex solutions
+      double const rpart = -b / 2 * a;
+      double const ipart = std::sqrt(-radicant);
+      return {{rpart, ipart}, {rpart, -ipart}};
+    }
+
+    if (radicant < epsilon) { // just one real solution
+      return {{double(-b / 2 * a), 0}};
+    }
+
+    // two real solutions, use Citardauq formula and actively avoid cancellation in the
+    // denominator
+
+    const long double x1 =
+        2 * c / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+    return {{double(x1), 0}, {double(c / (a * x1)), 0}};
+  }
+
+  inline std::vector<double> solve_quadratic_real(long double a, long double b,
+                                                  long double c, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quadratic: a={} b={} c={}", a, b, c);
+
+    if (std::abs(a) < epsilon) { return solve_linear_real(b, c); }
+    if (std::abs(c) < epsilon) {
+      std::vector<double> lin_result = solve_linear_real(a, b);
+      lin_result.push_back(0.);
+      return lin_result;
+    }
+
+    // long double const radicant = std::pow(b, 2) - a * c * 4;
+    // Thanks!
+    // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+    long double w = a * 4 * c;
+    long double e = std::fma(a * 4, c, -w);
+    long double f = std::fma(b, b, -w);
+    long double radicant = f + e;
+
+    CORSIKA_LOG_TRACE("radicant={} {} ", radicant, b * b - a * c * 4);
+
+    if (std::abs(radicant) < epsilon) { // just one real solution
+      return {double(-b / (2 * a))};
+    }
+
+    if (radicant < 0) { return {}; }
+
+    // two real solutions, use Citardauq formula and actively avoid cancellation in the
+    // denominator
+
+    long double const x1 =
+        c * 2 / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+    CORSIKA_LOG_TRACE("quadratic x1={} x2={}", double(x1), double(c / (a * x1)));
+
+    return {double(x1), double(c / (a * x1))};
+  }
+} // namespace corsika

corsika/detail/framework/utility/QuarticSolver.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/utility/CubicSolver.hpp>
+#include <cmath>
+
+namespace corsika {
+
+  namespace andre {
+
+    inline std::vector<double> solve_quartic_real(long double a, long double b,
+                                                  long double c, long double d,
+                                                  long double e, double const epsilon) {
+
+      if (std::abs(a) < epsilon) { return solve_cubic_real(b, c, d, e, epsilon); }
+
+      b /= a;
+      c /= a;
+      d /= a;
+      e /= a;
+
+      long double a3 = -c;
+      long double b3 = b * d - 4. * e;
+      long double c3 = -b * b * e - d * d + 4. * c * e;
+
+      // cubic resolvent
+      // y^3 − c*y^2 + (bd−4e)*y − b^2*e−d^2+4*c*e = 0
+
+      std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
+      if (!x3.size()) {
+        return {}; // no solution, numeric problem (LCOV_EXCL_LINE)
+      }
+      long double y = x3[0]; // there is always at least one solution
+      // The essence - choosing Y with maximal absolute value.
+      if (x3.size() == 3) {
+        if (std::abs(x3[1]) > std::abs(y)) y = x3[1];
+        if (std::abs(x3[2]) > std::abs(y)) y = x3[2];
+      }
+
+      long double q1, q2, p1, p2;
+      // h1+h2 = y && h1*h2 = e  <=>  h^2 -y*h + e = 0    (h === q)
+
+      long double Det = y * y - 4 * e;
+      CORSIKA_LOG_TRACE("Det={}", Det);
+      if (std::abs(Det) < epsilon) // in other words - D==0
+      {
+        q1 = q2 = y * 0.5;
+        // g1+g2 = b && g1+g2 = c-y   <=>   g^2 - b*g + c-y = 0    (p === g)
+        Det = b * b - 4 * (c - y);
+        CORSIKA_LOG_TRACE("Det={}", Det);
+        if (std::abs(Det) < epsilon) { // in other words - D==0
+          p1 = p2 = b * 0.5;
+        } else {
+          if (Det < 0) return {};
+          long double sqDet = sqrt(Det);
+          p1 = (b + sqDet) * 0.5;
+          p2 = (b - sqDet) * 0.5;
+        }
+      } else {
+        if (Det < 0) return {};
+        long double sqDet1 = sqrt(Det);
+        q1 = (y + sqDet1) * 0.5;
+        q2 = (y - sqDet1) * 0.5;
+        // g1+g2 = b && g1*h2 + g2*h1 = c       ( && g === p )  Krammer
+        p1 = (b * q1 - d) / (q1 - q2);
+        p2 = (d - b * q2) / (q1 - q2);
+      }
+
+      // solving quadratic eqs.  x^2 + p1*x + q1 = 0
+      //                         x^2 + p2*x + q2 = 0
+
+      std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
+      std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
+      if (quad2.size() > 0) {
+        for (auto const val : quad2) quad1.push_back(val);
+      }
+      return quad1;
+    }
+  } // namespace andre
+
+  inline std::vector<double> solve_quartic_depressed_real(long double p, long double q,
+                                                          long double r,
+                                                          double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f},  epsilon={}", p, q, r,
+                      epsilon);
+
+    long double const p2 = static_pow<2>(p);
+    long double const q2 = static_pow<2>(q);
+
+    std::vector<double> const resolve_cubic =
+        solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+    CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(),
+                      fmt::join(resolve_cubic, ", "));
+
+    if (!resolve_cubic.size()) return {};
+
+    long double m = 0;
+    for (auto const& v : resolve_cubic) {
+      CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+                                             v * (p2 / 4 - r) - q2 / 8));
+      if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
+    }
+    CORSIKA_LOG_TRACE("check m={}", m);
+    if (m == 0) { return {0}; }
+    if (m < 0) {
+      // this is a rare numerical instability
+      // first: try analytical solution, second: discard (curved->straight tracking)
+      std::vector<double> const resolve_cubic_analytic =
+          andre::solve_cubic_real_analytic(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+      CORSIKA_LOG_TRACE("andre::resolve_cubic_analytic: N={}, m=[{}]",
+                        resolve_cubic_analytic.size(),
+                        fmt::join(resolve_cubic_analytic, ", "));
+
+      if (!resolve_cubic_analytic.size()) return {};
+
+      for (auto const& v : resolve_cubic_analytic) {
+        CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+                                               v * (p2 / 4 - r) - q2 / 8));
+        if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
+      }
+      CORSIKA_LOG_TRACE("check m={}", m);
+      if (m == 0) { return {0}; }
+      if (m < 0) {
+        return {}; // now we are out of options, cannot solve: curved->straight tracking
+      }
+    }
+    long double const quad_term1 = p / 2 + m;
+    long double const quad_term2 = std::sqrt(2 * m);
+    long double const quad_term3 = q / (2 * quad_term2);
+
+    std::vector<double> z_quad1 =
+        solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, 1e-5);
+    std::vector<double> z_quad2 =
+        solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, 1e-5);
+    for (auto const& z : z_quad2) z_quad1.push_back(z);
+    return z_quad1;
+  }
+
+  inline std::vector<double> solve_quartic_real(long double a, long double b,
+                                                long double c, long double d,
+                                                long double e, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quartic: a={:f}, b={:f}, c={:f}, d={:f}, e={:f}, epsilon={}", a, b,
+                      c, d, e, epsilon);
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_cubic_real(b, c, d, e, epsilon);
+    }
+
+    if ((std::abs(a - 1) < epsilon) &&
+        (std::abs(b) < epsilon)) { // this is a depressed quartic
+      return solve_quartic_depressed_real(c, d, e, epsilon);
+    }
+
+    long double const b2 = static_pow<2>(b);
+    long double const b3 = static_pow<3>(b);
+    long double const b4 = static_pow<4>(b);
+    long double const a2 = static_pow<2>(a);
+    long double const a3 = static_pow<3>(a);
+    long double const a4 = static_pow<4>(a);
+
+    long double const p = (c * a * 8 - b2 * 3) / (a4 * 8);
+    long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8);
+    long double const r =
+        (-b4 * 3 + e * a3 * 256 - b * d * a2 * 64 + b2 * c * a * 16) / (a4 * 256);
+
+    std::vector<double> zs = solve_quartic_depressed_real(p, q, r, epsilon);
+    CORSIKA_LOG_TRACE("quartic: solve_depressed={}, b/4a={}", fmt::join(zs, ", "),
+                      b / (4 * a));
+    for (auto& z : zs) { z -= b / (4 * a); }
+    CORSIKA_LOG_TRACE("quartic: solve_quartic_real returns={}", fmt::join(zs, ", "));
+    return zs;
+  }
+} // namespace corsika

corsika/detail/framework/utility/SaveBoostHistogram.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 <cnpy.hpp>
+
+#include <boost/histogram.hpp>
+#include <boost/filesystem.hpp> // can be changed to std::filesystem if compiler supports it
+
+#include <functional>
+#include <memory>
+#include <numeric>
+#include <utility>
+#include <vector>
+#include <string>
+
+namespace corsika {
+
+  template <class Axes, class Storage>
+  inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
+                        std::string const& filename, bool overwrite) {
+    if (boost::filesystem::exists(filename)) {
+      if (overwrite) {
+        boost::filesystem::remove(filename);
+      } else {
+        using namespace std::literals;
+        throw std::runtime_error(
+            ("save_hist(): "s + filename + " already exists"s).c_str());
+      }
+    }
+
+    unsigned const rank = h.rank();
+
+    std::vector<size_t> axes_dims;
+    axes_dims.reserve(rank);
+
+    auto overflow = std::make_unique<bool[]>(rank);
+    auto underflow = std::make_unique<bool[]>(rank);
+
+    std::vector<char> axis_types;
+    axis_types.reserve(rank);
+
+    for (unsigned int i = 0; i < rank; ++i) {
+      auto const& ax = h.axis(i);
+      unsigned const has_underflow = (ax.options() & 0x01) ? 1 : 0;
+      unsigned const has_overflow = (ax.options() & 0x02) ? 1 : 0;
+
+      underflow[i] = has_underflow;
+      overflow[i] = has_overflow;
+
+      axes_dims.emplace_back(ax.size() + has_underflow + has_overflow);
+
+      if (ax.continuous()) {
+        axis_types.push_back('c');
+        std::vector<double> ax_edges;
+        ax_edges.reserve(ax.size() + 1);
+
+        for (int j = 0; j < ax.size(); ++j) { ax_edges.push_back(ax.bin(j).lower()); }
+
+        ax_edges.push_back(ax.bin(ax.size() - 1).upper());
+
+        cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
+                       ax_edges.data(), {ax_edges.size()}, "a");
+      } else {
+        axis_types.push_back('d');
+        std::vector<int64_t> bins; // we assume that discrete axes have integer bins
+        bins.reserve(ax.size());
+
+        for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
+
+        cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
+                       {bins.size()}, "a");
+      }
+    }
+
+    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, "a");
+
+    auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
+                                                unsigned{1}, std::multiplies<>());
+    auto temp = std::make_unique<float[]>(prod_axis_size);
+
+    assert(prod_axis_size == h.size());
+
+    // reduce multi-dim. to 1-dim, row-major (i.e., last axis index is contiguous in
+    // memory) take special care of underflow bins which have -1 as index and thus need
+    // to be shifted by +1
+
+    for (auto&& x : indexed(h, boost::histogram::coverage::all)) {
+      int p = 0;
+
+      for (unsigned axis_index = 0; axis_index < rank; ++axis_index) {
+        int const offset_underflow = (h.axis(axis_index).options() & 0x01) ? 1 : 0;
+        auto k = x.index(axis_index) + offset_underflow;
+        p = k + p * axes_dims.at(axis_index);
+      }
+
+      temp[p] = *x;
+    }
+
+    cnpy::npz_save(filename, "data", temp.get(), axes_dims, "a");
+    // In Python this array can directly be assigned to a histogram view if that
+    // histogram has its axes correspondingly: hist.view(flow=True)[:] = file['data']
+  } // namespace corsika
+
+} // namespace corsika

corsika/detail/media/BaseExponential.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/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline BaseExponential<TDerived>::BaseExponential(Point const& point,
+                                                    LengthType const referenceHeight,
+                                                    MassDensityType const rho0,
+                                                    LengthType const lambda)
+      : rho0_(rho0)
+      , lambda_(lambda)
+      , invLambda_(1 / lambda)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {}
+
+  template <typename TDerived>
+  inline auto const& BaseExponential<TDerived>::getImplementation() const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+      LengthType const height) const {
+    return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+  }
+
+  template <typename TDerived>
+  inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj, DirectionVector const& axis) const {
+    LengthType const length = traj.getLength();
+    if (length == LengthType::zero()) { return GrammageType::zero(); }
+
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis);
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
+
+    if (uDotA == 0) {
+      return length * rhoStart;
+    } else {
+      return rhoStart * (lambda_ / uDotA) * expm1(uDotA * length * invLambda_);
+    }
+  }
+
+  template <typename TDerived>
+  inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage,
+      DirectionVector const& axis) const {
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis);
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
+
+    if (uDotA == 0) {
+      return grammage / rhoStart;
+    } else {
+      auto const logArg = grammage * invLambda_ * uDotA / rhoStart;
+      if (logArg > -1) {
+        return lambda_ / uDotA * log1p(logArg);
+      } else {
+        return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
+               meter;
+      }
+    }
+  }
+
+} // namespace corsika

corsika/detail/media/BaseTabular.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/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+
+#include <exception>
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline BaseTabular<TDerived>::BaseTabular(
+      Point const& point, LengthType const referenceHeight,
+      std::function<MassDensityType(LengthType)> const& rho, unsigned int const nBins,
+      LengthType const deltaHeight)
+      : nBins_(nBins)
+      , deltaHeight_(deltaHeight)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {
+    density_.resize(nBins_);
+    for (unsigned int bin = 0; bin < nBins; ++bin) {
+      density_[bin] = rho(deltaHeight_ * bin);
+      CORSIKA_LOG_DEBUG("new tabulated atm bin={} h={} rho={}", bin, deltaHeight_ * bin,
+                        density_[bin]);
+    }
+  }
+
+  template <typename TDerived>
+  inline auto const& BaseTabular<TDerived>::getImplementation() const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline MassDensityType BaseTabular<TDerived>::getMassDensity(
+      LengthType const height) const {
+    double const fbin = (height - referenceHeight_) / deltaHeight_;
+    int const bin = int(fbin);
+    if (bin < 0) return MassDensityType::zero();
+    if (bin >= int(nBins_ - 1)) {
+      CORSIKA_LOG_ERROR(
+          "invalid height {} (corrected {}) in BaseTabular atmosphere. Min 0, max {}. If "
+          "max is too low: increase!",
+          height, height - referenceHeight_, nBins_ * deltaHeight_);
+      throw std::runtime_error("invalid height");
+    }
+    return density_[bin] + (fbin - bin) * (density_[bin + 1] - density_[bin]);
+  }
+
+  template <typename TDerived>
+  inline GrammageType BaseTabular<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+
+    Point pCurr = traj.getPosition(0);
+    DirectionVector dCurr = traj.getDirection(0);
+    LengthType height1 = (traj.getPosition(0) - point_).getNorm() - referenceHeight_;
+    LengthType height2 = (traj.getPosition(1) - point_).getNorm() - referenceHeight_;
+
+    LengthType const fullLength = traj.getLength(1);
+
+    int sign = +1; // normal direction
+    if (height1 > height2) {
+      std::swap(height1, height2);
+      pCurr = traj.getPosition(1);
+      dCurr = traj.getDirection(1);
+      sign = -1; // inverted direction
+    }
+
+    double const fbin1 = height1 / deltaHeight_;
+    unsigned int const bin1 = int(fbin1);
+
+    double const fbin2 = height2 / deltaHeight_;
+    unsigned int const bin2 = int(fbin2);
+
+    if (fbin1 == fbin2) { return GrammageType::zero(); }
+
+    if (bin1 >= nBins_ - 1 || bin2 >= nBins_ - 1) {
+      CORSIKA_LOG_ERROR("invalid height {} {} in BaseTabular atmosphere integration",
+                        height1, height2);
+      throw std::runtime_error("invalid height");
+    }
+
+    // interpolated start/end densities
+    MassDensityType const rho1 = getMassDensity(height1 + referenceHeight_);
+    MassDensityType const rho2 = getMassDensity(height2 + referenceHeight_);
+
+    // within first bin
+    if (bin1 == bin2) { return fullLength * (rho2 + rho1) / 2; }
+
+    // inclination of trajectory (local)
+    DirectionVector axis((pCurr - point_).normalized()); // to gravity center
+    double cosTheta = axis.dot(dCurr);
+
+    // distance to next height bin
+    unsigned int bin = bin1;
+    LengthType dD = (deltaHeight_ * (bin + 1) - height1) / cosTheta * sign;
+    LengthType distance = dD;
+
+    GrammageType X = dD * (rho1 + density_[bin + 1]) / 2;
+    double frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength);
+    pCurr = traj.getPosition(frac);
+    dCurr = traj.getDirection(frac);
+
+    for (++bin; bin < bin2; ++bin) {
+      // inclination of trajectory
+      axis = (pCurr - point_).normalized();
+      cosTheta = axis.dot(dCurr);
+      // distance to next height bin
+      dD = deltaHeight_ / cosTheta * sign;
+      distance += dD;
+      GrammageType const dX = dD * (density_[bin] + density_[bin + 1]) / 2;
+      X += dX;
+      frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength);
+      pCurr = traj.getPosition(frac);
+      dCurr = traj.getDirection(frac);
+    }
+
+    // inclination of trajectory
+    axis = ((pCurr - point_).normalized());
+    cosTheta = axis.dot(dCurr);
+    // distance to next height bin
+    dD = (height2 - deltaHeight_ * bin2) / cosTheta * sign;
+    X += dD * (rho2 + density_[bin2]) / 2;
+    distance += dD;
+    return X;
+  }
+
+  template <typename TDerived>
+  inline LengthType BaseTabular<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+
+    if (grammage < GrammageType::zero()) {
+      CORSIKA_LOG_ERROR("cannot integrate negative grammage");
+      throw std::runtime_error("negative grammage error");
+    }
+    LengthType const height = (traj.getPosition(0) - point_).getNorm() - referenceHeight_;
+
+    double const fbin = height / deltaHeight_;
+    int bin = int(fbin);
+
+    if (bin >= int(nBins_ - 1)) {
+      CORSIKA_LOG_ERROR("invalid height {} in BaseTabular atmosphere integration",
+                        height);
+      throw std::runtime_error("invalid height");
+    }
+
+    // interpolated start/end densities
+    MassDensityType const rho = getMassDensity(height + referenceHeight_);
+
+    // inclination of trajectory
+    Point pCurr = traj.getPosition(0);
+    DirectionVector dCurr = traj.getDirection(0);
+    DirectionVector axis((pCurr - point_).normalized());
+    double cosTheta = axis.dot(dCurr);
+    int sign = +1; // height increasing along traj
+    if (cosTheta < 0) {
+      cosTheta = -cosTheta; // absolute value only
+      sign = -1;            // height decreasing along traj
+    }
+
+    // height -> distance
+    LengthType const deltaDistance = deltaHeight_ / cosTheta;
+
+    // start with 0 g/cm2
+    GrammageType X(GrammageType::zero());
+    LengthType distance(LengthType::zero());
+
+    // within first bin
+    distance =
+        (sign > 0 ? deltaDistance * (bin + 1 - fbin) : deltaDistance * (fbin - bin));
+    GrammageType binGrammage = (sign > 0 ? distance * (rho + density_[bin + 1]) / 2
+                                         : distance * (rho + density_[bin]) / 2);
+    if (X + binGrammage > grammage) {
+      double const binFraction = (grammage - X) / binGrammage;
+      return distance * binFraction;
+    }
+    X += binGrammage;
+
+    // the following bins (along trajectory)
+    for (bin += sign; bin < int(nBins_) && bin >= 0; bin += sign) {
+
+      binGrammage = deltaDistance * (density_[bin] + density_[bin + 1]) / 2;
+      if (X + binGrammage > grammage) {
+        double const binFraction = (grammage - X) / binGrammage;
+        return distance + deltaDistance * binFraction;
+      }
+      X += binGrammage;
+      distance += deltaDistance;
+    }
+    return std::numeric_limits<double>::infinity() * meter;
+  }
+
+} // namespace corsika

corsika/detail/media/CORSIKA7Atmospheres.inl

+/*
+ * (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
+
+namespace corsika {
+
+  template <typename TEnvironmentInterface, template <typename> typename TExtraEnv,
+            typename TEnvironment, typename... TArgs>
+  void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
+                                Point const& center, TArgs... args) {
+
+    // construct the atmosphere builder
+    auto builder = make_layered_spherical_atmosphere_builder<
+        TEnvironmentInterface, TExtraEnv>::create(center, constants::EarthRadius::Mean,
+                                                  std::forward<TArgs>(args)...);
+
+    builder.setNuclearComposition(standardAirComposition);
+
+    // add the standard atmosphere layers
+    auto const params = atmosphereParameterList[static_cast<uint8_t>(atmId)];
+    for (int i = 0; i < 4; ++i) {
+      builder.addExponentialLayer(params[i].offset, params[i].scaleHeight,
+                                  params[i].altitude);
+    }
+    builder.addLinearLayer(params[4].offset, params[4].scaleHeight, params[4].altitude);
+
+    // and assemble the environment
+    builder.assemble(env);
+  }
+
+} // namespace corsika

corsika/detail/media/Environment.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/media/Environment.hpp>
+
+namespace corsika {
+
+  template <typename IEnvironmentModel>
+  inline Environment<IEnvironmentModel>::Environment()
+      : coordinateSystem_{get_root_CoordinateSystem()}
+      , universe_(std::make_unique<BaseNodeType>(
+            std::make_unique<Universe>(coordinateSystem_))) {}
+
+  template <typename IEnvironmentModel>
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr&
+  Environment<IEnvironmentModel>::getUniverse() {
+    return universe_;
+  }
+
+  template <typename IEnvironmentModel>
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr const&
+  Environment<IEnvironmentModel>::getUniverse() const {
+    return universe_;
+  }
+
+  template <typename IEnvironmentModel>
+  inline CoordinateSystemPtr const& Environment<IEnvironmentModel>::getCoordinateSystem()
+      const {
+    return coordinateSystem_;
+  }
+
+  // factory method for creation of VolumeTreeNodes
+  template <typename IEnvironmentModel>
+  template <class TVolumeType, typename... TVolumeArgs>
+  std::unique_ptr<VolumeTreeNode<IEnvironmentModel> > inline Environment<
+      IEnvironmentModel>::createNode(TVolumeArgs&&... args) {
+    static_assert(std::is_base_of_v<IVolume, TVolumeType>,
+                  "unusable type provided, needs to be derived from "
+                  "\"Volume\"");
+
+    return std::make_unique<BaseNodeType>(
+        std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
+  }
+
+  template <typename IEnvironmentModel>
+  std::set<Code> const get_all_elements_in_universe(
+      Environment<IEnvironmentModel> const& env) {
+
+    auto const& universe = *(env.getUniverse());
+    auto const allElementsInUniverse = std::invoke([&]() {
+      std::set<Code> allElementsInUniverse;
+      auto collectElements = [&](auto& vtn) {
+        if (vtn.hasModelProperties()) {
+          auto const& comp =
+              vtn.getModelProperties().getNuclearComposition().getComponents();
+          for (auto const c : comp) allElementsInUniverse.insert(c);
+        }
+      };
+      universe.walk(collectElements);
+      return allElementsInUniverse;
+    });
+    return allElementsInUniverse;
+  }
+
+} // namespace corsika

corsika/detail/media/ExponentialRefractiveIndex.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/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
+      double const n0, InverseLengthType const lambda, Point const center,
+      LengthType const radius, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , n0_(n0)
+      , lambda_(lambda)
+      , center_(center)
+      , radius_(radius) {}
+
+  template <typename T>
+  inline double ExponentialRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
+  }
+
+} // namespace corsika