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"
 read(*,*) a, b,c, d 

  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
Zrzut ekranu 2022-05-07 202621

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"
   read(*,*) a, b, c, d 

   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