'============================================================================
'   G.H. George                                       Created:   1992 SEP 28
'                               BIN-EXPN.BAS         Modified:   1998 MAY 01
'   to find the first few terms of the binomial expansion of  (a + bx)^n .
'============================================================================

'   Procedures:
'       Header:   print welcome message
'     ValuesIn:   obtain values of coefficients from the user
'         Guts:   calculate the binomial expansion
'   ResultsOut:   print the binomial expansion
DECLARE SUB Header ()
DECLARE SUB ValuesIn (a AS LONG, b AS LONG, Num AS LONG, Den AS LONG)
DECLARE SUB Guts (a AS LONG, b AS LONG, Num AS LONG, Den AS LONG, n() AS LONG, d() AS LONG)
DECLARE FUNCTION NextPrime& (Prime AS LONG)
DECLARE SUB ReduceFraction (Num&, Den&)
DECLARE SUB ClearArray (Array() AS LONG)
DECLARE SUB ResultsOut (a AS LONG, b AS LONG, Num AS LONG, Den AS LONG, n() AS LONG, d() AS LONG)

'   Global Constants:
'       False:  logical false (= 0)
'        True:  logical true (= -1)
DEFINT F, T
CONST False = 0, True = NOT False

'   Variables:
'       a, b:   integer coefficients in  (a + bx)^n
'   Num, Den:   numerator and denominator of the exponent  n  in (a + bx)^n
'        n():   array of numerators of binomial coefficients
'        d():   array of denominators of binomial coefficients
DIM a AS LONG, b AS LONG, Num AS LONG, Den AS LONG
DIM n(1 TO 20) AS LONG, d(1 TO 20) AS LONG

CALL Header
CALL ValuesIn(a, b, Num, Den)
CALL Guts(a, b, Num, Den, n(), d())
CALL ResultsOut(a, b, Num, Den, n(), d())

END

DEFSNG F, T
SUB ClearArray (Array() AS LONG)
'============================================================================
'   Set all entries in the long-integer array parameter  Array()  to zero.
'============================================================================

'   Parameter: [as above]
'   No procedures.

'   Local variable:
'       i:  loop counter = array subscript

DIM i AS INTEGER

FOR i = LBOUND(Array) TO UBOUND(Array)   ' all valid subscripts for Array().
    LET Array(i) = 0&
NEXT i

END SUB

SUB Guts (a AS LONG, b AS LONG, Num AS LONG, Den AS LONG, n() AS LONG, d() AS LONG)
'============================================================================
'  Find the first few terms in the binomial expansion (up to x^6).
'============================================================================

'   Procedures:
'       ClearArray:  set all elements of an array to zero.
'   ReduceFraction:  reduce a fraction to its lowest terms.

'   Parameters:
'       a, b:   integer coefficients in  (a + bx)^n
'   Num, Den:   numerator and denominator of the exponent  n  in (a + bx)^n
'        n():   array of numerators of binomial coefficients
'        d():   array of denominators of binomial coefficients

'   Local variables:
'       Up:   [leads to] numerator of the current term in the expansion
'     Down:   [leads to] denominator of the current term in the expansion
'     Term:   loop counter = index number of the current term
'    NewUp:   numerator of factor  (n + 1 - Term)  in the current term
DIM Term AS INTEGER, NewUp AS INTEGER
DIM Up AS LONG, Down AS LONG

'                  {                           2                    3       }
'           n    n {       n bx     n(n-1) (bx)     n(n-1)(n-2) (bx)        }
'   (a + bx)  = a  { 1  +  - --  +  ------ (--)  +  ----------- (--)  +  ...}
'                  {       1  a       2x1  ( a)        3x2x1    ( a)        }

'   First (x^1) term:   ---^^^^

CALL ClearArray(n())
CALL ClearArray(d())
LET Up = Num * b
LET Down = Den * a
CALL ReduceFraction(Up, Down)
LET n(1) = Up
LET d(1) = Down

'   Second to sixth terms:

LET Term = 2
LET NewUp = Num - Den

DO WHILE NewUp <> 0 AND Term < 7
    LET Up = Up * NewUp * b
    LET Down = Down * Term * Den * a
    CALL ReduceFraction(Up, Down)
    LET n(Term) = Up
    LET d(Term) = Down
    LET Term = Term + 1
    LET NewUp = NewUp - Den
LOOP

END SUB

SUB Header
'============================================================================
'  Print instructions to the user.
'============================================================================

'   No procedures, parameters or local variables.

