diff --git a/rng/philox/philox.hpp b/rng/philox/philox.hpp
index 51cdbe2..01cfcae 100644
--- a/rng/philox/philox.hpp
+++ b/rng/philox/philox.hpp
@@ -38,8 +38,46 @@ namespace detail {
     {
         static constexpr std::size_t value = Is * 2;
     };
+
+    template <class _CharT, class _Traits>
+    class save_stream_flags
+    {
+        typedef ::std::basic_ios<_CharT, _Traits> stream_type;
+
+    public:
+        save_stream_flags(const save_stream_flags&) = delete;
+        save_stream_flags&
+        operator=(const save_stream_flags&) = delete;
+
+        explicit save_stream_flags(stream_type& stream)
+            : stream_(stream), fmtflags_(stream.flags()), fill_(stream.fill())
+        {
+        }
+        ~save_stream_flags()
+        {
+            stream_.flags(fmtflags_);
+            stream_.fill(fill_);
+        }
+
+    private:
+        typename stream_type::fmtflags fmtflags_;
+        stream_type& stream_;
+        _CharT fill_;
+    };
 } // namespace detail
 
+template <typename UIntType, std::size_t w, std::size_t n, std::size_t r, UIntType... consts>
+struct philox_engine;
+
+template<class CharT, class Traits, typename UIntType_, std::size_t w_, std::size_t n_, std::size_t r_, UIntType_... consts_>
+std::basic_ostream<CharT, Traits>&
+operator<<(std::basic_ostream<CharT, Traits>& os, const philox_engine<UIntType_, w_, n_, r_, consts_...>& x);
+
+template<class CharT, class Traits, typename UIntType_, std::size_t w_, std::size_t n_, std::size_t r_, UIntType_... consts_>
+std::basic_istream<CharT, Traits>&
+operator>>(std::basic_istream<CharT, Traits>& is, philox_engine<UIntType_, w_, n_, r_, consts_...>& x);
+
+
 template <typename UIntType, std::size_t w, std::size_t n, std::size_t r, UIntType... consts>
 struct philox_engine
 {
@@ -73,7 +111,7 @@ struct philox_engine
     static constexpr std::array<result_type, array_size> multipliers = extract_elements(even_indices_sequence{});
     static constexpr std::array<result_type, array_size> round_consts = extract_elements(odd_indices_sequence{});
     static constexpr result_type min() { return 0; }
-    static constexpr result_type max() { return max_impl(); }
+    static constexpr result_type max() { return result_mask; }
     static constexpr result_type default_seed = 20111115u;
     // constructors and seeding functions
     philox_engine() : philox_engine(default_seed) {}
@@ -171,13 +209,13 @@ struct philox_engine
     }
 
     // inserters and extractors
-    template<class charT, class traits>
-    friend std::basic_ostream<charT, traits>&
-    operator<<(std::basic_ostream<charT, traits>& os, const philox_engine& x);
+    template<class CharT, class Traits, typename UIntType_, std::size_t w_, std::size_t n_, std::size_t r_, UIntType_... consts_>
+    friend std::basic_ostream<CharT, Traits>&
+    operator<<(std::basic_ostream<CharT, Traits>& os, const philox_engine<UIntType_, w_, n_, r_, consts_...>& x);
 
-    template<class charT, class traits>
-    friend std::basic_istream<charT, traits>&
-    operator>>(std::basic_istream<charT, traits>& is, philox_engine& x);
+    template<class CharT, class Traits, typename UIntType_, std::size_t w_, std::size_t n_, std::size_t r_, UIntType_... consts_>
+    friend std::basic_istream<CharT, Traits>&
+    operator>>(std::basic_istream<CharT, Traits>& is, philox_engine<UIntType_, w_, n_, r_, consts_...>& x);
 
 
 private: // utilities
@@ -194,8 +232,8 @@ struct philox_engine
     using counter_type = std::tuple_element_t<get_log_index(w), uint_types>;
     using promotion_type = std::tuple_element_t<get_log_index(w), promotion_types>;
 
-    static constexpr counter_type counter_mask = ~counter_type(0) >> (sizeof(counter_type) * CHAR_BIT - w);
-    static constexpr result_type result_mask = ~result_type(0) >> (sizeof(result_type) * CHAR_BIT - w);
+    static constexpr counter_type counter_mask = static_cast<counter_type>(~counter_type(0)) >> (std::numeric_limits<counter_type>::digits - w);
+    static constexpr result_type result_mask = static_cast<result_type>(~result_type(0)) >> (std::numeric_limits<result_type>::digits - w);
 
 
 private: // functions
