PolynomRealRoots.h

Go to the documentation of this file.

/* Copyright (C) 2020 Facility for Antiproton and Ion Research in Europe,
   Darmstadt SPDX-License-Identifier: GPL-3.0-only
   Authors: Pierre-Alain Loizeau [committer] */
 
#include <cmath>
 
//*************************************************************************
float polinom(int n, float x, float* k)
{
    float s = 1;
    for (int i = n - 1; i >= 0; i--)
        s = s * x + k[i];
    return s;
}   // polinom
 
float dihot(int degree, float edgeNegativ, float edgePositiv, float* kf)
{
    for (;;) {
        float x = 0.5 * (edgeNegativ + edgePositiv);
        if (std::abs(x - edgeNegativ) < 1e-3 || std::abs(x - edgePositiv) < 1e-3)
            return x;
        if (polinom(degree, x, kf) < 0)
            edgeNegativ = x;
        else
            edgePositiv = x;
    }
}   // dihot
 
void stepUp(int level, float** A, float** B, int* currentRootsCount)
{
 
    float major = 0;
    for (int i = 0; i < level; i++) {
        float s = fabs(A[level][i]);
        if (s > major)
            major = s;
    }
    major += 1.0;
 
    currentRootsCount[level] = 0;
 
    for (int i = 0; i <= currentRootsCount[level - 1]; i++) {
        int signLeft, signRight;
        float edgeLeft, edgeRight;
        float edgeNegativ, edgePositiv;
 
        if (i == 0)
            edgeLeft = -major;
        else
            edgeLeft = B[level - 1][i - 1];
 
        float rb = polinom(level, edgeLeft, A[level]);
 
        if (rb == 0) {
            B[level][currentRootsCount[level]] = edgeLeft;
            currentRootsCount[level]++;
            continue;
        }
 
        if (rb > 0)
            signLeft = 1;
        else
            signLeft = -1;
 
        if (i == currentRootsCount[level - 1])
            edgeRight = major;
        else
            edgeRight = B[level - 1][i];
 
        rb = polinom(level, edgeRight, A[level]);
 
        if (rb == 0) {
            B[level][currentRootsCount[level]] = edgeRight;
            currentRootsCount[level]++;
            continue;
        }
 
        if (rb > 0)
            signRight = 1;
        else
            signRight = -1;
        if (signLeft == signRight)
            continue;
 
        if (signLeft < 0) {
            edgeNegativ = edgeLeft;
            edgePositiv = edgeRight;
        } else {
            edgeNegativ = edgeRight;
            edgePositiv = edgeLeft;
        }
 
        B[level][currentRootsCount[level]] = dihot(level, edgeNegativ, edgePositiv, A[level]);
        currentRootsCount[level]++;
    }
    return;
}   // stepUp
 
void polynomRealRoots(float* rootsArray, int n, float* kf_, int& rootsCount)
{
 
    float* kf = new float[n + 1];
 
    float** A = new float*[n + 1];
    float** B = new float*[n + 1];
 
    int* currentRootsCount = new int[n + 1];
 
    for (int i = 1; i <= n; i++) {
        A[i] = new float[i];
        B[i] = new float[i];
    }
 
    for (int i = 0; i <= n; i++)
        kf[i] = kf_[n - i];
 
    for (int i = 0; i < n; i++)
        A[n][i] = kf[i] / kf[n];
 
    for (int i1 = n, i = n - 1; i > 0; i1 = i, i--) {
        for (int j1 = i, j = i - 1; j >= 0; j1 = j, j--) {
            A[i][j] = A[i1][j1] * j1 / i1;
        }
    }
 
    currentRootsCount[1] = 1;
    B[1][0] = -A[1][0];
 
    for (int i = 2; i <= n; i++)
        stepUp(i, A, B, currentRootsCount);
 
    rootsCount = currentRootsCount[n];
    for (int i = 0; i < rootsCount; i++)
        rootsArray[i] = B[n][i];
 
    for (int i = 1; i <= n; i++) {
        delete[] A[i];
        delete[] B[i];
    }
    delete[] A;
    delete[] B;
    delete[] currentRootsCount;
    delete[] kf;
}   // polynomRealRoots
 
//*************************************************************************