Lab 8: Threaded sine function

sine.cpp

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <unistd.h>
#include <pthread.h>
#define EPS 1e-16L

using namespace std;
typedef long double ld;

void* sinl_function(void*);
int threads;
ld *sum, x2t;

ld sinl_threaded(ld x) {

	// initialize variables
	pthread_t *pthreads = new pthread_t[threads];
	ld x2 = x * x;
	sum = new ld[threads];
	x2t = pow(x2, threads); // (x^2)^threads

	// create threads
	for (int i = 0; i < threads; ++i) {
		sum[i] = x;
		int error = pthread_create(pthreads + i, NULL, sinl_function, reinterpret_cast<void*>(i));
		if (error) printf("Thread %d returned error code(%d)\n", i, error);
		x *= -x2 / (2*i + 2) / (2*i + 3); // base jumps by two
	}

	// join threads
	for (int i = 0; i < threads; ++i)
		pthread_join(pthreads[i], NULL);

	// add values
	ld total = 0;
	for (int i = 0; i < threads; ++i)
		total += sum[i];

	// clean up
	delete[] pthreads, sum;
	return total;

}

int main(int argc, char *args[]) {
	if (argc < 3)
		printf("You are missing some arguments.\nFormat is \"%s N X\", where N and X are the number of threads and the real number for the sin() function respectively.\n", args[0]);
	else {
		int N; double X;
		sscanf(args[1], "%d", &N);
		sscanf(args[2], "%lf", &X); // scan as double, process as long double
		::threads = N;
		printf("cmath  sin(%lf): %.16lf\n", X, (double) sinl(X));
		printf("thread sin(%lf): %.16lf\n", X, (double) sinl_threaded(X));
	}
}

void* sinl_function(void* args) {
	
	// get the real values of the arguments
	int i = *(int*) &args;
	ld& sum = ::sum[i];

	// make sign alternate if number of threads is odd
	int sign = threads % 2 == 0 ? 1 : -1;
	
	// current sum is x^(2*i + 1) / (2*i + 1)!
	// iterate values after sum
	for (ld term = sum;; sum += term, i += threads) {

		// x^(2*i + 1) / (2*i + 1)! * factor = x^(2*(i + t) + 1) / (2*(i + t) + 1)!
		// factor = x^(2*(i + t) + 1 - 2*i - 1) / (2*i + 2) / (2*i + 3) / ... / (2*i + 2*t + 1)
		// factor = x^(2*t) / (2*i + 2) / (2*i + 3) / ... / (2*i + 2*t + 1)
		ld factor = x2t;
		for (int j = 2*i + 2; j <= 2*(i + threads) + 1; ++j)
			factor /= j;
		
		// check if term is okay
		term *= factor * sign;
		if (abs(term) < EPS) break;
		// printf("%d: x^%d / %d!\n", cur, sign * (2*i + 1), 2*i + 1);
		// sleep(1);
	}
}