# Cubic function problem

while writing the program, I ran into an error that I can’t fix, the problem is on lines 35 and 63. I would ask for help.

``````subroutine cubic(a,b,c,d,x)

real,INTENT(IN) :: a,b,c,d
complex,INTENT(INOUT) :: x(*)
real,parameter:: pi = 4.0d0*atan(1.0d0)
real :: p,q,DD,phi,temp1,y1,y2,y3,temp2,u,v,y2r,y2i

p = (3.0*c/a-(b/a)**2.0)/3.0
q = (2.0*(b/a)**3.0-9.0*b*c/a/a + 27.0*d/a)/27.0

DD = (p/3.0)**3.0 + (q/2.0)**2.0

if(DD < 0.0)then

phi=acos(-q/2.0/sqrt(abs(p)**3/27.0))
temp1=2.0*sqrt(abs(p)/3.0)
y1=temp1*cos(phi/3.0)
y2=-temp1*cos((phi+pi)/3.0)
y3=-temp1*cos((phi-pi)/3.0)
else

temp1=-q/2.0+sqrt(DD)
temp2=-q/2.0-sqrt(DD)
u=abs(temp1)**(1.0/3.0)
v=abs(temp2)**(1.0/3.0)
if(temp1 < 0.0) u=-u
if(temp2 < 0.0) v=-v
y1=u+v
y2r=-(u+v)/2.0
y2i=(u-v)*sqrt(3.0)/2.0
end if

temp1=b/a/3.0
y1=y1-temp1
y2=y2-temp1
y3=y3-temp1
y2r=y2r-temp1

if(DD < 0.0) then
x(1)=cmplx(y1,0.0)
x(2)=cmplx(y2,0.0)
x(3)=cmplx(y3,0.0)
else if(DD > 0.0) then
x(1)=cmplx(y1,0.0)
x(2)=cmplx(y2r,y2i)
x(3)=cmplx(y2r,-y2i)
else
x(1)=cmplx(y1,0.0)
x(2)=cmplx(y2r,0.0)
x(3)=cmplx(y2r,0.0)
endif
end subroutine cubic

program cubicfunction
implicit none

real :: a,b,c,d
complex :: x(1:3)

print*,"Cubic function: ax^3 + bx^2 + cx + d = 0"

call cubic(a,b,c,d,x)
print *, "x1 =",x(1)
print *, "x2 =",x(2)
print *, "x3 =",x(3)

end program cubicfunction
``````

when i try to activate the program i get this error Well, if `dd > 0`, then `y2` is not set to a value when you get to line 35.

Edit: you should also use integer exponents when possible. Oh, and `y3` is not set if `dd > 0`.

How should i set this in dd > 0?

You need to re-arrange the code. Something lik this untested modification.

``````program cubicfunction

use iso_c_binding, only : knd => c_float

implicit none

real(knd) :: a, b, c, d
complex(knd) :: x(3)

print *, "cubic function: ax^3 + bx^2 + cx + d = 0"

call cubic(a,b,c,d,x)

print *, "x1 =", x(1)
print *, "x2 =", x(2)
print *, "x3 =", x(3)

contains

subroutine cubic(a,b,c,d,x)

real(knd), intent(in) :: a, b, c, d
complex(knd), intent(inout) :: x(:)

real(knd), parameter:: pi = 4 * atan(1._knd)
real(knd) p, q, dd, phi, temp1, y1, y2, y3, temp2, u, v, y2r, y2i, ba

interface
function cbrt(x) bind(c,name='cbrt') result(r)
import knd
real(knd) r
real(knd), intent(in) :: x
end function cbrt
end interface

ba = b / a
p = (3 * c / a - ba**2) /3
q = (2 * ba**3 - 9 * ba * c / a + 27 * d/ a) / 27

dd = (p / 3)**3 + (q / 2)**2

if (dd < 0) then
phi = acos(- q / 2 / sqrt(abs(p)**3 / 27))
temp1 = 2 * sqrt(abs(p) / 3)
y1 = temp1 * cos(phi / 3)
y2 = -temp1 * cos((phi + pi) / 3)
y3 = -temp1 * cos((phi - pi) / 3)
else
temp1 = -q / 2 + sqrt(dd)
temp2 = -q / 2 - sqrt(dd)
u = cbrt(abs(temp1))
v = cbrt(abs(temp2))
if (temp1 < 0) u = -u
if (temp2 < 0) v = -v
y1 = u + v
y2r =-(u + v) / 2
y2i = (u - v) * sqrt(3._knd) / 2
end if

temp1 = ba / 3

if (dd < 0) then
x(1) = y1 - temp1
x(2) = y2 - temp1
x(3) = y3 - temp1
else if (dd > 0) then
x(1) = y1 - temp1
x(2)=cmplx(y2r - temp1, y2i, knd)
x(3)=cmplx(y2r - temp1,-y2i, knd)
else
x(1) = y1 - temp1
x(2) = y2r - temp1
x(3) = y2r - temp1
end if
end subroutine cubic

end program cubicfunction
``````

i will try and thank you

If you calculate pi in double precision you should store it and all your other real and complex variables that way, e.g. with

integer,parameter:: dp = kind(1d0)

and then turn all the real and complex declarations into real(dp) and complex(dp).

1 Like