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
 
//*************************************************************************