Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
Francesco Biscani
@bluescarni
@andrewcorrigan I managed to squeeze in some debug time in the afternoon. I think the issue you reported is a bug in the polynomial integration function, I have opened a PR with a tentative solution.
Andrew Corrigan
@andrewcorrigan
thanks! looks good
I'm testing my application and a lot of things are working the same as Piranha, however for certain cases, I see errors like the below. Would you have any hints on what to do about this?
# 2 | :0 | decltype(series_default_in_place_addsub_algorithm<true, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&&&, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&&>.second)::type obake::detail::series_default_in_place_addsub_impl<true, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void> >(obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&&&, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&&)
# 1 | :0 | void obake::detail::series_sym_extender<obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&>(obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&, obake::series<obake::polynomials::packed_monomial<unsigned long long, void>, mppp::rational<1ul>, obake::polynomials::tag, void>&&&, boost::container::flat_map<unsigned long, boost::container::flat_set<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::less<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >, void>, std::__1::less<unsigned long>, void> const&)
# 0 | :0 | obake::polynomials::packed_monomial<unsigned long long, std::__1::enable_if<is_k_packable_v<unsigned long long>, void>::type> obake::polynomials::key_merge_symbols<unsigned long long>(obake::polynomials::packed_monomial<unsigned long long, std::__1::enable_if<is_k_packable_v<unsigned long long>, void>::type> const&, boost::container::flat_map<unsigned long, boost::container::flat_set<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::less<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >, void>, std::__1::less<unsigned long>, void> const&, boost::container::flat_set<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::less<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >, void> const&)

Exception type   : std::overflow_error
Exception message: Invalid size specified in the constructor of a Kronecker packer for the type 'unsigned long long': the maximum possible size is 21, but a size of 22 was specified instead
Andrew Corrigan
@andrewcorrigan
Would I need d_packed_monomial instead? I'm not clear from the examples what reasonable template arguments would be. If this is the more appropriate data structure, how would I adapt the unpacking code above? Hope you're enjoying the holiday. Any help is appreciated.
Francesco Biscani
@bluescarni
Yes a d_packed_monomial does not have limits in the number of variables it can represent. The first template parameter, T, is the same as in packed_monomial. The second template parameter is the (approximate) number of bits you want to devote to each exponent. So for instance if you use NBits==8 and T an unsigned integral type, then each exponent will have a range of circa [0, 256] (it will be less than that in reality, but same order of magnitude).
A d_packed_monomial is kind of like a vector of packed_monomials
so the packing/unpacking is a bit more complicated
brb
Francesco Biscani
@bluescarni
From a d_packed_monomial, you can access the underlying vector of Ts via the _container() function. You can then proceed to unpack each T value as if it was the value in a packed_monomial. For instance:
using d_type = d_packed_monomial<unsigned long, 8>;
// This is the number of exponents packed in each unsigned long stored in a d_type.
// This is a compile-time constant.
constexpr auto psize = d_type::psize;

// Create a monomial.
d_type d{....};

