'*****************************************************************
'*   Program to demonstrate Bisection & Quartile subroutines     *
'* ------------------------------------------------------------- *
'* References:                                                   *
'*                                                               *
'*  Bisection: BASIC Scientific Subroutines, Vol. II             *
'*             By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
'*                                                               *
'*  Quartile:  After an algorithm provided By Namir Shammas.     *
'*                                                               *
'*                           Basic Release By J-P Moreau, Paris. *
'*                                   (www.jpmoreau.fr)           *
'* ------------------------------------------------------------- *
'* Example:  Find a real root of f(x)=(x+1)^5                    *
'*                                                               *
'* Sample run:                                                   *
'*                                                               *
'* What is the initial range (X0,X1):                            *
'*    X0 = -5                                                    *
'*    X1 = 0                                                     *
'* Convergence criterion: 1e-6                                   *
'*                                                               *
'*                                                               *
'* Bisection: The calculated zero is X =  -.9999996423721313     *
'*                                                               *
'* The associated Y value is Y = 5.850012206076229D-33           *
'*                                                               *
'* The number of steps was: 23                                   *
'*                                                               *
'*                                                               *
'* Quartile: The calculated zero is X =  -1.000000089406967      *
'*                                                               *
'* The associated Y value is Y = -5.712902544996317D-36          *
'*                                                               *
'* The number of steps was: 12                                   *
'*                                                               *
'*****************************************************************
'Note: The quartile method is slightly faster and more accurate
'      than the classical Bisection Method.
'-----------------------------------------------------------------
'PROGRAM TQUART
DEFDBL A-H, O-Z
DEFINT I-N

' real*8  a,b,root, e,x,x0,x1
' integer m

  CLS
  PRINT
  PRINT " What is the initial range (X0,X1):"
  INPUT "   X0 = ", x0
  INPUT "   X1 = ", x1
  INPUT " Convergence criterion: ", e
  PRINT

  a = x0: b = x1

  GOSUB 1000  'call Bisect(x0, x1, e, x, m)

  IF m <> 0 THEN
    PRINT
    PRINT " Bisection: The calculated zero is X = "; x
    PRINT
    GOSUB 500 'calculate Y(x)
    PRINT " The associated Y value is Y = "; Y
    PRINT
    PRINT " The number of steps was: "; m
  END IF

  GOSUB 2000  'call Quartile(a, b, e, root, istep)

  IF m <> 0 THEN
    PRINT
    PRINT " Quartile: The calculated zero is X = "; root
    PRINT
    x = root: GOSUB 500
    PRINT " The associated Y value is Y = "; Y
    PRINT
    PRINT " The number of steps was: "; istep
  END IF
  PRINT

END

'***********************************************
500 'function Y(x)
  Y = (x + 1) ^ 5
RETURN
'***********************************************


'***********************************************
'*        Bisection method subroutine          *
'* ------------------------------------------- *
'* This routine iteratively seeks the zero of  *
'* function Y(x) using the method of interval  *
'* halving until the interval is less than e   *
'* in width. It is assumed that the function   *
'* Y(x) is available from a function routine.  *
'* ------------------------------------------- *
'* Input values: range (x0,x1), and e.         *
'* Output values: root x, Y(x) and number of   *
'* steps m.                                    *
'***********************************************
1000 'subroutine Bisect(x0,x1,e,x,m)
'Label: 10
  m = 0
  x = x0: GOSUB 500: Y1 = Y
  x = x1: GOSUB 500
  IF Y1 * Y > 0# THEN
    PRINT " Bisection: No root found in given interval."
    RETURN
  ELSE
10  x = x0: GOSUB 500: y0 = Y
    x = (x0 + x1) / 2: GOSUB 500: yy = Y
    m = m + 1
    IF yy * y0 = 0# THEN RETURN
    IF yy * y0 < 0# THEN x1 = x
    IF yy * y0 > 0# THEN x0 = x
    IF ABS(x1 - x0) > e THEN GOTO 10
  END IF
RETURN

'***********************************************
'         Quartile method subroutine           *
' -------------------------------------------- *
' This routine iteratively seeks the zero of   *
' function Y(x) using the method of interval   *
' quarting until the interval is less than tol *
' in width. It is assumed that the function    *
' Y(x) is available from a function routine.   *
' -------------------------------------------- *
' Input values: range (a, b), and tol.         *
' Output values: found root and number of used *
' steps.                                       *
'***********************************************
2000 'subroutine Quartile(a, b, tol, root, istep)
  
' integer istep
' double xm, co

  tol = e
  co = .25# 'can be 0.3
  istep = 0
  x = a: GOSUB 500: Y1 = Y
  x = b: GOSUB 500
  IF Y1 * Y > 0# THEN
    PRINT " Quartile: No root found in given interval."
    RETURN
  ELSE
20  x = a: GOSUB 500: Y1 = Y
    x = b: GOSUB 500
    IF ABS(Y1) < ABS(Y) THEN
      xm = a + co * (b - a)
    ELSE
      xm = a + (1# - co) * (b - a)
    END IF

    x = xm: GOSUB 500: Y1 = Y
    x = a: GOSUB 500
    IF Y1 * Y > 0# THEN
      a = xm
    ELSE
      b = xm
    END IF
    istep = istep + 1
    IF ABS(a - b) > tol THEN GOTO 20
    root = (a + b) / 2#
  END IF
RETURN

'End of file tquart.bas