home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / MAGAZINE / MISC / QBNWS203.ZIP / CUBERT.ZIP / CUBERT.BAS
Encoding:
BASIC Source File  |  1991-01-04  |  3.3 KB  |  124 lines
  1. DECLARE SUB CUBERT (C#(), RT#(), FLAG%)
  2. ' Cubic equation solve routine.
  3. ' Illustrates usage of subroutine CUBERT
  4. start:
  5. CLS 0
  6. LOCATE 10
  7. DIM A#(3), X#(3)
  8.  
  9. INPUT "Enter coefficient of cubic term"; A#(3)
  10. INPUT "Enter coefficient of quadratic term"; A#(2)
  11. INPUT "Enter coefficient of linear term"; A#(1)
  12. INPUT "Enter constant term"; A#(0)
  13. PRINT "The roots of the equation - "
  14. CALL CUBERT(A#(), X#(), F%)
  15. IF F% = 2 THEN          ' 1 real and conjugate pair
  16.    X1# = R# - AAS# / 3
  17.    PRINT "x1 = "; X#(1)
  18.    PRINT "x2 = "; X#(2); " + "; X#(3); " i "
  19.    PRINT "x3 = "; X#(2); " - "; X#(3); " i "
  20. ELSE              ' 3 real and distinct roots
  21.    PRINT "x1 = "; X#(1)
  22.    PRINT "x2 = "; X#(2)
  23.    PRINT "x3 = "; X#(3)
  24. END IF
  25. P# = A#(3)
  26. FOR I% = 1 TO 3
  27.   P# = P# * X#(1) + A#(3 - I%)
  28. NEXT I%
  29. PRINT
  30. ' Display value of equation at root as check
  31. PRINT "Value of equation at root1 =: "; P#
  32. END
  33.  
  34. SUB CUBERT (C#(), RT#(), FLAG%) STATIC
  35.  
  36. ' Double Precision Cubic equation solve routine.
  37. ' Solves by Cardan's method.
  38. ' Written by Richard G. Jones 3/7/88
  39. '
  40. ' Call with:
  41. '           C#() = 4 coefficients in elements 0 - 3
  42. '           RT#() = array with at least 4 elements.
  43. ' Returns:
  44. '           FLAG% = status; one of the following:
  45. '
  46. '           = 1:   3 real and distinct roots in RT#(1) - RT#(3) (r1 - r3)
  47. '           = 2:   roots are given as follows
  48. '                  x1 = RT#(1)
  49. '                  x2 = RT#(2) + RT#(3) * i
  50. '                  x3 = RT#(2) - RT#(3) * i   where i = SquareRoot(-1)
  51. '
  52. '          note: Possible root multiplicities in case 2 could occur.
  53. '          if r3 = 0 and r1 <> r2 there is a double root.
  54. '          if r3 = 0 and r1 = r2  there is a triple root
  55.  
  56. CONST PI# = 3.141592653589793#, BIG# = 1D+40
  57.  
  58. ' assign coef's to locals to eliminate array indexing
  59.  
  60. A3# = C#(3)
  61.  
  62. ' calculate coefficients of Cardan's reduced cubic (p & q):
  63. '  u^3 + p * u + q = 0
  64.  
  65. AAS# = C#(2) / A3#: BS# = C#(1) / A3#
  66. CS# = C#(0) / A3#
  67. P# = BS# - AAS# ^ 2# / 3#
  68. Q# = 2# * (AAS# / 3#) ^ 3# - BS# * AAS# / 3# + CS#
  69.  
  70. ' calculate discriminant D.
  71. ' If D > 0 there are 1 real and a pair of complex conjugate roots.
  72. '    D = 0 all roots are real and at least 2 are equal.
  73. '    D < 0 all roots are real and distinct.
  74.  
  75.  
  76. D# = (P# / 3#) ^ 3# + (Q# / 2#) ^ 2#
  77. U# = -Q# / 2#
  78.  
  79. IF D# >= 0# THEN
  80. '  case of 1 real root and a pair of complex conjugates or multiple roots
  81.  
  82.    A33# = (U# + SQR(D#))
  83.    B33# = (U# - SQR(D#))
  84.    IF A33# < 0# THEN
  85.      A# = -(ABS(A33#) ^ (1# / 3#))
  86.    ELSE
  87.      A# = A33# ^ (1# / 3#)
  88.    END IF
  89.    IF B33# < 0# THEN
  90.      B# = -(ABS(B33#) ^ (1# / 3#))
  91.    ELSE
  92.      B# = B33# ^ (1# / 3#)
  93.    END IF
  94.    R# = A# + B#
  95.    S# = (A# - B#) * 3# ^ .5#
  96.    RT#(1) = R# - AAS# / 3           ' real root
  97.    RT#(2) = -(R# / 2# + AAS# / 3#)  ' real component of conjugates
  98.    RT#(3) = ABS(S#)                 ' imaginary component
  99.    FLAG% = 2
  100.  
  101. ELSE
  102. '  case of 3 real and distinct roots.
  103.  
  104.    MAG# = SQR(U# ^ 2# - D#)
  105.    DG# = SQR(-D#)
  106.    IF ABS(DG#) < ABS(U#) * BIG# THEN
  107.       THETA# = ATN(ABS(DG# / U#))
  108.       IF U# < 0 THEN THETA# = PI# - THETA#
  109.    ELSE
  110.       THETA# = PI# / 2#
  111.    END IF
  112.  
  113.    R# = MAG# ^ (1# / 3#) * COS(THETA# / 3#)
  114.    S# = MAG# ^ (1# / 3#) * SIN(THETA# / 3#) * 3# ^ .5#
  115.    C# = AAS# / 3#
  116.    RT#(1) = 2# * R# - C#      ' return 3 roots
  117.    RT#(2) = -R# + S# - C#
  118.    RT#(3) = -R# - S# - C#
  119.    FLAG% = 1
  120.  
  121. END IF
  122.  
  123. END SUB
  124.