How to calculate roots of a higher order polynomial in C

[john josef]

I want to calculate roots of higher order polynomial (5th order). The equation has complex coefficients. I have a python program which does that (calculates l1 to l5). I want to know if it could be done using C. I used gsl library but in the HPC that I am working doesn't provide support to gsl libraries. Also, it was able to calculate the roots of equations with real coefficients only.

Other way to is store values of l1 to l5 from python code and then use it in the C program which has exactly similar structure. But the files generated are really huge. I think it would increase the overall time of execution of the C program.

Please help me out here. The python program is below:
https://bit.ly/37BQcq2

from mpi4py import MPI
comm = MPI.COMM_WORLD
proc_id = comm.Get_rank()
n_procs = comm.Get_size()

import time
import numpy as np
from numpy.lib.scimath import sqrt
import scipy.special as sp

delta = 10.0;
lund = 1e04;
alpha = 1e-7;
g = 10.0;
eta = 1.0;
VA = lund*eta/delta;
k0 = sqrt(6)/(1.0*delta);
Omega = VA*k0*10;
ratio = 500000.0;
beta = 1e05;#VA*VA*k0*k0/(g*alpha*ratio*ratio);

t0 = VA*VA/(1.0*g*alpha*beta*eta)
t11 = 4*Omega*Omega/(1.0*eta*VA**2*k0**4);
teta = 1.0/(1.0*eta*k0*k0);

grid_point = 10000;

ks_max = 100.0/(1.0*delta);
kz_max = 100.0/(1.0*delta);
ks_min = 1e-10;
kz_min = 1e-10;    
     
kz_div = 1000001
ks_div = 10001
ks_inc = (ks_max-ks_min)/(1.0*ks_div)#increment in ks
kz_inc = (kz_max-kz_min)/(1.0*kz_div)#increment in kz

x=np.loadtxt("xh1.txt");
z=np.loadtxt("zh2.txt");

start = time.time()

work_size = (kz_div) // n_procs
extra_work = (kz_div) % n_procs
my_work = work_size + (1 if proc_id<extra_work else 0)
l_start = work_size * proc_id + (proc_id if proc_id<extra_work else extra_work)
l_end = l_start + my_work

for l in range(l_start, l_end):
       
       if (l == 0 or l == kz_div):
           p = 1
       elif (l % 2 != 0):
           p = 4
       else:
           p = 2

       kz = kz_min + l * kz_inc
       
       for i in range(int(ks_div)):
           
           if (i == 0 or i == ks_div):
               q = 1
           elif (i % 2 != 0):
               q = 4
           else:
               q = 2
   
           ks = ks_min + i * ks_inc        
     
           k = sqrt(ks**2+kz**2);    
   
           omgM = VA*kz;
           omge = eta*k**2;
           omgA = sqrt(g * alpha*beta)*ks/k;
           omgO = Omega*kz/k;
           a1 = 1j*1;
           a2 = 2*omge;
           a3 = -1j*(omgA**2+omge**2+2*omgM**2+4*omgO**2);
           a4 = -2*omge*(omgA**2+omgM**2+4*omgO**2);
           a5 = 1j*(omgM**4+omgA**2*(omge**2+omgM**2)+4*omge**2*omgO**2);
           a6 = omgA**2*omge*omgM**2;

           eqn =[a1,a2,a3,a4,a5,a6];
           roots = np.sort_complex(np.roots(eqn));
           
           l2 = roots[0];
           l4 = roots[1];
           l1 = roots[4];
           l3 = roots[3];
           l5 = roots[2];

Thank you

[Campbell Ritchie]

Welcome to the Ranch

Of course you can calculate the roots in many languages. Don't however try translating language to language. Write down the algorithm and implement it in C anew. Please explain all abbreviations including gsl and HPC.

Moving to our C forum.

[Stephan van Hulst]

Welcome to CodeRanch!

That's an absolutely terribly written application. I hope you didn't pay money for the training program that produced that code.

Anyway, if you can't rely on libraries to do the heavy lifting for you, you must implement an iterative root finding algorithm. Look at the pointers here: https://en.wikipedia.org/wiki/Root-finding_algorithms#Roots_of_polynomials

Once you understand the algorithm, you can write it down in absolutely any programming language. If you need help on specific parts, tell us what you're stuck on.