// Print out all exponents.
for (const auto &n : d._container()) {
  k_unpacker<unsigned long> ku(n, psize);

  unsigned long tmp;
  for (auto j = 0u; j < psize; ++j) {
    ku >> tmp;
    std::cout << tmp << '\n';
  }
}
Note that like this it will print a number of exponents multiple of psize. That is, if psize is 8, then the code above will print either 8, 16, 24, 32, etc... exponents.
Even if in reality the actual number of exponents is not necessarily a multiple of 8.
So in actual code you need a bit more book-keeping on how to stop the inner for loop, but hopefully you get the gist.
Andrew Corrigan
@andrewcorrigan
thank you
where would I get the actual number of exponents from? would the polynomial containing the monomial give me the right size, as was the case with packed_monomial?
Francesco Biscani
@bluescarni
yes exactly. you can see an example of this in the function that computes the degree of a d_packed_monomial:
https://github.com/bluescarni/obake/blob/master/include/obake/polynomials/d_packed_monomial.hpp#L902
Andrew Corrigan
@andrewcorrigan
got it, thanks!
Francesco Biscani
@bluescarni
no problem! when you have time, I'd be interested to hear how obake is working for you wrt piranha, both in terms of API and in terms of performance.
Andrew Corrigan
@andrewcorrigan
so far with packed_monomial, without overflow, I saw an immediate 2x speed-up
Francesco Biscani
@bluescarni
that's great to hear
let me know if you find weak spots wrt piranha. obake was built to be faster but there could always be instances where it ends up being slower for some reason.
are you normally using polynomials over the rationals?
Andrew Corrigan
@andrewcorrigan
In terms of API, overall it's great and it's wonderful to have such a high-quality, open-source library available. You've helped me with my main struggles switching over. A few minor differences I noticed:
  1. in Piranha, I could divide polynomials as long as the denominator was actually just a rational number. In Obake, it seems to require an explicit conversion to a rational number to compile. Maybe that's better than a run-time error, but it was a difference that I noticed.
  2. In Piranha, I could pass a string into the constructor, like Polynomial{"x"}, Polynomial{"y"}. In Obake, I think that compiles, but didn't work. Instead, I had to change those to make_polynomials<Polynomial>("x", "y"); I like the make_polynomials function, but would prefer the string-constructor fail to compile rather than working in an unexpected way.
  3. I like that in Obake the range-based for-loop to iterate over the monomials and coefficients of a polynomial returns monomial as a monomial type, whereas in Piranha the monomial was still a polynomial type (I thought that was confusing).
Once I get d_packed_monomial I can measure performance for less trivial cases
yes, always over rationals
Francesco Biscani
@bluescarni
Regarding 2, the fact that the string constructor is still there but it does not operate as before is an unfortunate consequence of the fact that you can construct rationals from string. In piranha the poly ctor from string had a special meaning, while in obake the poly ctors other than copy/move simply forward the construction to the coefficient type.
Regarding the use of rationals as polynomial coefficients, obake does not implement yet an optimization that piranha had. piranha used to convert all rational coefficients in a poly multiplication to a common denominator, thus transforming a rational poly multiplication into an integer poly multiplication. This essentially means avoiding a lot of GCD computations which are not needed. obake will also have this optimisation, but it does not at the moment.
Andrew Corrigan
@andrewcorrigan
wow, so it can get even faster?
Francesco Biscani
@bluescarni
it should, but only if the polys are long enough I believe. In my tests I always used rather large polynomials and there the difference was very visible.
it might be possible that for shorter polynomials the difference is not that great, but it's something that needs to be studied/tuned.
Andrew Corrigan
@andrewcorrigan
large is important for me too
Francesco Biscani
@bluescarni
are your rational coefficients large in terms of bit width?
Andrew Corrigan
@andrewcorrigan
I don't know to be honest
how does that relate to digits?
sorry, I'm a bit ignorant of this stuff
Francesco Biscani
@bluescarni
no problem, just take the number of base 10 digits and multiply by ~3.3 (that is log(10)/log(2))
Andrew Corrigan
@andrewcorrigan
is that like if I could store the numerator/denominator using 32-bit / 64-bit / 128-bit integers?
Francesco Biscani
@bluescarni
yes
Andrew Corrigan
@andrewcorrigan
etc.
ok
Francesco Biscani
@bluescarni
I just wanted to get a feeling of whether your numerators/denominators are "small" (i.e., < 2**64) or not
Andrew Corrigan
@andrewcorrigan
I'm checking
yeah, they're small
well, at least at first, we perform many operations, so I imagine they grow quite a bit
Francesco Biscani
@bluescarni
ok, I see. I am pretty sure you'll see a speedup after I implement that optimisation, even though it might not be easy to say by how much exactly.
Andrew Corrigan
@andrewcorrigan
honestly, Piranha was already blazing fast (compared to Sympy), but I just implemented the d_packed_monomial , it's working great, and for a much larger problem it's close to 2x faster
Francesco Biscani
@bluescarni
great!
Andrew Corrigan
@andrewcorrigan
Here's an example of some the errors that -Werror triggers when the TBB headers are included:
[ 69%] Building CXX object /obake/CMakeFiles/obake.dir/src/series.cpp.o
In file included from /obake/src/series.cpp:14:
In file included from /obake/include/obake/series.hpp:41:
In file included from /third_party/tbb/include/tbb/blocked_range.h:20:
/third_party/tbb/include/tbb/tbb_stddef.h:350:16: error: use of old-style cast [-Werror,-Wold-style-cast]
    return 0==((uintptr_t)pointer & (alignment-1));
               ^          ~~~~~~~