CLS
PRINT "This program will evaluate, with rational coefficients, the first few"
PRINT "terms in the binomial expansion of"
PRINT
PRINT "                                           n"
PRINT "                           f(x) = (a + b x) "
PRINT
PRINT "where the parameters  a  and  b  are integers and the index  n  is"
PRINT "entered as a rational number.": PRINT
PRINT "Note:   fractions are reduced to their lowest terms only if there is"
PRINT "no prime factor larger than 1000 common to numerator and denominator."
PRINT

END SUB

FUNCTION NextPrime& (Number AS LONG)
'=============================================================================
'  **  User defined FUNCTION  **     (c) 1991, 1998  **>>  G.H.  George  <<**
'                    - given input "Number", find the next prime number.
'  Copied and modified from LINSYS.BAS .
'=============================================================================

'   Procedures:
'       <none>

'   Parameters:
'     Number:  incoming number, the next prime above which is to be found

'   Local Variables:
'      Found:  logical variable (= true when nature of candidate is found)
'       Cand:  candidate next prime number (initially = Number)
'     Factor:  factor by which to try to divide the candidate prime number
' LastFactor:  last factor to try for this candidate prime number

DIM Cand  AS LONG, Factor AS INTEGER, Found AS LONG

LET Cand = Number

DO
   LET Cand = Cand + 1
   LET LastFactor = SQR(Cand)
   LET Factor = 1
   LET Found = False

   DO
      LET Factor = Factor + 1
      IF Factor > LastFactor THEN LET Found = True      ' <-- Prime detected.
      IF Cand MOD Factor = 0 THEN LET Found = True  ' <-- "Cand" is composite.
   LOOP UNTIL Found

LOOP UNTIL Factor > LastFactor                 ' Exit only when prime found.

LET NextPrime& = Cand

END FUNCTION

SUB ReduceFraction (Num AS LONG, Den AS LONG)
'=============================================================================
' ** SUBprogram **                   (c) 1991, 1998  **>>  G.H.  George  <<**
' cancels the fraction (Num/Den) down to lowest terms.
' (copied and modified from LINSYS.BAS).
'                         {modified not to seek common prime factors > 1000 }
'=============================================================================

'   Procedures:
'  NextPrime:  returns the next prime number after the argument

'   Parameters:
'   Num, Den:  numerator and denominator of the fraction to be reduced

'   Local Variables:
'    Proceed:  logical variable (= true if further reduction may be needed)
'      Prime:  prime number; used to test Num & Den for common factors

DIM Proceed AS INTEGER, Prime AS LONG

LET Proceed = True                                ' no special cases.

IF Den = 0 THEN
   LET Proceed = False                    ' Avoid zero denominator entirely.
ELSE
   IF Num = 0 THEN                        ' Check for zero numerator.
      LET Den = 1
      LET Proceed = False                 '  and no need to reduce further.
   END IF

   IF Den < 0 THEN                        ' Make denominator positive.
      LET Num = -Num                      ' - but further reduction may
      LET Den = -Den                      '   be needed.
   END IF

   IF ABS(Num) = 1 OR Den = 1 THEN
      LET Proceed = False                 ' No reduction needed for unit
   END IF                                 '  numerator or denominator.

END IF

IF Proceed THEN                           ' Proceed with main routine ONLY
   LET Prime = 2                         '  when no special cases detected.

   DO
   
      DO
         IF Num / Prime = INT(Num / Prime) THEN
            LET Proceed = False           ' Test divisibility of numerator.
         ELSE
            LET Proceed = True
         END IF
'
         IF Den / Prime <> INT(Den / Prime) THEN
            LET Proceed = True            ' Test divisibility of denominator.
         END IF
'
         IF NOT Proceed THEN              ' Divide out by common factor found.
            LET Num = Num / Prime
            LET Den = Den / Prime
         END IF
    
      LOOP UNTIL Proceed                  ' Try same prime again if it's a CF.
   
      LET Prime = NextPrime&(Prime)          '  else go to next prime number.

   LOOP UNTIL Prime > ABS(Num) OR Prime > Den OR Prime > 1000
                                  ' until lowest terms (or prime too big).

END IF

END SUB

SUB ResultsOut (a AS LONG, b AS LONG, Num AS LONG, Den AS LONG, n() AS LONG, d() AS LONG)
'============================================================================
'  Print results out.      - note: temporarily disabled after  x^6 .
'============================================================================

'   No procedures.

