Skip to content

Commit

Permalink
Using integration to evaluate asin.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ZCG-coder
ZCG-coder committed Aug 11, 2024
1 parent c2eb37c commit 1ccb389
Show file tree
Hide file tree
Showing 3 changed files with 35 additions and 59 deletions.
61 changes: 16 additions & 45 deletions src/calc/trig/trig.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,7 +325,6 @@ namespace steppable::__internals::arithmetic
// 2 x
result = subtract(static_cast<std::string>(constants::PI_OVER_2), result, 0);
}
result = roundOff(result, decimals);

// Convert the result as needed.
switch (mode)
Expand All @@ -342,7 +341,7 @@ namespace steppable::__internals::arithmetic

if (_x.front() == '-')
result = "-" + result;
return result;
return roundOff(result, decimals);
}

std::string asin(const std::string& x, const int decimals, const int mode)
Expand All @@ -355,49 +354,21 @@ namespace steppable::__internals::arithmetic
if (compare(abs(x, 0), "0", 0) == "2")
return "0";

auto x2 = power(x, "2", 0);
auto x3 = multiply(x2, x, 0);
auto x5 = multiply(x3, x2, 0);
auto x7 = multiply(x5, x2, 0);
std::string onePlusX;
std::string oneMinusX;
std::string sqrtOnePlusX;
std::string sqrtOneMinusX;
std::string result;

// For x <= 0.5, use Taylor series.
// 3 5 7
// x 3x 5x
// arcsin(x) = x + --- + ---- + -----
// 6 40 112
if (compare(abs(x, 0), "0.5", 0) != "1")
{
result = add(x, divide(x3, "6", 0), 0);
result = add(result, divide(multiply(x5, "3", 0), "40", 0), 0);
result = add(result, divide(multiply(x7, "5", 0), "112", 0), 0);
goto out; // NOLINT(cppcoreguidelines-avoid-goto)
}
// / x
// | 1
// | ---------------- dy
// | /------|
// | / 2
// | \/ 1 - y
// / 0
auto integrand = [&](const std::string& y) {
auto y2 = power(y, "2", 0);
auto oneMinusY2 = subtract("1", y2, 0);
auto denominator = root(oneMinusY2, "2", decimals + 1);
return divide("1", denominator, 0, decimals + 1);
};
auto result = calculus::romberg(integrand, "0", x, 10, decimals + 2);

// Othman, S. B.; Bagul, Y. J. An Innovative Method for Approximating Arcsine Function. Preprints 2022,
// 2022070388. https://doi.org/10.20944/preprints202207.0388.v1
//
// /-----| /-----| 5 3
// arcsin(x) = 2 * 0.510774109 * (\/ 1 + x − \/ 1 − x ) + (0.239x − 0.138x + 0.005x)
//
// /-----| /-----| 5 3
// = 1.021548218 * (\/ 1 + x − \/ 1 − x ) + (0.239x − 0.138x + 0.005x)

onePlusX = add("1", x, 0);
oneMinusX = subtract("1", x, 0);
sqrtOnePlusX = root(onePlusX, "2", static_cast<long>(decimals) * 2);
sqrtOneMinusX = root(oneMinusX, "2", static_cast<long>(decimals) * 2);

result = multiply("1.021548218", subtract(sqrtOnePlusX, sqrtOneMinusX, 0), 0);
result = add(result,
add(multiply("0.239", x5, 0), subtract(multiply("0.138", x3, 0), multiply("0.005", x, 0), 0), 0),
0);

