Reputation:
Making c/++ function to compute sin(x)
Referancing by the formula;
But i don't know where to stop infinity sum, and i don't want libraries for quick-answer. I want to make it. I made factorial function well. (fact() function)
double sin(int x) {
int tp=0;
for(int k=1; ;k++)
{
tp+=pow(-1,k)*pow(x,2*k+1)/fact(2*k+1);
}
return tp;
}
Upvotes: 2
Views: 304
Reputation:
I ran the loop by 8 times, it's approaching to true value if you improve loop timing. Thanks @Sreeram TP
double sin(int x) {
int tp=0;
for(int k=1;k<8 ;k++)
{
tp+=pow(-1,k)*pow(x,2*k+1)/fact(2*k+1);
}
return tp;
}
Upvotes: 0
Reputation: 15045
Before we go any further, you must use double
for return types and accumulation, otherwise the result will be (heavily) truncated.
You should stop when adding the change term does not change the current accumulated summation value, i.e. there is underflow. This occurs when:
The magnitude of the change term is less than the value of the lowest binary digit of the accumulator.
To check this, it is useful to extract the exponent of the floating point number. This SO post provides a useful function: frexp
.
We can then check that the exponent of the change term is smaller than the exponent of the accumulator minus the number of bits in the mantissa of the double
type (52 bits). If you need to know more, search for IEEE floating point format.
Another point to improve numerical accuracy and speed: let's examine the relationship between successive terms:
So instead of calculating each term from scratch (using the power and factorial functions), we can calculate each term iteratively as we go.
Sample code:
// best to use an existing C-library define for this, if exists
#define DBL_MANTISSA_BIT 52
double mysin_new(double x)
{
double term = x;
double sum = term;
for (int n = 1; ; n++) {
// next term in the series
term *= - (x / (2*n)) * (x / (2*n+1));
// check if the exponent is small enough for us to stop
int e_term, e_sum;
frexp(term, &e_term);
frexp(sum, &e_sum);
if (e_term <= e_sum - DBL_MANTISSA_BIT) break;
sum += term;
}
return sum;
}
Some numerical tests:
angle = 0.157080
library : 0.15643446504023086896f, 0011111111000100000001100000101101100111101010000101001101110101b
mysin_old: 0.15643446504023086896f, 0011111111000100000001100000101101100111101010000101001101110101b
mysin_new: 0.15643446504023089672f, 0011111111000100000001100000101101100111101010000101001101110110b
angle = 0.314159
library : 0.30901699437494739575f, 0011111111010011110001101110111100110111001011111110100101001111b
mysin_old: 0.30901699437494745126f, 0011111111010011110001101110111100110111001011111110100101010000b
mysin_new: 0.30901699437494739575f, 0011111111010011110001101110111100110111001011111110100101001111b
angle = 0.471239
library : 0.45399049973954674897f, 0011111111011101000011100010111000101011010001001101111000000000b
mysin_old: 0.45399049973954674897f, 0011111111011101000011100010111000101011010001001101111000000000b
mysin_new: 0.45399049973954674897f, 0011111111011101000011100010111000101011010001001101111000000000b
angle = 0.628319
library : 0.58778525229247313710f, 0011111111100010110011110010001100000100011101010101101001011110b
mysin_old: 0.58778525229247313710f, 0011111111100010110011110010001100000100011101010101101001011110b
mysin_new: 0.58778525229247313710f, 0011111111100010110011110010001100000100011101010101101001011110b
angle = 0.785398
library : 0.70710678118654746172f, 0011111111100110101000001001111001100110011111110011101111001100b
mysin_old: 0.70710678118654746172f, 0011111111100110101000001001111001100110011111110011101111001100b
mysin_new: 0.70710678118654746172f, 0011111111100110101000001001111001100110011111110011101111001100b
angle = 0.942478
library : 0.80901699437494745126f, 0011111111101001111000110111011110011011100101111111010010101000b
mysin_old: 0.80901699437494734024f, 0011111111101001111000110111011110011011100101111111010010100111b
mysin_new: 0.80901699437494734024f, 0011111111101001111000110111011110011011100101111111010010100111b
angle = 1.099557
library : 0.89100652418836778779f, 0011111111101100100000110010000000011101001111010010110001101100b
mysin_old: 0.89100652418836767676f, 0011111111101100100000110010000000011101001111010010110001101011b
mysin_new: 0.89100652418836767676f, 0011111111101100100000110010000000011101001111010010110001101011b
angle = 1.256637
library : 0.95105651629515353118f, 0011111111101110011011110000111000010011010001000101010011111111b
mysin_old: 0.95105651629515353118f, 0011111111101110011011110000111000010011010001000101010011111111b
mysin_new: 0.95105651629515353118f, 0011111111101110011011110000111000010011010001000101010011111111b
angle = 1.413717
library : 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
mysin_old: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
mysin_new: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
angle = 1.570796
library : 1.00000000000000000000f, 0011111111110000000000000000000000000000000000000000000000000000b
mysin_old: 1.00000000000000022204f, 0011111111110000000000000000000000000000000000000000000000000001b
mysin_new: 1.00000000000000000000f, 0011111111110000000000000000000000000000000000000000000000000000b
angle = 1.727876
library : 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
mysin_old: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
mysin_new: 0.98768834059513777035f, 0011111111101111100110110010010010010100001011111110010001011100b
angle = 1.884956
library : 0.95105651629515364220f, 0011111111101110011011110000111000010011010001000101010100000000b
mysin_old: 0.95105651629515364220f, 0011111111101110011011110000111000010011010001000101010100000000b
mysin_new: 0.95105651629515375323f, 0011111111101110011011110000111000010011010001000101010100000001b
angle = 2.042035
library : 0.89100652418836789881f, 0011111111101100100000110010000000011101001111010010110001101101b
mysin_old: 0.89100652418836800983f, 0011111111101100100000110010000000011101001111010010110001101110b
mysin_new: 0.89100652418836800983f, 0011111111101100100000110010000000011101001111010010110001101110b
angle = 2.199115
library : 0.80901699437494745126f, 0011111111101001111000110111011110011011100101111111010010101000b
mysin_old: 0.80901699437494756229f, 0011111111101001111000110111011110011011100101111111010010101001b
mysin_new: 0.80901699437494756229f, 0011111111101001111000110111011110011011100101111111010010101001b
angle = 2.356194
library : 0.70710678118654757274f, 0011111111100110101000001001111001100110011111110011101111001101b
mysin_old: 0.70710678118654768376f, 0011111111100110101000001001111001100110011111110011101111001110b
mysin_new: 0.70710678118654757274f, 0011111111100110101000001001111001100110011111110011101111001101b
angle = 2.513274
library : 0.58778525229247324813f, 0011111111100010110011110010001100000100011101010101101001011111b
mysin_old: 0.58778525229247324813f, 0011111111100010110011110010001100000100011101010101101001011111b
mysin_new: 0.58778525229247324813f, 0011111111100010110011110010001100000100011101010101101001011111b
angle = 2.670354
library : 0.45399049973954685999f, 0011111111011101000011100010111000101011010001001101111000000010b
mysin_old: 0.45399049973954685999f, 0011111111011101000011100010111000101011010001001101111000000010b
mysin_new: 0.45399049973954691550f, 0011111111011101000011100010111000101011010001001101111000000011b
angle = 2.827433
library : 0.30901699437494750677f, 0011111111010011110001101110111100110111001011111110100101010001b
mysin_old: 0.30901699437494745126f, 0011111111010011110001101110111100110111001011111110100101010000b
mysin_new: 0.30901699437494745126f, 0011111111010011110001101110111100110111001011111110100101010000b
angle = 2.984513
library : 0.15643446504023097998f, 0011111111000100000001100000101101100111101010000101001101111001b
mysin_old: 0.15643446504023095223f, 0011111111000100000001100000101101100111101010000101001101111000b
mysin_new: 0.15643446504023095223f, 0011111111000100000001100000101101100111101010000101001101111000b
angle = 3.141593
library : 0.00000000000000012246f, 0011110010100001101001100000000000000000000000000000000000000000b
mysin_old: 0.00000000000000023566f, 0011110010110000111110110010110111010110111000110001010101001000b
mysin_new: 0.00000000000000023566f, 0011110010110000111110110010110111010110111000110001010101001001b
The new adaptive method wins out some of the times compared to the original, draws most of them, and loses a number of others. (Note it is assumed that the library method sin
is perfectly accurate upto the precision of double
)
Upvotes: 0
Reputation: 340055
Posted for reference, here's my version that uses progressive calculation of the various components of the series expansion, instead of calling out to an external O(n)
factorial function:
double mysin(double x)
{
int sign = 1;
int iteration = 1;
double xn = x;
double factorial = 1.0;
double result = 0.0;
double delta;
while (iteration < 31) {
delta = xn / factorial;
result += sign * delta;
xn *= x * x;
sign = -sign;
factorial *= ++iteration;
factorial *= ++iteration;
}
return result;
}
From my tests there doesn't appear to be any significant advantage in accuracy (or performance) by jumping out early if fabs(delta)
becomes too small, unless the initial value of x
is also very small.
Upvotes: 0
Reputation: 73444
You should stop when the value computed (tp
) is not changed much. Check for the difference of the newly computed value with the old one, and if the absolute difference (use fabs()
for that) is less than epsilon (you define that), stop iterating.
When you stop, you will end up with an approximation of the result of sin.
Moreover, int tp=0;
doesn't seem correct to me. You should initialize it with the formula result, when k
is 0. And more importantly, it should be of type double
, since int
is for integers.
Furthermore, not only must the value x
be in radians, it must also typically be in the range -π to +π.
Upvotes: 2
Reputation: 3416
What you need is an epsilon value, a value that you say when I am more accurate than that number we want to stop. It would look something like if(abs(prev_ans - cur_ans) < epsilon) break;
In math words the n
and n+1
partial sum of the series vary by so little that continuing the series wouldn't produce any more desired accuracy
Upvotes: 0