5

I have attempted to implement the sine function in C, yet I am getting weird results. Here are the three functions I am using to calculate sine:

#define PI 3.14159265358979323846
#define DEPTH 16

double sine(long double);
long double pow(long double, unsigned int);
unsigned int fact(unsigned int);

double sine(long double x) {
    long double i_x = x *= PI/180;
    int n = 3, d = 0, sign = -1;

    // fails past 67 degrees
    for (; d < DEPTH; n += 2, d++, sign *= -1) {
        x += pow(i_x, n) / fact(n) * sign;
    }

    return x;
}

long double pow(long double base, unsigned int exp) {
    double answer = 1;
    while (exp) {
        answer *= base;
        exp--;
    }
    return answer;
}

unsigned int fact(unsigned int n) {
    unsigned int answer = 1;
    while (n > 1) {
        answer *= n--;
    }
    return answer;
}

To test it I have been testing it against the built in sine function as follows:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

main() {
    for (int i = 0; i <= 180; i++) {
        printf("sin(%i) = %lf, %lf\n", i, sine(i), sin(i*3.14159265358979323846/180));
    }

    exit(EXIT_SUCCESS);
}

Up through 67 degrees, it calculates the same as the built in function. Though, as it increases past 67 it typically strays gets further and further from the actual value.

Here is an example output:

>> sin(100) = 0.987711, 0.984808
>> sin(101) = 0.986885, 0.981627
>> sin(102) = 0.987056, 0.978148
>> sin(103) = 0.988830, 0.974370
>> sin(104) = 0.993060, 0.970296
>> sin(105) = 1.000948, 0.965926
>> sin(106) = 1.014169, 0.961262
>> sin(107) = 1.035052, 0.956305
>> sin(108) = 1.066807, 0.951057
>> sin(109) = 1.113846, 0.945519
>> sin(110) = 1.182194, 0.939693
>> sin(111) = 1.280047, 0.933580
>> sin(112) = 1.418502, 0.927184
>> sin(113) = 1.612527, 0.920505
>> sin(114) = 1.882224, 0.913545
>> sin(115) = 2.254492, 0.906308
>> sin(116) = 2.765192, 0.898794
>> sin(117) = 3.461969, 0.891007
              ...
>> sin(180) = 8431648.192239, 0.000000

Does anybody know why this is happening? I am using Visual Studio 2017 on Windows 7, if that provides any useful information.

ndim
  • 35,870
  • 12
  • 47
  • 57