out:
switch (mode)
{
case 1:
Expand Down Expand Up @@ -486,7 +457,7 @@ int main(int _argc, const char* _argv[])
Utf8CodePage _;
ProgramArgs program(_argc, _argv);
program.addPosArg('c', $("trig", "47dcf91b-847c-48f0-9889-f5ce1b6831e3"), false);
program.addPosArg('n', $("trig", "bcd0a3e9-3d89-4921-94b3-d7533d60911f"), false);
program.addPosArg('n', $("trig", "bcd0a3e9-3d89-4921-94b3-d7533d60911f"));
program.addKeywordArg("mode", 0, $("trig", "03fdd1f2-6ea5-49d4-ac3f-27f01f04a518"));
program.addKeywordArg("decimals", 5, $("trig", "d1df3b60-dac1-496c-99bb-ba763dc551df"));
program.addSwitch("profile", false, $("trig", "162adb13-c4b2-4418-b3df-edb6f9355d64"));
Expand Down
29 changes: 17 additions & 12 deletions src/calculus/nInt/nInt.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
**************************************************************************************************/

#include "fn/basicArithm.hpp"
#include "output.hpp"
#include "rounding.hpp"

#include <cmath>
Expand All @@ -41,13 +42,13 @@ namespace steppable::__internals::calculus
const int decimals)
{
auto acc = "0." + std::string(decimals - 1, '0') + "1";
auto Rp = std::vector(max_steps, "0"s);
auto Rc = std::vector(max_steps, "0"s);
auto previous = std::vector(max_steps, "0"s);
auto current = std::vector(max_steps, "0"s);

auto h = subtract(b, a, 0);
auto fAB = add(f(a), f(b), 0);
auto halfH = multiply(h, "0.5", 0);
Rp.front() = multiply(halfH, fAB, 0);
previous.front() = multiply(halfH, fAB, 0);

for (int i = 1; i < max_steps; i++)
{
Expand All @@ -59,26 +60,30 @@ namespace steppable::__internals::calculus
auto d = multiply(std::to_string((2 * j) - 1), h, 0);
c = add(c, f(add(a, d, 0)), 0);
}
Rc.front() = add(multiply(h, c, 0), multiply("0.5", Rp.front(), 0), 0);
current.front() = add(multiply(h, c, 0), multiply("0.5", previous.front(), 0), 0);

for (int j = 1; j < (i + 1); j++)
{
long double n_k = pow(4, j);
auto one = multiply(std::to_string(n_k), Rc.at(j - 1), 0);
auto top = subtract(one, Rp.at(j - 1), 0);
Rc.at(j) = divide(top, std::to_string(n_k - 1), 0, decimals + 1);
auto one = multiply(std::to_string(n_k), current.at(j - 1), 0);
auto top = subtract(one, previous.at(j - 1), 0);
current.at(j) = divide(top, std::to_string(n_k - 1), 0, decimals + 1);

if (i > 1 and compare(abs(subtract(Rp.at(i - 1), Rc.at(i), 0), 0), acc, 0) == "0")
return roundOff(Rc.at(i), decimals);
if (i > 1 and compare(abs(subtract(previous.at(i - 1), current.at(i), 0), 0), acc, 0) == "0")
return roundOff(current.at(i), decimals);
}

std::ranges::swap(Rp, Rc);
std::ranges::swap(previous, current);
}

return roundOff(Rp.at(max_steps - 1), decimals);
return roundOff(previous.at(max_steps - 1), decimals);
}
} // namespace steppable::__internals::calculus

#ifndef NO_MAIN
int main() {}
int main()
{
steppable::output::error("calculus::nInt"s, "This program should not be run directly."s);
return 1;
}
#endif
4 changes: 2 additions & 2 deletions tests/testTrig.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,13 @@ _.assertIsEqual(tan("90", 2, 1), "Infinity");
SECTION_END()

SECTION(Test arc cosine)
_.assertIsEqual(acos("0.5", 2, 1), "60");
_.assertIsEqual(acos("0.5", 2, 1), "60.00");
// Zero check test
_.assertIsEqual(acos("1", 2, 1), "0");
SECTION_END()

SECTION(Test arc sine)
_.assertIsEqual(asin("0.5", 2, 1), "30.00");
_.assertIsEqual(asin("0.5", 0, 1), "30");
// Zero check test
_.assertIsEqual(asin("0", 2, 1), "0");
SECTION_END()
Expand Down

0 comments on commit 1ccb389

Please sign in to comment.