/third_party/tbb/include/tbb/tbb_stddef.h:476:33: error: use of old-style cast [-Werror,-Wold-style-cast]
    static const size_t value = (size_t)((sizeof(size_t)==sizeof(u)) ? u : ull);
                                ^       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /obake/src/series.cpp:14:
In file included from /obake/include/obake/series.hpp:42:
In file included from /third_party/tbb/include/tbb/parallel_for.h:21:
In file included from /third_party/tbb/include/tbb/task.h:21:
In file included from /third_party/tbb/include/tbb/tbb_machine.h:235:
In file included from /third_party/tbb/include/tbb/machine/linux_intel64.h:24:
/third_party/tbb/include/tbb/machine/gcc_ia32_common.h:29:12: error: implicit conversion changes signedness: 'uintptr_t' (aka 'unsigned long') to 'intptr_t' (aka 'long') [-Werror,-Wsign-conversion]
    return j;
    ~~~~~~ ^
In file included from /obake/src/series.cpp:14:
In file included from /obake/include/obake/series.hpp:42:
In file included from /third_party/tbb/include/tbb/parallel_for.h:21:
In file included from /third_party/tbb/include/tbb/task.h:21:
In file included from /third_party/tbb/include/tbb/tbb_machine.h:235:
/third_party/tbb/include/tbb/machine/linux_intel64.h:70:1: error: use of old-style cast [-Werror,-Wold-style-cast]
__TBB_MACHINE_DEFINE_ATOMICS(1,int8_t,"")
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:44:49: note: expanded from macro '__TBB_MACHINE_DEFINE_ATOMICS'
                          : "=a"(result), "=m"(*(volatile T*)ptr)                    \
                                                ^            ~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:70:1: error: use of old-style cast [-Werror,-Wold-style-cast]
__TBB_MACHINE_DEFINE_ATOMICS(1,int8_t,"")
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:45:62: note: expanded from macro '__TBB_MACHINE_DEFINE_ATOMICS'
                          : "q"(value), "0"(comparand), "m"(*(volatile T*)ptr)       \
                                                             ^            ~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:70:1: error: use of old-style cast [-Werror,-Wold-style-cast]
__TBB_MACHINE_DEFINE_ATOMICS(1,int8_t,"")
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:54:48: note: expanded from macro '__TBB_MACHINE_DEFINE_ATOMICS'
                          : "=r"(result),"=m"(*(volatile T*)ptr)                     \
                                               ^            ~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:70:1: error: use of old-style cast [-Werror,-Wold-style-cast]
__TBB_MACHINE_DEFINE_ATOMICS(1,int8_t,"")
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:55:47: note: expanded from macro '__TBB_MACHINE_DEFINE_ATOMICS'
                          : "0"(addend), "m"(*(volatile T*)ptr)                      \
                                              ^            ~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:70:1: error: use of old-style cast [-Werror,-Wold-style-cast]
__TBB_MACHINE_DEFINE_ATOMICS(1,int8_t,"")
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/third_party/tbb/include/tbb/machine/linux_intel64.h:64:48: note: expanded from macro '__TBB_MACHINE_DEFINE_ATOMICS'
                          : "=r"(result),"=m"(*(volatile T*)ptr)                     \
                                               ^            ~~~
/third_party/tb
This is using Clang 10 on Linux. Something similar happened with Xcode on OS X 10.14.
Francesco Biscani
@bluescarni
could you run make VERBOSE=1 and check if the TBB headers are included via -I or -isystem?
and what version of TBB is this?