Robbie
  • 1,733
  • 2
  • 13
  • 20
  • What happens between 80 and 100: does it increase to 1, then decrease again? –  Jul 23 '17 at 03:37
  • `printf` `%lf` is pointless. It's either equivalent to `%f` (C99 and later) or undefined behavior (C89). – melpomene Jul 23 '17 at 03:38
  • 3
    Why are you redefining the standard `pow` function? – melpomene Jul 23 '17 at 03:38
  • 6
    Factorial 12 (12!) is the largest factorial that can be represented in a 32-bit integer; 20! is the largest that can be represented in a 64-bit integer. One plausible source of trouble is your `fact()` function. – Jonathan Leffler Jul 23 '17 at 03:40
  • What is `DEPTH`? Your iteration may be too short, that far out from 0 (but given the above comment about the limits of your factorial function, it may be too long). –  Jul 23 '17 at 03:41
  • It does not, it continues to increase @Evert – Robbie Jul 23 '17 at 03:44
  • Does that mean it's not (anywhere near) 1 at 90 degrees? –  Jul 23 '17 at 03:44
  • I was not intending to, I was attempting to write my own power function @melpomene – Robbie Jul 23 '17 at 03:44
  • 3
    The standard defines `pow()` and declares it in ``. Use a different name for safety. – Jonathan Leffler Jul 23 '17 at 03:46
  • It is still very close. It is accurate up to the hundreds place. Also, DEPTH is defined as 16 (its included in the code I posted) @Evert – Robbie Jul 23 '17 at 03:48
  • It appears that decreasing the DEPTH to around 5 improves accuracy 100% up until 100 degrees. Past there, the accuracy begins to decrease very slightly, but it is very close. For example, it says that sine(180) = -0.000445, which is very close to 0. – Robbie Jul 23 '17 at 03:53
  • 2
    Please note that using `long double` gives you no extra precision as Microsoft's C Runtime Library maps `long double` to `double` and hence only 64-bit double precision floating point as outlined in IEEE 754 standard. – m0h4mm4d Jul 23 '17 at 04:07
  • 1
    2 major problems: [factional overflow](https://stackoverflow.com/questions/45261259/implementation-of-sine-function-in-c-not-working#comment77485990_45261259) and lack of range reduction. For ever larger values of `x` , Taylor's series fail unless more and more terms used. Use `remquo(x,45.0, &quo)` to bring into +/- 45 degree range and then apply trig identities and the series solutions. [`remquo()` example](https://stackoverflow.com/a/31525208/2410359) – chux - Reinstate Monica Jul 24 '17 at 02:05

5 Answers5

7

Each time your for loop progresses, n is increased by 2 and hence for DEPTH = 16, near the end of loop you have to calculate factorials of numbers as big as 30 and you are using unsigned int that can only store values as big as 2^32 = 4294967296 ~= 12! and this causes overflow in your factorial function which in turn gives you the wrong factorial.

Even if you used long double for it and I already stated in my comments that long double in MSCRT is mapped to double (Reference) you'd still see some anomalies probably at larger angles because although double can store values as big as 1.8E+308 but it loses its granularity at 2^53 = 9007199254740992 ~= 18! (i.e. 2^53 + 1 stored as a double is equal to 2^53). So once you go up in angles, the effect of this behavior becomes larger and larger to the point that it is noticeable in the 6 decimal precision that you are using with printf().

Although you are on the right track, you should use a bignum library like GMP or libcrypto. They can perform these calculations without the loss of precision.

BTW, since your are developing on Windows 7 that means you are either using x86 or x86-64. On these platforms, x87 is capable of performing extended precision (as per 754 standard) operations with 80 bits but I am unaware of compiler intrinsics that can give you that capability without resorting to assembly code.

I also would like to direct your attention to range reduction techniques. Although I still recommend using bignum libs, if you are good between 0 and 90 degrees (0 and 45 if I'm to be more strict), you can compute the sine() of all other angles just by simple trigonometric identities.

UPDATE:

Actually I'm gonna correct myself about using doubles in factorial calculations. After writing a simple program I verified that when I usedouble to store factorials, they are correct even if I go upper than 18. After giving it some thought I realized that in the case of factorials, the situation with double's granularity is a little bit more complex. I'll give you an example to make it clear:

19! = 19 * 18 * ... * 2 * 1

in this number 18, 16, 14, ... , 2 are all multiples of 2 and since a multiplication by 2 is equivalent to a shift to the left in binary representation, all lower bits in 19! are already 0 and hence when double's rounding kicks in for integers greater than 2^53, these factorials are unaffected. You can compute the number of least significant zeroes in the binary representation of 19! by counting the number of 2's which is 16. (for 20!, it is 18)

I'm gonna go up to 1.8e+308 and check if all the factorials are unaffected or not. I'll update you with the results.

UPDATE 2:

If we use doubles to hold factorials, they are affected by rounding from 23! onward. It can be easily shown, because 2^74 < 23! < 2^75 which means that at least 75 bits of precision is required to represent it, but since 23! has 19 least significant bits with the value of 0, so it needs 75 - 19 = 56 which is larger than 53 bits provided by double.

For 22!, it is 51 bits (you can calculate it yourself).

m0h4mm4d
  • 400
  • 4
  • 12
4

There are multiple issues in your code:

  • You redefine the standard function pow() with a different prototype. This may cause problems when you link the program as en executable. Use a different anme such as pow_int.

  • You should define the pow_int and fact functions as static before the sine function. It may allow for better optimisation at compile time.

  • Indeed fact is limited by the range of type unsigned int which is much less than the precision of type long double. Factorials beyond 12 have an incorrect value, causing a loss of precision.

  • You could actually compute the terms incrementally, saving a lot of computations and avoiding potential loss of precision.

  • The prototype for main() without arguments is int main(void)

  • The computation of PI/180 is performed as double, which is less precise than long double. You should write the expression as x = x * PI / 180;

  • DEPTH should be increased to improve the precision. At least 4 more terms bring a substantial improvement.

  • You should apply a range reduction: taking advantage of the sine function symmetric and periodic nature, computation can be performed with fewer terms on x modulo 90 or even 45 degrees.

Here is a modified version:

#include <stdio.h>
#include <math.h>

#define PI_L   3.14159265358979323846264338327950288L
#define PI     3.14159265358979323846264338327950288
#define DEPTH  24

double sine(long double x) {
    long double res, term, x2, t1;
    int phase;

    x = remquol(x, 90, &phase);
    if (phase & 1)
        x = 90 - x;

    x = x * PI_L / 180; // convert x to radians
    x2 = x * x;         // pre-compute x^2

    // compute the sine series: x - x^3/3! + x^5/5! ...
    res = term = x;   // the first term is x
    for (int n = 1; n < DEPTH; n += 4) {
        // to reduce precision loss, compute 2 terms for each iteration
        t1 = term * x2 / ((n + 1) * (n + 2));
        term = t1 * x2 / ((n + 3) * (n + 4));
        // update the result with the difference of the terms
        res += term - t1;
    }
    if (phase & 2)
        res = -res;

    return (double)res;
}

int main(void) {
    printf("deg            sin                  sine         delta\n\n");
    for (int i = 0; i <= 360; i += 10) {
        double s1 = sin(i * PI / 180);
        double s2 = sine(i);
        printf("%3i  %20.17f  %20.17f  %g\n", i, s1, s2, s2 - s1);
    }
    return 0;
}

The output is:

deg            sin                  sine         delta

  0   0.00000000000000000   0.00000000000000000  0
 10   0.17364817766693033   0.17364817766693036  2.77556e-17
 20   0.34202014332566871   0.34202014332566871  0
 30   0.49999999999999994   0.50000000000000000  5.55112e-17
 40   0.64278760968653925   0.64278760968653936  1.11022e-16
 50   0.76604444311897801   0.76604444311897801  0
 60   0.86602540378443860   0.86602540378443860  0
 70   0.93969262078590832   0.93969262078590843  1.11022e-16
 80   0.98480775301220802   0.98480775301220802  0
 90   1.00000000000000000   1.00000000000000000  0
100   0.98480775301220802   0.98480775301220802  0
110   0.93969262078590843   0.93969262078590843  0
120   0.86602540378443882   0.86602540378443860  -2.22045e-16
130   0.76604444311897812   0.76604444311897801  -1.11022e-16
140   0.64278760968653947   0.64278760968653936  -1.11022e-16
150   0.49999999999999994   0.50000000000000000  5.55112e-17
160   0.34202014332566888   0.34202014332566871  -1.66533e-16
170   0.17364817766693025   0.17364817766693036  1.11022e-16
180   0.00000000000000012  -0.00000000000000000  -1.22465e-16
190  -0.17364817766693047  -0.17364817766693036  1.11022e-16
200  -0.34202014332566866  -0.34202014332566871  -5.55112e-17
210  -0.50000000000000011  -0.50000000000000000  1.11022e-16
220  -0.64278760968653925  -0.64278760968653936  -1.11022e-16
230  -0.76604444311897790  -0.76604444311897801  -1.11022e-16
240  -0.86602540378443837  -0.86602540378443860  -2.22045e-16
250  -0.93969262078590821  -0.93969262078590843  -2.22045e-16
260  -0.98480775301220802  -0.98480775301220802  0
270  -1.00000000000000000  -1.00000000000000000  0
280  -0.98480775301220813  -0.98480775301220802  1.11022e-16
290  -0.93969262078590854  -0.93969262078590843  1.11022e-16
300  -0.86602540378443860  -0.86602540378443860  0
310  -0.76604444311897812  -0.76604444311897801  1.11022e-16
320  -0.64278760968653958  -0.64278760968653936  2.22045e-16
330  -0.50000000000000044  -0.50000000000000000  4.44089e-16
340  -0.34202014332566855  -0.34202014332566871  -1.66533e-16
350  -0.17364817766693127  -0.17364817766693036  9.15934e-16
360  -0.00000000000000024   0.00000000000000000  2.44929e-16

As can be seen above, the sine() function seems more precise than the standard sin function on my system: sin(180 * M_PI / 128) should be 0 precisely. Similarly, sin(150 * M_PI / 128) should be 0.5 exactly.

chqrlie
  • 131,814
  • 10
  • 121
  • 189
  • Impressive implementation. If I may ask though, why do you have a hoisted declaration immediately above the definition? – Patrick Roberts Jul 24 '17 at 14:41
  • 1
    @PatrickRoberts: good question: `clang -Weverything` issues a warning for functions with external linkage for which it has not seen a prototype. The purpose of this warning it the ensure the header file with the declaration is duly included in the file with the implementation and thus consistency between the prototype and the definition can be verified. It is not needed here, so I shall remove it. – chqrlie Jul 24 '17 at 15:02
3

Your way of polynomial series evaluation is numerically unstable. Try horner's method which is more stable than power calculations.

arbitUser1401
  • 575
  • 2
  • 8
  • 25
2

Your problem is here:

for (; d < DEPTH; n += 2, d++, sign *= -1) {
    x += pow(i_x, n) / fact(n) * sign;
}

You are using d < DEPTH in error when it should be n < DEPTH, d is irrelevant to your computations within the loop. The following should work -- although I have not compiled to test.

for (; n < DEPTH; n += 2, sign *= -1) {
    x += pow(i_x, n) / fact(n) * sign;
}

note: a DEPTH of 12 (e.g. Taylor Series expansion with terms 1, 3, 5, ... 11) is sufficient for 3e-10 error -- 3 ten billionths at 60-degrees. (though error increases as angle increases between 0-360, a DEPTH of 20 will keep error less than 1.0e-8 over the entire range.)

Enabling compiler warnings would have caught the unused d in sine.

Here is an example of code with the changes (note: Gnu provides a constant M_PI for PI):

#include <stdio.h>
#include <stdint.h>
#include <math.h>

#define DEPTH 16

/* n factorial */
uint64_t nfact (int n)
{
    if (n <= 0) return 1;

    uint64_t s = n;

    while (--n)
        s *= n;

    return s;
}

/* y ^ x */
double powerd (const double y, const int x)
{
    if (!x) return 1;

    double r = y;

    for (int i = 1; i < x; i++)
        r *= y;

    return r;
}

double sine (double deg)
{
    double rad = deg * M_PI / 180.0,
        x = rad;
    int sign = -1;

    for (int n = 3; n < DEPTH; n += 2, sign *= -1)
        x += sign * powerd (rad, n) / nfact (n);

    return x;
}

int main (void) {

    printf (" deg       sin            sine\n\n");
    for (int i = 0; i < 180; i++)
        printf ("%3d    %11.8f    %11.8f\n", i, sin (i * M_PI / 180.0), sine (i));

    return 0;
}

Example Use/Output

$ ./bin/sine
 deg       sin            sine

  0     0.00000000     0.00000000
  1     0.01745241     0.01745241
  2     0.03489950     0.03489950
  3     0.05233596     0.05233596
  4     0.06975647     0.06975647
  5     0.08715574     0.08715574
  6     0.10452846     0.10452846
  7     0.12186934     0.12186934
  8     0.13917310     0.13917310
  9     0.15643447     0.15643447
 10     0.17364818     0.17364818
 11     0.19080900     0.19080900
 12     0.20791169     0.20791169
 13     0.22495105     0.22495105
 14     0.24192190     0.24192190
 15     0.25881905     0.25881905
 16     0.27563736     0.27563736
 17     0.29237170     0.29237170
 18     0.30901699     0.30901699
 19     0.32556815     0.32556815
 20     0.34202014     0.34202014
 21     0.35836795     0.35836795
 22     0.37460659     0.37460659
 23     0.39073113     0.39073113
 24     0.40673664     0.40673664
 25     0.42261826     0.42261826
 26     0.43837115     0.43837115
 27     0.45399050     0.45399050
 28     0.46947156     0.46947156
 29     0.48480962     0.48480962
 30     0.50000000     0.50000000
 31     0.51503807     0.51503807
 32     0.52991926     0.52991926
 33     0.54463904     0.54463904
 34     0.55919290     0.55919290
 35     0.57357644     0.57357644
 36     0.58778525     0.58778525
 37     0.60181502     0.60181502
 38     0.61566148     0.61566148
 39     0.62932039     0.62932039
 40     0.64278761     0.64278761
 41     0.65605903     0.65605903
 42     0.66913061     0.66913061
 43     0.68199836     0.68199836
 44     0.69465837     0.69465837
 45     0.70710678     0.70710678
 46     0.71933980     0.71933980
 47     0.73135370     0.73135370
 48     0.74314483     0.74314483
 49     0.75470958     0.75470958
 50     0.76604444     0.76604444
 51     0.77714596     0.77714596
 52     0.78801075     0.78801075
 53     0.79863551     0.79863551
 54     0.80901699     0.80901699
 55     0.81915204     0.81915204
 56     0.82903757     0.82903757
 57     0.83867057     0.83867057
 58     0.84804810     0.84804810
 59     0.85716730     0.85716730
 60     0.86602540     0.86602540
 61     0.87461971     0.87461971
 62     0.88294759     0.88294759
 63     0.89100652     0.89100652
 64     0.89879405     0.89879405
 65     0.90630779     0.90630779
 66     0.91354546     0.91354546
 67     0.92050485     0.92050485
 68     0.92718385     0.92718385
 69     0.93358043     0.93358043
 70     0.93969262     0.93969262
 71     0.94551858     0.94551858
 72     0.95105652     0.95105652
 73     0.95630476     0.95630476
 74     0.96126170     0.96126170
 75     0.96592583     0.96592583
 76     0.97029573     0.97029573
 77     0.97437006     0.97437006
 78     0.97814760     0.97814760
 79     0.98162718     0.98162718
 80     0.98480775     0.98480775
 81     0.98768834     0.98768834
 82     0.99026807     0.99026807
 83     0.99254615     0.99254615
 84     0.99452190     0.99452190
 85     0.99619470     0.99619470
 86     0.99756405     0.99756405
 87     0.99862953     0.99862953
 88     0.99939083     0.99939083
 89     0.99984770     0.99984770
 90     1.00000000     1.00000000
 91     0.99984770     0.99984770
 92     0.99939083     0.99939083
 93     0.99862953     0.99862953
 94     0.99756405     0.99756405
 95     0.99619470     0.99619470
 96     0.99452190     0.99452190
 97     0.99254615     0.99254615
 98     0.99026807     0.99026807
 99     0.98768834     0.98768834
100     0.98480775     0.98480775
101     0.98162718     0.98162718
102     0.97814760     0.97814760
103     0.97437006     0.97437006
104     0.97029573     0.97029573
105     0.96592583     0.96592583
106     0.96126170     0.96126170
107     0.95630476     0.95630476
108     0.95105652     0.95105652
109     0.94551858     0.94551858
110     0.93969262     0.93969262
111     0.93358043     0.93358043
112     0.92718385     0.92718385
113     0.92050485     0.92050485
114     0.91354546     0.91354546
115     0.90630779     0.90630779
116     0.89879405     0.89879405
117     0.89100652     0.89100652
118     0.88294759     0.88294759
119     0.87461971     0.87461971
120     0.86602540     0.86602540
121     0.85716730     0.85716730
122     0.84804810     0.84804810
123     0.83867057     0.83867057
124     0.82903757     0.82903757
125     0.81915204     0.81915204
126     0.80901699     0.80901699
127     0.79863551     0.79863551
128     0.78801075     0.78801075
129     0.77714596     0.77714596
130     0.76604444     0.76604444
131     0.75470958     0.75470958
132     0.74314483     0.74314482
133     0.73135370     0.73135370
134     0.71933980     0.71933980
135     0.70710678     0.70710678
136     0.69465837     0.69465836
137     0.68199836     0.68199835
138     0.66913061     0.66913060
139     0.65605903     0.65605902
140     0.64278761     0.64278760
141     0.62932039     0.62932038
142     0.61566148     0.61566146
143     0.60181502     0.60181501
144     0.58778525     0.58778523
145     0.57357644     0.57357642
146     0.55919290     0.55919288
147     0.54463904     0.54463901
148     0.52991926     0.52991924
149     0.51503807     0.51503804
150     0.50000000     0.49999996
151     0.48480962     0.48480958
152     0.46947156     0.46947152
153     0.45399050     0.45399045
154     0.43837115     0.43837109
155     0.42261826     0.42261820
156     0.40673664     0.40673657
157     0.39073113     0.39073105
158     0.37460659     0.37460651
159     0.35836795     0.35836786
160     0.34202014     0.34202004
161     0.32556815     0.32556804
162     0.30901699     0.30901686
163     0.29237170     0.29237156
164     0.27563736     0.27563720
165     0.25881905     0.25881887
166     0.24192190     0.24192170
167     0.22495105     0.22495084
168     0.20791169     0.20791145
169     0.19080900     0.19080873
170     0.17364818     0.17364788
171     0.15643447     0.15643414
172     0.13917310     0.13917274
173     0.12186934     0.12186895
174     0.10452846     0.10452803
175     0.08715574     0.08715526
176     0.06975647     0.06975595
177     0.05233596     0.05233537
178     0.03489950     0.03489886
179     0.01745241     0.01745170

Error Check based on DEPTH

In response to the comment regarding computing the error, you investigate the error associated with Taylor-Series expansions for both sin and cos base on the number of terms by varying DEPTH and setting an max error of EMAX 1.0e-8 using something similar to the following for the range of 0-360 (or 0-2PI),

#define DEPTH 20
#define EMAX 1.0e-8
...
/* sine as above */
...
/* cos with taylor series expansion to n = DEPTH */
long double cose (const long double deg)
{
    long double rad = deg * M_PI / 180.0,
        x = 1.0;
    int sign = -1;

    for (int n = 2; n < DEPTH; n += 2, sign *= -1)
        x += sign * powerd (rad, n) / nfact (n);

    return x;
}

int main (void) {

    for (int i = 0; i < 180; i++) {
        long double sinlibc = sin (i * M_PI / 180.0),
            coslibc = cos (i * M_PI / 180.0),
            sints = sine (i),
            costs = cose (i),
            serr = fabs (sinlibc - sints),
            cerr = fabs (coslibc - costs);

        if (serr > EMAX)
            fprintf (stderr, "sine error exceeds limit of %e\n"
            "%3d    %11.8Lf    %11.8Lf    %Le\n",
            EMAX, i, sinlibc, sints, serr);

        if (cerr > EMAX)
            fprintf (stderr, "cose error exceeds limit of %e\n"
            "%3d    %11.8Lf    %11.8Lf    %Le\n",
            EMAX, i, coslibc, costs, cerr);
    }

    return 0;
}

If you check, you will find that for anything less than DEPTH 20 (10 terms in each expansion), error will exceed 1.0e-8 for higher angles. Surprisingly, the expansions are very accurate over the first quadrant for values of DEPTH as low as 12 (6-terms).

Addemdum - Improved Taylor-Series Accuracy Using 0-90 & Quadrants

In the normal Taylor-Series expansion, error grows as angle grows. And... because some just can't not tinker, I wanted to further compare accuracy between the libc sin/cos and the Taylor-Series if computations were limited to 0-90 degrees and the remainder of the period from 90-360 were handled by quadrant (2, 3 & 4) mirroring of results from 0-90. It works -- marvelously.

For example, the results of handing only angles 0-90 and bracketing angles between 90 - 180, 180 - 270 and 270 - 360 with an initial angle % 360 produces results comparable to the libc math lib functions. The maximum error between the libc and 8 & 10 term Taylor-Series expansions were, respectably:

Max Error from libc sin/cos

With TSLIM 16

sine_ts max err at :  90.00 deg  -- 6.023182e-12
cose_ts max err at : 270.00 deg  -- 6.513370e-11

With TSLIM 20

sine_ts max err at : 357.00 deg  -- 5.342948e-16
cose_ts max err at : 270.00 deg  -- 3.557149e-15

(with a large number of angles showing no difference at all)

The tweaked versions of sine and cose with Taylor-Series were as follows:

double sine (const double deg)
{
    double fp = deg - (int64_t)deg, /* save fractional part of deg */
        qdeg = (int64_t)deg % 360,  /* get equivalent 0-359 deg angle */
        rad, sine_deg;              /* radians, sine_deg */
    int pos_quad = 1,               /* positive quadrant flag 1,2  */
        sign = -1;                  /* taylor series term sign */

    qdeg += fp;                     /* add fractional part back to angle */

    /* get equivalent 0-90 degree angle, set pos_quad flag */
    if (90 < qdeg && qdeg <= 180)           /* in 2nd quadrant */
        qdeg = 180 - qdeg;
    else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
        qdeg = qdeg - 180;
        pos_quad = 0;
    }
    else if (270 < qdeg && qdeg <= 360) {   /* in 4th quadrant */
        qdeg = 360 - qdeg;
        pos_quad = 0;
    }

    rad = qdeg * M_PI / 180.0;      /* convert to radians */
    sine_deg = rad;                 /* save copy for computation */

    /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
    for (int n = 3; n < TSLIM; n += 2, sign *= -1) {
        double p = rad;
        uint64_t f = n;

        for (int i = 1; i < n; i++)     /* pow */
            p *= rad;

        for (int i = 1; i < n; i++)     /* nfact */
            f *= i;

        sine_deg += sign * p / f;       /* Taylor-series term */
    }

    return pos_quad ? sine_deg : -sine_deg;
}

and for cos

double cose (const double deg)
{
    double fp = deg - (int64_t)deg, /* save fractional part of deg */
        qdeg = (int64_t)deg % 360,  /* get equivalent 0-359 deg angle */
        rad, cose_deg = 1.0;        /* radians, cose_deg */
    int pos_quad = 1,               /* positive quadrant flag 1,4  */
        sign = -1;                  /* taylor series term sign */

    qdeg += fp;                     /* add fractional part back to angle */

    /* get equivalent 0-90 degree angle, set pos_quad flag */
    if (90 < qdeg && qdeg <= 180) {         /* in 2nd quadrant */
        qdeg = 180 - qdeg;
        pos_quad = 0;
    }
    else if (180 < qdeg && qdeg <= 270) {   /* in 3rd quadrant */
        qdeg = qdeg - 180;
        pos_quad = 0;
    }
    else if (270 < qdeg && qdeg <= 360)     /* in 4th quadrant */
        qdeg = 360 - qdeg;

    rad = qdeg * M_PI / 180.0;      /* convert to radians */

    /* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
    for (int n = 2; n < TSLIM; n += 2, sign *= -1) {
        double p = rad;
        uint64_t f = n;

        for (int i = 1; i < n; i++)     /* pow */
            p *= rad;

        for (int i = 1; i < n; i++)     /* nfact */
            f *= i;

        cose_deg += sign * p / f;       /* Taylor-series term */
    }

    return pos_quad ? cose_deg : -cose_deg;
}