'   Parameters:
'       a, b:   integer coefficients in  (a + bx)^n
'   Num, Den:   numerator and denominator of the exponent  n  in (a + bx)^n
'        n():   array of numerators of binomial coefficients
'        d():   array of denominators of binomial coefficients

'   Local variables:
'   u$, c$, l$:  format strings for upper, centre and lower printed lines
'          Row:  loop counter = row of output
'         Term:  keeps track of terms (= (first power of x) on current row)
'            i:  inner loop counter
DIM Row AS INTEGER, Term AS INTEGER, i AS INTEGER

'   First line:  original expression
CLS
PRINT "The first few terms of the binomial expansion of"
PRINT TAB(24); "("; a; " + "; b; "x) ^ ("; Num; "/ "; Den; ")"
PRINT "are:"
PRINT

'   Second line: binomial expansion up to the x^2 term
LET u$ = "       (#####/#####) {     #######     #######  2"
LET c$ = "#######              { 1 + ------- x + ------- x + "
LET l$ = "                     {     #######     #######   "

PRINT USING u$; Num; Den; n(1); n(2)
PRINT USING c$; a
PRINT USING l$; d(1); d(2)
PRINT

'   All remaining lines (four terms at a time)
LET u$ = "  ##########  ##"
LET c$ = "+ ---------- x  "
LET l$ = "  ##########    "

FOR Row = 1 TO 1
    LET Term = 4 * Row - 1
   
    FOR i = Term TO Term + 3
        PRINT USING u$; n(i); i;
    NEXT i
   
    PRINT "   }"                    ' <-- <strip this when  x^7  restored>
   
    FOR i = Term TO Term + 3
        PRINT c$;
    NEXT i
   
    PRINT "   }"                    ' <-- <strip this when  x^7  restored>
   
    FOR i = Term TO Term + 3
        PRINT USING l$; d(i);
    NEXT i
   
    PRINT "   }"                    ' <-- <strip this when  x^7  restored>
    PRINT
NEXT Row

'FOR i = 19 TO 20                   ' <-- <for terms beyond  x^6 >
'    PRINT USING u$; n(i); i;
'NEXT i
'PRINT
'PRINT c$; c$; "}"
'FOR i = 19 TO 20
'    PRINT USING l$; d(i);
'NEXT i

END SUB

SUB ValuesIn (a AS LONG, b AS LONG, Num AS LONG, Den AS LONG)
'============================================================================
'   Ask user for parameter values in   (a + bx)^(Num/Den) .
'============================================================================

'   Procedures:
'   ReduceFraction:  reduce a fraction to its lowest terms.

'   Parameters:
'       a, b:   integer coefficients in  (a + bx)^n
'   Num, Den:   numerator and denominator of the exponent  n  in (a + bx)^n

'   Local variable:
'       Temp:   single precision temporary variable for various user inputs

DIM Temp AS SINGLE
                                ' -- CHR$(30)  forces cursor UP one line.
' a  must be a non-zero integer, < 10,000,000
DO
    PRINT CHR$(30); SPACE$(70)
    PRINT CHR$(30);
    INPUT "Enter a non-zero INTEGER value for  a :", Temp
LOOP UNTIL INT(Temp) = Temp AND Temp <> 0 AND ABS(Temp) < 1E+07

LET a = CLNG(Temp)
PRINT

' b  must be an integer, < 10,000,000
DO
    PRINT CHR$(30); SPACE$(70)
    PRINT CHR$(30);
    INPUT "Enter an INTEGER value for  b :", Temp
LOOP UNTIL INT(Temp) = Temp AND ABS(Temp) < 1E+07

LET b = CLNG(Temp)
PRINT

' Num  must be an integer, < 10,000,000
DO
    PRINT CHR$(30); SPACE$(70)
    PRINT CHR$(30);
    INPUT "Enter the numerator of the exponent  n (INTEGER): ", Temp
LOOP UNTIL Temp = INT(Temp) AND ABS(Temp) < 1E+07

LET Num = CLNG(Temp)
PRINT

' Den  must be a non-zero integer, < 10,000,000
DO
    PRINT CHR$(30); SPACE$(70)
    PRINT CHR$(30);
    INPUT "Enter the denominator of  n (POSITIVE INTEGER): ", Temp
LOOP UNTIL Temp = INT(Temp) AND Temp <> 0 AND ABS(Temp) < 1E+07

LET Den = CLNG(Temp)
PRINT

CALL ReduceFraction(Num, Den)

END SUB