diff --git a/engine/CMakeLists.txt b/engine/CMakeLists.txt index 7e47911..72abac3 100644 --- a/engine/CMakeLists.txt +++ b/engine/CMakeLists.txt @@ -72,7 +72,9 @@ add_library( ${PARSER_DIR}/Scanner.cu ${PARSER_DIR}/Parser.cu) -add_executable(cusymint ${SRC_DIR}/main.cu) +add_executable(cusymint + ${SRC_DIR}/Performance.cu + ${SRC_DIR}/main.cu) add_executable( testcusymint diff --git a/engine/src/Evaluation/ComputationHistory.cu b/engine/src/Evaluation/ComputationHistory.cu index e53ac6c..242c5fe 100644 --- a/engine/src/Evaluation/ComputationHistory.cu +++ b/engine/src/Evaluation/ComputationHistory.cu @@ -234,6 +234,22 @@ namespace Sym { fmt::print("}}\n"); } + size_t ComputationStep::get_memory_occupance_in_bytes() const { + size_t size = 0; + for (const auto& expr : expression_tree) { + size += expr.size(); + } + return sizeof(ComputationStep) + size * sizeof(Symbol); + } + + size_t ComputationHistory::get_memory_occupance_in_bytes() const { + size_t size = 0; + for (const auto& step : computation_steps) { + size += step.get_memory_occupance_in_bytes(); + } + return sizeof(ComputationHistory) + size; + } + void ComputationHistory::complete() { if (!is_solution_at_end_of_steps()) { Util::crash("ComputationHistory must have solution at the end to complete"); @@ -290,19 +306,24 @@ namespace Sym { else { size_t trans_idx = 1; for (const auto& trans : step.get_operations(*prev_step)) { - result.push_back(fmt::format( - R"(\quad \text{{{}{} {}:}}\:{}:)", trans_idx, get_ordinal_suffix(trans_idx), - Globalization::INTEGRAL, trans->get_description())); + result.push_back(fmt::format(R"(\quad \text{{{}{} {}:}}\:{}:)", trans_idx, + get_ordinal_suffix(trans_idx), + Globalization::INTEGRAL, + trans->get_description())); ++trans_idx; } } - result.push_back(fmt::format(R"(=\qquad {})", expression.data()->to_tex())); + result.push_back(fmt::format(R"(=\qquad {})", expression_str)); } prev_str = expression_str; prev_step = &step; } + if (!result.empty()) { + result[result.size() - 1] += " + C"; + } + return result; } diff --git a/engine/src/Evaluation/ComputationHistory.cuh b/engine/src/Evaluation/ComputationHistory.cuh index 01b283a..474c533 100644 --- a/engine/src/Evaluation/ComputationHistory.cuh +++ b/engine/src/Evaluation/ComputationHistory.cuh @@ -33,6 +33,7 @@ namespace Sym { std::vector get_expression() const; TransformationList get_operations(const ComputationStep& previous_step) const; void print_step() const; + size_t get_memory_occupance_in_bytes() const; }; using ComputationStepCollection = std::list; @@ -58,6 +59,8 @@ namespace Sym { void print_history() const; const ComputationStepCollection& get_steps() const { return computation_steps; } + + size_t get_memory_occupance_in_bytes() const; }; } diff --git a/engine/src/Evaluation/Integrator.cu b/engine/src/Evaluation/Integrator.cu index e33b3f1..5998392 100644 --- a/engine/src/Evaluation/Integrator.cu +++ b/engine/src/Evaluation/Integrator.cu @@ -4,9 +4,9 @@ #include #include +#include "ComputationHistory.cuh" #include "IntegratorKernels.cuh" #include "StaticFunctions.cuh" -#include "ComputationHistory.cuh" #include "Evaluation/Heuristic/Heuristic.cuh" #include "Symbol/SubexpressionCandidate.cuh" @@ -188,7 +188,8 @@ namespace Sym { scan_array_2.zero_mem(); cudaDeviceSynchronize(); - help_space.reoffset_like(integrals.iterator(), HELP_SPACE_MULTIPLIER, Heuristic::COUNT); + help_space.reoffset_like(integrals.iterator(), + HELP_SPACE_MULTIPLIER, Heuristic::COUNT); Kernel::check_heuristics_applicability<<>>( integrals, expressions, help_space, scan_array_1, scan_array_2); diff --git a/engine/src/Evaluation/Integrator.cuh b/engine/src/Evaluation/Integrator.cuh index ed3077a..a99c74d 100644 --- a/engine/src/Evaluation/Integrator.cuh +++ b/engine/src/Evaluation/Integrator.cuh @@ -238,6 +238,131 @@ namespace Sym { return std::nullopt; } + + size_t memory_usage_for_integral(const std::vector& integral, + size_t& total_memory) { + size_t initial_usage; + size_t max_usage; + size_t usage; + + size_t initial_size = + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + integrals_swap.symbols_capacity() + + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size()); + size_t max_shown_size = initial_size; + + cudaMemGetInfo(&initial_usage, &total_memory); + max_usage = initial_usage; + + expressions.load_from_vector({single_integral_vacancy()}); + integrals.load_from_vector({first_expression_candidate(integral)}); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + integrals_swap.symbols_capacity() + + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + for (size_t i = 0;; ++i) { + simplify_integrals(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + check_for_known_integrals(); + apply_known_integrals(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + if (is_original_expression_solved()) { + break; + } + + remove_unnecessary_candidates(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + check_heuristics_applicability(); + apply_heuristics(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + + if (has_original_expression_failed()) { + break; + } + + remove_failed_candidates(); + + cudaMemGetInfo(&usage, nullptr); + max_usage = std::min(usage, max_usage); + + max_shown_size = std::max( + max_shown_size, + sizeof(Symbol) * + (expressions.symbols_capacity() + integrals.symbols_capacity() + + expressions_swap.symbols_capacity() + + integrals_swap.symbols_capacity() + help_space.symbols_capacity()) + + sizeof(size_t) * (scan_array_1.size() + scan_array_2.size()) + + sizeof(EvaluationStatus) * + (evaluation_statuses_1.size() + evaluation_statuses_2.size())); + } + printf("Declared bytes:%lu\n", max_shown_size); + return initial_usage - max_usage; + } }; } diff --git a/engine/src/Performance.cu b/engine/src/Performance.cu new file mode 100644 index 0000000..2d356fb --- /dev/null +++ b/engine/src/Performance.cu @@ -0,0 +1,313 @@ +#include "Performance.cuh" +#include "Symbol/Integral.cuh" +#include "Symbol/SymbolType.cuh" +#include +#include +#include +#include + +namespace Performance { + namespace { + std::string exec_and_read_output(const std::string cmd) { + // stolen from SolverProcessManager.cu + std::array buffer{}; + std::string result; + std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); + if (!pipe) { + throw std::runtime_error("popen() failed!"); + } + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + result += buffer.data(); + } + return result; + } + + void capitalize_substring(std::string& string, const std::string& what) { + if (string.size() < what.size()) { + return; + } + for (size_t i = 0; i <= string.size() - what.size(); ++i) { + if (string.substr(i, what.size()) == what && + (i == 0 || !std::isalpha(string[i - 1])) && + (i == string.size() - what.size() || !std::isalpha(string[i + what.size()]))) { + string[i] = string[i] - 'a' + 'A'; + } + } + } + + void replace(std::string& string, const std::string& from, const std::string& to) { + if (string.size() < from.size()) { + return; + } + for (ssize_t i = 0; i <= string.size() - from.size(); ++i) { + if (string.substr(i, from.size()) == from) { + string = string.substr(0, i) + to + string.substr(i + from.size()); + } + } + } + + std::string make_mathematica_command(const std::string& integral) { + std::string expression = integral; + capitalize_substring(expression, "e"); + capitalize_substring(expression, "pi"); + capitalize_substring(expression, "sqrt"); + + return fmt::format( + R"(wolframscript -code 'expr=ToExpression["{}",TraditionalForm];Integrate[expr,x];{{t,b}}=AbsoluteTiming[Integrate[expr,x]];t')", + expression); + } + + std::string make_sympy_command(const std::string& integral) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "^", "**"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + + return fmt::format( + R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");integrate({},x);print(utilities.timeutils.timed(lambda:integrate({},x))[1])')", + expression, expression); + } + + std::string make_matlab_command(const std::string& integral) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "E", "exp(sym(1))"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + + return fmt::format( + R"(matlab -batch 'syms x z;f=@()int({},x);f();fprintf("%f\n",timeit(f));')", + expression); + } + + std::string make_mathematica_command_batch(const std::vector& integrals) { + std::string result = "wolframscript -code \'"; + for (const auto& integral : integrals) { + std::string expression = integral; + capitalize_substring(expression, "e"); + capitalize_substring(expression, "pi"); + capitalize_substring(expression, "sqrt"); + result += fmt::format("{{t,b}}=AbsoluteTiming[Integrate[ToExpression[\"{}\"," + "TraditionalForm],x]];Print[t];ClearSystemCache[];", + expression); + } + + return result + "\'"; + } + + std::string make_sympy_command_batch(const std::vector& integrals) { + std::string result = R"(python3 -c 'from sympy import *;x=Symbol("x");t=Symbol("t");)"; + for (const auto& integral : integrals) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "^", "**"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + result += fmt::format( + "print(utilities.timeutils.timed(lambda:integrate({},x))[1]);", expression); + } + + return result + "\'"; + } + + std::string make_matlab_command_batch(const std::vector& integrals) { + std::string result = R"(matlab -batch 'syms x t;)"; + for (const auto& integral : integrals) { + std::string expression = integral; + capitalize_substring(expression, "e"); + replace(expression, "E", "exp(sym(1))"); + replace(expression, "arctan", "atan"); + replace(expression, "arccos", "acos"); + replace(expression, "arcsin", "asin"); + replace(expression, "sgn", "sign"); + replace(expression, "ln", "log"); + + // result += fmt::format(R"(f=@()int({},x);fprintf("%f\n",timeit(f));)", + // expression); + result += fmt::format(R"(tic;int({},x);t=toc;fprintf("%f\n",t);)", expression); + } + + return result + "\'"; + } + } + + void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result, const std::string& matlab_result) {} + + void print_human_readable_results(const std::string& integral_str, + const double& cusymint_seconds, bool cusymint_success, + const std::string& mathematica_result, + const std::string& sympy_result, + const std::string& matlab_result) { + printf("%s:\n" + " Cusymint time: %f%s\n" + " Mathematica time: %s\n" + " Sympy time: %s\n" + " Matlab time: %s\n\n", + integral_str.c_str(), cusymint_seconds, cusymint_success ? "" : " (failure)", + mathematica_result.c_str(), sympy_result.c_str(), matlab_result.c_str()); + } + + void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result, const std::string& matlab_result) { + printf("%s;%s;%f;%s;%s;%s\n", integral_str.c_str(), cusymint_success ? "TRUE" : "FALSE", + cusymint_seconds, mathematica_result.c_str(), sympy_result.c_str(), + matlab_result.c_str()); + } + + void print_csv_headers() { + printf("integral;solved_by_cusymint;cusymint_time[s];mathematica_time[s];sympy_time[s];" + "matlab_time[s]\n"); + } + + void test_with_other_solutions(const std::string& integral_str, PrintRoutine print_results) { + const auto integral = Sym::integral(Parser::parse_function(integral_str)); + + Sym::Integrator integrator; + + const clock_t start = clock(); + const auto result = integrator.solve_integral(integral); + const clock_t end = clock(); + + const double cusymint_seconds = static_cast(end - start) / CLOCKS_PER_SEC; + + auto mathematica_result = exec_and_read_output(make_mathematica_command(integral_str)); + auto sympy_result = exec_and_read_output(make_sympy_command(integral[1].to_string())); + auto matlab_result = exec_and_read_output(make_matlab_command(integral[1].to_string())); + + print_results(integral_str, cusymint_seconds, result.has_value(), mathematica_result, + sympy_result, matlab_result); + } + + void test_cuda_and_print_commands(const std::vector& integrals) { + printf("Computing on CUDA...\n"); + + for (int i = 0; i < integrals.size(); ++i) { + Sym::Integrator integrator; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + const clock_t start = clock(); + const auto result = integrator.solve_integral(integral); + const clock_t end = clock(); + + fmt::print("{}{}\n", static_cast(end - start) / CLOCKS_PER_SEC, + result.has_value() ? "" : " (failure)"); + } + + printf("%s\n", make_mathematica_command_batch(integrals).c_str()); + printf("%s\n", make_sympy_command_batch(integrals).c_str()); + printf("%s\n", make_matlab_command_batch(integrals).c_str()); + } + + void test_memory_occupance(const std::vector& integrals) { + printf("Memory occupance:\n"); + for (int i = 0; i < integrals.size(); ++i) { + printf(" Integral %s:\n", integrals[i].c_str()); + + size_t initial_mem; + size_t after_init; + + cudaMemGetInfo(&initial_mem, nullptr); + + Sym::Integrator integrator; + + cudaMemGetInfo(&after_init, nullptr); + initial_mem -= after_init; + + size_t total; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + + const size_t usage = integrator.memory_usage_for_integral(integral, total); + + printf(" %10lu/%10lu B used (%f%%)\n", usage + initial_mem, total, + static_cast(usage + initial_mem) / total * 100); + } + } + + void test_history_cpu_memory_occupance(const std::vector& integrals) { + printf("Memory occupance for history:\n"); + for (int i = 0; i < integrals.size(); ++i) { + printf(" Integral %s:\n", integrals[i].c_str()); + + Sym::Integrator integrator; + Sym::ComputationHistory history; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + + const auto result = integrator.solve_integral_with_history(integral, history); + + printf(" %10luB used%s\n", history.get_memory_occupance_in_bytes(), result.has_value() ? "" : " (failure)"); + } + } + + void test_performance(const std::vector& integrals, PrintRoutine print_results) { + std::vector cusymint_seconds_vector(integrals.size()); + std::vector cusymint_results_vector(integrals.size()); + std::vector mathematica_results_vector(integrals.size()); + std::vector sympy_results_vector(integrals.size()); + std::vector matlab_results_vector(integrals.size()); + + printf("Computing on CUDA...\n"); + + for (int i = 0; i < integrals.size(); ++i) { + Sym::Integrator integrator; + const auto integral = Sym::integral(Parser::parse_function(integrals[i])); + const clock_t start = clock(); + const auto result = integrator.solve_integral(integral); + const clock_t end = clock(); + + cusymint_seconds_vector[i] = static_cast(end - start) / CLOCKS_PER_SEC; + cusymint_results_vector[i] = result.has_value(); + } + + printf("Computing on Mathematica...\n"); + auto mathematica_result = exec_and_read_output(make_mathematica_command_batch(integrals)); + // Warning: computing na integral which requires many substitutions on SymPy may hang your + // computer and is extremely slow! + printf("Computing on SymPy...\n"); + auto sympy_result = exec_and_read_output(make_sympy_command_batch(integrals)); + printf("Computing on MATLAB...\n"); + auto matlab_result = exec_and_read_output(make_matlab_command_batch(integrals)); + + size_t mc_idx = 0; + size_t sp_idx = 0; + size_t ml_idx = 0; + for (int i = 0; i < integrals.size(); ++i) { + const auto mc_inc = mathematica_result.find('\n', mc_idx); + if (mc_inc != std::string::npos) { + mathematica_results_vector[i] = mathematica_result.substr(mc_idx, mc_inc - mc_idx); + mc_idx = mc_inc + 1; + } + + const auto sp_inc = sympy_result.find('\n', sp_idx); + if (sp_inc != std::string::npos) { + sympy_results_vector[i] = sympy_result.substr(sp_idx, sp_inc - sp_idx); + sp_idx = sp_inc + 1; + } + + const auto ml_inc = matlab_result.find('\n', ml_idx); + if (ml_inc != std::string::npos) { + matlab_results_vector[i] = matlab_result.substr(ml_idx, ml_inc - ml_idx); + ml_idx = ml_inc + 1; + } + } + + if (print_results == print_csv_results) { + print_csv_headers(); + } + + for (int i = 0; i < integrals.size(); ++i) { + print_results(integrals[i], cusymint_seconds_vector[i], cusymint_results_vector[i], + mathematica_results_vector[i], sympy_results_vector[i], + matlab_results_vector[i]); + } + } + +} \ No newline at end of file diff --git a/engine/src/Performance.cuh b/engine/src/Performance.cuh new file mode 100644 index 0000000..85f5a84 --- /dev/null +++ b/engine/src/Performance.cuh @@ -0,0 +1,51 @@ +#ifndef PERFORMACNE_CUH +#define PERFORMANCE_CUH + +#include +#include +#include +#include +#include +#include + +#include "Evaluation/Integrator.cuh" +#include "Parser/Parser.cuh" +#include "Symbol/Macros.cuh" +#include "Symbol/Symbol.cuh" + +namespace Performance { + + using PrintRoutine = void (*)(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result, + const std::string& matlab_result); + + void do_not_print_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result, const std::string& matlab_result); + void print_human_readable_results(const std::string& integral_str, + const double& cusymint_seconds, bool cusymint_success, + const std::string& mathematica_result, + const std::string& sympy_result, + const std::string& matlab_result); + + void print_csv_results(const std::string& integral_str, const double& cusymint_seconds, + bool cusymint_success, const std::string& mathematica_result, + const std::string& sympy_result, const std::string& matlab_result); + + void print_csv_headers(); + + void test_with_other_solutions(const std::string& integral_str, + PrintRoutine print_results = print_human_readable_results); + + void test_memory_occupance(const std::vector& integrals); + + void test_history_cpu_memory_occupance(const std::vector& integrals); + + void test_cuda_and_print_commands(const std::vector& integrals); + + void test_performance(const std::vector& integrals, + PrintRoutine print_results = print_human_readable_results); +} + +#endif \ No newline at end of file diff --git a/engine/src/main.cu b/engine/src/main.cu index d5b9a56..c5aedb2 100644 --- a/engine/src/main.cu +++ b/engine/src/main.cu @@ -1,8 +1,10 @@ +#include #include - +#include #include #include +#include #include #include @@ -10,6 +12,7 @@ #include "Evaluation/Integrator.cuh" #include "Evaluation/StaticFunctions.cuh" +#include "Performance.cuh" #include "Symbol/Constants.cuh" #include "Symbol/ExpressionArray.cuh" #include "Symbol/Integral.cuh" @@ -43,32 +46,172 @@ std::string e_tower(size_t n) { return res; } -int main() { - if constexpr (Consts::DEBUG) { - fmt::print("Running in debug mode\n"); +std::vector generate_e_towers(size_t end) { + std::vector result(end + 1); + result[0] = "1"; + if (end > 0) { + result[1] = "e^x"; } + for (size_t i = 2; i <= end; ++i) { + std::string e_power = "*"; + for (size_t j = 1; j <= i; ++j) { + e_power += "e^"; + } + e_power += "x"; + result[i] = result[i - 1] + e_power; + } + return result; +} - Sym::Static::init_functions(); - - - const auto integral = Sym::integral(Parser::parse_function("tg(x)^2")); - - fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); - - Sym::Integrator integrator; - Sym::ComputationHistory history; - const auto solution = integrator.solve_integral_with_history(integral, history); +std::vector generate_sin_towers(size_t end) { + std::vector result(end + 1); + result[0] = "cos(x)"; + for (size_t i = 1; i <= end; ++i) { + std::string trig = "*cos("; + for (size_t j = 0; j < i; ++j) { + trig += "sin("; + } + trig += "x"; + for (size_t j = 0; j <= i; ++j) { + trig += ")"; + } + result[i] = result[i - 1] + trig; + } + return result; +} - if (solution.has_value()) { - fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); +std::vector generate_geometric_sums(size_t start, size_t end, size_t step = 1) { + std::vector result; + for (size_t i = start; i <= end; i += step) { + result.push_back(fmt::format("(x^{}-1)/(x-1)", i)); + } + return result; +} - fmt::print("\nComputation steps:\n\n"); +std::vector generate_random_polynomials(size_t min_rank, size_t max_rank, size_t count, + size_t seed) { + constexpr size_t MAX_COEFF = 50; + std::vector result(count); + srand(seed); + for (size_t i = 0; i < count; ++i) { + size_t rank = (rand() % (max_rank - min_rank + 1)) + min_rank; + result[i] = fmt::format("{}", rand() % (MAX_COEFF + 1)); + for (size_t j = 1; j <= rank; ++j) { + result[i] += fmt::format("{}{}*x^{}", (rand() % 2 == 0) ? "+" : "-", + rand() % (MAX_COEFF + 1), j); + } + } + return result; +} - for (const auto& step: history.get_tex_history()) { - fmt::print(" {}\\\\\n\n", step); +std::vector generate_trig_polynomials(size_t start, size_t end, size_t step) { + std::vector result; + for (size_t rank = start; rank <= end; rank += step) { + std::string str = "cos(x)"; + for (size_t j = 1; j <= rank; ++j) { + str += fmt::format("+sin(x)^{}*cos(x)", j); } + result.push_back(str); } - else { - fmt::print("No solution found\n"); + return result; +} + +int main() { + if constexpr (Consts::DEBUG) { + fmt::print("Running in debug mode\n"); } + + Sym::Static::init_functions(); + + fmt::print("First two integrals are for warmup!\n"); + + std::vector integrals = { + "x+sin(x)", "e^x", + // "x^4-5*x+3", + // "(1-x^3)^2", + // "x^3*t^2", + // "x^(2/3)", + // "1/x^(2/3)", + // "(x^3-3*x^2+1)/x^5", + // "(x^(2/3)-sqrt(x))/x^(1/3)", + // "2*e^x+3*5^x", + // "x*cos(x)", + // "x*e^x", + // "x^2*e^x", + // "x*e^(-x)", + // "x^2*cos(x)", + // "sqrt(x)*ln(x)", + // "ln(x)/x^4", + // "(2*x-1)^20", + // "cos(3*x-1)", + // "e^x/(3+e^x)", + // "tan(x)", + // "tan(x)/cos(x)^2", + // "cos(x)*e^sin(x)", + // "cos(x)/sqrt(1+sin(x))", + // "6^(1-x)", + // "(sqrt(x)-2)^2/x^2", + // "1/(sin(x)^2*cos(x)^2)", + // "x*arctan(x)", + // "abs(x)", + // "(abs(1-x^2)+1+x^2)/2", + // "tan(x)^2", + // "(2*x-3)^10", + // "sin(x)^5*cos(x)", + // "ln(x)", + // "x*cos(x)", + // "x^2*e^(1-x)", + // "sin(x)^5", + // "sin(x)^4*cos(x)^3", + // "e^x*(x+1)*ln(x)", + // "(1-x^200)/(1-x)", + // "(1-sin(x)^20)/(1-sin(x))*cos(x)", + // e_tower(70), + }; + + // add huge sum + // integrals.push_back(fmt::format("{}", fmt::join(integrals, "+"))); + + // fmt::print("{}", fmt::join(generate_sin_towers(5), ",")); + + // add e towers + // const auto e_towers = generate_e_towers(10); + // integrals.insert(integrals.end(), e_towers.begin(), e_towers.end()); + + // add sin towers + // const auto sin_towers = generate_sin_towers(8); + // integrals.insert(integrals.end(), sin_towers.begin(), sin_towers.end()); + + // add geometric sums + const auto geo_sums = generate_geometric_sums(0, 200, 5); + //integrals.insert(integrals.end(), geo_sums.begin(), geo_sums.end()); + + // add random polynomials + const auto polynomials = generate_trig_polynomials(0, 100, 5); + integrals.insert(integrals.end(), polynomials.begin(), polynomials.end()); + + // Performance::test_performance(integrals, Performance::print_csv_results); + Performance::test_memory_occupance(integrals); + // Performance::test_history_cpu_memory_occupance(integrals); + // const auto integral = Sym::integral(Parser::parse_function(integrals.back())); + + // fmt::print("Trying to solve an integral: {}\n", integral.data()->to_tex()); + + // Sym::Integrator integrator; + // Sym::ComputationHistory history; + // const auto solution = + // integrator.solve_integral_with_history(integral, history); + + // if (solution.has_value()) { + // fmt::print("Success! Solution:\n{} + C\n", solution.value().data()->to_tex()); + + // fmt::print("\nComputation steps:\n\n"); + + // for (const auto& step: history.get_tex_history()) { + // fmt::print(" {}\\\\\n\n", step); + // } + // } + // else { + // fmt::print("No solution found\n"); + // } } diff --git a/engine/test/AdvancedIntegral.cu b/engine/test/AdvancedIntegral.cu index 7968fb0..854f7cb 100644 --- a/engine/test/AdvancedIntegral.cu +++ b/engine/test/AdvancedIntegral.cu @@ -1,5 +1,7 @@ #include "IntegralCommons.cuh" +#include "Evaluation/Integrator.cuh" + #define ADVANCED_INTEGRAL_TEST(_name, _integral, _expected_result) \ INTEGRAL_TEST(AdvancedIntegral, _name, _integral, _expected_result) @@ -19,24 +21,24 @@ namespace Test { ADVANCED_INTEGRAL_TEST(CosineByParts, "x cos(x)", "x sin(x)+cos(x)") ADVANCED_INTEGRAL_TEST(ExponentByParts, "x e^x", "x e^x-e^x") ADVANCED_INTEGRAL_TEST(ExponentByParts2, "x^2e^x", "e^x(x^2-2x+2)") - //ADVANCED_INTEGRAL_TEST(NegativeExponentByParts, "x e^(-x)", "-e^(-x)(x+1)") // linear substitution for nested-vars-only required + ADVANCED_INTEGRAL_TEST(NegativeExponentByParts, "x e^(-x)", "-e^(-x)(x+1)") ADVANCED_INTEGRAL_TEST(CosineByParts2, "x^2cos(x)", "x^2sin(x)+2x cos(x)-2sin(x)") ADVANCED_INTEGRAL_TEST(LogarithmTimesRoot, "sqrt(x)ln(x)", "2/3x^(3/2)(ln(x)-2/3)") ADVANCED_INTEGRAL_TEST(LogarithmDividedByPower, "ln(x)/x^4", "-1/(3x^3)(ln(x)+1/3)") ADVANCED_INTEGRAL_TEST(VeryLongPolynomial, "(2x-1)^20", "1/42(2x-1)^21+1/42") ADVANCED_INTEGRAL_TEST(CosineOfLinearFunction, "cos(3x-1)", "sin(3x-1)/3") - //ADVANCED_INTEGRAL_TEST(RationalWithSubstitution, "2x/(5x^2+1)", "1/5ln(5x^2+1)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(RationalWithSubstitution, "2x/(5x^2+1)", "1/5ln(5x^2+1)") // substitution t=x^2+a required ADVANCED_INTEGRAL_TEST(SimpleRationalWithExponent, "e^x/(3+e^x)", "ln(3+e^x)") - //ADVANCED_INTEGRAL_TEST(ExponentWithSquare, "x*e^(-x^2)", "-1/2e^(-x^2)") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(ExponentOfInverse, "e^(1/x)/x^2", "-e^(1/x)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(ExponentWithSquare, "x*e^(-x^2)", "-1/2e^(-x^2)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(ExponentOfInverse, "e^(1/x)/x^2", "-e^(1/x)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(Tangent, "tg(x)", "-ln(cos(x))") ADVANCED_INTEGRAL_TEST(TangentDividedByCosSquared, "tg(x)/cos^2(x)", "tg^2(x)/2") ADVANCED_INTEGRAL_TEST(CosTimesExpSine, "cos(x)e^sin(x)", "e^sin(x)") ADVANCED_INTEGRAL_TEST(CosineDividedByRootOfLinSine, "cos(x)/sqrt(1+sin(x))", "2sqrt(1+sin(x))") - //ADVANCED_INTEGRAL_TEST(TangentDerivativeWithSubstitution, "x^3/cos^2(x^4)", "tg(x^4)/4") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(TangentDerivativeWithSubstitution, "x^3/cos^2(x^4)", "tg(x^4)/4") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(ExponentOfLinear, "6^(1-x)", "-6^(1-x)/ln(6)") - //ADVANCED_INTEGRAL_TEST(ArctanSubtitution, "1/(1+x^2)/arctan(x)", "ln(arctan(x))") // substitution t=arctan(x) required + // ADVANCED_INTEGRAL_TEST(ArctanSubtitution, "1/(1+x^2)/arctan(x)", "ln(arctan(x))") // substitution t=arctan(x) required // task 4 contains examples which require nested integrals (not implemeneted) @@ -45,39 +47,44 @@ namespace Test { ADVANCED_INTEGRAL_TEST(SimpleRationalWithSquare, "(sqrt(x)-2)^2/x^2", "-4/x+8/sqrt(x)+ln(x)") ADVANCED_INTEGRAL_TEST(SineCosineSquaredInDenominator, "1/(sin^2(x)cos^2(x))", "-1/tg(x)+tg(x)") - //ADVANCED_INTEGRAL_TEST(LogarithmDividedByX, "ln^5(x)/x", "ln^6(x)/6") // substitution t=ln(x) or cyclic integrals required - //ADVANCED_INTEGRAL_TEST(ArctangentWithX3Substitution, "x^3/(x^8+1)", "arctg(x^4)/4") // substitution t=x^4 required - //ADVANCED_INTEGRAL_TEST(LogarithmMultipliedByX, "x^n*ln(x)", "x^(n+1)ln(x)/(n+1)-x^(n+1)/(n+1)^2") // non-numeric powers integration by parts required + // ADVANCED_INTEGRAL_TEST(LogarithmDividedByX, "ln^5(x)/x", "ln^6(x)/6") // substitution t=ln(x) or cyclic integrals required + // ADVANCED_INTEGRAL_TEST(ArctangentWithX3Substitution, "x^3/(x^8+1)", "arctg(x^4)/4") // substitution t=x^4 required + // ADVANCED_INTEGRAL_TEST(LogarithmMultipliedByX, "x^n*ln(x)", "x^(n+1)ln(x)/(n+1)-x^(n+1)/(n+1)^2") // non-numeric powers integration by parts required ADVANCED_INTEGRAL_TEST(ArctangentMultipliedByX, "x*arctg(x)", "x^2/2*arctg(x)-x/2+arctg(x)/2") - //ADVANCED_INTEGRAL_TEST(Arcsine, "arcsin(x)", "x*arcsin(x)+2*sqrt(1-x^2)") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(LogarithmWithX2Substitution, "x*ln(x^2+1)", "(x^2+1)ln(x^2+1)-x^2-1") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(CyclicIntegral1, "e^x*sin(2x)", "e^x*sin(2x)/5 - e^x*2cos(2x)/5") // cyclic integrals required - //ADVANCED_INTEGRAL_TEST(CyclicIntegral2, "sin(ln(x))", "-x/2(cos(ln(x))-sin(ln(x)))") // cyclic integrals required - //ADVANCED_INTEGRAL_TEST(CotangentDividedByLogOfSine, "ctg(x)/ln(sin(x))", "ln(ln(sin(x)))") // substitution t=ln(x) required - //ADVANCED_INTEGRAL_TEST(InvertedCosH, "1/(e^x+e^-x)", "arctg(e^x)") // extract_function required + // ADVANCED_INTEGRAL_TEST(Arcsine, "arcsin(x)", "x*arcsin(x)+2*sqrt(1-x^2)") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(LogarithmWithX2Substitution, "x*ln(x^2+1)", "(x^2+1)ln(x^2+1)-x^2-1") // substitution t=x^2+a required + // ADVANCED_INTEGRAL_TEST(CyclicIntegral1, "e^x*sin(2x)", "e^x*sin(2x)/5 - e^x*2cos(2x)/5") // cyclic integrals required + // ADVANCED_INTEGRAL_TEST(CyclicIntegral2, "sin(ln(x))", "-x/2(cos(ln(x))-sin(ln(x)))") // cyclic integrals required + // ADVANCED_INTEGRAL_TEST(CotangentDividedByLogOfSine, "ctg(x)/ln(sin(x))", "ln(ln(sin(x)))") // substitution t=ln(x) required + // ADVANCED_INTEGRAL_TEST(InvertedCosH, "1/(e^x+e^-x)", "arctg(e^x)") // extract_function required ADVANCED_INTEGRAL_TEST(Absolute, "abs(x)", "x*abs(x)/2") // abs required - ADVANCED_INTEGRAL_TEST(MaxOfOneAndSquare, "(abs(1-x^2)+1+x^2)/2", "(sgn(1-x^2)(x-x^3/3)+x+x^3/3)/2") // abs required + ADVANCED_INTEGRAL_TEST(MaxOfOneAndSquare, "(abs(1-x^2)+1+x^2)/2", + "(sgn(1-x^2)(x-x^3/3)+x+x^3/3)/2") // abs required // tasks 5-7 involve rational integrals // https://pages.mini.pw.edu.pl/~dembinskaa/www/?download=Inf_I_PowtKol3_2022-2023.pdf, task 11 ADVANCED_INTEGRAL_TEST(TangentSquared, "tg^2(x)", "tg(x)-arctg(tg(x))") // may sometimes fail ADVANCED_INTEGRAL_TEST(LongPolynomial, "(2x-3)^10", "1/22(2x-3)^11+3^11/22") - //ADVANCED_INTEGRAL_TEST(SquareRootInDenominator, "1/(2+sqrt(x))", "2sqrt(x)-4ln(sqrt(x)+2)") // substitution t=sqrt(x)+a required - //ADVANCED_INTEGRAL_TEST(InvertedCosH2, "1/cosh(x)", "2arctg(e^x)") // extract_funcion required + // ADVANCED_INTEGRAL_TEST(SquareRootInDenominator, "1/(2+sqrt(x))", "2sqrt(x)-4ln(sqrt(x)+2)") // substitution t=sqrt(x)+a required + // ADVANCED_INTEGRAL_TEST(InvertedCosH2, "1/cosh(x)", "2arctg(e^x)") // extract_funcion required ADVANCED_INTEGRAL_TEST(Sine5Cosine, "sin^5(x)cos(x)", "sin^6(x)/6") - //ADVANCED_INTEGRAL_TEST(ExpressionWithSquareInDenominator, "x/(x^2-1)^(3/2)", "-1/sqrt(x^2-1)") // substitution t=ax^b+c required - //ADVANCED_INTEGRAL_TEST(ArcsineWithDerivative, "arcsin^2(x)/sqrt(1-x^2)", "arcsin^3(x)/3") // substitution t=arcsin(x) required - //ADVANCED_INTEGRAL_TEST(RootInDenominator, "1/(x^(1/3)+1)", "3/2*x^(2/3)-3x^(1/3)+3ln(x^(1/3)+1)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) + // ADVANCED_INTEGRAL_TEST(ExpressionWithSquareInDenominator, "x/(x^2-1)^(3/2)", "-1/sqrt(x^2-1)") // substitution t=ax^b+c required + // ADVANCED_INTEGRAL_TEST(ArcsineWithDerivative, "arcsin^2(x)/sqrt(1-x^2)", "arcsin^3(x)/3") // substitution t=arcsin(x) required + // ADVANCED_INTEGRAL_TEST(RootInDenominator, "1/(x^(1/3)+1)", "3/2*x^(2/3)-3x^(1/3)+3ln(x^(1/3)+1)") // substitution t=ax^b+c required (dx = 1/b((t-c)/a)^((1-b)/b)) ADVANCED_INTEGRAL_TEST(Logarithm, "ln(x)", "x ln(x) - x") ADVANCED_INTEGRAL_TEST(XTimesCosine, "x*cos(x)", "x*sin(x)+cos(x)") - //ADVANCED_INTEGRAL_TEST(XTimesExponential, "x^2e^(1-x)", "-(x^2+2x+2)e^(1-x)") // linear substitution for nested-vars-only required - //ADVANCED_INTEGRAL_TEST(SquareTimesExponential, "8x^2e^(4-x^3)", "-8/3e^(4-x^3)") // substitution t=x^2+a required + ADVANCED_INTEGRAL_TEST(XTimesExponential, "x^2e^(1-x)", "-(x^2+2x+2)e^(1-x)") + // ADVANCED_INTEGRAL_TEST(SquareTimesExponential, "8x^2e^(4-x^3)", "-8/3e^(4-x^3)") // substitution t=x^2+a required - //ADVANCED_INTEGRAL_TEST(SineSquared, "sin^2(x)", "x/2-sin(2x)/4") // rational integrals required + // ADVANCED_INTEGRAL_TEST(SineSquared, "sin^2(x)", "x/2-sin(2x)/4") // rational integrals required ADVANCED_INTEGRAL_TEST(Sine5, "sin^5(x)", "-cos^5(x)/5+2/3cos^3(x)-cos(x)") ADVANCED_INTEGRAL_TEST(Sine4Cos3, "sin^4(x)cos^3(x)", "sin^5(x)/5-sin^7(x)/7") - //ADVANCED_INTEGRAL_TEST(Sine4Cos2, "sin^4(x)cos^2(x)", "1/6sin^5(x)cos(x)+1/24sin^3(x)cos(x)-1/32sin(x)cos(x)+1/16x") // rational integrals required - //ADVANCED_INTEGRAL_TEST(Cos3XCos5X, "cos(3x)cos(5x)", "1/16sin(8x)+1/4sin(2x)") // more trigonometric identitied required + // ADVANCED_INTEGRAL_TEST(Sine4Cos2, "sin^4(x)cos^2(x)", "1/6sin^5(x)cos(x)+1/24sin^3(x)cos(x)-1/32sin(x)cos(x)+1/16x") // rational integrals required + // ADVANCED_INTEGRAL_TEST(Cos3XCos5X, "cos(3x)cos(5x)", "1/16sin(8x)+1/4sin(2x)") // more trigonometric identitied required + + // own tests + ADVANCED_INTEGRAL_TEST(ETower, "e^x*e^e^x*e^e^e^x*e^e^e^e^x*e^e^e^e^e^x*e^e^e^e^e^e^x", + "e^e^e^e^e^e^x") } \ No newline at end of file