Blog Entry © Thursday and Friday, March 12 – 13, 2026, Naive Power Series Computations of the Cosine and Sine Trigonometric Functions and Polynomial Fit to the Trig Functions Using Least Squares

// TrigLSPolyFit.cpp (c) Thursday, March 12, 2026
// by James Pate Williams, Jr.
// Least Squares Curve Fitting to the Trigonometric
// functions Cosine and Sine

#include <iomanip>
#include <iostream>
#include <vector>
#include "PolyLSFit.h"

static const int HornerN = 50;
double even_factors[HornerN], odd_factors[HornerN + 1];

static double factorial(int n)
{
	double nfactorial = 1.0;

	for (int i = 2; i <= n; i++)
		nfactorial *= i;
	return nfactorial;
}

static void computeEvenFactors()
{
	int one = 1;

	for (int n = 0; n < 50; n += 2)
	{
		even_factors[n] = one / factorial(n);
		one = -one;
	}
}

static void computeOddFactors()
{
	int one = 1;

	for (int n = 1; n < 51; n += 2)
	{
		odd_factors[n] = one / factorial(n);
		one = -one;
	}
}

static double fx_cos(double x, int n)
{
	double sum = even_factors[n];

	for (int i = n - 1; i >= 0; i--)
	{
		sum = sum * x + even_factors[i];
	}

	return sum;
}

static double fx_sin(double x, int n)
{
	double sum = odd_factors[n];

	for (int i = n - 1; i >= 0; i--)
	{
		sum = sum * x + odd_factors[i];
	}

	return sum;
}

static const int M = 33, N = 100, N21 = N + N + 1;
std::vector<double> cosX(N21), cosY(N21), cosCoeff;
std::vector<double> sinX(N21), sinY(N21), sinCoeff;
std::vector<double> pX(M), pY(M);

static void LSPolyFit(int HornerN, int degree)
{
	double pi = 4.0 * atan(1.0);
	double deltaX = 2.0 * pi / N;
	double X = -pi;

	for (int i = 0; i < N21; i++)
	{
		cosX[i] = X;
		cosY[i] = fx_cos(X, HornerN);
		sinX[i] = X;
		sinY[i] = fx_sin(X, HornerN + 1);
		X += deltaX;
	}

	PolyLSFit::leastSquares(cosX, cosY, cosCoeff, degree, N21);
	PolyLSFit::leastSquares(sinX, sinY, sinCoeff, degree, N21);

	std::cout << std::setw(6) << "x";
	std::cout << std::setw(17) << "fit cos";
	std::cout << std::setw(17) << "C++ cos";
	std::cout << std::setw(17) << "% Error";
	std::cout << std::endl;

	deltaX = 2.0 * pi / (M - 1);
	X = -pi;

	for (int i = 0; i < M; i++)
	{
		double Y = cosCoeff[degree];

		for (int j = degree - 1; j >= 0; j--)
		{
			Y = Y * X + cosCoeff[j];
		}

		pX[i] = X;
		pY[i] = Y;

		std::cout << std::fixed << std::setprecision(3);
		std::cout << std::setw(6);
		std::cout << X << ' ';

		if (Y >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);
		
		std::cout << fabs(Y) << ' ';
		double cx = cos(X);

		if (cx >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);
		
		std::cout << fabs(cx) << ' ';
		double pcDifference = 0.0;

		if (fabs(cx) >= 1.0e-12)
			pcDifference = 100.0 * fabs(Y - cx) / fabs(cx);
		else
			pcDifference = -1.0;

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		if (pcDifference != -1.0)
			std::cout << pcDifference << std::endl;
		else
			std::cout << "not a number" << std::endl;

		X += deltaX;
	}

	std::cout << std::endl;

	std::cout << std::setw(6) << "x";
	std::cout << std::setw(17) << "fit sin";
	std::cout << std::setw(17) << "C++ sin";
	std::cout << std::setw(17) << "% Error";
	std::cout << std::endl;

	deltaX = 2.0 * pi / (M - 1);
	X = -pi;

	for (int i = 0; i < M; i++)
	{
		double Y = sinCoeff[degree];

		for (int j = degree - 1; j >= 0; j--)
		{
			Y = Y * X + sinCoeff[j];
		}

		pX[i] = X;
		pY[i] = Y;

		std::cout << std::fixed << std::setprecision(3);
		std::cout << std::setw(6);
		std::cout << X << ' ';

		if (Y >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		std::cout << fabs(Y) << ' ';
		double sx = sin(X);

		if (sx >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		std::cout << fabs(sx) << ' ';
		double pcDifference = 0.0;

		if (fabs(sx) >= 1.0e-12)
			pcDifference = 100.0 * fabs(Y - sx) / fabs(sx);
		else
			pcDifference = -1.0;

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		if (pcDifference != -1.0)
			std::cout << pcDifference << std::endl;
		else
			std::cout << "not a number" << std::endl;

		X += deltaX;
	}

	std::cout << std::endl;
}

int main()
{
	computeEvenFactors();
	computeOddFactors();
	std::cout << "# terms 16" << std::endl;
	LSPolyFit(HornerN, 16);
	std::cout << "# terms 20" << std::endl;
	LSPolyFit(HornerN, 20);
	std::cout << "# terms 24" << std::endl;
	LSPolyFit(HornerN, 24);
	return 0;
}