Rabbit-trail end found...

David C. Rankin
  • 81,885
  • 6
  • 58
  • 85
  • David, would you please elaborate on how you reached that `3e-10` number? – m0h4mm4d Jul 23 '17 at 12:52
  • 1
    Sure, that was just a quick check at 60 degrees of `fabs (sin (60 * M_PI/180.0) - sine (60))` with `sine` above and a `DEPTH` of `12`. The error gets progressively larger as angle increases (which would make breaking the computation of `sine` from `0-360` in to quadrants much like is done for `atan` advantageous). A `DEPTH` of `20` will keep error less than `1.0e-8` over the entire range `0-360`. Use of `double` or `long double` is irrelevant`. – David C. Rankin Jul 24 '17 at 01:18
  • Thanks for the clarification but since the library implementation is not the true value of `sin` (it's just a different approximation and probably a better one) I'd say calling it `difference` is better than `error`, cause it might imply that `3e-10` is a theoretically calculated difference from the true value of `sin`. Anyway, thanks again for making it clear to us. – m0h4mm4d Jul 24 '17 at 06:02
  • Hey, good catch. The error reporting is only as good as the difference between what libc `sin` thinks the **sine** of the angle is and what the manual Taylor series with `DEPTH / 2` terms thinks it is. In the computer world it's all approximation when it comes down to it. With a pencil and paper and the *unit circle* we can get exact closed form representations, With computers, its just the product of all the combined errors inherent to the process along the way `:)` – David C. Rankin Jul 24 '17 at 06:40
0

Changing the angle range in main to -90 to 90 will still cover the whole sine range. but as the Taylor serie is starting from zero the DEPTH value can be reduced to 7. As earlier mentioned, making the fact function 64 bits unsigned will fix the 67 degree problem.