@@ -274,13 +312,6 @@ struct philox_engine
         state_i = n - 1;
     }
 
-    static constexpr result_type max_impl()
-    {
-        return w == std::numeric_limits<result_type>::digits
-            ? std::numeric_limits<result_type>::digits - 1
-            : (result_type(1) << w) - 1;
-    }
-
 public: // state
     std::array<counter_type, n> x;
     std::array<result_type, array_size> k;
@@ -289,6 +320,58 @@ struct philox_engine
     std::uint32_t state_i;
 };
 
+template<class CharT, class Traits, typename UIntType, std::size_t w, std::size_t n, std::size_t r, UIntType... consts>
+std::basic_ostream<CharT, Traits>&
+operator<<(std::basic_ostream<CharT, Traits>& os, const philox_engine<UIntType, w, n, r, consts...>& engine) {
+
+    detail::save_stream_flags<CharT, Traits> flags(os);
+
+    os.setf(std::ios_base::dec | std::ios_base::left);
+    CharT sp = os.widen(' ');
+    os.fill(sp);
+
+    for (std::size_t i = 0; i < n / 2; ++i) {
+        os << engine.k[i] << sp;
+    }
+    for (std::size_t i = 0; i < n; ++i) {
+        os << (UIntType)engine.x[i] << sp;
+    }
+    for (std::size_t i = 0; i < n; ++i) {
+        os << engine.y[i] << sp;
+    }
+    os << engine.state_i;
+    return os;
+}
+
+template<class CharT, class Traits, typename UIntType, std::size_t w, std::size_t n, std::size_t r, UIntType... consts>
+std::basic_istream<CharT, Traits>&
+operator>>(std::basic_istream<CharT, Traits>& is, philox_engine<UIntType, w, n, r, consts...>& engine) {
+    detail::save_stream_flags<CharT, Traits> flags(is);
+
+    is.setf(std::ios_base::dec | std::ios_base::skipws);
+
+    // need a check for the different types?
+    UIntType tmp[5 * n / 2 + 1];
+
+    for (std::size_t i = 0; i < 5 * n / 2 + 1; ++i) {
+        is >> tmp[i];
+    }
+    if(!is.fail()) {
+        std::size_t j = 0;
+        for (std::size_t i = 0; i < n / 2; ++i) {
+            engine.k[i] = tmp[j++];
+        }
+        for (std::size_t i = 0; i < n; ++i) {
+            engine.x[i] = tmp[j++];
+        }
+        for (std::size_t i = 0; i < n; ++i) {
+            engine.y[i] = tmp[j++];
+        }
+        engine.state_i = tmp[j]; // do we change state_i or set it to n - 1 ?
+    }
+    return is;
+}
+
 using philox4x32 = philox_engine<std::uint_fast32_t, 32, 4, 10, 0xCD9E8D57, 0x9E3779B9, 0xD2511F53, 0xBB67AE85>;
 using philox4x64 = philox_engine<std::uint_fast64_t, 64, 4, 10, 0xCA5A826395121157, 0x9E3779B97F4A7C15, 0xD2E7470EE14C6C93, 0xBB67AE8584CAA73B>;
 
