Hello, I am trying to use Piranha to compute some symbolic polynomial expressions. Roughly speaking, I have a loops that perform some operations and add the result to the global polynomial. It works fine with a very small problem (6 nested loops over 2 values), but when I increase th size of the problem (size 4), the program crashes in a Piranha operation. I have it using both piranha::polynomial<piranha::integer,piranha::k_monomial> and piranha::polynomial<piranha::integer,piranha::monomial<int>>. Using the former, it throws an std::bad_alloc exception.
Using the latter, I have the following error message:
terminate called after throwing an instance of 'std::invalid_argument'
what():
Function name : encode
Location : /usr/local/include/piranha/kronecker_array.hpp, line 281
Exception type : std::invalid_argument
Exception message: size of vector to be encoded is too large
I can also provide the Valgrind trace.
Is Piranha intended to be used with large polynomials?
This error:
terminate called after throwing an instance of 'std::invalid_argument'
what():
Function name : encode
Location : /usr/local/include/piranha/kronecker_array.hpp, line 281
Exception type : std::invalid_argument
Exception message: size of vector to be encoded is too large
indicates that there are too many exponents for the selected monomial representation (which should be k_monomial
). The other message is more generic, but I suspect it might also have to do with the total number of exponents.
@bluescarni Thank you very much for your reply. Sorry about the flu. Don't worry about the reply time, I totally understand.
My polynomials are multivariate and they all have a coefficient of 1. Actually, polynomials might not be the best idea. I want to compute some symbolic values for tensors. I am using symbols for the elements of my tensors, and I am adding them into a polynomial. I have reproduced the error with a simpler example (my actual code is computing something more interesting ;) ):
#define N 8
int main( int argc, char** argv ) {
p_type poly;
for( int i = 0 ; i < N ; i++ ) {
for( int j = 0 ; j < N ; j++ ) {
for( int k = 0 ; k < N ; k++ ) {
std::string name = "T_" + std::to_string( i ) + std::to_string( j ) + std::to_string( k );
p_type p{name};
poly += p;
// std::cout << poly << std::endl;
}
}
}
return EXIT_SUCCESS;
}
But again, since my polynomials are multivariant, this might not be the most efficien way to do it with Piranha.
I have tried two different definitions of p_type:
using p_type = piranha::polynomial<piranha::integer,piranha::monomial<int>>; // throws a bad_alloc exception
using p_type = piranha::polynomial<piranha::integer,piranha::k_monomial>; // throws an invalid_argument exception
Couldn't find a chat room, for questions about Obake. Asking here:
Is checking for a monomial that is 1 best accomplished via:
return m.get_value() == 0;
Or is this using internal implementation details and subject to break?
More generally, could you please offer some pointers on accessing the symbols and exponents of a monomial in Obake? With your help, in Piranha, I accessed a monomial (as a polynomial):
const auto& container = polynomial._container();
const auto& args = polynomial.get_symbol_set();
const auto& monomial = container.begin()->m_key;
for(auto i : range(monomial.size()))
{
std::cout << "exponent = " << monomial[i] << std::endl;
std::cout << "symbol = " << args[i].get_name() << std::endl;
}
What would the equivalent be in Obake directly on the monomial type (since the monomials of a polynomial are monomials in Obake)?
obake has a specialised primitive for detecting monomials which are equivalent to one. It is called key_is_one()
and you can see some examples of its usage here:
https://github.com/bluescarni/obake/blob/master/test/polynomials_d_packed_monomial_00.cpp#L229
Like all monomial-related primitives, you need to pass also the symbol set as last argument. E.g.,
using monomial_t = packed_monomial<unsigned long>;
symbol_set s{"x", "y", "z"};
monomial_t m{0, 0, 0};
assert(key_is_one(m, s));
p
is a polynomial object):for (const auto &[m, cf] : p) {
// In here "m" is a const reference to the monomial,
// "cf" a const reference to its coefficient.
}
packed_monomial
, you need first to extract the packed exponents via the get_value()
member function, and then you can unpack this value using facilities from the k_packing.hpp
header.
using monomial_t = packed_monomial<unsigned long>;
monomial_t m{1, 2, 3};
k_unpacker<unsigned long> ku(m.get_value(), 3);
unsigned long tmp;
ku >> tmp;
std::cout << "The first exponent is: " << tmp << '\n';
ku >> tmp;
std::cout << "The second exponent is: " << tmp << '\n';
ku >> tmp;
std::cout << "The third exponent is: " << tmp << '\n';
p.get_symbol_set()
? Does the unpacker report exponents in the same order as p.get_symbol_set()
?
In terms of builds, one thing wasn't clear to me: I see that you have a very nice modern CMake build infrastructure for Obake, mp++ and other dependencies. In my projects, I try to include all dependencies in one parent project, and then build them all hierarchically via add_subdirectory. I've been modifying your CMake file as follows:
# mp++.
# NOTE: put the minimum version in a variable
# so that we can re-use it in the config-file package
# machinery below.
#set (_OBAKE_MIN_MPPP_VERSION 0.17)
#find_package(mp++ REQUIRED)
#if(${mp++_VERSION} VERSION_LESS ${_OBAKE_MIN_MPPP_VERSION})
# message(FATAL_ERROR "The minimum mp++ version required by obake is ${_OBAKE_MIN_MPPP_VERSION}, but version ${mp++_VERSION} was found instead.")
#endif()
add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/mppp")
Is there a way I could have a parent project that includes both Obake and mp++ via add_subdirectory
, where I instead leave your Obake CMakeLists.txt unmodified and find_package(mp++ REQUIRED)
detects that mp++
is being configured/built as part of the same project?
find_package()
ultimately does is to create an mp++::mp++
imported target. Perhaps you can try to write:if(NOT TARGET mp++::mp++)
set (_OBAKE_MIN_MPPP_VERSION 0.17)
find_package(mp++ REQUIRED)
if(${mp++_VERSION} VERSION_LESS ${_OBAKE_MIN_MPPP_VERSION})
message(FATAL_ERROR "The minimum mp++ version required by obake is ${_OBAKE_MIN_MPPP_VERSION}, but version ${mp++_VERSION} was found instead.")
endif()
endif()
# 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
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).
d_packed_monomial
is kind of like a vector of packed_monomial
s