Math Function Solver & Optimizer

The program will approximate all the roots and local optima of a given one-dimensional function over a specific domain using the method described below.

The root-finding method is based on a loop that finds initial values for classical algorithms used to locate roots and local optima. The loop will use values represented as 32-bit floating-point numbers (float) and convert them to 64-bit floats (double) using the command:
xv_double = (double) xv;

Basic Idea:
We define a 32-bit unsigned integer variable (e.g., intvalue) and increment it in integer intervals. The increment intdiff is determined using a method explained below (after the declaration of general_solver_optimizer). After each increment, we convert intvalue into a float xv using:
memcpy((void *)&xv, (void *)&intvalue, sizeof(float));

Then we convert xv to 64-bit using:
xv_double = (double) xv;
Throughout the loop, we maintain variables that store the two previous values of xv_double (e.g., prev_prev_xv_double, prev_xv_double), the current function value at xv_double (e.g., current_value), and the function value at the previous point (prev_value).

You will also need arrays of double that store root intervals and optima intervals.

Main Loop Structure:
The main loop (initialized beforehand) proceeds as follows:
prev_prev_xv_double = prev_xv_double;
prev_xv_double = xv_double;

prev_prev_value = prev_value;
prev_value = current_value;
Calculate the current xv_double and its function value current_value.

Then perform the following checks:

If current_value and prev_value have opposite signs → a root interval is found.

If current_value is smaller than a given epsilon → a single-point root interval is found.

If the function was decreasing, and:
prev_prev_value > prev_value;
prev_value < current_value;
then a minimum interval is found.

If the function was increasing, and:
prev_prev_value < prev_value;
prev_value > current_value;
then a maximum interval is found.

End of the Loop:
Run a classical root-finding algorithm on each root interval (e.g., bisection or Steffensen on the midpoint).

Run a classical optimization algorithm on each optima interval. For maxima, apply the method to the function multiplied by -1. Possible methods: parabolic interpolation, golden section, or quadratic fit (QF).

After finding the optima, check if the absolute function value at each optimum is not smaller than a given epsilon — this helps find roots of functions that merely touch the x-axis without crossing. Add those to the roots array if not already found.

Technical Considerations:
When finding optima, if the difference between current_value and prev_value is smaller than 1e-11, it’s better to ignore this measurement. One way is to treat it as neither rising nor falling.

Once a root is found, ignore nearby points that represent the same solution. Resume checking only after encountering points that aren’t considered close approximations.

You may want to explicitly handle the case where 0.0 is a root.

You'll need to define a series of data structures to implement all these ideas.

Function Declaration:
extern void general_solver_optimizer(FUN_PTR objfun, int num_of_checks,
double roots_arr[], int *nr,
double optimum_arr[], int direction_arr[],
int *no, double solution_epsilon, double low_range, double high_range);

The integer intervals (intdiff) are calculated by:
nbits = 32 - (int) round(log(num_of_checks) / log(2));
intdiff = (int) round(exp(log(2.0) * nbits));

Example 1 – Main Program:
int main() {
double roots_arr[MAX_SOLUTIONS_NUMBER];
double optimum_arr[MAX_SOLUTIONS_NUMBER];
int direction_arr[MAX_SOLUTIONS_NUMBER];

int nr, no, i;
double solution_epsilon;
int num_of_checks;

solution_epsilon = 1.0e-6;
num_of_checks = 1048576;

general_solver_optimizer(f, num_of_checks, roots_arr, &nr,
optimum_arr, direction_arr, &no, solution_epsilon, -11.0, 10.0);

printf("\n\n %d roots:\n\n", nr);
for (i = 0; i < nr; i++)
printf("\nx* = %lg, function value = %lg\n", roots_arr[i], f(roots_arr[i]));

printf("\n\n %d optimum points:\n\n", no);
for (i = 0; i < no; i++)
if (direction_arr[i] == 1)
printf("\nx* = %lg: maximum, function value = %lg\n", optimum_arr[i], f(optimum_arr[i]));
else
printf("\nx* = %lg: minimum, function value = %lg\n", optimum_arr[i], f(optimum_arr[i]));
}

Example 1 – Objective Function:
double f(double x) {
double temp;
temp = sin(3.14159265 * x);
return temp;
}

Sample Output:
20 roots:

x* = -10, function value = 3.58979e-008
x* = -9, function value = -3.23081e-008
...
x* = 9, function value = 3.23081e-008

21 optimum points:

x* = -10.5016: minimum, function value = -0.999988
x* = -9.50156: maximum, function value = 0.999988
...
x* = 9.5: minimum, function value = -1


Example 2 – Objective Function:
double f(double x) {
double temp;
temp = x*x*x - 9.0 * x*x + 26*x - 24.0;
return temp;
}

Sample Output:
3 roots:

x* = 2, function value = 0
x* = 3, function value = 0
x* = 4, function value = 0

2 optimum points:

x* = 2.42265: maximum, function value = 0.3849
x* = 3.57735: minimum, function value = -0.3849



Example 3 – Objective Function:
double f(double x) {
double temp;
temp = (x - 3.0) * (0.25*x*x*x*x - 3.0*x*x*x + 13.0*x*x - 24.0*x);
return temp;
}

Sample Output:
2 roots:

x* = -0.483609, function value = -4.88869e-007
x* = 6.48361, function value = 4.55032e-010

3 optimum points:

x* = 2: minimum, function value = -31
x* = 3: maximum, function value = -30.75
x* = 4: minimum, function value = -31