diff --git a/rng/philox/test.cpp b/rng/philox/test.cpp
new file mode 100644
index 0000000..ba0337e
--- /dev/null
+++ b/rng/philox/test.cpp
@@ -0,0 +1,308 @@
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <limits.h>
+#include <random>
+
+#include "philox.hpp"
+
+using philox2x32_w5 = std::philox_engine<std::uint_fast32_t, 5, 2, 10, 0xD256D193, 0x0>;
+using philox2x32_w30 = std::philox_engine<std::uint_fast32_t, 30, 2, 10, 0xD256D193, 0x0>;
+using philox2x64_w15 = std::philox_engine<std::uint_fast64_t, 15, 2, 10, 0xD2B74407B1CE6E93, 0x0>;
+using philox2x64_w49 = std::philox_engine<std::uint_fast64_t, 49, 2, 10, 0xD2B74407B1CE6E93, 0x0>;
+
+// Test the conformance of the implementation with the ISO C++ standard
+template <typename Engine>
+void conformance_test() {
+    Engine engine;
+    for(int i = 0; i < 9999; i++) {
+        engine();
+    }
+    typename Engine::result_type reference;
+    if(std::is_same_v<Engine, std::philox4x32>) {
+        reference = 1955073260;
+    }
+    else {
+        reference = 3409172418970261260;
+    }
+    if(engine() == reference) {
+        std::cout << __PRETTY_FUNCTION__ << " passed" << std::endl;
+    } else {
+        std::cout << __PRETTY_FUNCTION__ << " failed" << std::endl;
+    }
+}
+
+// Test public API
+template <typename Engine>
+void api_test() {
+    {
+        Engine engine;
+        engine.seed();
+    }
+    {
+        Engine engine(1);
+        engine.seed(1);
+    }
+    {
+        std::seed_seq s;
+        Engine engine(s);
+        engine.seed(s);
+    }
+    {
+        Engine engine;
+        Engine engine2;
+        if(!(engine == engine2) || (engine != engine2)) {
+            std::cout << __PRETTY_FUNCTION__ << " failed !=, == for the same engines" << std::endl;
+            return;
+        }
+        engine2.seed(42);
+        if((engine == engine2) || !(engine != engine2)) {
+            std::cout << __PRETTY_FUNCTION__ << " failed !=, == for the different engines" << std::endl;
+            return;
+        }
+    }
+    {
+        std::ostringstream os;
+        Engine engine;
+        os << engine << std::endl;
+        Engine engine2;
+        engine2();
+        std::istringstream in(os.str());
+        in >> engine2;
+        if(engine != engine2) {
+            std::cout << __PRETTY_FUNCTION__ << " failed for >> << operators" << std::endl;
+            return;
+        }
+    }
+    {
+        Engine engine;
+        engine.min();
+        engine.max();
+    }
+    std::cout << __PRETTY_FUNCTION__ << " passed" << std::endl;
+}
+
+template <typename Engine>
+void seed_test() {
+    for(int i = 1; i < 5; i++) { // make sure that the state is reset properly for all idx positions
+        Engine engine;
+        typename Engine::result_type res;
+        for(int j = 0; j < i - 1; j++) {
+            engine();
+        }
+        res = engine();
+        engine.seed();
+        for(int j = 0; j < i - 1; j++) {
+            engine();
+        }
+        if(res != engine()) {
+            std::cout << __PRETTY_FUNCTION__ << " failed while generating " << i  << " elements" << std::endl;
+        }
+    }
+    std::cout << __PRETTY_FUNCTION__ << " passed" << std::endl;
+}
+
+template <typename Engine>
+void discard_test() {
+    {
+        constexpr size_t n = 10; // arbitrary length we want to check
+        typename Engine::result_type reference[n];
+        Engine engine;
+        for(int i = 0; i < n; i++) {
+            reference[i] = engine();
+        }
+        for(int i = 0; i < n; i++) {
+            engine.seed();
+            engine.discard(i);
+            for(size_t j = i; j < n; j++) {
+                if(reference[j] != engine()) {
+                    std::cout << __PRETTY_FUNCTION__ << " failed with error in element " << j << " discard " << i << std::endl;
+                    break;
+                }
+            }
+        }
+        std::cout << __PRETTY_FUNCTION__ << " passed step 1 discard from the intial state" << std::endl;
+
+        for(int i = 1; i < n; i++) {
+            for(int j = 1; j < i; j++) {
+                engine.seed();
+                for(size_t k = 0; k < i - j; k++) {
+                    engine();
+                }
+                engine.discard(j);
+                if(reference[i] != engine()) {
+                    std::cout << __PRETTY_FUNCTION__ << " failed on step " << i << " " << j << std::endl;
+                    break;
+                }
+            }
+        }
+        std::cout << __PRETTY_FUNCTION__ << " passed step 2 discard after generation" << std::endl;
+    }
+}
+
+template <typename Engine>
+void set_counter_conformance_test() {
+    Engine engine;
+    std::array<typename Engine::result_type, Engine::word_count> counter;
+    for(int i = 0; i < Engine::word_count - 1; i++) {
+        counter[i] = 0;
+    }
+    
+    counter[Engine::word_count - 1] = 2499; // to get 10'000 element
+    engine.set_counter(counter);
+
+    for(int i = 0; i < Engine::word_count - 1; i++) {
+        engine();
+    }
+
+    typename Engine::result_type reference;
+    if(std::is_same_v<Engine, std::philox4x32>) {
+        reference = 1955073260;
+    }
+    else {
+        reference = 3409172418970261260;
+    }
+    if(engine() == reference) {
+        std::cout << __PRETTY_FUNCTION__ << " passed" << std::endl;
+    } else {
+        std::cout << __PRETTY_FUNCTION__ << " failed" << std::endl;
+    }
+}
+
+template <typename Engine>
+void skip_test() {
+    using T = typename Engine::result_type;
+    for(T i = 1; i <= Engine::word_count + 1; i++) {
+        Engine engine1;
+        std::array<T, Engine::word_count> counter = {0};
+        counter[Engine::word_count - 1] = i / Engine::word_count;
+        engine1.set_counter(counter);
+        for(T j = 0; j < i % Engine::word_count; j++) {
+            engine1();
+        }
+
+        Engine engine2;
+        engine2.discard(i);
+
+        if(engine1() != engine2()) {
+            std::cout << __PRETTY_FUNCTION__ << " failed for " << i << " skip" << std::endl;
+            return;
+        }
+    }
+    std::cout << __PRETTY_FUNCTION__ << " passed" << std::endl;
+}
+
+template <typename Engine>
+void counter_overflow_test() {
+    using T = typename Engine::result_type;
+    Engine engine1;
+    std::array<T, Engine::word_count> counter;
+    for(int i = 0; i < Engine::word_count; i++) {
+        counter[i] = std::numeric_limits<T>::max();
+    }
+
+    engine1.set_counter(counter);
+    for(int i = 0; i < Engine::word_count; i++) {
+        engine1();
+    } // all counters overflowed == start from 0 0 0 0
+
+    Engine engine2;
+
+    if(engine1() == engine2()) {
+        std::cout << __PRETTY_FUNCTION__ << " passed" << std::endl;
+    } else {
+        std::cout << __PRETTY_FUNCTION__ << " failed" << std::endl;
+    }
+}
+
+template <typename Engine>
+void discard_overflow_test() {
+    using T = typename Engine::result_type;
+    for (int overflow_position = 0; overflow_position < Engine::word_count - 1; overflow_position++) {
+        Engine engine1;
+        std::array<T, Engine::word_count> counter = {0};
+
+        int raw_counter_position = (Engine::word_count - overflow_position - 2) % Engine::word_count;
+        std::cout << "Testing discard overflow for position " << raw_counter_position << std::endl;
+        counter[raw_counter_position] = 1;
+
+        engine1.set_counter(counter);
+
+        Engine engine2;
+
+        std::array<T, Engine::word_count> counter2 = {0};
+        for (int i = Engine::word_count - overflow_position - 1; i < Engine::word_count - 1; i++) {
+            counter2[i] = std::numeric_limits<T>::max();
+        }
+
+        engine2.set_counter(counter2);
+
+        for (int i = 0; i < Engine::word_count; i++) {
+            engine2();
+        }
+
+        for (int i = 0; i < Engine::word_count; i++) {
+            engine2.discard(engine2.max());
+        }
+
+        if (engine1() == engine2()) {
+            std::cout << __PRETTY_FUNCTION__ << " passed for overflow_position " << overflow_position << std::endl;
+        }
+        else {
+            std::cout << __PRETTY_FUNCTION__ << " failed for overflow_position " << overflow_position << std::endl;
+            break;
+        }
+    }
+}
+
+int main() {
+    conformance_test<std::philox4x32>();
+    conformance_test<std::philox4x64>();
+
+    set_counter_conformance_test<std::philox4x32>();
+    set_counter_conformance_test<std::philox4x64>();
+
+    api_test<std::philox4x32>();
+    api_test<std::philox4x64>();
+    api_test<philox2x32_w5>();
+    api_test<philox2x32_w30>();
+    api_test<philox2x64_w15>();
+    api_test<philox2x64_w49>();
+
+    seed_test<std::philox4x32>();
+    seed_test<std::philox4x64>();
+    seed_test<philox2x32_w5>();
+    seed_test<philox2x32_w30>();
+    seed_test<philox2x64_w15>();
+    seed_test<philox2x64_w49>();
+
+    discard_test<std::philox4x32>();
+    discard_test<std::philox4x64>();
+    discard_test<philox2x32_w5>();
+    discard_test<philox2x32_w30>();
+    discard_test<philox2x64_w15>();
+    discard_test<philox2x64_w49>();
+
+    skip_test<std::philox4x32>();
+    skip_test<std::philox4x64>();
+    skip_test<philox2x32_w5>();
+    skip_test<philox2x32_w30>();
+    skip_test<philox2x64_w15>();
+    skip_test<philox2x64_w49>();
+
+    counter_overflow_test<std::philox4x32>();
+    counter_overflow_test<std::philox4x64>();
+    counter_overflow_test<philox2x32_w5>();
+    counter_overflow_test<philox2x32_w30>();
+    counter_overflow_test<philox2x64_w15>();
+    counter_overflow_test<philox2x64_w49>();
+
+    discard_overflow_test<std::philox4x32>();
+    discard_overflow_test<std::philox4x64>();
+    discard_overflow_test<philox2x32_w5>();
+    discard_overflow_test<philox2x32_w30>();
+    discard_overflow_test<philox2x64_w15>();
+    discard_overflow_test<philox2x64_w49>();
+
+    return 0;
+}