This week's book giveaway is in the Cloud/Virtualization forum.
We're giving away four copies of Mastering Corda: Blockchain for Java Developers and have Jamiel Sheikh on-line!
See this thread for details.
Win a copy of Mastering Corda: Blockchain for Java Developers this week in the Cloud/Virtualization forum!
  • Post Reply Bookmark Topic Watch Topic
  • New Topic
programming forums Java Mobile Certification Databases Caching Books Engineering Micro Controllers OS Languages Paradigms IDEs Build Tools Frameworks Application Servers Open Source This Site Careers Other all forums
this forum made possible by our volunteer staff, including ...
Marshals:
  • Campbell Ritchie
  • Paul Clapham
  • Ron McLeod
  • Bear Bibeault
  • Liutauras Vilda
Sheriffs:
  • Jeanne Boyarsky
  • Tim Cooke
  • Junilu Lacar
Saloon Keepers:
  • Tim Moores
  • Tim Holloway
  • Stephan van Hulst
  • Jj Roberts
  • Carey Brown
Bartenders:
  • salvin francis
  • Frits Walraven
  • Piet Souris

How to calculate roots of a higher order polynomial in C

 
Greenhorn
Posts: 1
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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
 
Marshal
Posts: 71682
312
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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.
 
Saloon Keeper
Posts: 12608
273
  • Mark post as helpful
  • send pies
    Number of slices to send:
    Optional 'thank-you' note:
  • Quote
  • Report post to moderator
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.
 
Don't get me started about those stupid light bulbs.
reply
    Bookmark Topic Watch Topic
  • New Topic