Skip to content

Commit 1ccb389

Browse files
committed
Using integration to evaluate asin.
1 parent c2eb37c commit 1ccb389

File tree

3 files changed

+35
-59
lines changed

3 files changed

+35
-59
lines changed

src/calc/trig/trig.cpp

Lines changed: 16 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -325,7 +325,6 @@ namespace steppable::__internals::arithmetic
325325
// 2 x
326326
result = subtract(static_cast<std::string>(constants::PI_OVER_2), result, 0);
327327
}
328-
result = roundOff(result, decimals);
329328

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

343342
if (_x.front() == '-')
344343
result = "-" + result;
345-
return result;
344+
return roundOff(result, decimals);
346345
}
347346

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

358-
auto x2 = power(x, "2", 0);
359-
auto x3 = multiply(x2, x, 0);
360-
auto x5 = multiply(x3, x2, 0);
361-
auto x7 = multiply(x5, x2, 0);
362-
std::string onePlusX;
363-
std::string oneMinusX;
364-
std::string sqrtOnePlusX;
365-
std::string sqrtOneMinusX;
366-
std::string result;
367-
368-
// For x <= 0.5, use Taylor series.
369-
// 3 5 7
370-
// x 3x 5x
371-
// arcsin(x) = x + --- + ---- + -----
372-
// 6 40 112
373-
if (compare(abs(x, 0), "0.5", 0) != "1")
374-
{
375-
result = add(x, divide(x3, "6", 0), 0);
376-
result = add(result, divide(multiply(x5, "3", 0), "40", 0), 0);
377-
result = add(result, divide(multiply(x7, "5", 0), "112", 0), 0);
378-
goto out; // NOLINT(cppcoreguidelines-avoid-goto)
379-
}
357+
// / x
358+
// | 1
359+
// | ---------------- dy
360+
// | /------|
361+
// | / 2
362+
// | \/ 1 - y
363+
// / 0
364+
auto integrand = [&](const std::string& y) {
365+
auto y2 = power(y, "2", 0);
366+
auto oneMinusY2 = subtract("1", y2, 0);
367+
auto denominator = root(oneMinusY2, "2", decimals + 1);
368+
return divide("1", denominator, 0, decimals + 1);
369+
};
370+
auto result = calculus::romberg(integrand, "0", x, 10, decimals + 2);
380371

381-
// Othman, S. B.; Bagul, Y. J. An Innovative Method for Approximating Arcsine Function. Preprints 2022,
382-
// 2022070388. https://doi.org/10.20944/preprints202207.0388.v1
383-
//
384-
// /-----| /-----| 5 3
385-
// arcsin(x) = 2 * 0.510774109 * (\/ 1 + x − \/ 1 − x ) + (0.239x − 0.138x + 0.005x)
386-
//
387-
// /-----| /-----| 5 3
388-
// = 1.021548218 * (\/ 1 + x − \/ 1 − x ) + (0.239x − 0.138x + 0.005x)
389-
390-
onePlusX = add("1", x, 0);
391-
oneMinusX = subtract("1", x, 0);
392-
sqrtOnePlusX = root(onePlusX, "2", static_cast<long>(decimals) * 2);
393-
sqrtOneMinusX = root(oneMinusX, "2", static_cast<long>(decimals) * 2);
394-
395-
result = multiply("1.021548218", subtract(sqrtOnePlusX, sqrtOneMinusX, 0), 0);
396-
result = add(result,
397-
add(multiply("0.239", x5, 0), subtract(multiply("0.138", x3, 0), multiply("0.005", x, 0), 0), 0),
398-
0);
399-
400-
out:
401372
switch (mode)
402373
{
403374
case 1:
@@ -486,7 +457,7 @@ int main(int _argc, const char* _argv[])
486457
Utf8CodePage _;
487458
ProgramArgs program(_argc, _argv);
488459
program.addPosArg('c', $("trig", "47dcf91b-847c-48f0-9889-f5ce1b6831e3"), false);
489-
program.addPosArg('n', $("trig", "bcd0a3e9-3d89-4921-94b3-d7533d60911f"), false);
460+
program.addPosArg('n', $("trig", "bcd0a3e9-3d89-4921-94b3-d7533d60911f"));
490461
program.addKeywordArg("mode", 0, $("trig", "03fdd1f2-6ea5-49d4-ac3f-27f01f04a518"));
491462
program.addKeywordArg("decimals", 5, $("trig", "d1df3b60-dac1-496c-99bb-ba763dc551df"));
492463
program.addSwitch("profile", false, $("trig", "162adb13-c4b2-4418-b3df-edb6f9355d64"));

src/calculus/nInt/nInt.cpp

Lines changed: 17 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
**************************************************************************************************/
2222

2323
#include "fn/basicArithm.hpp"
24+
#include "output.hpp"
2425
#include "rounding.hpp"
2526

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

4748
auto h = subtract(b, a, 0);
4849
auto fAB = add(f(a), f(b), 0);
4950
auto halfH = multiply(h, "0.5", 0);
50-
Rp.front() = multiply(halfH, fAB, 0);
51+
previous.front() = multiply(halfH, fAB, 0);
5152

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

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

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

75-
std::ranges::swap(Rp, Rc);
76+
std::ranges::swap(previous, current);
7677
}
7778

78-
return roundOff(Rp.at(max_steps - 1), decimals);
79+
return roundOff(previous.at(max_steps - 1), decimals);
7980
}
8081
} // namespace steppable::__internals::calculus
8182

8283
#ifndef NO_MAIN
83-
int main() {}
84+
int main()
85+
{
86+
steppable::output::error("calculus::nInt"s, "This program should not be run directly."s);
87+
return 1;
88+
}
8489
#endif

tests/testTrig.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -45,13 +45,13 @@ _.assertIsEqual(tan("90", 2, 1), "Infinity");
4545
SECTION_END()
4646

4747
SECTION(Test arc cosine)
48-
_.assertIsEqual(acos("0.5", 2, 1), "60");
48+
_.assertIsEqual(acos("0.5", 2, 1), "60.00");
4949
// Zero check test
5050
_.assertIsEqual(acos("1", 2, 1), "0");
5151
SECTION_END()
5252

5353
SECTION(Test arc sine)
54-
_.assertIsEqual(asin("0.5", 2, 1), "30.00");
54+
_.assertIsEqual(asin("0.5", 0, 1), "30");
5555
// Zero check test
5656
_.assertIsEqual(asin("0", 2, 1), "0");
5757
SECTION_END()

0 commit comments

Comments
 (0)