1SECTION : 1.4                (OCTOBER 1974)                      PAGE 1
  
  
  
 AUTHOR: H.J.J. TE RIELE.
  
  
 INSTITUTE: MATHEMATICAL CENTRE.
  
  
 RECEIVED: 740125; REVISED: 740514;
  
  
 BRIEF DESCRIPTION:
  
     THIS SECTION CONTAINS A SET OF FIVE PROCEDURES FOR
     THE BASIC ARITHMETIC OPERATIONS WITH LONG INTEGERS:
     LNG INT ADD EXACTLY COMPUTES THE SUM OF TWO NONNEGATIVE INTEGERS.
     LNG INT SUBTRACT EXACTLY COMPUTES THE DIFFERENCE OF TWO NONNEGATIVE
     INTEGERS.
     LNG INT MULT EXACTLY COMPUTES THE PRODUCT OF TWO NONNEGATIVE
     INTEGERS.
     LNG INT DIVIDE EXACTLY COMPUTES THE QUOTIENT WITH REMAINDER OF TWO
     NONNEGATIVE INTEGERS.
     LNG INT POWER EXACTLY COMPUTES U**POWER,
     WHERE U IS A NONNEGATIVE LONG INTEGER AND POWER IS THE
     POSITIVE (SINGLE-LENGTH) EXPONENT.
  
  
 KEYWORDS:
  
     LONG INTEGER ARITHMETIC,
     ADDITION,
     SUBTRACTION,
     MULTIPLICATION,
     DIVISION WITH REMAINDER,
     EXPONENTIATION.
  
  
1SECTION : 1.4                (OCTOBER 1974)                      PAGE 2
  
  
  
 SUBSECTION : LNG INT ADD.
  
  
 CALLING SEQUENCE :
  
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" LNG INT ADD(U,V,SUM);
     "INTEGER" "ARRAY" U,V,SUM;
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     U,V,SUM: <ARRAY IDENTIFIER>;
             "INTEGER" "ARRAY" U[0:U[0]], V[0:V[0]],
                               SUM[0:MAX(U[0],V[0])+1];
             BEFORE THE CALL OF LNG INT ADD, U AND V MUST
             CONTAIN THE LONG INTEGERS TO BE ADDED;
             AFTER THE CALL, SUM CONTAINS THE MULTI-LENGTH
             SUM OF U AND V, WHILE U AND V REMAIN UNCHANGED.
  
  
 PROCEDURES USED : NONE.
  
  
 REQUIRED CENTRAL MEMORY :
  
     EXECUTION FIELD LENGTH : 7.
  
  
 RUNNING TIME :
  
     WE GIVE A FORMULA FOR THE RUNNING TIME  IN MILLISECONDS ON THE
     CD CYBER 73-28 COMPUTER; THE RELATIVE PRECISION OF THE
     COEFFICIENTS IS AT MOST ONE OR TWO DIGITS:
     .10*MAX(U[0],V[0]) + .06*MIN(U[0],V[0]) + .56.
  
  
 LANGUAGE : ALGOL 60.
  
  
 METHOD AND PERFORMANCE : SEE LNG INT POWER (THIS SECTION).
  
  
 EXAMPLE OF USE : SEE LNG INT POWER (THIS SECTION).
  
  
1SECTION : 1.4                (OCTOBER 1974)                      PAGE 3
  
  
  
 SUBSECTION : LNG INT SUBTRACT.
  
  
 CALLING SEQUENCE :
  
     THE HEADING OF THE PROCEDURE READS :
     "PROCEDURE" LNG INT SUBTRACT (U,V,DIFFERENCE);
     "INTEGER" "ARRAY" U,V,DIFFERENCE;
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     U,V,DIFFERENCE: <ARRAY IDENTIFIER>;
             "INTEGER" "ARRAY" U[0:U[0]],V[0:V[0]],DIFFERENCE[0:U[0]];
             BEFORE THE CALL OF LNG INT SUBTRACT, U AND V MUST
             CONTAIN THE LONG INTEGERS TO BE SUBTRACTED(V FROM U);
             AFTER THE CALL, DIFFERENCE CONTAINS THE MULTI-LENGTH
             DIFFERENCE U-V; IF U<V THEN DIFFERENCE[0]=0
             IS DELIVERED; U AND V REMAIN UNCHANGED.
  
  
 PROCEDURES USED : NONE.
  
  
 REQUIRED CENTRAL MEMORY:
  
     EXECUTION FIELD LENGTH : 7.
  
  
 RUNNING TIME :
  
     WE GIVE A FORMULA FOR THE RUNNING TIME  IN MILLISECONDS ON THE
     CD CYBER 73-28 COMPUTER; THE RELATIVE PRECISION OF THE
     COEFFICIENTS IS AT MOST ONE OR TWO DIGITS:
     .10*U[0] + .06*V[0] + .64.
  
  
 LANGUAGE : ALGOL 60.
  
  
 METHOD AND PERFORMANCE : SEE LNG INT POWER (THIS SECTION).
  
  
 EXAMPLE OF USE : SEE LNG INT POWER (THIS SECTION).
  
  
1SECTION : 1.4                (OCTOBER 1974)                      PAGE 4
  
  
  
 SUBSECTION : LNG INT MULT.
  
  
 CALLING SEQUENCE :
  
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" LNG INT MULT(U,V,PRODUCT);
     "INTEGER" "ARRAY" U,V,PRODUCT;
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     U,V,PRODUCT: <ARRAY IDENTIFIER>;
             "INTEGER" "ARRAY" U[0:U[0]], V[0:V[0]],
                               PRODUCT[0:U[0]+V[0]];
             BEFORE THE CALL OF LNG INT MULT, U AND V MUST
             CONTAIN THE LONG INTEGERS TO BE MULTIPLIED;
             AFTER THE CALL, PRODUCT CONTAINS THE MULTI-LENGTH
             PRODUCT OF U AND V, WHILE U AND V REMAIN UNCHANGED.
  
  
 PROCEDURES USED : NONE.
  
  
 REQUIRED CENTRAL MEMORY  :
  
     EXECUTION FIELD LENGTH : 7.
  
  
 RUNNING TIME :
  
     WE GIVE A FORMULA FOR THE RUNNING TIME  IN MILLISECONDS ON THE
     CD CYBER 73-28 COMPUTER; THE RELATIVE PRECISION OF THE
     COEFFICIENTS IS AT MOST ONE OR TWO DIGITS:
     .18*U[0]*V[0] + .15*U[0] + .06*V[0] + .46.
  
  
 LANGUAGE : ALGOL 60.
  
  
 METHOD AND PERFORMANCE : SEE LNG INT POWER (THIS SECTION).
  
  
 EXAMPLE OF USE : SEE LNG INT POWER (THIS SECTION).
  
  
1SECTION : 1.4                (MARCH 1977)                        PAGE 5
  
  
  
 SUBSECTION : LNG INT DIVIDE.
  
  
 CALLING SEQUENCE :
  
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" LNG INT DIVIDE(U,V,QUOTIENT,REMAINDER); "VALUE" U;
     "INTEGER" "ARRAY" U,V,QUOTIENT,REMAINDER;
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     U,V,QUOTIENT,REMAINDER: <ARRAY IDENTIFIER>;
             "INTEGER" "ARRAY" U[0:U[0]], V[0:V[0]],
             QUOTIENT[0:U[0]-V[0]+1], REMAINDER[0:V[0]];
             BEFORE THE CALL OF LNG INT DIVIDE, U MUST CONTAIN THE
             DIVIDEND, V THE DIVISOR(V ^= 0);
             AFTER THE CALL, THE RESULTS OF THE LONG DIVISION
             OF U BY V (I.E. U//V AND U-U//V) ARE STORED INTO
             QUOTIENT AND REMAINDER; U AND V REMAIN UNCHANGED.
  
  
 PROCEDURES USED : NONE.
  
  
 REQUIRED CENTRAL MEMORY :
  
     11 + U[0] + (IF V[0]=1 OR U[0]<V[0] THEN 0 ELSE V[0]+1),
  
  
 RUNNING TIME :
  
     WE GIVE A FORMULA FOR THE RUNNING TIME  IN MILLISECONDS ON THE
     CD CYBER 73-28 COMPUTER; THE RELATIVE PRECISION OF THE
     COEFFICIENTS IS AT MOST ONE OR TWO DIGITS:
     IF V[0]=1 THEN (.34*U[0] + .67) ELSE IF V[1] >= 5 000 000 THEN
                         (.26*DIFF*V[0] + .57*DIFF + .10*V[0] + 1.8)
               ELSE  (.27*DIFF*V[0] + .66*DIFF + .66*V[0] + 2.0)
     (HERE DIFF=U[0]-V[0]+1, I.E. THE NUMBER OF EXECUTIONS
     OF THE STATEMENT, IN WHICH DIVISION OF A (V[0]+1)-PLACE
     NUMBER BY A V[0]-PLACE NUMBER IS PERFORMED).
  
  
 LANGUAGE : ALGOL 60.
  
  
 METHOD AND PERFORMANCE : SEE LNG INT POWER (THIS SECTION).
  
  
 EXAMPLE OF USE : SEE LNG INT POWER (THIS SECTION).
  
  
1SECTION : 1.4                (OCTOBER 1974)                      PAGE 6
  
  
  
 SUBSECTION : LNG INT POWER.
  
  
 CALLING SEQUENCE :
  
     THE HEADING OF THE PROCEDURE READS:
     "PROCEDURE" LNG INT POWER(U,EXPONENT,RESULT);
     "VALUE" EXPONENT; "INTEGER" EXPONENT; "INTEGER" "ARRAY" U,RESULT;
  
     THE MEANING OF THE FORMAL PARAMETERS IS:
     EXPONENT: <ARITHMETIC EXPRESSION>;
             THE (POSITIVE) POWER TO WHICH THE LONG INTEGER U
             WILL BE RAISED;
     U,RESULT: <ARRAY IDENTIFIER>;
             "INTEGER" "ARRAY" U[0:U[0]], RESULT[0:U[0]*EXPONENT];
             BEFORE THE CALL OF LNG INT POWER, U MUST CONTAIN
             THE LONG INTEGER WHICH HAS TO BE RAISED TO THE
             POWER EXPONENT;
             AFTER THE CALL, RESULT CONTAINS THE VALUE OF THE
             LONG INTEGER U**EXPONENT; U REMAINS UNCHANGED.
  
  
 PROCEDURES USED :
  
     LNG INT MULT = CP31202.
  
  
 REQUIRED CENTRAL MEMORY :
  
     EXECUTION FIELD LENGTH : 4 + 3 * (U[0] * EXPONENT + 1).
  
  
 RUNNING TIME :
  
     FOR THIS PROCEDURE THE TIME FORMULA IS A COMPLICATED FUNCTION OF
     U[0], EXPONENT AND  THE NUMBER OF ONES IN THE BINARY REPRESENTATION
     OF EXPONENT, BUT ROUGHLY THE TIME IS  OF THE ORDER :
     (U(0)*EXPONENT)**2).
     TWO TESTCASES :
     EXPONENT    TIME(IN SEC.) FOR:
                   U[0]=1     U[0]=2
         20           .04        .10
         40           .13        .34
        100           .68       1.94
        300          5.48      16.6
        500         16.8       51.0
  
  
1SECTION : 1.4                (OCTOBER 1974)                      PAGE 7
  
  
  
 LANGUAGE : ALGOL 60.
  
  
 METHOD AND PERFORMANCE:
  
     DEFINITION:
     A LONG INTEGER OF LENGTH N, OR AN N-PLACE INTEGER
     (N>0) IS ANY NONNEGATIVE INTEGER LESS THAN BASE**N,
     AND GREATER THAN OR EQUAL TO BASE**(N-1), WHERE BASE
     IS THE (POSITIVE) RADIX OF THE POSITIONAL NOTATION, IN
     WHICH THE INTEGERS ARE EXPRESSED.
  
     ALL FIVE PROCEDURES USE THE BASE 10 000 000; THIS IS THE LARGEST
     POWER OF 10, THE SQUARE OF WHICH CAN BE REPRESENTED EXACTLY
     ON THE CD CYBER 73-28 COMPUTER. IF ONE WANTS TO USE THE
     PROCEDURES WITH ANOTHER VALUE OF THE BASE, SAY R (NOT
     NECESSARILY A POWER OF 10), THEN IN THE SOURCE TEXTS OF THE
     PROCEDURES THE NUMBER 10 000 000 HAS TO BE REPLACED BY R
     (8 TIMES IN LNG INT ADD,
      2 TIMES IN LNG INT SUBTRACT,
      2 TIMES IN LNG INT MULT AND
     16 TIMES IN LNG INT DIVIDE).
     MOREOVER, IN LNG INT DIVIDE THE NUMBER 9 999 999 HAS TO BE
     REPLACED BY THE NUMBER R - 1.
  
     IF A[1], A[2], ..., A[N] ARE THE N "DIGITS" OF THE LONG
     INTEGER M OF LENGTH N (A[1] ^= 0), THEN
  
     M=((...(A[1]*BASE+A[2])*BASE+...+A[N-2])*BASE+A[N-1])*BASE+A[N].
  
     ACCORDINGLY, A LONG INTEGER M OF LENGTH N ALWAYS WILL
     BE STORED INTO A CORRESPONDING "INTEGER" "ARRAY" A,
     THE LENGTH N WILL BE STORED INTO THE ARRAY ELEMENT A[0].
  
     FOR THE METHOD OF THE PROCEDURES LNG INT ADD, LNG INT SUBTRACT,
     LNG INT MULT AND LNG INT DIVIDE, SEE [1; PP.229-248];
     PROCEDURE LNG INT POWER USES THE BINARY METHOD FOR
     EXPONENTIATION (SEE [1; PP.398-401]).
  
  
 REFERENCES:
  
     [1]. DONALD E. KNUTH.
          THE ART OF COMPUTER PROGRAMMING, VOLUME 2/
          SEMINUMERICAL ALGORITHMS.
          ADDISON-WESLEY PUBLISHING COMPANY, 1969.
  
1SECTION : 1.4                (MARCH 1977)                        PAGE 8
  
  
  
 EXAMPLE OF USE:
  
     "BEGIN"
        "PROCEDURE" LNG INT ADD(U,V,SUM); "CODE" 31200;
        "PROCEDURE" LNG INT SUBTRACT(U,V,DIFFERENCE); "CODE" 31201;
        "PROCEDURE" LNG INT MULT(U,V,PRODUCT); "CODE" 31202;
        "PROCEDURE" LNG INT DIVIDE(U,V,QUOTIENT,REMAINDER);"CODE" 31203;
        "PROCEDURE" LNG INT POWER (U,EXPONENT,RESULT); "CODE" 31204;
  
        "PROCEDURE" OUT(A); "INTEGER" "ARRAY" A;
        "BEGIN" "INTEGER" I,L; L:= A[0];
           OUTPUT(61,"("B6ZD,(B7D)")",(A[I],I:= 1:L));
           OUTPUT(61,"("/")")
        "END" OUT;
        "INTEGER" "ARRAY" U,V,R1,R2[0:100];
  
        U[0]:=5; U[1]:=333; U[2]:=U[3]:=U[4]:=U[5]:=7 000 000; OUT(U);
        V[0]:=2; V[1]:=4 444; V[2]:=4 444 444; OUT(V);
  
        LNG INT ADD(U,V,R1); OUT(R1);
        LNG INT SUBTRACT(U,V,R1); OUT(R1);
        LNG INT MULT(U,V,R1); OUT(R1);
        LNG INT DIVIDE(U,V,R1,R2); OUT(R1); OUT(R2);
        LNG INT POWER(V,5,R1); OUT(R1)
     "END"
  
     DELIVERS:
         333 7000000 7000000 7000000 7000000
        4444 4444444
         333 7000000 7000000 7004445 1444444
         333 7000000 7000000 6995556 2555556
     1483111 1114073 9114221 9114221 9111110 8000000
      750825 0001650 0826575
         734 0700700
       17341 5299149 6553709 6327185 8964586 9972395 8069589 6628224
  
  
1SECTION : 1.4                (OCTOBER 1974)                      PAGE 9
  
  
  
 SOURCE TEXT(S):
        
 "CODE" 31200;    
 "PROCEDURE" LNG INT ADD(U,V,SUM); "INTEGER""ARRAY" U,V,SUM;        
 "BEGIN""INTEGER" LU,LV,DIFF,CARRY,I,T,MAX;     
    LU:=U[0]; LV:=V[0];     
    "IF" LU >= LV "THEN"    
    "BEGIN" MAX:=LU; DIFF:=LU - LV + 1; CARRY:=0;         
       "FOR" I:=LU "STEP" -1 "UNTIL" DIFF "DO"  
       "BEGIN" T:=U[I] + V[I-DIFF+1] + CARRY;   
          CARRY:="IF" T<10 000 000 "THEN" 0 "ELSE" 1;     
          SUM[I]:=T - CARRY * 10 000 000        
       "END";     
       "FOR" I:=DIFF - 1 "STEP" -1 "UNTIL" 1 "DO"         
       "BEGIN" T:=U[I] + CARRY;       
          CARRY:="IF" T<10 000 000 "THEN" 0 "ELSE" 1;     
          SUM[I]:=T - CARRY * 10 000 000        
       "END"      
    "END" "ELSE"  
    "BEGIN" MAX:=LV; DIFF:=LV - LU + 1; CARRY:=0;         
       "FOR" I:=LV "STEP" -1 "UNTIL" DIFF "DO"  
       "BEGIN" T:=V[I] + U[I-DIFF+1] + CARRY;   
          CARRY:="IF" T<10 000 000 "THEN" 0 "ELSE" 1;     
          SUM[I]:=T - CARRY * 10 000 000        
       "END";     
       "FOR" I:=DIFF - 1 "STEP" -1 "UNTIL" 1 "DO"         
       "BEGIN" T:=V[I] + CARRY;       
          CARRY:="IF" T<10 000 000 "THEN" 0 "ELSE" 1;     
          SUM[I]:=T - CARRY * 10 000 000        
       "END"      
    "END";        
    "IF" CARRY=1 "THEN"     
    "BEGIN" "FOR" I:=MAX "STEP" -1 "UNTIL" 1 "DO"         
       SUM[I+1]:=SUM[I]; SUM[1]:=1; MAX:=MAX + 1
    "END";        
    SUM[0]:=MAX   
 "END" LNG INT ADD;         
         "EOP"    
        
1SECTION : 1.4                (OCTOBER 1974)                     PAGE 10      
        
        
        
 "CODE" 31201;    
 "PROCEDURE" LNG INT SUBTRACT(U,V,DIFFERENCE);  
 "INTEGER""ARRAY" U,V,DIFFERENCE;     
 "BEGIN""INTEGER" LU,LV,DIFF,I,T,J,CARRY;       
    LU:=U[0]; LV:=V[0];     
    "IF" LU<LV "OR" LU=LV "AND" U[1]<V[1] "THEN" DIFFERENCE[0]:=0 "ELSE"      
    "BEGIN" DIFF:=LU - LV + 1; CARRY:=0;        
       "FOR" I:=LU "STEP" -1 "UNTIL" DIFF "DO"  
       "BEGIN" T:=U[I] - V[I-DIFF+1] + CARRY;   
          CARRY:="IF" T<0 "THEN" -1 "ELSE" 0;   
          DIFFERENCE[I]:=T - CARRY * 10 000 000 
       "END";     
       "FOR" I:=DIFF - 1 "STEP" -1 "UNTIL" 1 "DO"         
       "BEGIN" T:=U[I] + CARRY; CARRY:="IF" T<0 "THEN" -1 "ELSE" 0; 
          DIFFERENCE[I]:=T - CARRY * 10 000 000 
       "END";     
       "IF" CARRY=-1 "THEN" 
       "BEGIN" DIFFERENCE[0]:=0; "GOTO" READY "END";      
       I:=0; J:=LU;         
       "FOR" I:=I+1 "WHILE" DIFFERENCE[I]=0 "AND" J>1 "DO" J:=J-1;  
       DIFFERENCE[0]:=J;    
       "IF" J<LU "THEN"     
       "FOR" I:=1 "STEP" 1 "UNTIL" J "DO"       
       DIFFERENCE[I]:=DIFFERENCE[LU+I-J]        
    "END";        
    READY:        
 "END" LNG INT SUBTRACT;    
         "EOP"    
        
 "CODE" 31202;    
 "PROCEDURE" LNG INT MULT(U,V,PRODUCT);         
 "INTEGER""ARRAY" U,V,PRODUCT;        
 "BEGIN""INTEGER" LU,LV,LUV,I,J,CARRY,T;        
    LU:=U[0]; LV:=V[0]; LUV:=LU + LV; 
    "FOR" I:=LU + 1 "STEP" 1 "UNTIL" LUV "DO" PRODUCT[I]:=0;        
    "FOR" J:=LU "STEP" -1 "UNTIL" 1 "DO"        
    "BEGIN" CARRY:=0;       
       "FOR" I:=LV "STEP" -1 "UNTIL" 1 "DO"     
       "BEGIN" T:=U[J] * V[I] + PRODUCT[J+I] + CARRY;     
          CARRY:=T//10 000 000; PRODUCT[J+I]:=T - CARRY * 10 000 000
       "END"; PRODUCT[J]:=CARRY       
    "END";        
    "IF" PRODUCT[1]=0 "THEN"
    "BEGIN" "FOR" I:=2 "STEP" 1 "UNTIL" LUV "DO"
       PRODUCT[I-1]:=PRODUCT[I]; LUV:=LUV - 1   
    "END"; PRODUCT[0]:=LUV  
 "END" LNG INT MULT;        
         "EOP"    
1SECTION : 1.4                (OCTOBER 1974)                     PAGE 11      
        
        
        
 "CODE" 31203;    
 "PROCEDURE" LNG INT DIVIDE(U,V,QUOTIENT,REMAINDER); "VALUE" U;     
 "INTEGER""ARRAY" U,V,QUOTIENT,REMAINDER;       
 "BEGIN""INTEGER" LU,LV,V1,DIFF,I,T,SCALE,D,Q1,J,CARRY;   
    LU:=U[0]; LV:=V[0]; V1:=V[1]; DIFF:=LU - LV;
        
    "IF" LV=1 "THEN"        
    "BEGIN" CARRY:=0;       
       "FOR" I:=1 "STEP" 1 "UNTIL" LU "DO"      
       "BEGIN" T:=CARRY * 10 000 000 + U[I]; QUOTIENT[I]:=T//V1;    
          CARRY:=T - QUOTIENT[I] * V1 
       "END"; REMAINDER[0]:=1; REMAINDER[1]:=CARRY;       
       "IF" QUOTIENT[1]=0 "THEN"      
       "BEGIN" "FOR" I:=2 "STEP" 1 "UNTIL" LU "DO"        
          QUOTIENT[I-1]:=QUOTIENT[I]; 
          QUOTIENT[0]:=LU - ("IF" LU=1 "THEN" 0 "ELSE" 1) 
       "END" "ELSE" QUOTIENT[0]:=LU   
    "END" LV=1    
    "ELSE"        
        
    "IF" LU<LV "THEN"       
    "BEGIN" QUOTIENT[0]:=1; QUOTIENT[1]:=0;     
       "FOR" I:=0 "STEP" 1 "UNTIL" LU "DO" REMAINDER[I]:=U[I]       
    "END" LU<LV   
    "ELSE"        
        
    "BEGIN""INTEGER""ARRAY" A[0:LV];  
       SCALE:=10 000 000//(V1+1);     
       "IF" SCALE>1 "THEN"  
       "BEGIN""COMMENT" NORMALIZE U; CARRY:=0;  
          "FOR" I:=LU "STEP" -1 "UNTIL" 1 "DO"  
          "BEGIN" T:=SCALE * U[I] + CARRY; CARRY:=T//10 000 000;    
             U[I]:=T - CARRY * 10 000 000       
          "END"; U[0]:=CARRY;         
          "COMMENT" NORMALIZE V; CARRY:=0;      
          "FOR" I:=LV "STEP" -1 "UNTIL" 1 "DO"  
          "BEGIN" T:=SCALE * V[I] + CARRY; CARRY:=T//10 000 000;    
             V[I]:=T - CARRY * 10 000 000       
          "END"; V1:=V[1]   
       "END" NORMALIZATION  
       "ELSE" U[0]:=0;                                         "COMMENT"      
        
1SECTION : 1.4                (OCTOBER 1974)                     PAGE 12      
                                                                 ;  
        
        
       "COMMENT" COMPUTE QUOTIENT AND REMAINDER;
       "FOR" I:=0 "STEP" 1 "UNTIL" DIFF "DO"    
       "BEGIN" D:=U[I] * 10 000 000 + U[I+1];   
          Q1:="IF" U[I]=V1 "THEN" 9 999 999 "ELSE" D//V1; 
          "IF" V[2] * Q1 > (D-Q1*V1) * 10 000 000 + U[I+2] "THEN"   
          "BEGIN" Q1:=Q1 - 1;         
             "IF" V[2]*Q1>(D-Q1*V1)*10 000 000+U[I+2] "THEN" Q1:=Q1-1         
          "END";  
        
          "COMMENT" U[I:I+LV]:=U[I:I+LV] - Q1 * V[1:LV];  
          CARRY:=0;         
          "FOR" J:=LV "STEP" -1 "UNTIL" 1 "DO"  
          "BEGIN" T:=Q1 * V[J] + CARRY; CARRY:=T//10 000 000;       
             A[J]:=T - CARRY * 10 000 000       
          "END"; A[0]:=CARRY;         
          CARRY:=0;         
          "FOR" J:=LV "STEP" -1 "UNTIL" 0 "DO"  
          "BEGIN" T:=U[I+J] - A[J] + CARRY; CARRY:="IF" T<0 "THEN" -1         
             "ELSE" 0; U[I+J]:=T - CARRY * 10 000 000     
          "END";  
        
          "COMMENT" IF CARRY=-1 THEN Q1 IS ONE TOO LARGE, 
          AND V MUST BE ADDED BACK TO U[I:I+LV];
          "IF" CARRY=-1 "THEN"        
          "BEGIN" Q1:=Q1 - 1; CARRY:=0;         
             "FOR" J:=LV "STEP" -1 "UNTIL" 1 "DO"         
             "BEGIN" T:=U[I+J] + V[J] + CARRY; CARRY:="IF" T<10 000 000       
                "THEN" 0 "ELSE" 1; U[I+J]:=T - CARRY * 10 000 000   
             "END"
          "END"; QUOTIENT[I]:=Q1      
       "END" I;   
        
       "COMMENT" CORRECT STORAGE OF QUOTIENT;   
       "IF" QUOTIENT[0] ^= 0 "THEN"   
       "BEGIN" "FOR" I:=DIFF "STEP" -1 "UNTIL" 0 "DO"     
          QUOTIENT[I+1]:=QUOTIENT[I]; QUOTIENT[0]:=DIFF + 1         
       "END" "ELSE"         
       "IF" QUOTIENT[1] ^= 0 "THEN" QUOTIENT[0]:=DIFF "ELSE"        
       "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" DIFF - 1 "DO"  
          QUOTIENT[I]:=QUOTIENT[I+1]; QUOTIENT[0]:=DIFF - 1         
       "END";                                                  "COMMENT"      
        
1SECTION : 1.4                (OCTOBER 1974)                     PAGE 13      
                                                                 ;  
        
        
       "COMMENT" REMAINDER := U[DIFF+1:LU]//SCALE;        
       "IF" SCALE>1 "THEN"  
       "BEGIN" CARRY:=0;    
          "FOR" I:=1 "STEP" 1 "UNTIL" LV "DO"   
          "BEGIN" T:=CARRY * 10 000 000 + U[DIFF+I];      
             REMAINDER[I]:=T//SCALE; CARRY:=T - REMAINDER[I] * SCALE
          "END"   
       "END" "ELSE"         
       "FOR" I:=1 "STEP" 1 "UNTIL" LV "DO" REMAINDER[I]:=U[DIFF+I]; 
        
       "COMMENT" CORRECT STORAGE OF REMAINDER;  
       I:=0; J:=LV;         
       "FOR" I:=I+1 "WHILE" REMAINDER[I]=0 "AND" J>1 "DO" J:=J-1;   
       REMAINDER[0]:=J;     
       "IF" J<LV "THEN"     
       "FOR" I:=1 "STEP" 1 "UNTIL" J "DO"       
       REMAINDER[I]:=REMAINDER[LV + I - J];     
        
       "COMMENT" UNNORMALIZE THE DIVISOR V;     
       "IF" SCALE>1 "THEN"  
       "BEGIN" CARRY:=0;    
          "FOR" I:=1 "STEP" 1 "UNTIL" LV "DO"   
          "BEGIN" T:=CARRY * 10 000 000 + V[I]; 
             V[I]:=T//SCALE; CARRY:=T - V[I] * SCALE      
          "END"   
       "END"      
    "END"         
 "END" LNG INT DIVIDE;      
         "EOP"    
        
 "CODE" 31204;    
 "PROCEDURE" LNG INT POWER(U,EXPONENT,RESULT);  
 "VALUE" EXPONENT; "INTEGER" EXPONENT; "INTEGER""ARRAY" U,RESULT;   
 "BEGIN""INTEGER" MAX,I,N;  
    "PROCEDURE" LNG INT MULT(U,V,PRODUCT); "CODE" 31202;  
    MAX:=U[0] * EXPONENT;   
    "BEGIN""INTEGER""ARRAY" Y,Z,H[0:MAX];       
       "COMMENT" Y:=1, Z:=U;
       Y[0]:=Y[1]:=1;       
       "FOR" I:=U[0] "STEP" -1 "UNTIL" 0 "DO" Z[I]:=U[I]; 
    HALVE: N:=EXPONENT//2; "IF" N+N=EXPONENT "THEN" "GOTO" SQUARE Z;
       LNG INT MULT(Y,Z,H); 
       "FOR" I:=H[0] "STEP" -1 "UNTIL" 0 "DO" Y[I]:=H[I]; 
       "IF" N=0 "THEN" "GOTO" READY;  
    SQUARE Z: LNG INT MULT(Z,Z,H);    
       "FOR" I:=H[0] "STEP" -1 "UNTIL" 0 "DO" Z[I]:=H[I]; 
       EXPONENT:=N; "GOTO"HALVE;      
    READY: "FOR" I:=Y[0] "STEP" -1 "UNTIL" 0 "DO" RESULT[I]:=Y[I]   
    "END"         
 "END" LNG INT POWER;       
         "EOP"    
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 1      
        
        
        
 AUTHOR:     C.G. VAN DER LAAN.       
        
        
 INSTITUTE:  RIJKSUNIVERSITEIT GRONINGEN.       
        
        
 RECEIVED:   740701.        
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS THE PROCEDURES:      
     SINSER FOR EVALUATING A SINE SERIES;       
     COSSER FOR EVALUTING A COSINE SERIES;      
     FOUSER,FOUSER1,FOUSER2 FOR EVALUATING A FOURIER SERIES         
     (IN FOUSER THE SERIES IS RESTRICTED TO A SERIES WITH SINE      
      COEFFICIENTS EQUAL TO COSINE COEFFICIENTS);         
     COMFOUSER,COMFOUSER1,COMFOUSER2 FOR EVALUATING A COMPLEX FOURIER         
     SERIES       
     (IN COMFOUSER THE SERIES IS RESTRICTED TO A SERIES WITH REAL   
     COEFFICIENTS).         
        
        
 KEYWORDS:        
        
     FINITE FOURIER SERIES EVALUATION,
     TRIGONOMETRIC POLYNOMIAL EVALUATION,       
     GOERTZEL,WATT,CLENSHAW,REINSCH ALGORITHM,  
     LINEAR THREE-TERM INHOMOGENEOUS RECURRENCE RELATION. 
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 2      
        
        
        
 SUBSECTION : SINSER.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "REAL""PROCEDURE"SINSER(N,THETA,B);        
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"B;      
        
     SINSER:= THE VALUE OF THE SINE SERIES      
              B[1]*SIN(THETA)+...+B[N]*SIN(N*THETA).      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE NUMBER OF TERMS IN THE SINE SERIES;         
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE SINE SERIES;      
     B:      <ARRAY IDENTIFIER>;      
             "ARRAY"B[1:N]; 
             ENTRY: THE COEFFICIENTS OF THE SINE SERIES.  
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER: N MULTIPLICATIONS; 3N ADDITIONS;    
                3 SINE/COSINE EVALUATIONS).     
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 3      
        
        
        
 SUBSECTION : COSSER.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "REAL""PROCEDURE"COSSER(N,THETA,A);        
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A;      
        
     COSSER:= THE VALUE OF THE COSINE SERIES    
              A[0]+A[1]*COS(THETA)+...+A[N]*COS(N*THETA). 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE TRIGONOMETRIC POLYNOMIAL.     
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE COSINE SERIES.    
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"A[0:N]; 
             ENTRY: THE COEFFICIENTS OF THE COSINE SERIES.
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: PROPORTIONAL TO N      
              (IN FIRST ORDER: N MULTIPLICATIONS; 3N ADDITIONS;     
               2 COSINE/SINE EVALUATIONS).      
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 4      
        
        
        
 SUBSECTION : FOUSER.       
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "REAL""PROCEDURE"FOUSER (N,THETA,A);       
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A;      
        
     FOUSER := THE VALUE OF THE FOURIER SERIES  
               A[0]+A[1]*(COS(THETA)+SIN(THETA))+...+A[N]*(COS(N*THETA)       
               +SIN(N*THETA)).        
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE TRIGONOMETRIC POLYNOMIAL;     
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE FOURIER SERIES;   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"A[0:N]; 
             ENTRY: THE COEFFICIENTS OF THE (FINITE) FOURIER SERIES.
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER: N MULTIPLICATIONS; 3N ADDITIONS;    
                3 COSINE/SINE EVALUATIONS).     
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 5      
        
        
        
 SUBSECTION : FOUSER1.      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "REAL""PROCEDURE"FOUSER1(N,THETA,A,B);     
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A,B;    
        
     FOUSER1:= THE VALUE OF THE FOURIER SERIES  
               A[0]+A[1]*COS(THETA)+B[1]*SIN(THETA)+...   
               +A[N]*COS(N*THETA)+B[N]*SIN(N*THETA).      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE TRIGONOMETRIC POLYNOMIAL;     
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE FOURIER SERIES;   
     A,B:    <ARRAY IDENTIFIER>;      
             "ARRAY"A[0:N],B[1:N];    
             ENTRY: THE COEFFICIENTS OF THE (FINITE) FOURIER SERIES,
                    WITH A[K] COEFFICIENT OF COS(K*THETA), (K=0,...,N)        
                    AND  B[K] COEFFICIENT OF SIN(K*THETA), (K=1,...,N).       
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER: 4N MULTIPLICATIONS; 4N ADDITIONS;   
                2 COSINE/SINE EVALUATIONS).     
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 6      
        
        
        
 SUBSECTION : FOUSER2.      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "REAL""PROCEDURE"FOUSER2(N,THETA,A,B);     
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A,B;    
        
     FOUSER2:= THE VALUE OF THE FOURIER SERIES  
               A[0]+A[1]*COS(THETA)+B[1]*SIN(THETA)+...   
               +A[N]*COS(N*THETA)+B[N]*SIN(N*THETA).      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE TRIGONOMETRIC POLYNOMIAL;     
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE FOURIER SERIES;   
     A,B:    <ARRAY IDENTIFIER>;      
             "ARRAY"A[0:N],B[1:N];    
             ENTRY: THE COEFFICIENTS OF THE (FINITE) FOURIER SERIES,
                    WITH A[K] COEFFICIENT OF COS(K*THETA), (K=0,...,N)        
                    AND  B[K] COEFFICIENT OF SIN(K*THETA), (K=1,...,N).       
        
        
 PROCEDURES USED: SINSER = CP31090,   
                  COSSER = CP31091.   
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER: 2N MULTIPLICATIONS; 6N ADDITIONS;   
                6 COSINE/SINE EVALUATIONS).     
        
        
 LANGUAGE : ALGOL 60.       
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 7      
        
        
        
 SUBSECTION : COMFOUSER.    
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE"COMFOUSER (N,THETA,A,RR,RI);    
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA,RR,RI;"ARRAY"A;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE POLYNOMIAL IN EXP(I*THETA);   
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE FOURIER SERIES;   
     A:      <ARRAY IDENTIFIER>;      
             "ARRAY"A[0:N]; 
             ENTRY: THE REAL COEFFICIENTS A[K] (K=0,...,N) IN THE SERIES      
                    FN(THETA)=A[0]+A[1]*EXP(I*THETA)+...+A[N]*EXP(I*
                    THETA)**N, MUST BE GIVEN IN ARRAY A;  
     RR,RI:  <VARIABLE>;    
             EXIT:  THE REAL PART AND THE IMAGINARY PART OF FN(THETA)         
                    ARE DELIVERED IN RR AND RI, RESPECTIVELY.       
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER:  N MULTIPLICATIONS; 3N ADDITIONS;   
                3 COSINE/SINE EVALUATIONS).     
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 8      
        
        
        
 SUBSECTION : COMFOUSER1.   
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE"COMFOUSER1(N,THETA,AR,AI,RR,RI);
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA,RR,RI;"ARRAY"AR,AI;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE POLYNOMIAL IN EXP(I*THETA);   
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE FOURIER SERIES;   
     AR,AI:  <ARRAY IDENTIFIER>;      
             "ARRAY"AR,AI[0:N];       
             ENTRY: THE REAL PART AND THE IMAGINARY PART OF THE COMPLEX       
                    COEFFICIENTS C[K] (K=0,...,N) IN THE SERIES     
                    FN(THETA)=C[0]+C[1]*EXP(I*THETA)+...+C[N]*EXP(I*
                              THETA)**N         
                    MUST BE GIVEN IN ARRAY AR AND AI, RESPECTIVELY; 
     RR,RI:  <VARIABLE>;    
             EXIT:  THE REAL PART AND THE IMAGINARY PART OF FN(THETA)         
                    ARE DELIVERED IN RR AND RI, RESPECTIVELY.       
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER: 4N MULTIPLICATIONS; 4N ADDITIONS;   
                2 COSINE/SINE EVALUATIONS).     
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORMANCE : SEE COMFOUSER2 (THIS SECTION).  
        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                      PAGE 9      
        
        
        
 SUBSECTION : COMFOUSER2.   
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THE PROCEDURE READS:        
     "PROCEDURE"COMFOUSER2(N,THETA,AR,AI,RR,RI);
     "VALUE"N,THETA;"INTEGER"N;"REAL"THETA,RR,RI;"ARRAY"AR,AI;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY: THE DEGREE OF THE POLYNOMIAL IN EXP(I*THETA);   
     THETA:  <ARITHMETIC EXPRESSION>; 
             ENTRY: THE ARGUMENT OF THE FOURIER SERIES;   
     AR,AI:  <ARRAY IDENTIFIER>;      
             "ARRAY"AR,AI[0:N];       
             ENTRY: THE REAL PART AND THE IMAGINARY PART OF THE COMPLEX       
                    COEFFICIENTS C[K] (K=0,...,N) IN THE SERIES     
                    FN(THETA)=C[0]+C[1]*EXP(I*THETA)+...+C[N]*EXP(I*
                              THETA)**N         
                    MUST BE GIVEN IN ARRAY AR AND AI, RESPECTIVELY; 
     RR,RI:  <VARIABLE>;    
             EXIT:  THE REAL PART AND THE IMAGINARY PART OF FN(THETA)         
                    ARE DELIVERED IN RR AND RI, RESPECTIVELY.       
        
        
 PROCEDURES USED: COMFOUSER= CP31095. 
        
        
 RUNNING TIME: PROPORTIONAL TO N      
               (IN FIRST ORDER: 2N MULTIPLICATIONS; 6N ADDITIONS;   
                6 COSINE/SINE EVALUATIONS).     
        
        
 LANGUAGE: ALGOL 60.        
        
        
 METHOD AND PERFORANCE:     
        
     FOR THE EVALUATION OF A FINITE FOURIER SERIES        
     (=TRIGONOMETRIC POLYNOMIAL OF DEGREE N SEE POLYA AND SZEGOE, 1971,       
     P. 76)       
     FN(THETA)=A[0]+A[1]*COS(THETA)+B[1]*SIN(THETA)+...+  
               A[N]*COS(N*THETA)+B[N]*SIN(N*THETA),       
     TWO ALGORITHMS ARE USED:         
1SECTION : 2.2.3.1            (OCTOBER 1974)                     PAGE 10      
        
        
        
     1.  HORNER SCHEME      
         LET C[K]=A[K]+I*B[K], K=0,...,N        
         AND Z=EXP(-I*THETA)
         THEN     
             FN(THETA)=RE(C[0]+C[1]*Z+...+C[N]*Z**N).     
         THE ALGORITHM IS GIVEN BY:   
             P:=C[N]        
             P:=P*Z+C[K], K=N-1,...,0 
             FN(THETA):=RE(P).        
         (FOUSER1)
     2.  A COMBINATION OF THE CLENSHAW ALGORITHM (SEE GENTLEMAN(1969,II)      
         , VAN DER LAAN, LUKE(1969, P.327-329) OR STOER(1972, P.62,63))       
         AND THE MODIFICATION OF REINSCH (SEE REINSCH(1967), VAN DER
         LAAN, STOER(1972, P.64,65)). 
         (SINSER,COSSER,FOUSER,FOUSER2)         
     A MODIFICATION OF THE IDEA OF NEWBERY IS NOT IMPLEMENTED BECAUSE         
     OF THE INTRODUCTION OF SINE (COSINE) TERMS IN A COSINE (SINE)  
     SERIES AND THE THE INEFFICIENCY OF THE ALGORITHM (SEE VAN DER  
     LAAN OR NEWBERY(1973)).
        
     FOR THE EVALUATION OF A FINITE COMPLEX FOURIER SERIES
     FN(THETA)=AR[0]+I*AI[0]+(AR[1]+I*AI[1])*EXP(I*THETA)+...       
               +(AR[N]+I*AI[N])*EXP(I*THETA)**N,
     TWO ALGORITHMS, IN REAL ARITHMETIC, ARE USED:        
     1. HORNER SCHEME       
        LET C[K]=AR[K]+I*AI[K], K=0,...,N       
        AND Z=EXP(I*THETA)  
        THEN      
            FN(THETA)=C[0]+C[1]*Z+...+C[N]*Z**N.
        THE ALGORITHM IS GIVEN BY     
            P:=C[N]         
            P:=P*Z+C[K], K=N-1,N-2,...,0        
            FN(THETA):=P.   
        (COMFOUSER1)        
     2. A COMBINATION OF THE CLENSHAW ALGORITHM AND THE MODIFICATION OF       
        REINSCH.  
        LET CAR=AR[0]+AR[1]*COS(THETA)+...+AR[N]*COS(N*THETA),      
            SAI=      AI[1]*SIN(THETA)+...+AI[N]*SIN(N*THETA),      
            SAR=      AR[1]*SIN(THETA)+...+AR[N]*SIN(N*THETA),      
            CAI=AI[0]+AI[1]*COS(THETA)+...+AI[N]*COS(N*THETA)       
        THEN FN(THETA)=CAR-SAI+I*(SAR+CAI).     
        (COMFOUSER,COMFOUSER2)        
     THE HORNER SCHEME IS IMPLEMENTED BECAUSE OF THE SIMPLICITY OF  
     THE ALGORITHM (ALTHOUGH THIS ALGORITHM IS LESS EFFICIENT THAN THE        
     GOERTZEL/WATT/CLENSHAW/REINSCH ALGORITHM) AND THE STABLE NATURE
     OF ORTHOGONAL TRANSFORMATIONS.   
     A COMBINATION OF THE ALGORITHM OF GOERTZEL/WATT/CLENSHAW AND THE         
     MODIFICATION OF REINSCH IS IMPLEMENTED BECAUSE OF THE EFFICIENCY         
     OF THE GWC ALGORITHM AND THE STABILITY OF THE MODIFICATION OF  
     REINSCH, ESPECIALLY FOR SMALL VALUES OF THE ARGUMENT (MOD. PI).
     AN UPPER BOUND FOR THE ERROR GROWTH IS GIVEN BY A LINEAR FUNCTION        
     OF THE DEGREE FOR BOTH (IMPLEMENTED) ALGORITHMS (SEE VAN DER LAAN).      
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                     PAGE 11      
        
        
        
 REFERENCES:      
        
     GENTLEMAN,W.M.(1969):  
     AN ERROR ANALYSIS OF GOERTZEL'S(WATT'S) METHOD FOR COMPUTING   
     FOURIER COEFFICIENTS.  
     COMP.J.,VOL.12,P.160-165.        
        
     LAAN,C.G.VAN DER(TO APPEAR):     
     ORTHOGONAL POLYNOMIALS IN NUMERICAL ANALYSIS 1.      
     ERROR ANALYSIS OF LINEAR TWO-TERM AND THREE-TERM RECURRENCE    
     RELATIONS.   
        
     LUKE,Y.L.(1969):       
     THE SPECIAL FUNCTIONS AND THEIR APPROXIMATIONS.VOL.1.
     ACADEMIC PRESS.        
        
     NEWBERY,A.C.R.(1973):  
     ERROR ANALYSIS FOR FOURIER SERIES EVALUATION.        
     MATH.COMP.,VOL.26,P.923-924.     
        
     POLYA,G. AND G.SZEGOE(1971):     
     AUFGABEN UND LEHRSAETZE AUS DER ANALYSIS II.         
     HEIDELBERGER TASCHENBUECHER 74. SPRINGER.  
        
     REINSCH,C.(1967):      
     A NOTE ON TRIGONOMETRIC INTERPOLATION.     
     BERICHT NR. 6709.      
     ABTEILUNG MATHEMATIK DER TECHNISCHEN UNIVERSITAET MUENCHEN.    
        
     STOER,J.(1972):        
     EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK 1.
     HEIDELBERGER TASCHENBUECHER 105. SPRINGER. 
        
        
 EXAMPLE OF USE:  
        
     THE FOURIER SERIES .5+COS(THETA)+SIN(THETA)
     IS EVALUATED FOR THE ARGUMENTS 0,PI/2,PI, BY MEANS OF FOUSER   
        
     "BEGIN""REAL"THETA,PI;"ARRAY"A[0:1];       
        "REAL""PROCEDURE"FOUSER(N,THETA,A  );"CODE" 31092;
        PI:=ARCTAN(1)*4;A[0]:=.5;A[1]:=1;       
        "FOR"THETA:=0,PI/2,PI"DO"     
        OUTPUT(61,"("/,B-D.DD")",FOUSER(1,THETA,A))       
     "END"        
        
      1.50        
      1.50        
     -0.50        
        
1SECTION : 2.2.3.1            (OCTOBER 1974)                     PAGE 12      
        
        
        
 SOURCE TEXTS:    
        
 "CODE" 31090;    
 "REAL""PROCEDURE"SINSER(N,THETA,B);  
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"B;
 "BEGIN""INTEGER"K;"REAL"C,CC,LAMBDA,H,DUN,UN,UN1;        
    C:=COS(THETA);
    "IF"C<-.5"THEN"         
    "BEGIN"LAMBDA:= 4*COS(THETA/2)**2;UN:=DUN:=0;         
           "FOR"K:=N"STEP"-1"UNTIL"1"DO"        
           "BEGIN"DUN:=LAMBDA*UN-DUN+B[K] ;     
                  UN:=DUN-UN;         
           "END"  
    "END""ELSE""IF"C> .5"THEN"        
    "BEGIN"LAMBDA:=-4*SIN(THETA/2)**2;UN:=DUN:=0;         
           "FOR"K:=N"STEP"-1"UNTIL"1"DO"        
           "BEGIN"DUN:=LAMBDA*UN+DUN+B[K] ;     
                  UN:=DUN+UN;         
           "END"  
    "END""ELSE"   
    "BEGIN"CC:=C+C;UN:=UN1:=0;        
           "FOR"K:=N"STEP"-1"UNTIL"1"DO"        
           "BEGIN"H:=CC*UN-UN1+B[K]; UN1 := UN; UN := H; "END"      
    "END";        
    SINSER:=UN*SIN(THETA)   
 "END"SINSER;     
         "EOP"    
        
 "CODE" 31091;    
 "REAL""PROCEDURE"COSSER(N,THETA,A);  
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A;
 "BEGIN""INTEGER"K;"REAL"C,CC,LAMBDA,H,DUN,UN,UN1;        
    C:=COS(THETA);
    "IF"C<-.5"THEN"         
    "BEGIN"LAMBDA:= 4*COS(THETA/2)**2;UN:=DUN:=0;         
           "FOR"K:=N"STEP"-1"UNTIL"0"DO"        
           "BEGIN"UN:=DUN-UN;         
                  DUN:=LAMBDA*UN-DUN+A[K]       
           "END";COSSER:=DUN-LAMBDA/2*UN        
    "END""ELSE""IF"C> .5"THEN"        
    "BEGIN"LAMBDA:=-4*SIN(THETA/2)**2;UN:=DUN:=0;         
           "FOR"K:=N"STEP"-1"UNTIL"0"DO"        
           "BEGIN"UN:=DUN+UN;         
                  DUN:=LAMBDA*UN+DUN+A[K]       
           "END";COSSER:=DUN-LAMBDA/2*UN        
    "END""ELSE"   
    "BEGIN"CC:=C+C;UN:=UN1:=0;        
           "FOR"K:=N"STEP"-1"UNTIL"1"DO"        
           "BEGIN"H:=CC*UN-UN1+A[K];  
                  UN1:=UN;UN:=H       
           "END";COSSER:=A[0]+UN*C-UN1
    "END"         
 "END"COSSER;     
         "EOP"    
1SECTION : 2.2.3.1            (OCTOBER 1974)                     PAGE 13      
        
        
        
 "CODE" 31092;    
 "REAL""PROCEDURE"FOUSER (N,THETA,A); 
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A;
 "BEGIN""INTEGER"K;"REAL"C,CC,LAMBDA,H,DUN,UN,UN1,C2,S2;  
    C:=COS(THETA);
    "IF"C<-.5"THEN"         
    "BEGIN"C2:=COS(THETA/2);LAMBDA:=4*C2**2;UN:=DUN:=0;   
           "FOR"K:=N"STEP"-1"UNTIL"0"DO"        
           "BEGIN"UN:=DUN-UN;         
                  DUN:=LAMBDA*UN-DUN+A[K]       
           "END";FOUSER :=DUN+2*C2*(SIN(THETA/2)-C2)*UN   
    "END""ELSE""IF"C> .5"THEN"        
    "BEGIN"S2:=SIN(THETA/2);LAMBDA:=-4*S2*S2;UN:=DUN:=0;  
           "FOR"K:=N"STEP"-1"UNTIL"0"DO"        
           "BEGIN"UN:=DUN+UN;         
                  DUN:=LAMBDA*UN+DUN+A[K]       
           "END";FOUSER :=DUN+2*S2*(S2+COS(THETA/2))*UN   
    "END""ELSE"   
    "BEGIN"CC:=C+C;UN:=UN1:=0;        
           "FOR"K:=N"STEP"-1"UNTIL"1"DO"        
           "BEGIN"H:=CC*UN-UN1+A[K];  
                  UN1:=UN;UN:=H       
           "END";FOUSER :=A[0]-UN1+(C+SIN(THETA))*UN      
    "END"         
 "END"FOUSER;     
         "EOP"    
        
 "CODE" 31093;    
 "REAL""PROCEDURE" FOUSER1(N,THETA,A,B);        
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A,B;        
 "BEGIN""INTEGER"I;"REAL"R,S,H,CO,SI; 
    R:=S:=0;CO:=COS(THETA);SI:=SIN(THETA);      
    "FOR"I:=N"STEP"-1"UNTIL"1"DO"     
    "BEGIN" H:=CO*R+SI*S+A[I];        
            S:=CO*S-SI*R+B[I];        
            R:=H  
    "END";FOUSER1:=CO*R+SI*S+A[0]     
 "END"FOUSER1;    
         "EOP"    
        
 "CODE" 31094;    
 "REAL""PROCEDURE"FOUSER2(N,THETA,A,B);         
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA;"ARRAY"A,B;        
 "BEGIN"
    "REAL""PROCEDURE"SINSER(N,THETA,B);"CODE" 31090;      
    "REAL""PROCEDURE"COSSER(N,THETA,A);"CODE" 31091;      
    FOUSER2:=COSSER(N,THETA,A)+SINSER(N,THETA,B);         
 "END"FOUSER2;    
         "EOP"    
1SECTION : 2.2.3.1            (OCTOBER 1974)                     PAGE 14      
        
        
        
 "CODE" 31095;    
 "PROCEDURE"COMFOUSER(N,THETA,A,RR,RI);         
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA,RR,RI;"ARRAY"A;    
 "BEGIN""INTEGER"K;"REAL"C,CC,LAMBDA,H,DUN,UN,UN1;        
    C:=COS(THETA);
    "IF"C<-.5"THEN"         
    "BEGIN"LAMBDA:= 4*COS(THETA/2)**2;UN:=DUN:=0;         
           "FOR"K:=N"STEP"-1"UNTIL"0"DO"        
           "BEGIN"UN:=DUN-UN;         
                  DUN:=LAMBDA*UN-DUN+A[K]       
           "END";RR    :=DUN-LAMBDA/2*UN        
    "END""ELSE""IF"C> .5"THEN"        
    "BEGIN"LAMBDA:=-4*SIN(THETA/2)**2;UN:=DUN:=0;         
           "FOR"K:=N"STEP"-1"UNTIL"0"DO"        
           "BEGIN"UN:=DUN+UN;         
                  DUN:=LAMBDA*UN+DUN+A[K]       
           "END";RR    :=DUN-LAMBDA/2*UN        
    "END""ELSE"   
    "BEGIN"CC:=C+C;UN:=UN1:=0;        
           "FOR"K:=N"STEP"-1"UNTIL"1"DO"        
           "BEGIN"H:=CC*UN-UN1+A[K];  
                  UN1:=UN;UN:=H       
           "END";RR    :=A[0]+UN*C-UN1
    "END";RI:=UN*SIN(THETA) 
 "END"COMFOUSER;  
         "EOP"    
        
 "CODE" 31096;    
 "PROCEDURE"COMFOUSER1(N,THETA,AR,AI,RR,RI);    
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA,RR,RI;"ARRAY"AR,AI;
 "BEGIN""INTEGER"K;"REAL"H,HR,HI,CO,SI;         
    HR:=HI:=0;CO:=COS(THETA);SI:=SIN(THETA);    
    "FOR"K:=N"STEP"-1"UNTIL"1"DO"     
    "BEGIN"H:=CO*HR-SI*HI+AR[K];      
          HI:=CO*HI+SI*HR+AI[K];      
          HR:=H   
    "END";        
    RR:=CO*HR-SI*HI+AR[0];  
    RI:=CO*HI+SI*HR+AI[0]   
 "END"COMFOUSER1; 
         "EOP"    
        
 "CODE" 31097;    
 "PROCEDURE"COMFOUSER2(N,THETA,AR,AI,RR,RI);    
 "VALUE"N,THETA;"INTEGER"N;"REAL"THETA,RR,RI;"ARRAY"AR,AI;
 "BEGIN""REAL"CAR,CAI,SAR,SAI;        
    "PROCEDURE"COMFOUSER(N,THETA,A,RR,RI);"CODE" 31095;   
    COMFOUSER(N,THETA,AR,CAR,SAR);    
    COMFOUSER(N,THETA,AI,CAI,SAI);    
    RR:=CAR-SAI;  
    RI:=CAI+SAR   
 "END"COMFOUSER2; 
         "EOP"    
1SECTION : 2.4.3              (OCTOBER 1974)                      PAGE 1      
        
        
 AUTHOR : C.G. VAN DER LAAN.
        
        
 INSTITUTE : RIJKSUNIVERSITEIT GRONINGEN.       
        
        
 RECEIVED : 740131.         
        
        
 BRIEF DESCRIPTION :        
        
     INTCHS COMPUTES THE INDEFINITE INTEGRAL OF A GIVEN CHEBYSHEV   
     SERIES.      
        
        
 KEYWORDS :       
        
     INDEFINITE INTEGRATION,
     CHEBYSHEV SERIES.      
        
        
 CALLING SEQUENCE :         
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE"INTCHS(N,A,B);        
     "VALUE"N;"INTEGER"N;"ARRAY"A,B;  
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     N  : <ARITHMETIC EXPRESSION>;    
          ENTRY:  
          THE DEGREE OF THE POLYNOMIAL REPRESENTED BY THE CHEBYSHEV 
          SERIES; 
     A,B: <ARRAY IDENTIFIER>;         
          "ARRAY" A[0:N],B[1:N+1];    
          ENTRY:  
          THE COEFFICIENTS OF THE CHEBYSHEV SERIES,A[0]+A[1]*T1(X)+...+       
          +A[N]*TN(X),SHOULD BE GIVEN IN ARRAY A.         
          EXIT:   
          THE COEFFICIENTS OF THE INTEGRAL CHEBYSHEV SERIES,        
          B[1]*T1(X)+...+B[N+1]*TN+1(X), ARE DELIVERED IN ARRAY B.  
          (T1(X),...TN+1(X) DENOTE CHEBYSHEV POLYNOMIALS OF THE FIRST         
          KIND,OF DEGREE 1,...N+1,RESPECTIVELY).
        
        
1SECTION : 2.4.3              (OCTOBER 1974)                      PAGE 2      
        
        
        
 METHOD AND PERFORMANCE :   
        
     FOR A DESCRIPTION OF THE ALGORITHM SEE AMONG OTHERS :
     CLENSHAW,1962,P.11,OR FOX AND PARKER,1968,P.59.      
        
        
 REFERENCES :     
        
     BROUCKE,R.(1973):      
     TEN SUBROUTINES FOR THE MANIPULATION OF CHEBYSHEV SERIES.      
     ALGORITHM 446.(FORTRAN).         
     COMM.ACM,VOL.16,1,P.254-256.     
        
     CLENSHAW,C.W.(1962):   
     CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS.         
     MATH.TAB.NAT.PHYS.LAB. 5,LONDON. 
     H.M. STATIONARY OFFICE.
        
     FOX,L.&I.B.PARKER(1968):         
     CHEBYSHEV POLYNOMIALS IN NUMERICAL ANALYSIS.         
     OXFORD UNIVERSITY PRESS.         
        
        
 EXAMPLE OF USE : 
        
     AS A FORMAL TEST OF THE PROCEDURE INTCHS THE CHEBYSHEV SERIES :
     1+1/2*T1(X)+1/5*T2(X)+1/10*T3(X) IS TRANSFORMED INTO ITS INTEGRAL.       
        
     "BEGIN""ARRAY"A[0:3],B[1:4];     
        "PROCEDURE"INTCHS(N,A,B);"CODE"31248;   
        A[0]:=1;A[1]:=.5;A[2]:=.2;A[3]:="-1;    
        INTCHS(3,A,B);      
        OUTPUT(61,"("/,4(BZ.4D)")",B[1],B[2],B[3],B[4]);  
     "END"        
        
        
   .9000  .1000  .0333  .0125         
1SECTION : 2.4.3              (OCTOBER 1974)                      PAGE 3      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE"31248;     
 "PROCEDURE"INTCHS(N,A,B);  
 "VALUE"N;"INTEGER"N;"ARRAY"A,B;      
 "COMMENT"        
     INTCHS DELIVERS THE COEFFICIENTS B[I],I=1,...N+1, OF THE INTEGRAL        
     CHEBYSHEV SERIES B[1]*T1(X)+...+B[N]*TN(X)+B[N+1]*TN+1(X).     
     THESE COEFFICIENTS ARE OBTAINED BY MEANS OF INDEFINITE INTEGRATION       
     OF THE CHEBYSHEV SERIES A[0]+A[1]*T1(X)+...+A[N]*TN(X).        
     T1(X),...TN+1(X) DENOTE CHEBYSHEV POLYNOMIALS OF THE FIRST KIND,         
     OF DEGREE 1,...N+1,RESPECTIVELY; 
 "IF"N=0"THEN"B[1]:=A[0]    
 "ELSE""IF"N=1"THEN""BEGIN"B[2]:=A[1]/4;B[1]:=A[0]"END"   
       "ELSE""BEGIN""INTEGER"I;"REAL"H,L,DUM;   
          H:=A[N];DUM:=A[N-1];B[N+1]:=H/((N+1)*2);B[N]:=DUM/(N*2);  
          "FOR"I:=N-1"STEP"-1"UNTIL"2"DO"       
          "BEGIN"L:=A[I-1];B[I]:=(L-H)/(2*I);H:=DUM;DUM:=L
          "END";B[1]:=A[0]-H/2        
       "END"INTCHS;         
         "EOP"    
1SECTION : 3.1.1.2.1.3        (OCTOBER 1974)                      PAGE 1      
        
        
        
 CONTRIBUTOR :    J. KOK.   
        
        
 INSTITUTE   :    MATHEMATICAL CENTRE.
        
        
 RECEIVED    :    740617.   
        
        
 BRIEF DESCRIPTION    :     
        
     THIS SECTION CONTAINS ONE PROCEDURE,       
     LSQINV, FOR THE CALCULATION OF THE INVERSE OF THE MATRIX  S'S, 
     WHERE  S  IS THE COEFFICIENT MATRIX OF A LINEAR LEAST SQUARES  
     PROBLEM.     
        
        
 KEYWORDS   :     
        
     INVERSE MATRIX,        
     LEAST SQUARES PROBLEM. 
        
        
 CALLING SEQUENCE    :      
        
     THE HEADING OF THE PROCEDURE READS :       
     "PROCEDURE" LSQINV(A, M, AID, C); "VALUE" M; "INTEGER" M;      
     "ARRAY" A, AID; "INTEGER""ARRAY" C;        
        
     THE MEANING OF THE FORMAL PARAMETERS IS   :
        
         A   :   <ARRAY IDENTIFIER>;  
                 "ARRAY" A[1 : M, 1 : M];       
                 ENTRY :  IN THE UPPER TRIANGLE OF A (THE ELEMENTS  
                 A[I,J] WITH  1 <= I < J <= M) THE SUPERDIAGONAL    
                 ELEMENTS SHOULD BE GIVEN OF THE UPPERTRIANGULAR MATRIX       
                 THAT IS PRODUCED BY THE HOUSEHOLDER TRIANGULARIZATION        
                 IN A CALL OF THE PROCEDURE  LSQORTDEC  (SECTION    
                 3.1.1.2.1.1.) WITH A NORMAL EXIT (AUX[3] = M).     
                 SEE ALSO THE MEANING OF THE PARAMETER  AID;        
                 EXIT :  THE UPPER TRIANGLE OF THE (SYMMETRIC) INVERSE        
                 MATRIX IS DELIVERED IN THE UPPERTRIANGULAR ELEMENTS OF       
                 THE ARRAY A (A[I,J] FOR  1 <= I <= J <= M);        
         M   :   <ARITHMETIC EXPRESSION>;       
                 NUMBER OF COLUMNS OF THE MATRIX OF THE LINEAR LEAST
                 SQUARES PROBLEM;     
         AID :   <ARRAY IDENTIFIER>;  
                 "ARRAY" AID[1 : M];  
                 ENTRY :  AID CONTAINS THE DIAGONAL ELEMENTS OF THE 
                 UPPERTRIANGULAR MATRIX THAT IS PRODUCED BY LSQORTDEC;        
         C   :   <ARRAY IDENTIFIER>;  
                 "INTEGER""ARRAY" C[1 : M];     
                 ENTRY :  C CONTAINS THE PIVOTAL INDICES AS PRODUCED BY       
                 A CALL OF LSQORTDEC. 
        
        
1SECTION : 3.1.1.2.1.3        (OCTOBER 1974)                      PAGE 2      
        
        
        
 PROCEDURES USED    :       
        
     CHLINV2   = CP34400,   
     ICHCOL    = CP34031,   
     ICHROW    = CP34032,   
     ICHROWCOL = CP34033.   
        
        
 REQUIRED CENTRAL MEMORY    :         
        
     EXECUTION FIELD LENGTH :  A REAL ARRAY OF M ELEMENTS IS        
     DECLARED (IN THE CALL OF CHLINV2).         
        
        
 RUNNING TIME    :  PROPORTIONAL TO M CUBED.    
        
        
 LANGUAGE    :  ALGOL 60.   
        
        
 METHOD AND PERFORMANCE  :  
        
     LSQINV  SHOULD  BE CALLED  AFTER  A SUCCESSFUL  CALL  OF  LSQORTDEC      
     (SECTION  3.1.1.2.1.1.).  LSQINV CAN BE USED FOR THE CALCULATION OF      
     THE COVARIANCE MATRIX OF A LINEAR LEAST SQUARES PROBLEM.       
     LET  S  BE THE MATRIX  OF THE LEAST SQUARES  SYSTEM  WITH  PERMUTED      
     COLUMNS  AND  Q * R  THE  HOUSEHOLDER DECOMPOSITION OF  S. THEN THE      
     INVERSE  OF   S'S  ALSO IS  THE INVERSE  OF  R'R.  SINCE  R  IS THE      
     CHOLESKY MATRIX OF  R'R,  THE INVERSE MATRIX IS COMPUTED  IN A CALL      
     OF  CHLINV2  (SECTION 3.1.1.1.1.2.4.).  AFTERWARDS  THE  COVARIANCE      
     MATRIX IS OBTAINED BY INTERCHANGES  OF THE COLUMNS AND ROWS  OF THE      
     INVERSE MATRIX.        
        
        
 EXAMPLE OF USE  :
        
     THE FOLLOWING PROGRAM COMPUTES THE INVERSE  T  OF  S'S, WHERE  S IS      
     THE COEFFICIENT  MATRIX  OF THE SYSTEM  IN THE EXAMPLE  OF  USE  OF      
     LSQORTDEC AND  LSQDGLINV (SECTION 3.1.1.2.1.1.). THE DIAGONAL OF  S      
     CAN BE COMPARED WITH THE RESULT OF  LSQDGLINV. TO CHECK THE ANSWERS      
     S' * (S * T)  IS PRINTED.        
        
1SECTION : 3.1.1.2.1.3        (OCTOBER 1974)                      PAGE 3      
        
        
        
     "BEGIN""COMMENT" JKOK, 740530, EXAMPLE OF USE OF  LSQORTDEC AND
         LSQINV;  
        
         "ARRAY" A, C[1 : 5, 1 : 2], AID[1 : 2], T[1 : 2, 1 : 2],   
         AUX[2 : 5];        
         "INTEGER""ARRAY" PIV[1 : 2]; 
         "INTEGER" I, J; "REAL" H;    
        
         "REAL""PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013;  
         "REAL""PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014;  
         "PROCEDURE" LSQORTDEC(A, N, M, AUX, AID, CI); "CODE" 34134;
         "PROCEDURE" LSQINV(A, M, AID, CI); "CODE" 34136; 
        
         AUX[2]:= "-12; I:= J:= 1;    
         "FOR" H:= - 2, - 1, 1, 2, 1, 1, 1, 1, 1, 2 "DO"  
         "BEGIN" A[I,J]:= C[I,J]:= H; "IF" I < 5 "THEN" I:= I + 1 "ELSE"      
             "BEGIN" I:= 1; J:= J + 1 "END"     
         "END";   
        
         LSQORTDEC(A, 5, 2, AUX, AID, PIV); "IF" AUX[3] = 2 "THEN"  
         "BEGIN" LSQINV(A, 2, AID, PIV);        
             T[1,1]:= A[1,1]; T[2,2]:= A[2,2]; T[2,1]:= T[1,2]:= A[1,2];      
             "FOR" J:= 1, 2 "DO""FOR" I:= 1 "STEP" 1 "UNTIL" 5 "DO" 
             A[I,J]:= MATMAT(1, 2, I, J, C, T); 
             OUTPUT(61, "("/4B, "(" AUX[2, 3, 5] = ")",   
                 -.4D"+DD5B, 3ZD5B, +.4D"+DD/,  
                 2(/4B, 30S, /, 2(/4B, 2(2B+.8D"+DD)), /) ")",      
                 AUX[2], AUX[3], AUX[5],        
                 "(" INVERSE :")",  ((T[I,J], J:= 1 : 2), I:= 1 : 2),         
                 "(" CHECK : S' * (S * T) :")", 
                 ((TAMMAT(1, 5, I, J, C, A), J:= 1 : 2), I:= 1 : 2) )         
         "END"    
     "END"        
         "EOP"    END OF PROGRAM      
        
     OUTPUT  :    
        
      AUX[2, 3, 5] =  .1000"-11        2     +.3317"+01   
        
      INVERSE :   
        
       +.95238095"-01  -.23809524"-01 
       -.23809524"-01  +.13095238"+00 
        
      CHECK : S' * (S * T) :
        
       +.10000000"+01  +.17763568"-14 
       +.00000000"+00  +.10000000"+01 
        
1SECTION : 3.1.1.2.1.3        (OCTOBER 1974)                      PAGE 4      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 34136;    
     "PROCEDURE" LSQINV(A, M, AID, C); "VALUE" M; "INTEGER" M;      
     "ARRAY" A, AID; "INTEGER""ARRAY" C;        
     "BEGIN""INTEGER" I, CI;
        "REAL" W; 
        
        "PROCEDURE" CHLINV2(A, N); "CODE" 34400;
        "PROCEDURE" ICHCOL(L, U, I, J, A); "CODE" 34031;  
        "PROCEDURE" ICHROW(L, U, I, J, A); "CODE" 34032;  
        "PROCEDURE" ICHROWCOL(L, U, I, J, A); "CODE" 34033;         
        
        "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" A[I,I]:= AID[I];        
        CHLINV2(A, M);      
        "FOR" I:= M "STEP" - 1 "UNTIL" 1 "DO"   
        "BEGIN" CI:= C[I]; "IF" CI ^= I "THEN"  
           "BEGIN" ICHCOL(1, I - 1, I, CI, A); ICHROWCOL(I + 1, CI - 1,       
              I, CI, A); ICHROW(CI + 1, M, I, CI, A);     
              W:= A[I,I]; A[I,I]:= A[CI,CI]; A[CI,CI]:= W 
           "END"  
        "END"     
     "END" LSQINV;
         "EOP"    
1SECTION : 0.0               (JULY 1979)                         PAGE 0       
1SECTION : 0.0               (JULY 1979)                         PAGE 0       
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  1      
        
        
 AUTHORS:        M. BAKKER, P.J. HARINGHUIZEN AND C.G. VAN DER LAAN.
        
        
 CONTRIBUTORS:   M. BAKKER, P.J. HARINGHUIZEN,  
                 C.G. VAN DER LAAN AND M. VOORINTHOLT.    
        
        
 INSTITUTE:      MATHEMATICAL CENTRE AND RIJKSUNIVERSITEIT GRONINGEN. 
        
        
 RECEIVED:       780601.    
        
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS FIVE PROCEDURES FOR CALCULATING ZEROS OF 
     ORTHOGONAL POLYNOMIALS WHICH ARE GIVEN BY THE COEFFICIENTS OF  
     THEIR RECURRENCE RELATION:       
     ALLZERORTPOL : CALCULATES ALL ZEROS,       
     LUPZERORTPOL : CALCULATES A NUMBER OF ADJACENT UPPER OR LOWER ZEROS,     
     SELZERORTPOL : CALCULATES A NUMBER OF ADJACENT ZEROS. IT IS    
                    EFFICIENT TO USE ALLZERORTPOL IF MORE THAN 50   
                    PERCENT OF EXTREME ZEROS OR MORE THAN 25 PERCENT OF       
                    SELECTED ZEROS ARE WANTED.  
     ALLJACZER    : CALCULATES THE ZEROS OF THE N-TH JACOBIAN POLYNOMIAL.     
     ALLLAGZER    : CALCULATES THE ZEROS OF THE N-TH LAGUERRE POLYNOMIAL.     
        
        
 KEYWORDS:
        
     ZEROS,       
     ORTHOGONAL POLYNOMIALS,
     CHRISTOFFEL ABSCISSAS. 
        
        
 REFERENCES:      
        
     ABRAMOWITZ, M. AND I.A. STEGUN (1964):     
     HANDBOOK OF MATHEMATICAL FUNCTIONS.
     DOVER PUBLICATIONS INC.
        
     GOLUB, G.H. AND J.H. WELSCH (1969):
     CALCULATION OF GAUSS QUADRATURE RULES.     
     MATH. COMP. VOL. 23, P.221-230.  
        
     LANCZOS, C. (1957):    
     APPLIED ANALYSIS.      
     PRENTICE HALL. 
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  2      
        
        
        
     STOER, J. (1972):      
     EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK 1.
     HEIDELBERG TASCHENBUECHER 105, SPRINGER.   
        
     WILKINSON,J AND REINSCH,C. :     
     HANDBOOK OF AUTOMATIC COMPUTATION. VOL. 2. 
     LINEAR ALGEBRA 
     HEIDELBERG (1971).     
        
        
 SUBSECTION:         ALLZERORTPOL.    
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
     "PROCEDURE" ALLZERORTPOL (N, B, C, ZER, EM); 
     "VALUE" N; "INTEGER" N; "ARRAY" B, C, ZER, EM;       
     "CODE" 31362;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY:  THE DEGREE OF THE ORTHOGONAL POLYNOMIAL OF WHICH 
                     THE ZEROS ARE TO BE CALCULATED;      
     B, C:   <ARRAY IDENTIFIER>;      
             "ARRAY" B, C [0:N-1];    
             ENTRY:  THE ELEMENTS B[I] AND C[I], I = 0, 1, ... , N-1, 
                     CONTAIN THE COEFFICIENTS OF THE RECURRENCE RELATION      
                     P[I+1](X) = (X - B[I]) * P[I](X) - C[I] * P[I-1](X)      
                     I = 0, 1, ... , N-1; ASSUMED IS C[0]=0, WHILE THE
                     CONTENTS OF THE ARRAYS ARE PRESERVED;
     ZER:    <ARRAY IDENTIFIER>;      
             "ARRAY" ZER[1:N];
             EXIT:   THE ZEROS OF THE N-TH DEGREE ORTHOGONAL POLYNOMIAL;      
                     (B MAY BE USED FOR ZER, BUT THEN THIS RECURRENCE 
                      COEFFICIENTS ARE OVERWRITTEN BY THE ZEROS AND 
                      THE ORIGINAL CONTENTS OF B ARE NOT PRESERVED.)
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:5]; 
             ENTRY:  EM[0]:  THE MACHINE PRECISION;       
                     EM[2]:  THE RELATIVE TOLERANCE OF THE ZEROS;   
                     EM[4]:  THE MAXIMUM ALLOWED NUMBER OF ITERATIONS,
                             E.G. 5 * N;
             EXIT:   EM[1]:  THE MAXIMUM OF ABS(B[0])+1,C[I]+ABS(B[I])+1,     
                             (I=1,...N-2),C[N-1]+ABS(B[N-1]);       
                     EM[3]:  INFORMATION CONCERNING THE PROCESS USED; 
                             I.E. THE MAXIMUM ABSOLUTE VALUE OF THE 
                             CODIAGONAL ELEMENTS NEGLECTED, 
                             (SEE ALSO SECTION 3.3.1.1.1.); 
                     EM[5]:  THE NUMBER OF ITERATIONS PERFORMED.    
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  3      
        
        
        
 PROCEDURES USED: 
        
     QRIVALSYMTRI    =  CP34160,      
     DUPVEC          =  CP31030.      
        
        
 METHOD AND PERFORMANCE : SEE SELZERORTPOL (THIS SECTION).
        
        
 EXAMPLE OF USE:  
        
     AS A FORMAL TEST OF THE PROCEDURE WE CALCULATE THE ZEROS OF    
     THE CHEBYSHEV POLYNOMIAL (OF THE FIRST KIND) OF THE THIRD DEGREE.
     THE RECURRENCE COEFFICIENTS ARE: 
     B[I] = 0, I = 0, 1, ....;
     C[0] = 0, C[1] = .5, C[I] = .25 , I = 2, 3, .....    
     (IT IS RECOMMENDED TO STORE THE ELEMENTS OF THE ARRAYS B AND C IN
      REVERSED ORDER IF THESE ELEMENTS ARE STRONGLY INCREASING).    
        
 "BEGIN" "ARRAY" B, C[0:3], ZER[1:3], EM[0:5];  
     "PROCEDURE" ALLZERORTPOL(N,B,C,ZER,EM);    
     "VALUE" N; "INTEGER" N; "ARRAY" B,C,ZER,EM;
     "CODE" 31362;
     EM[0]:= EM[2]:= "-14; EM[4]:=15; 
     B[2]:=B[1]:=B[0]:=0;   
     C[0]:= 0; C[1]:= .5; C[2]:= .25; 
     ALLZERORTPOL (3, B, C, ZER, EM); 
     OUTPUT(61,"(""("THE THREE ZEROS:")",/,3(/ZD5B,N),2/, 
               "("EM[1]:")",5BD.2D"+2D, /,      
               "("EM[3]:")",5BD.2D"+2D, /,"("EM[5]:")",5ZD")",      
               1,ZER[1],2,ZER[2],3,ZER[3],EM[1],EM[3],EM[5])
 "END"  
        
        
 THE THREE ZEROS: 
        
  1     -8.6602540378444"-001 
  2     +8.6602540378444"-001 
  3     -1.0000000000000"-014 
        
 EM[1]:     1.50"+00
 EM[3]:     7.07"-15
 EM[5]:     1     
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  4      
        
        
        
 SUBSECTION:         LUPZERORTPOL.    
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
     "PROCEDURE" LUPZERORTPOL (N, M, B, C, ZER, EM);      
     "VALUE" N, M; "INTEGER" N, M; "ARRAY" B, C, ZER, EM; 
     "CODE" 31363;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             ENTRY:  THE DEGREE OF THE ORTHOGONAL POLYNOMIAL OF WHICH 
                     THE ZEROS ARE TO BE CALCULATED;      
     M:      <ARITHMETIC EXPRESSION>; 
             ENTRY:  THE NUMBER OF ZEROS TO BE CALCULATED;
     B, C:   <ARRAY IDENTIFIER>;      
             "ARRAY" B, C [0:N-1];    
             ENTRY:  THE ELEMENTS B[I] AND C[I], I = 0, 1, ... , N-1, 
                     CONTAIN THE COEFFICIENTS OF THE RECURRENCE RELATION      
                     P[I+1](X) = (X - B[I]) * P[I](X) - C[I] * P[I-1](X)      
                     I = 0, 1, ... , N-1;       
                     ASSUMED IS C[0]=0, WHILE THE CONTENTS
                     OF THE ARRAYS ARE NOT PRESERVED;     
     ZER:    <ARRAY IDENTIFIER>;      
             "ARRAY" ZER[1:M];
             EXIT:   THE M LOWEST ZEROS ARE DELIVERED;    
                     IF HOWEVER THE ARRAY B[0:N-1] CONTAINED THE OPPOSIT      
                     VALUES OF THE CORRESPONDING RECURRENCE COEFFICIENTS      
                     THEN THE OPPOSITE VALUES OF THE M UPPER ZEROS  
                     ARE DELIVERED.   
                     IN EITHER CASE, ZER[I]<ZER[I+1], I = 1,...,M-1;
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:6]; 
             ENTRY:  EM[0]:  THE MACHINE PRECISION.       
                     EM[2]:  THE RELATIVE TOLERANCE OF THE ZEROS;   
                     EM[4]:  THE MAXIMUM ALLOWED NUMBER OF ITERATIONS,
                             E.G. 15 * M;       
                     EM[6]:  IF ALL ZEROS ARE KNOWN TO BE POSITIVE  
                             THEN 1 ELSE 0;     
             EXIT:   EM[1]:  THE MAXIMUM OF ABS(B[0]) + 1 , 
                             C[I]+ABS(B[I])+1,(I=1,...N-2), 
                             C[N-1]+ABS(B[N-1]);
                     EM[3]:  INFORMATION CONCERNING THE PROCESS USED, 
                             I.E. THE MAXIMUM ABSOLUTE VALUE OF THE 
                             THEORETICAL ERRORS OF THE ZEROS
                             (SEE WILKINSON AND REINSCH, 1971, P.263);
                     EM[5]:  THE NUMBER OF ITERATIONS PERFORMED.    
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  5      
        
        
        
 PROCEDURES USED: 
        
     DUPVEC          =  CP31030,      
     INFNRMVEC       =  CP31061.      
        
        
 METHOD AND PERFORMANCE : SEE SELZERORTPOL (THIS SECTION).
        
        
 EXAMPLE OF USE:  
        
     WE CALCULATE THE TWO LOWER AND THE TWO UPPER ZEROS OF THE      
     LAGUERRE POLYNOMIAL OF THE THIRD DEGREE.   
     THE RECURRENCE COEFFICIENTS ARE OBTAINED FROM [1], P.782:      
     B[I] = - A2I / A3I = 2I + 1;     
     C[I] = A4I / (A3I * A3(I-1)) * A1(I-1) = I * I, I = 0, 1, 2.   
     (IT IS RECOMMENDED TO STORE THE ELEMENTS OF THE ARRAYS B AND C IN
      REVERSED ORDER IF THESE ELEMENTS ARE STRONGLY DECREASING).    
        
     "BEGIN" "ARRAY" B, C[0:3], ZER[1:2], EM[0:6];
     "INTEGER" I; 
     "PROCEDURE"LUPZERORTPOL(N,M,B,C,ZER,EM);   
     "VALUE" N,M; "INTEGER" N,M; "ARRAY" B,C,ZER,EM;      
     "CODE" 31363;
     EM[0]:= EM[2]:= "-14; EM[4]:= 45; EM[6]:= 1; 
     "FOR" I:= 0, 1, 2 "DO" 
     "BEGIN" B[I]:= 2 * I + 1; C[I]:= I * I "END";
     LUPZERORTPOL (3, 2, B, C, ZER, EM);
     OUTPUT(61,"(""("THE TWO LOWER ZEROS:")",/,2(/ZD5B,N),2/,       
               "("EM[1]:")",5BD.2D"+2D, /,      
               "("EM[3]:")",5BD.2D"+2D, /,"("EM[5]:")",5ZD")",      
               1,ZER[1],2,ZER[2],EM[1],EM[3],EM[5]);      
     EM[6]:= 0;   
     "FOR" I:= 0, 1, 2 "DO" 
     "BEGIN" B[I]:= - 2 * I - 1; C[I]:= I * I "END";      
     LUPZERORTPOL (3, 2, B, C, ZER, EM);
     OUTPUT(61,"("3/,"("THE TWO UPPER ZEROS:")",/,2(/ZD5B,N),2/,    
               "("EM[1]:")",5BD.2D"+2D, /,      
               "("EM[3]:")",5BD.2D"+2D, /,"("EM[5]:")",5ZD")",      
               1,-ZER[1],2,-ZER[2],EM[1],EM[3],EM[5]);    
     "END"
        
        
     THE TWO LOWER ZEROS:   
        
      1     +4.1577455678348"-001     
      2     +2.2942803602791"+000     
        
     EM[1]:     9.00"+00    
     EM[3]:     5.72"-16    
     EM[5]:    12 
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  6      
        
        
        
     THE TWO UPPER ZEROS:   
        
      1     +6.2899450829375"+000     
      2     +2.2942803602791"+000     
        
     EM[1]:     9.00"+00    
     EM[3]:     4.70"-20    
     EM[5]:    14 
        
        
 SUBSECTION:         SELZERORTPOL.    
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
     "PROCEDURE" SELZERORTPOL (N, N1, N2, B, C, ZER, EM); 
     "VALUE" N, N1, N2; "INTEGER" N, N1, N2; "ARRAY" B, C, ZER, EM; 
     "CODE" 31364;
        
     THE MEANING OF THE FORMAL PARAMETERS IS :  
     N,N1,N2:<ARITHMETIC EXPRESSION>; 
             ENTRY:  N IS THE DEGREE OF THE ORTHOGONAL POLYNOMIAL OF
                     WHICH THE N1-ST UP TO AND INCLUDING N2-ND ZEROS ARE      
                     TO BE CALCULATED(ZER[N1]>=ZER[N2]);  
     B, C:   <ARRAY IDENTIFIER>;      
             "ARRAY" B,C [0 : N-1];   
             ENTRY:  THE ELEMENTS B[I] AND C[I], I = 0, 1, ... , N-1, 
                     CONTAIN THE COEFFICIENTS OF THE RECURRENCE RELATION      
                     P[I+1](X) = (X - B[I]) * P[I](X) - C[I] * P[I-1](X)      
                     I = 0, 1, ... , N-1,       
                     ASSUMED IS C[0]=0, WHILE THE CONTENTS
                     OF THE ARRAYS IS PRESERVED;
     ZER:    <ARRAY IDENTIFIER>;      
             "ARRAY" ZER [N1:N2];     
             EXIT:   THE N2-N1+1 CALCULATED ZEROS IN DECREASING ORDER.
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:5]; 
             ENTRY:  EM[0]:  THE MACHINE PRECISION.       
                     EM[2]:  THE RELATIVE TOLERANCE OF THE ZEROS;   
             EXIT:   EM[1]:  THE MAXIMUM OF ABS(B[0]) + 1,
                             C[I]+ABS(B[I])+1 (I=1,...N-2) AND      
                             C[N-1]+ABS(B[N-1]);
                     EM[5]:  THE NUMBER OF ITERATIONS PERFORMED.    
        
        
 PROCEDURES USED: 
        
    VALSYMTRI       =  CP34151.       
1SECTION : 3.6.2                  (DECEMBER 1978)                PAGE  7      
        
        
        
 METHOD AND PERFORMANCE:    
        
     THE ZEROS OF AN ORTHOGONAL POLYNOMIAL ARE THE EIGENVALUES OF A 
     SYMMETRIC TRIDIAGONAL MATRIX (SEE GOLUB AND WELSCH (1969),     
     LANCZOS (1957),P.375,376, STOER (1972),P.120).       
     THE ORTHOGONAL POLYNOMIAL IS DEFINED BY A LINEAR THREE-TERM    
     HOMOGENEOUS RECURRENCE RELATION. 
        
        
 EXAMPLE OF USE : 
        
     WE CALCULATE THE THIRD ZERO OF THE LEGENDRE POLYNOMIAL OF THE  
     FOURTH DEGREE. THE RECURRENCE COEFFICIENTS ARE OBTAINED FROM   
     ABRAMOWITZ AND STEGUN (1964),P.782:
     B[I] = 0, I = 0, 1, ....;
     C[I] = A4I / (A3I * A3(I-1)) * A1(I-1) = I * I / ( 4 * I * I - 1),       
             I = 0, 1, .....
     (IT IS RECOMMENDED TO STORE THE ELEMENTS OF THE ARRAYS B AND C IN
      REVERSED ORDER IF THESE ELEMENTS ARE STRONGLY DECREASING).    
        
     "BEGIN" "ARRAY" B, C[0:4], ZER[3:3], EM[0:5];
     "INTEGER" I; 
     "PROCEDURE" SELZERORTPOL (N,N1,N2,B,C,ZER,EM);       
     "VALUE" N,N1,N2; "INTEGER" N,N1,N2; "ARRAY" B,C,ZER,EM;
     "CODE" 31364;
     EM[0]:= EM[2]:= "-14;  
     "FOR" I:= 0, 1, 2, 3 "DO"
     "BEGIN" B[I]:= 0; C[I]:= I * I / (4 * I * I - 1) "END";
     SELZERORTPOL (4, 3, 3, B, C, ZER, EM);     
     OUTPUT(61,"(""("THE THIRD ZERO:")",2/,ZD5B,N,2/,     
               "("EM[1]:")",5BD.2D"+2D, /,      
               "("EM[5]:")",5ZD")",3,ZER[3],EM[1],EM[5])  
     "END"
        
     THE THIRD ZERO:
        
      3     -3.3998104358486"-001     
        
     EM[1]:     1.33"+00    
     EM[5]:    12 
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 8       
        
        
        
  SUBSECTION: ALL JAC ZER.  
        
        
  CALLING SEQUENCE: 
        
       THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
       "PROCEDURE" ALL JAC ZER(N, ALFA, BETA, ZER);       
       "VALUE" N, ALFA, BETA; 
       "INTEGER" N; "REAL" ALFA, BETA; "ARRAY" ZER;       
       "CODE" 31370;
        
       THE MEANING OF THE FORMAL PARAMETERS IS: 
       N:      <ARITHMETIC EXPRESSION>; 
               THE UPPER BOUND OF THE ARRAY ZER; N >= 1;  
       ALFA,BETA: 
               <ARITHMETIC EXPRESSION>; 
               THE PARAMETERS OF THE JACOBI POLYNOMIAL,   
               SEE ABRAMOWITZ AND STEGUN (1964);
               ALFA, BETA > - 1;      
       ZER:    <ARRAY IDENTIFIER>;    
               "ARRAY" ZER[1 : N];    
               EXIT: ZER[1], ..., ZER[N] ARE THE ZEROS OF THE N-TH  
                     JACOBI POLYNOMIAL WITH PARAMETERS ALFA AND BETA. 
        
        
  PROCEDURES USED:
        
       ALL ZER ORT POL = CP 31362.    
        
        
  REQUIRED CENTRAL MEMORY:  
        
       IF ALFA = BETA THEN TWO AUXILIARY ARRAYS OF N//2 REALS ARE   
       USED, OTHERWISE TWO AUXILIARY ARRAYS OF N REALS ARE DECLARED.
        
        
  METHOD AND PERFORMANCE:   
        
       THE JACOBI POLYNOMIALS ARE A SPECIAL CASE OF ORTHOGONAL      
       POLYNOMIALS (SEE ABRAMOWITZ AND STEGUN (1964)); ALL JAC ZER  
       COMPUTES THE COEFFICIENTS OF THE THREE-TERM RECURRENCE RELATION
       AND CALLS THE PROCEDURE ALL ZER ORT POL TO COMPUTE THE ZEROS;
       IF ALFA=BETA, THE POLYNOMIALS ARE ODD OR EVEN, HENCE ONLY THE
       POSITIVE ZEROS ARE CALCULATED; THIS IS DONE BY MEANS OF THE  
       FORMULAS   
        
             P(2*M, ALFA, ALFA, X) = C(M)*P(M, ALFA, -0.5, 2*X*X - 1),
        
         P(2*M - 1, ALFA, ALFA, X) = D(M)*P(M, ALFA, +0.5, 2*X*X - 1)*X       
        
(      (SEE ABRAMOWITZ AND STEGUN (1964), FORMULAS 22.5.20 - 22.5.27).
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 9       
        
        
        
  EXAMPLE OF USE: 
        
       THE PROGRAM
        
       "BEGIN" "ARRAY" X[1:3];
         "PROCEDURE" ALL JAC ZER(N, ALFA, BETA, ZER);     
         "VALUE" N,ALFA,BETA; "INTEGER" N;      
         "REAL" ALFA,BETA; "ARRAY" ZER; 
         "CODE" 31370;      
         ALL JAC ZER(3,-.5,-.5,X);    
         OUTPUT(61,"("3(4B-D.13D"-ZD)")",X[1],X[2],X[3])  
       "END"      
        
       DELIVERS THE FOLOWING RESULTS: 
        
       -8.6602540378444"-1    0.0000000000000" 0     8.6602540378444"-1       
        
        
  SUBSECTION: ALL LAG ZER.  
        
        
  CALLING SEQUENCE: 
        
       THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
        
       "PROCEDURE" ALL LAG ZER(N, ALFA, ZER);   
       "VALUE" N, ALFA;     
       "INTEGER" N; "REAL" ALFA; "ARRAY" ZER;   
       "CODE" 31371;
        
       THE MEANING OF THE FORMAL PARAMETERS IS: 
       N:      <ARITHMETIC EXPRESSION>; 
               THE UPPER BOUND OF THE ARRAY ZER; N >= 1;  
       ALFA:   <ARITHMETIC EXPRESSION>; 
               THE PARAMETER OF THE LAGUERRE POLYNOMIAL,  
               SEE ABRAMOWITZ AND STEGUN (1964); ALFA > -1; 
       ZER:    <ARRAY IDENTIFIER>;    
               "ARRAY" ZER[1 : N];    
               EXIT: ZER[1], ..., ZER[N] ARE THE ZEROS OF THE N-TH  
                     LAGUERRE POLYNOMIAL WITH PARAMETER ALFA.       
        
        
  PROCEDURES USED:
        
       ALL ZER ORT POL = CP 31362.    
        
        
  REQUIRED CENTRAL MEMORY:  
        
       TWO AUXILIARY ARRAYS OF N REALS ARE USED.
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 10      
        
        
        
  METHOD AND PERFORMANCE:   
        
       THE LAGUERRE POLYNOMIALS ARE A SPECIAL CASE OF ORTHOGONAL    
       POLYNOMIALS (SEE ABRAMOWITZ AND STEGUN (1964)); ALL LAG ZER  
       COMPUTES THE COEFFICIENTS OF THE THREE-TERM RECURRENCE RELATION
       AND CALLS THE PROCEDURE ALL ZER ORT POL TO COMPUTE THE ZEROS.
        
        
  EXAMPLE OF USE: 
        
       "BEGIN" "ARRAY" X[1:3];
         "PROCEDURE" ALL LAG ZER(N, ALFA, ZER); 
         "VALUE" N,ALFA; "INTEGER" N; "REAL" ALFA; "ARRAY" ZER;     
         "CODE" 31371;      
         ALL LAG ZER(3,-.5,X);
         OUTPUT(61,"("3(4B-D.13D"-ZD)")",X[1],X[2],X[3])  
       "END"      
        
       DELIVERS THE FOLOWING RESULTS: 
        
       5.5253437422633" 0.    1.7844927485432" 0     1.9016350919350" -1      
        
        
 SOURCE TEXT(S) : 
        
 "CODE"31362;     
 "PROCEDURE" ALLZERORTPOL (N, B, C, ZER, EM);   
 "VALUE" N; "INTEGER" N; "ARRAY" B, C, ZER, EM; 
 "BEGIN" "INTEGER"I;"REAL"NRM;"ARRAY"BB[1:N];   
     "INTEGER" "PROCEDURE" QRIVALSYMTRI (D, BB, N, EM);   
     "VALUE"N;"INTEGER"N;"ARRAY"D,BB,EM;
     "CODE" 34160;
     "PROCEDURE" DUPVEC (L, U, SHIFT, A, B);    
     "VALUE"L,U,SHIFT;"INTEGER"L,U,SHIFT;"ARRAY"A,B;      
     "CODE" 31030;
     "PROCEDURE" DUPCEV (L, U, SHIFT, A, B);    
     "VALUE"L,U,SHIFT;"INTEGER"L,U,SHIFT;"ARRAY"A,B;      
     "FOR" U:=U "STEP" -1 "UNTIL" L "DO" A[U]:=B[U+SHIFT];
     NRM:=ABS(B[0]);
     "FOR"I:=1"STEP"1"UNTIL"N-2"DO""IF"C[I]+ABS(B[I])>NRM"THEN"     
          NRM:=C[I]+ABS(B[I]);
      "IF"N>1"THEN"NRM:="IF"NRM+1>=C[N-1]+ABS(B[N-1])"THEN"NRM+1"ELSE"
                        C[N-1]+ABS(B[N-1]);     
      EM[1]:=NRM; DUPCEV(1,N,-1,ZER,B); 
      DUPVEC(1,N-1,0,BB,C);BB[N]:=0;  
      QRIVALSYMTRI(ZER,BB,N,EM)       
 "END" ALLZERORTPOL;
         "EOP"    
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 11      
        
        
        
 "CODE"31363;     
 "PROCEDURE" LUPZERORTPOL (N, M, B, C, ZER, EM);
 "VALUE" N, M; "INTEGER" N, M; "ARRAY" B, C, ZER, EM;     
 "BEGIN"
 "PROCEDURE" RATQR(N,M,POSDEF,DLAM,EPS)TRANS:(D,B2);      
   "VALUE" N,M,POSDEF,DLAM,EPS;       
   "INTEGER" N,M; 
   "BOOLEAN" POSDEF;
   "REAL" DLAM,EPS; 
   "ARRAY" D,B2;  
 "COMMENT" QR ALGORITHM FOR THE COMPUTATION OF THE LOWEST EIGENVALUES 
   OF A SYMMETRIC TRIDIAGONAL MATRIX. A RATIONAL VARIANT OF THE     
   QR TRANSFORMATION IS USED, CONSISTING OF TWO SUCCESSIVE QD STEPS 
   PER ITERATION. 
   A SHIFT OF THE SPECTRUM AFTER EACH ITERATION GIVES AN ACCELERATED
   RATE OF CONVERGENCE. A NEWTON CORRECTION,DERIVED FROM THE
   CHARACTERISTIC POLYNOMIAL,IS USED AS SHIFT.  
    RATQR IS IMPLEMENTED BY REINSCH AND BAUER, SEE WILKINSON AND REINSCH      
    ,1971, CONTR. II-6.     
   FORMATS: D,B2[1:N];      
 "BEGIN"
   "INTEGER" I,J,K,T; "REAL" DELTA,E,EP,ERR,P,Q,QP,R,S,TOT; 
   "REAL""PROCEDURE" INFNRMVEC(L,U,K,A);
   "VALUE" L,U; "INTEGER" L,U,K; "ARRAY" A;     
   "CODE" 31061;  
 "COMMENT" LOWER BOUND FOR EIGENVALUES FROM GERSHGORIN, INITIAL SHIFT;
   B2[1]:= ERR:= Q:= S:= 0; TOT:= D[1]; 
   "FOR" I:= N "STEP" -1   "UNTIL" 1 "DO"       
   "BEGIN"
     P:= Q; Q:= SQRT(B2[I]); E:= D[I]-P-Q;      
     "IF" E < TOT "THEN" TOT:= E      
   "END" I;       
   "IF" POSDEF & TOT < 0 "THEN" TOT:= 0 "ELSE"  
   "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" D[I]:= D[I]-TOT;   
   T:= 0; 
   "FOR" K:= 1 "STEP" 1 "UNTIL" M "DO"
   "BEGIN"
 NEXT QR TRANSFORMATION:  T:= T + 1;  
   TOT:= TOT + S; DELTA:= D[N]-S; I:= N;
   E:= ABS(EPS*TOT); "IF" DLAM < E "THEN" DLAM:= E;       
   "IF" DELTA < = DLAM "THEN" "GOTO" CONVERGENCE; 
   E:= B2[N]/DELTA; QP:= DELTA+E; P:= 1;
   "FOR" I:= N-1 "STEP" -1 "UNTIL" K "DO"       
   "BEGIN"
     Q:= D[I]-S-E; R:= Q/QP; P:= P*R+1; 
     EP:= E*R; D[I+1]:= QP+EP; DELTA:= Q-EP;    
     "IF" DELTA < = DLAM "THEN" "GOTO" CONVERGENCE;       
     E:= B2[I]/Q;QP:= DELTA+E; B2[I+1]:= QP*EP  
   "END" I;       
   D[K]:= QP; S:= QP/P;     
                                                               "COMMENT"      
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 12      
                                                                 ;  
        
        
   "IF" TOT+S > TOT "THEN" "GOTO" NEXT QR TRANSFORMATION; 
 "COMMENT" IRREGULAR END OF ITERATION,
   DEFLATE MINIMUM DIAGONAL ELEMENT;  
   S:= 0; I:= K; DELTA:= QP;
   "FOR" J:= K+1 "STEP" 1 "UNTIL" N "DO"
     "IF" D[J] < DELTA "THEN" 
     "BEGIN" I:= J; DELTA:= D[J] "END"; 
 CONVERGENCE:     
   "IF" I < N "THEN" B2[I+1]:= B2[I]*E/QP;      
     "FOR" J:= I-1 "STEP" -1 "UNTIL" K "DO"     
       "BEGIN" D[J+1]:= D[J]-S; B2[J+1]:= B2[J] "END" J;  
     D[K]:= TOT; B2[K]:= ERR:= ERR+ABS(DELTA)   
   "END" K;       
 EM[5]:=T;EM[3]:=INFNRMVEC(1,M,T,B2); 
 "END" RATQR;     
     "PROCEDURE" DUPCEV (L, U, SHIFT, A, B);    
     "VALUE"L,U,SHIFT;"INTEGER"L,U,SHIFT;"ARRAY"A,B;      
     "FOR" U:=U "STEP" -1 "UNTIL" L "DO" A[U]:=B[U+SHIFT];
     "PROCEDURE" DUPVEC (L, U, SHIFT, A, B);    
     "VALUE"L,U,SHIFT;"INTEGER"L,U,SHIFT;"ARRAY"A,B;      
     "CODE" 31030;
     "INTEGER" I;"REAL"NRM; 
     NRM:=ABS(B[0]);
     "FOR"I:=1"STEP"1"UNTIL"N-2"DO""IF"C[I]+ABS(B[I])>NRM"THEN"     
          NRM:=C[I]+ABS(B[I]);
      "IF"N>1"THEN"NRM:="IF"NRM+1>=C[N-1]+ABS(B[N-1])"THEN"NRM+1"ELSE"
                        C[N-1]+ABS(B[N-1]);     
      EM[1]:=NRM; 
     DUPCEV(1,N,-1,B,B);    
     DUPCEV(2,N,-1,C,C);    
     RATQR (N, M, EM[6] = 1, EM[2], EM[0], B, C); 
     DUPVEC (1, M, 0, ZER, B) 
 "END" LUPZERORTPOL;
         "EOP"    
        
        
 "CODE"31364;     
 "PROCEDURE" SELZERORTPOL (N, N1, N2, B, C, ZER, EM);     
 "VALUE" N, N1, N2; "INTEGER" N, N1, N2; "ARRAY" B, C, ZER, EM;     
 "BEGIN" "INTEGER"I;"REAL"NRM;"ARRAY"D[1:N];    
     "PROCEDURE" DUPCEV (L, U, SHIFT, A, B);    
     "VALUE"L,U,SHIFT;"INTEGER"L,U,SHIFT;"ARRAY"A,B;      
     "FOR" U:=U "STEP" -1 "UNTIL" L "DO" A[U]:=B[U+SHIFT];
     "PROCEDURE" VALSYMTRI (D, BB, N, N1, N2, VAL, EM);   
     "VALUE"N,N1,N2;"INTEGER"N,N1,N2;"ARRAY"D,BB,VAL,EM;  
     "CODE" 34151;
     NRM:=ABS(B[0]);
     "FOR"I:=N-2"STEP"-1"UNTIL"1"DO""IF"C[I]+ABS(B[I])>NRM"THEN"    
          NRM:=C[I]+ABS(B[I]);
                                                               "COMMENT"      
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 13      
                                                                 ;  
        
        
      "IF"N>1"THEN"NRM:="IF"NRM+1>=C[N-1]+ABS(B[N-1])"THEN"NRM+1"ELSE"
                        C[N-1]+ABS(B[N-1]);     
      EM[1]:=NRM; 
     DUPCEV(1,N,-1,D,B);    
     VALSYMTRI (D, C, N, N1, N2, ZER, EM);      
     EM[5]:=EM[3] 
 "END" SELZERORTPOL;
         "EOP"    
        
        
  "CODE" 31370;   
     "PROCEDURE" ALL JAC ZER(N, ALFA, BETA, ZER); 
     "VALUE" N, ALFA, BETA     ; "INTEGER" N;   
     "REAL" ALFA, BETA     ; "ARRAY" ZER;       
     "IF" ALFA = BETA "THEN"
     "BEGIN" "INTEGER" I, M;
         "ARRAY" A, B[0:N//2], EM[0:5]; 
         "REAL" MIN, GAMMA, SUM, ZERI;
         "REAL""PROCEDURE" ARREB;     
         "CODE" 30002;      
         "PROCEDURE" ALL ZER ORT POL(N, B, C, ZER, EM);   
         "VALUE" N; "INTEGER" N; "ARRAY" B,C,ZER,EM;      
         "CODE" 31362;      
         M:= N//2; "IF" N ^= 2*M "THEN" 
         "BEGIN" GAMMA:= + 0.5; ZER[M + 1]:= 0 "END"      
         "ELSE" GAMMA:= - 0.5;
         MIN:= 0.25 - ALFA*ALFA; SUM:= ALFA + GAMMA + 2;  
         A[0]:= (GAMMA - ALFA)/SUM; A[1]:= MIN/SUM/(SUM + 2);       
         B[1]:= 4*(1 + ALFA)*(1 + GAMMA)/SUM/SUM/(SUM + 1); 
         "FOR" I:= 2 "STEP" 1 "UNTIL" M - 1 "DO"
         "BEGIN" SUM:= I + I + ALFA + GAMMA;    
             A[I]:= MIN/SUM/(SUM + 2); SUM := SUM*SUM;    
             B[I]:= 4*I*(I + ALFA + GAMMA)*(I + ALFA)*(I + GAMMA)/  
             SUM/(SUM - 1)  
         "END";   
         EM[0]:=ARREB; EM[2]:="-10; EM[4]:= 6*M;
         ALL ZER ORT POL (M, A, B, ZER, EM);    
         "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"    
         "BEGIN" ZER[I]:= ZERI:= - SQRT((1 + ZER[I])/2);  
             ZER[N + 1 - I]:= - ZERI  
         "END"    
     "END" "ELSE" 
     "BEGIN" "INTEGER" I; "REAL" SUM, MIN;      
         "ARRAY" A, B[0:N], EM[0:5];  
         "REAL""PROCEDURE" ARREB;     
         "CODE" 30002;      
         "PROCEDURE" ALL ZER ORT POL(N, B, C, ZER, EM);   
         "VALUE" N; "INTEGER" N; "ARRAY" B,C,ZER,EM;      
         "CODE" 31362;      
                                                               "COMMENT"      
1SECTION : 3.6.2             (DECEMBER 1978)                     PAGE 14      
                                                                 ;  
        
        
         MIN:= (BETA - ALFA)*(BETA + ALFA);     
         SUM:= ALFA + BETA + 2; B[0]:= 0;       
         A[0]:= (BETA - ALFA)/SUM;    
         A[1]:= MIN/SUM/(SUM + 2);    
         B[1]:= 4*(1 + ALFA)*(1 + BETA)/SUM/SUM/(SUM + 1);
         "FOR" I:= 2 "STEP" 1 "UNTIL" N - 1 "DO"
         "BEGIN" SUM:= I + I + ALFA + BETA;     
             A[I]:= MIN/SUM/(SUM + 2); SUM:= SUM*SUM;     
             B[I]:= 4*I*(I + ALFA + BETA)*(I + ALFA)*(I + BETA)/    
             (SUM - 1)/SUM  
         "END";   
         EM[0]:=ARREB; EM[2]:= 1.0"-8; EM[4]:= 6*N;       
         ALL ZER ORT POL(N, A, B, ZER, EM)      
     "END" ALL JAC ZER;     
         "EOP"    
        
        
  "CODE" 31371;   
     "PROCEDURE"  ALL LAG ZER(N, ALFA, ZER);    
     "VALUE" N, ALFA     ; "INTEGER"N; "REAL" ALFA     ; "ARRAY" ZER; 
     "BEGIN" "INTEGER" I; "ARRAY" A, B[0:N], EM[0:5];     
         "REAL""PROCEDURE" ARREB;     
         "CODE" 30002;      
         "PROCEDURE" ALL ZER ORT POL(N, B, C, ZER, EM);   
         "VALUE" N; "INTEGER" N; "ARRAY" B,C,ZER,EM;      
         "CODE" 31362;      
         B[0]:= 0; A[N - 1]:= N + N + ALFA - 1; 
         "FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
         "BEGIN" A[I - 1]:= I + I + ALFA - 1;   
             B[I]:= I*(I + ALFA)      
         "END";   
        EM[0]:=ARREB; EM[2]:= "-10;EM[4]:= 6*N; 
         ALL ZER ORT POL(N, A, B, ZER, EM)      
     "END" ALL LAG ZER;     
         "EOP"    
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 1      
        
        
        
 AUTHOR: J.C.P.BUS.         
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 740218.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS PROCEDURES FOR CALCULATING THE DERIVATIVES OF      
     FUNCTIONS OF MORE VARIABLES, USING DIFFERENCE FORMULAS;        
     JACOBNNF   CALCULATES  THE  JACOBIAN  MATRIX  OF  AN  N-DIMENSIONAL      
         FUNCTION OF N VARIABLES USING FORWARD DIFFERENCES;         
     JACOBNMF   CALCULATES  THE  JACOBIAN  MATRIX  OF  AN  N-DIMENSIONAL      
         FUNCTION OF M VARIABLES USING FORWARD DIFFERENCES;         
     JACOBNBNDF  CALCULATES  THE  JACOBIAN  MATRIX  OF AN  N-DIMENSIONAL      
         FUNCTION OF  N  VARIABLES,  IF THIS JACOBIAN  IS KNOWN TO BE  A      
         BAND MATRIX AND HAVE TO BE STORED ROW-WISE IN A ONE-DIMENSIONAL      
         ARRAY.   
        
        
 KEYWORDS:        
        
     NUMERICAL DIFFERENTIATION,       
     FUNCTIONS OF MORE VARIABLES,     
     DIFFERENCE FORMULAS.   
        
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 2      
        
        
        
 SUBSECTION: JACOBNNF.      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS:       
     "PROCEDURE" JACOBNNF(N, X, F, JAC, I, DI, FUNCT);    
     "VALUE" N; "INTEGER" I, N; "REAL" DI; "ARRAY" X, F, JAC;       
     "PROCEDURE" FUNCT;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER  OF INDEPENDENT VARIABLES  AND  THE DIMENSION OF      
             THE FUNCTION;  
     X:      <ARRAY IDENTIFIER>;      
             "ARRAY" X[1:N];
             ENTRY: THE POINT AT WHICH THE JACOBIAN HAS TO BE CALCULATED      
     F:      <ARRAY IDENTIFIER>;      
             "ARRAY" F[1:N];
             ENTRY : THE VALUES OF THE FUNCTION-COMPONENTS AT THE POINT       
             GIVEN IN ARRAY X;        
     JAC:    <ARRAY IDENTIFIER>;      
             "ARRAY" JAC[1:N, 1:N];   
             EXIT : THE JACOBIAN MATRIX IN SUCH A WAY THAT THE PARTIAL        
             DERIVATIVE OF  F[I] TO  X[J] IS  GIVEN IN    
             JAC[I, J], I, J = 1, ..., N;       
     I:      <INTEGER VARIABLE>;      
             A JENSEN PARAMETER; DI MAY BE DEPENDENT OF I;
     DI:     <ARITHMETIC EXPRESSION>; 
             THE  PARTIAL DERIVATIVES  TO  X[I]  ARE  APPROXIMATED  WITH      
             FORWARD  DIFFERENCES,   USING  AN  INCREMENT  TO  THE  I-TH      
             VARIABLE THAT EQUALS THE VALUE OF DI, I = 1, ..., N;   
     FUNCT:  <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE SHOULD READ:   
             "PROCEDURE" FUNCT(N, X, F);        
             "VALUE" N; "INTEGER" N; "ARRAY" X, F;        
             THE MEANING OF THE FORMAL PARAMETERS IS:     
             N:  <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF INDEPENDENT VARIABLES OF THE FUNCTION  F;      
             X:  <ARRAY IDENTIFIER>;  
                 THE INDEPENDENT VARIABLES ARE GIVEN IN X[1:N];     
             F:  <ARRAY IDENTIFIER>;  
                 AFTER A CALL OF FUNCT THE FUNCTION COMPONENTS SHOULD BE      
                 GIVEN IN F[1:N].     
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 3      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 REQUIRED CENTRAL MEMORY :  
        
     EXECUTION FIELD LENGTH: JACOBNNF DECLARES ONE  AUXILIARY  ARRAY  OF      
                             ORDER N. 
        
        
 RUNNING TIME: PROPORTIONAL TO  N ** 2.         
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     JACOBNNF   CALCULATES  THE  JACOBIAN  MATRIX  OF  AN  N-DIMENSIONAL      
     FUNCTION OF N VARIABLES; THE ELEMENTS OF THIS MATRIX, WHICH ARE THE      
     PARTIAL DERIVATIVES OF THE FUNCTION,  ARE CALCULATED  USING FORWARD      
     DIFFERENCES WITH AN INCREMENT TO THE I-TH  VARIABLE OF  DI, (I = 1,      
     ..., N).     
        
        
 EXAMPLE OF USE:  
        
     LET F BE DEFINED BY:   
             F[1]= X[1] ** 3 + X[2],  
             F[2]= 10 * X[2];         
     THE JACOBIAN MATRIX  AT  THE POINT  (2, 1)  MAY  BE CALCULATED  AND      
     PRINTED BY THE FOLLOWING PROGRAM:
        
     "BEGIN"      
         "INTEGER" I; "ARRAY" JAC[1:2, 1:2], X, F[1:2];   
         "PROCEDURE" JACOBNNF(N, X, F, JAC, I, DI, FU); "CODE" 34437;         
         "PROCEDURE" F1(N, X, F); "VALUE" N; "INTEGER" N; "ARRAY" X, F;       
         "BEGIN" F[1]:= X[1] ** 3 + X[2]; F[2]:= X[2] * 10 "END" F1;
         X[1]:= 2; X[2]:= 1; F1(2, X, F);       
         JACOBNNF(2, X, F, JAC, I, "IF" I = 1 "THEN" "-6 "ELSE" 1, F1);       
         OUTPUT(71, "("*,4B,"("THE CALCULATED JACOBIAN IS:")",//,   
         2(4B,2(N),/)")", JAC[1, 1], JAC[1, 2], JAC[2, 1], JAC[2, 2])         
     "END"        
        
     RESULTS:     
        
     THE CALCULATED JACOBIAN IS:      
        
     +1.2000005938262"+001  +1.0000000000000"+000         
     +0.0000000000000"+000  +1.0000000000000"+001         
        
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 4      
        
        
        
 SUBSECTION: JACOBNMF.      
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS:       
     "PROCEDURE" JACOBNMF(N, M, X, F, JAC, I, DI, FUNCT); 
     "VALUE" N, M; "INTEGER" I, N, M; "REAL" DI; "ARRAY" X, F, JAC; 
     "PROCEDURE" FUNCT;     
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF FUNCTION COMPONENTS; 
     M:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF INDEPENDENT VARIABLES;         
     X:      <ARRAY IDENTIFIER>;      
             "ARRAY" X[1:M];
             ENTRY: THE POINT AT WHICH THE JACOBIAN HAS TO BE CALCULATED      
     F:      <ARRAY IDENTIFIER>;      
             "ARRAY" F[1:N];
             ENTRY: THE VALUES OF THE FUNCTION-COMPONENTS AT THE POINT        
             GIVEN IN ARRAY X;        
     JAC:    <ARRAY IDENTIFIER>;      
             "ARRAY" JAC[1:N, 1:M];   
             EXIT : THE JACOBIAN MATRIX IN SUCH A WAY THAT THE PARTIAL        
             DERIVATIVE OF  F[I] TO  X[J] IS  GIVEN IN    
             JAC[I, J], I = 1, ..., N, J = 1, ... M;      
     I:      <INTEGER VARIABLE>;      
             A JENSEN PARAMETER; DI MAY BE DEPENDENT OF I;
     DI:     <ARITHMETIC EXPRESSION>; 
             THE  PARTIAL DERIVATIVES  TO  X[I]  ARE  APPROXIMATED  WITH      
             FORWARD  DIFFERENCES,   USING  AN  INCREMENT  TO  THE  I-TH      
             VARIABLE THAT EQUALS THE VALUE OF DI, I = 1, ..., M;   
     FUNCT:  <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS :        
             "PROCEDURE" FUNCT(N, M, X, F);     
             "VALUE" N, M; "INTEGER" N, M; "ARRAY" X, F;  
             THE MEANING OF THE FORMAL PARAMETERS IS :    
             N:  <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF FUNCTION COMPONENTS;       
             M:  <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF INDEPENDENT VARIABLES OF THE FUNCTION  F;      
             X:  <ARRAY IDENTIFIER>;  
                 THE INDEPENDENT VARIABLES ARE GIVEN IN X[1:M];     
             F:  <ARRAY IDENTIFIER>;  
                 AFTER A CALL OF FUNCT THE FUNCTION COMPONENTS SHOULD BE      
                 GIVEN IN F[1:N].     
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 5      
        
        
        
 PROCEDURES USED: NONE.     
        
        
 EXECUTION FIELD LENGTH: JACOBNMF   DECLARES   ONE  AUXILIARY  ARRAY  OF      
                         ORDER N.     
        
        
 RUNNING TIME: PROPORTIONAL TO  N * M.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     JACOBNMF   CALCULATES  THE  JACOBIAN  MATRIX  OF  AN  N-DIMENSIONAL      
     FUNCTION OF M VARIABLES; THE ELEMENTS OF THIS MATRIX, WHICH ARE THE      
     PARTIAL DERIVATIVES OF THE FUNCTION,  ARE CALCULATED  USING FORWARD      
     DIFFERENCES  WITH AN INCREMENT TO THE I-TH  VARIABLE OF  DI,(I = 1,      
     ..., M).     
        
        
 EXAMPLE OF USE:  
        
     LET F BE DEFINED BY:   
             F[1]= X[1] ** 3 + X[2],  
             F[2]= 10 * X[2] + X[2] * X[1],     
             F[3]= X[1] * X[2];       
     THE  JACOBIAN MATRIX  AT  THE POINT  (2, 1)  MAY BE CALCULATED  AND      
     PRINTED BY THE FOLLOWING PROGRAM:
        
     "BEGIN"      
         "INTEGER" I; "ARRAY" JAC[1:3, 1:2], X[1:2], F[1:3];        
         "PROCEDURE" JACOBNMF(N, X, F, JAC, I, DI, FU); "CODE" 34438;         
         "PROCEDURE" F1(N, M, X, F); "VALUE" N, M; "INTEGER" N, M;  
         "ARRAY" X, F;      
         "BEGIN" F[1]:= X[1] ** 3 + X[2];       
             F[2]:= X[2] * 10 + X[2] * X[1] ** 2; F[3]:= X[1] * X[2]
         "END" F1;
         X[1]:= 2; X[2]:= 1; F1(3, 2, X, F);    
         JACOBNMF(3, 2, X, F, JAC, I, "IF" I=2 "THEN" 1 "ELSE" "-5, F1);      
         OUTPUT(71, "("*,4B,"("THE CALCULATED JACOBIAN IS:")",//,   
         3(4B,2(N),/)")", JAC[1, 1], JAC[1, 2], JAC[2, 1], JAC[2, 2],         
         JAC[3, 1], JAC[3, 2])        
     "END"        
        
     RESULTS:     
        
     THE CALCULATED JACOBIAN IS:      
        
     +1.2000060002038"+001  +1.0000000000000"+000         
     +4.0000100000270"+000  +1.4000000000000"+001         
     +1.0000000003174"+000  +2.0000000000000"+000         
        
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 6      
        
        
        
 SUBSECTION: JACOBNBNDF.    
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS :      
     "PROCEDURE" JACOBNBNDF(N, LW, RW, X, F, JAC, I, DI, FUNCT);    
     "VALUE" N, LW, RW; "INTEGER" N, I, LW, RW; "REAL" DI;
     "ARRAY" X, F, JAC; "PROCEDURE" FUNCT;      
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER  OF INDEPENDENT VARIABLES  AND  THE DIMENSION OF      
             THE FUNCTION;  
     LW:     <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF CODIAGONALS  TO THE LEFT OF THE MAIN DIAGONAL      
             OF THE JACOBIAN MATRIX, WHICH IS KNOWN TO BE A BAND MATRIX;      
     RW:     <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF CODIAGONALS TO THE RIGHT OF THE MAIN DIAGONAL      
             OF THE JACOBIAN MATRIX;  
     X:      <ARRAY IDENTIFIER>;      
             "ARRAY" X[1:N];
             ENTRY: THE POINT AT WHICH THE JACOBIAN HAS TO BE CALCULATED      
     F:      <ARRAY IDENTIFIER>;      
             "ARRAY" F[1:N];
             ENTRY: THE VALUES OF THE FUNCTION-COMPONENTS AT THE POINT        
             GIVEN IN ARRAY X;        
     JAC:    <ARRAY IDENTIFIER>;      
             "ARRAY" JAC [1 : (LW + RW ) * (N - 1) + N];  
             EXIT: THE JACOBIAN MATRIX  IN SUCH A WAY THAT THE (I, J)-TH      
             ELEMENT OF THE JACOBIAN, I.E. THE PARTIAL DERIVATIVE OF
             F[I] TO X[J], IS  GIVEN IN         
             JAC[(LW + RW) * (I - 1) + J], FOR I = 1, ..., N        
             J= MAX(1, I - LW), ..., MIN(N, I + RW);      
     I:      <INTEGER VARIABLE>;      
             A JENSEN PARAMETER; DI MAY BE DEPENDENT OF I;
     DI:     <ARITHMETIC EXPRESSION>; 
             THE  PARTIAL DERIVATIVES  TO  X[I]  ARE  APPROXIMATED  WITH      
             FORWARD  DIFFERENCES,   USING  AN  INCREMENT  TO  THE  I-TH      
             VARIABLE THAT EQUALS THE VALUE OF DI, I = 1, ..., N;   
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 7      
        
        
        
     FUNCT:  <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS :        
             "PROCEDURE" FUNCT(N, L, U, X, F);  
             "VALUE" N, L, U; "INTEGER" N, L, U; "ARRAY" X, F;      
             THE MEANING OF THE FORMAL PARAMETERS IS :    
             N:  <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF FUNCTION COMPONENTS;       
             L,U:<ARITHMETIC EXPRESSION>;       
                 LOWER  AND  UPPER  BOUND   OF  THE  FUNCTION  COMPONENT      
                 SUBSCRIPT; 
             X:  <ARRAY IDENTIFIER>;  
                 THE INDEPENDENT VARIABLES ARE GIVEN IN X[1:N];     
             F:  <ARRAY IDENTIFIER>;  
                 AFTER  A CALL  OF  FUNCT  THE FUNCTION COMPONENTS F[I],      
                 I = L, ..., U, SHOULD BE GIVEN IN F[L:U].
        
        
 PROCEDURES USED: NONE.     
        
        
 EXECUTION FIELD LENGTH: JACOBNMF   DECLARES   ONE  AUXILIARY  ARRAY  OF      
                         MAXIMUM ORDER LW + RW + 1;       
        
        
 RUNNING TIME: PROPORTIONAL TO  N * (LW + RW + 1).        
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     JACOBNBNDF   CALCULATES  THE  JACOBIAN MATRIX  OF AN  N-DIMENSIONAL      
     FUNCTION OF  N  VARIABLES,  IF THIS JACOBIAN IS KNOWN TO BE  A BAND      
     MATRIX AND HAVE TO BE STORED  ROW-WISE  IN A ONE-DIMENSIONAL ARRAY;      
     THE ELEMENTS OF THIS JACOBIAN MATRIX ARE CALCULATED  USING  FORWARD      
     DIFFERENCES, WITH AN INCREMENT TO THE I-TH VARIABLE OF DI,  (I = 1,      
     ..., N).     
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 8      
        
        
        
 EXAMPLE OF USE:  
        
     LET F BE DEFINED BY:   
     F[1] = (3 - 2 * X[1]) * X[1] + 1 - 2 * X[2],         
     F[I] = (3 - 2 * X[I]) * X[I] + 1 - X[I-1] - 2 * X[I+1], I= 2, 3, 4,      
     F[5] = 4 - 2 * X[5] - X[4];      
     THE TRIDIAGONAL JACOBIAN MATRIX AT THE POINT X, GIVEN BY X[I] = -1,      
     I = 1, ..., 5,  MAY  BE  CALCULATED  AND  PRINTED  BY THE FOLLOWING      
     PROGRAM:     
        
     "BEGIN"      
         "INTEGER" I; "ARRAY" JAC[1:13], X, F[1:5];       
         "PROCEDURE" JACOBNBNDF(N, L, R, X, F, J, I, D,G); "CODE" 34439;      
         "PROCEDURE" F1(N, L, U, X, F); "VALUE" N, L, U;  
         "INTEGER" N, L, U; "ARRAY" X, F;       
         "BEGIN" "INTEGER" I;         
             "FOR" I:= L "STEP" 1 "UNTIL" ("IF" U = 5 "THEN" 4 "ELSE" U)      
             "DO" 
             "BEGIN" F[I]:= (3 - 2 * X[I]) * X[I] + 1 - 2 * X[I + 1];         
                 "IF" I ^= 1 "THEN" F[I]:= F[I] - X[I - 1]
             "END";         
             "IF" U = 5 "THEN" F[5]:= 4 - X[4] - X[5] * 2 
         "END" F1;
        
         "PROCEDURE" LIST(ITEM); "PROCEDURE" ITEM;        
         "BEGIN" "INTEGER" I;         
             ITEM("("THE CALCULATED TRIDIAGONAL JACOBIAN IS:")");   
             "FOR" I:= 1 "STEP" 1 "UNTIL" 13 "DO" ITEM(JAC[I])      
          "END" LIST;       
        
         "PROCEDURE" LAYOUT;
         FORMAT("("//,4B,40S,//,4B,2(+.5D"+D2B),/,4B,3(+.5D"+D2B),/,
         16B,3(+.5D"+D2B),/,28B,3(+.5D"+D2B),/,40B,2(+.5D"+D2B),/")");        
         "FOR" I := 1 "STEP" 1 "UNTIL" 5 "DO" X[I]:= -1;  
         F1(5, 1, 5, X, F); 
         JACOBNBNDF(5, 1, 1, X, F, JAC, I, "IF" I = 5 "THEN" 1 "ELSE"         
         "-6, F1);
         OUTLIST(71, LAYOUT, LIST)    
     "END"        
        
     RESULTS:     
        
     THE CALCULATED TRIDIAGONAL JACOBIAN IS:    
        
     +.70000"+1  -.20000"+1 
     -.10000"+1  +.70000"+1  -.20000"+1         
                 -.10000"+1  +.70000"+1  -.20000"+1       
                             -.10000"+1  +.70000"+1  -.20000"+1     
                                         -.10000"+1  -.20000"+1     
        
        
        
1SECTION : 4.3.2.1            (OCTOBER 1974)                      PAGE 9      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 34437;    
 "PROCEDURE" JACOBNNF(N, X, F, JAC, I, DI, FUNCT); "VALUE" N;       
 "INTEGER" N, I; "REAL" DI; "ARRAY" X, F, JAC; "PROCEDURE" FUNCT;   
 "BEGIN" "INTEGER" J; "REAL" STEP, AID; "ARRAY" F1[1:N];  
     "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"        
     "BEGIN" STEP:= DI; AID:= X[I]; X[I]:= AID + STEP;    
         STEP:= 1 / STEP; FUNCT(N, X, F1);      
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         JAC[J,I]:= (F1[J] - F[J]) * STEP; X[I]:= AID     
     "END"        
 "END" JACOBNNF;  
         "EOP"    
        
 "CODE" 34438;    
 "PROCEDURE" JACOBNMF(N, M, X, F, JAC, I, DI, FUNCT); "VALUE" N, M; 
 "INTEGER" N, M, I; "REAL" DI; "ARRAY" X, F, JAC; "PROCEDURE" FUNCT;
 "BEGIN" "INTEGER" J; "REAL" STEP, AID; "ARRAY" F1[1:N];  
     "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"        
     "BEGIN" STEP:= DI; AID:= X[I]; X[I]:= AID + STEP;    
         STEP:= 1 / STEP; FUNCT(N, M, X, F1);   
         "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"    
         JAC[J,I]:= (F1[J] - F[J]) * STEP; X[I]:= AID     
     "END"        
 "END" JACOBNMF;  
         "EOP"    
        
 "CODE" 34439;    
     "PROCEDURE" JACOBNBNDF(N, LW, RW, X, F, JAC, I, DI, FUNCT);    
     "VALUE" N, LW, RW; "INTEGER" I, N, LW, RW; "REAL" DI;
     "ARRAY" X, F, JAC; "PROCEDURE" FUNCT;      
     "BEGIN" "INTEGER" J, K, L, U, T, B; "REAL" AID, STEPI;         
         L:= 1; U:= LW + 1; T:= RW + 1; B:= LW + RW;      
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" "ARRAY" F1[L:U];     
             STEPI:= DI; AID:= X[I]; X[I]:= AID + STEPI;  
             FUNCT(N, L, U, X, F1); X[I]:= AID; 
             K:= I + ("IF" I <= T "THEN" 0 "ELSE" I - T) * B;       
             "FOR" J:= L "STEP" 1 "UNTIL" U "DO"
             "BEGIN" JAC[K]:= (F1[J] - F[J]) / STEPI; K:=K + B "END";         
             "IF" I >= T "THEN" L:= L + 1;      
             "IF" U < N "THEN" U:= U + 1        
         "END"    
     "END" JACOBNBNDF;      
         "EOP"    
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 1      
        
        
        
 AUTHORS: J. C. P. BUS AND  T. J. DEKKER.       
        
        
 CONTRIBUTOR: J. C. P. BUS. 
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 750615.
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS TWO PROCEDURES  FOR FINDING A ZERO OF A GIVEN      
     FUNCTION IN A GIVEN INTERVAL;    
     ZEROIN    APPROXIMATES A ZERO  MAINLY  BY LINEAR INTERPOLATION  AND      
               EXTRAPOLATION;         
     ZEROINRAT APPROXIMATES  A  ZERO  BY  INTERPOLATION   WITH  RATIONAL      
               FUNCTIONS.   
     ZEROIN IS PREFERABLE FOR SIMPLE (I.E. CHEAPLY TO CALCULATE)    
     FUNCTIONS AND/OR WHEN NO HIGH PRECISION IS REQUIRED. ZEROINRAT IS        
     PREFERABLE FOR COMPLICATED (I.E. EXPENSIVE) FUNCTIONS WHEN A ZERO        
     IS REQUIRED IN RATHER HIGH PRECISION AND ALSO FOR FUNCTIONS HAVING       
     A POLE NEAR THE ZERO. WHEN THE ANALYTIC DERIVATIVE OF THE FUNCTION       
     IS EASILY OBTAINED, THEN ZEROINDER (SECTION 5.1.1.1.2) SHOULD BE         
     TAKEN INTO CONSIDERATION.        
        
        
 KEYWORDS:        
        
     ZERO SEARCHING,        
     ANALYTIC EQUATIONS,    
     SINGLE NON-LINEAR EQUATION.      
        
        
        
        
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 2      
        
        
        
 SUBSECTION: ZEROIN.        
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS :      
     "BOOLEAN" "PROCEDURE" ZEROIN(X, Y, FX, TOLX);        
     "REAL" X, Y, FX, TOLX; 
        
     ZEROIN  SEARCHES FOR A ZERO OF A REAL FUNCTION  F  DEFINED ON A
             CERTAIN INTERVAL J;      
             ZEROIN := "TRUE" WHEN A (SUFFICIENTLY SMALL) SUBINTERVAL OF      
             J CONTAINING A  ZERO OF  F HAS  BEEN  FOUND;  OTHERWISE,         
             ZEROIN := "FALSE";       
             LET A REAL FUNCTION T DEFINED ON J, DENOTE A TOLERANCE 
             FUNCTION DEFINING THE REQUIRED PRECISION OF THE ZERO;  
             (FOR INSTANCE  
                           T(X) = ABS(X) * RE + AE,       
             WHERE  RE  AND  AE  ARE THE REQUIRED  RELATIVE AND ABSOLUTE      
             PRECISION RESPECTIVELY); 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:      <REAL VARIABLE>;         
             A JENSEN VARIABLE;  THE ACTUAL PARAMETERS FOR  FX AND TOLX       
             (MAY) DEPEND ON THE ACTUAL PARAMETER FOR X;  
             ENTRY:  ONE  ENDPOINT OF INTERVAL  J  IN  WHICH  A  ZERO IS      
                     SEARCHED FOR;    
             EXIT:   A VALUE APPROXIMATING THE ZERO WITHIN THE TOLERANCE      
                     2 * T(X) WHEN ZEROIN HAS THE VALUE "TRUE", AND A         
                     PRESUMABLY WORTHLESS ARGUMENT VALUE OTHERWISE; 
     Y:      <REAL VARIABLE>;         
             ENTRY:  THE OTHER ENDPOINT OF INTERVAL J IN WHICH A ZERO IS      
                     SEARCHED FOR;  UPON ENTRY X < Y AS WELL AS Y < X IS      
                     ALLOWED;         
             EXIT:   THE  OTHER  STRADDLING  APPROXIMATION  OF THE ZERO,      
                     I.E.  UPON  EXIT  THE  VALUES  OF  Y AND X  SATISFY      
                     1. F(X) * F(Y) <= 0,  2. ABS(X - Y) <= 2 * T(X) AND      
                     3. ABS(F(X)) <= ABS(F(Y)) WHEN ZEROIN HAS THE  
                     VALUE "TRUE", AND  A  PRESUMABLY WORTHLESS ARGUMENT      
                     VALUE  SATISFYING  CONDITIONS  2  AND  3  BUT NOT 1      
                     OTHERWISE;       
     FX:     <ARITHMETIC EXPRESSION>; 
             DEFINES  FUNCTION F AS A FUNCTION  DEPENDING ON THE  ACTUAL      
             PARAMETER FOR X;         
     TOLX:   <ARITHMETIC EXPRESSION>; 
             DEFINES THE  TOLERANCE  FUNCTION T WHICH MAY  DEPEND ON THE      
             ACTUAL PARAMETER FOR X;  
             ONE SHOULD CHOOSE TOLX  POSITIVE AND NEVER SMALLER THAN THE      
             PRECISION OF THE  MACHINE'S  ARITHMETIC  AT X, I.E. IN THIS      
             ARITHMETIC  X + TOLX  AND  X - TOLX   SHOULD  ALWAYS  YIELD      
             VALUES  DISTINCT  FROM  X; OTHERWISE  THE PROCEDURE MAY GET      
             INTO A LOOP.   
        
        
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 3      
        
        
        
 PROCEDURES USED: DWARF = CP30003;    
        
        
 EXECUTION FIELD LENGTH: NO AUXILIARY ARRAYS ARE DECLARED IN ZEROIN.
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IS DESCRIBED IN DETAIL IN [1].       
     THE NUMBER OF EVALUATIONS OF FX AND TOLX IS AT MOST  
             4 * LOG(ABS(X - Y)) / TAU,         
     WHERE X AND Y ARE THE ARGUMENT VALUES GIVEN UPON ENTRY, LOG DENOTES      
     THE BASE 2 LOGARITHM AND TAU IS THE MINIMUM OF THE TOLERANCE   
     FUNCTION TOLX ON THE INITIAL INTERVAL. IF UPON ENTRY X AND Y   
     SATISFY F(X) * F(Y) <= 0, THEN CONVERGENCE IS GUARANTEED AND THE         
     ASYMPTOTIC ORDER OF CONVERGENCE IS 1.618 FOR SIMPLE ZEROES.    
        
        
 EXAMPLE OF USE:  
        
     THE ZERO OF THE FUNCTION   EXP(-X * 3) * (X - 1) + X ** 3,  IN  THE      
     INTERVAL [0, 1], MAY BE CALCULATED BY THE FOLLOWING PROGRAM:   
        
     "BEGIN" "REAL" X, Y;   
         "BOOLEAN" "PROCEDURE" ZEROIN(X, Y, FX, TOLX); "CODE" 34150;
         "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;    
         F:= EXP(-X * 3) * ( X - 1) + X ** 3;   
         X:= 0; Y:= 1;      
         "IF" ZEROIN(X, Y, F(X), ABS(X) * "-14 + "-14) "THEN"       
         OUTPUT(71, "("/4B,"("CALCULATED ZERO:")"B+.15D"+3D")", X)  
         "ELSE" OUTPUT(71, "("/4B,"("NO ZERO FOUND")"")") 
     "END"        
        
     RESULT:      
        
     CALCULATED ZERO: +.489702748548240"+000    
        
        
        
        
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 4      
        
        
        
 SUBSECTION: ZEROINRAT.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS:       
     "BOOLEAN" "PROCEDURE" ZEROINRAT(X, Y, FX, TOLX);     
     "REAL" X, Y, FX, TOLX; 
        
     ZEROINRAT SEARCHES FOR A ZERO OF A REAL FUNCTION  F DEFINED ON A         
             CERTAIN INTERVAL J;      
             ZEROINRAT := "TRUE" WHEN A (SUFFICIENTLY SMALL) SUBINTERVAL      
             OF J CONTAINING A  ZERO OF  F HAS  BEEN  FOUND;  OTHERWISE,      
             ZEROINRAT := "FALSE";    
             LET A REAL FUNCTION T DEFINED ON J, DENOTE A TOLERANCE 
             FUNCTION DEFINING THE REQUIRED PRECISION OF THE ZERO;  
             (FOR INSTANCE  
                           T(X) = ABS(X) * RE + AE,       
             WHERE  RE  AND  AE  ARE THE REQUIRED  RELATIVE AND ABSOLUTE      
             PRECISION RESPECTIVELY); 
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     X:      <REAL VARIABLE>;         
             A JENSEN VARIABLE;  THE ACTUAL PARAMETERS FOR  FX AND TOLX       
             (MAY) DEPEND ON THE ACTUAL PARAMETER FOR X;  
             ENTRY:  ONE  ENDPOINT OF INTERVAL  J  IN  WHICH  A  ZERO IS      
                     SEARCHED FOR;    
             EXIT:   A VALUE APPROXIMATING THE ZERO WITHIN THE TOLERANCE      
                     2 * T(X) WHEN ZEROINRAT HAS THE VALUE "TRUE", AND A      
                     PRESUMABLY WORTHLESS ARGUMENT VALUE OTHERWISE; 
     Y:      <REAL VARIABLE>;         
             ENTRY:  THE OTHER ENDPOINT OF INTERVAL J IN WHICH A ZERO IS      
                     SEARCHED FOR;  UPON ENTRY X < Y AS WELL AS Y < X IS      
                     ALLOWED;         
             EXIT:   THE  OTHER  STRADDLING  APPROXIMATION  OF THE ZERO,      
                     I.E.  UPON  EXIT  THE  VALUES  OF  Y AND X  SATISFY      
                     1. F(X) * F(Y) <= 0,  2. ABS(X - Y) <= 2 * T(X) AND      
                     3. ABS(F(X)) <= ABS(F(Y))  WHEN  ZEROINRAT  HAS THE      
                     VALUE "TRUE", AND  A  PRESUMABLY WORTHLESS ARGUMENT      
                     VALUE  SATISFYING  CONDITIONS  2  AND  3  BUT NOT 1      
                     OTHERWISE;       
     FX:     <ARITHMETIC EXPRESSION>; 
             DEFINES  FUNCTION F AS A FUNCTION  DEPENDING ON THE  ACTUAL      
             PARAMETER FOR X;         
     TOLX:   <ARITHMETIC EXPRESSION>; 
             DEFINES THE  TOLERANCE  FUNCTION T WHICH MAY  DEPEND ON THE      
             ACTUAL PARAMETER FOR X;  
             ONE SHOULD CHOOSE TOLX  POSITIVE AND NEVER SMALLER THAN THE      
             PRECISION OF THE  MACHINE'S  ARITHMETIC  AT X, I.E. IN THIS      
             ARITHMETIC  X + TOLX  AND  X - TOLX   SHOULD  ALWAYS  YIELD      
             VALUES  DISTINCT  FROM  X; OTHERWISE  THE PROCEDURE MAY GET      
             INTO A LOOP.   
        
        
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 5      
        
        
        
 PROCEDURES USED: DWARF = CP30003;    
        
        
 EXECUTION FIELD LENGTH: NO AUXILIARY ARRAYS ARE DECLARED IN ZEROINRAT.       
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IS DESCRIBED IN DETAIL IN [1].       
     THE NUMBER OF EVALUATIONS OF FX AND TOLX IS AT MOST  
             5 * LOG(ABS(X - Y)) / TAU,         
     WHERE X AND Y ARE THE ARGUMENT VALUES GIVEN UPON ENTRY, LOG DENOTES      
     THE BASE 2 LOGARITHM AND TAU IS THE MINIMUM OF THE TOLERANCE   
     FUNCTION TOLX ON THE INITIAL INTERVAL. IF UPON ENTRY X AND Y   
     SATISFY F(X) * F(Y) <= 0, THEN CONVERGENCE IS GUARANTEED AND THE         
     ASYMPTOTIC ORDER OF CONVERGENCE IS 1.839 FOR SIMPLE ZEROES.    
        
        
 EXAMPLE OF USE:  
        
     THE ZERO OF THE FUNCTION   EXP(-X * 3) * (X - 1) + X ** 3,  IN  THE      
     INTERVAL [0, 1], MAY BE CALCULATED BY THE FOLLOWING PROGRAM:   
        
     "BEGIN" "REAL" X, Y;   
         "BOOLEAN" "PROCEDURE" ZEROINRAT(X, Y, FX, TOLX); "CODE" 34436;       
         "REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;    
         F:= EXP(-X * 3) * ( X - 1) + X ** 3;   
         X:= 0; Y:= 1;      
         "IF" ZEROINRAT(X, Y, F(X), ABS(X) * "-14 + "-14) "THEN"    
         OUTPUT(71, "("/4B,"("CALCULATED ZERO:")"B+.15D"+3D")", X)  
         "ELSE" OUTPUT(71, "("/4B,"("NO ZERO FOUND")"")") 
     "END"        
        
     RESULT:      
        
     CALCULATED ZERO: +.489702748548240"+000    
        
        
        
 REFERENCES:      
        
     [1]: BUS, J.C.P. AND  DEKKER, T.J.,        
         TWO EFFICIENT ALGORITHMS WITH GUARANTEED CONVERGENCE FOR   
         FINDING A ZERO OF A FUNCTION.
         MATHEMATICAL CENTRE, REPORT NW 13/74, AMSTERDAM (1974),    
         (ALSO TO APPEAR IN TOMS 1975).         
        
        
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 6      
        
        
        
 SOURCE TEXT(S):  
0"CODE" 34150;    
 "BOOLEAN" "PROCEDURE" ZEROIN(X, Y, FX, TOLX);  
 "REAL" X, Y, FX, TOLX;     
 "BEGIN" "INTEGER" EXT;     
         "REAL" C, FC, B, FB, A, FA, D, FD, FDB, FDA, W, MB,        
         TOL, M, P, Q, DW;  
         DW:= DWARF; B:= X; FB:= FX; A:= X:= Y; FA:= FX;  
     INTERPOLATE: C:= A; FC:= FA; EXT:= 0;      
     EXTRAPOLATE: "IF" ABS(FC) < ABS(FB) "THEN" 
         "BEGIN" "IF" C ^= A "THEN" "BEGIN" D:= A; FD:= FA "END";   
             A:= B; FA:= FB; B:= X:= C; FB:= FC; C:= A; FC:= FA     
         "END" INTERCHANGE; 
         TOL:= TOLX; M:= (C + B) * 0.5; MB:= M - B;       
         "IF" ABS(MB) > TOL "THEN"    
         "BEGIN" "IF" EXT > 2 "THEN" W:= MB "ELSE"        
             "BEGIN" TOL:= TOL * SIGN(MB);      
                 P:= (B - A) * FB; "IF" EXT <= 1 "THEN"   
                 Q:= FA - FB "ELSE"   
                 "BEGIN" FDB:= (FD - FB) / (D - B);       
                     FDA:= (FD - FA) / (D - A); 
                     P:= FDA * P; Q:= FDB * FA - FDA * FB 
                 "END"; "IF" P < 0 "THEN"       
                 "BEGIN" P:= -P; Q:= -Q "END";  
                 W:= "IF" P < DW "OR" P <= Q * TOL "THEN" TOL "ELSE"
                 "IF" P < MB * Q "THEN" P / Q "ELSE" MB   
             "END"; D:= A; FD:= FA; A:= B; FA:= FB;       
             X:= B:= B + W; FB:= FX;  
             "IF" ("IF" FC >= 0 "THEN" FB >= 0 "ELSE" FB <= 0) "THEN"         
             "GOTO" INTERPOLATE "ELSE"
             "BEGIN" EXT:= "IF" W = MB "THEN" 0 "ELSE" EXT + 1;     
                 "GOTO" EXTRAPOLATE   
             "END"
         "END"; Y:= C;      
         ZEROIN:= "IF" FC >= 0 "THEN" FB <= 0 "ELSE" FB >= 0        
 "END" ZEROIN;    
         "EOP"    
1SECTION : 5.1.1.1.1         (OCTOBER 1975)                      PAGE 7      
        
        
        
 "CODE" 34436;    
 "BOOLEAN" "PROCEDURE" ZEROINRAT(X, Y, FX, TOLX);         
 "REAL" X, Y, FX, TOLX;     
 "BEGIN" "INTEGER" EXT; "BOOLEAN" FIRST;        
         "REAL" B, FB, A, FA, D, FD, C, FC, FDB, FDA, W,  
         MB, TOL, M, P, Q, DW;        
         "REAL" "PROCEDURE" DWARF; "CODE" 30003;
         DW:= DWARF; B:= X; FB:= FX; A:= X:= Y; FA:= FX; FIRST:= "TRUE";      
     INTERPOLATE: C:= A; FC:= FA; EXT:= 0;      
     EXTRAPOLATE: "IF" ABS(FC) < ABS(FB) "THEN" 
         "BEGIN" "IF" C ^= A "THEN" "BEGIN" D:= A; FD:= FA "END";   
             A:= B; FA:= FB; B:= X:= C; FB:= FC; C:= A; FC:= FA     
         "END" INTERCHANGE; 
         TOL:= TOLX; M:= (C + B) * .5; MB:= M - B;        
         "IF" ABS(MB) > TOL "THEN"    
         "BEGIN" "IF" EXT > 3 "THEN" W:= MB "ELSE"        
             "BEGIN" TOL:= TOL * SIGN(MB);      
                 P:= (B - A) * FB; "IF" FIRST "THEN"      
                 "BEGIN" Q:= FA - FB; FIRST:= "FALSE" "END" "ELSE"  
                 "BEGIN" FDB:= (FD - FB) / (D - B);       
                     FDA:= (FD - FA) / (D - A); 
                     P:= FDA * P; Q:= FDB * FA - FDA * FB 
                 "END"; "IF" P < 0 "THEN"       
                 "BEGIN" P:= -P; Q:= -Q "END";  
                 "IF" EXT = 3 "THEN" P:= P * 2; 
                 W:= "IF" P < DW "OR" P <= Q * TOL "THEN" TOL "ELSE"
                 "IF" P < MB * Q "THEN" P / Q "ELSE" MB   
             "END"; D:= A; FD:= FA; A:= B; FA:= FB;       
             X:= B:= B + W; FB:= FX;  
             "IF" ("IF" FC >= 0 "THEN" FB >= 0 "ELSE" FB <= 0) "THEN"         
             "GOTO" INTERPOLATE "ELSE"
             "BEGIN" EXT:= "IF" W = MB "THEN" 0 "ELSE" EXT + 1;     
                 "GOTO" EXTRAPOLATE   
             "END"
         "END"; Y:= C;      
         ZEROINRAT:= "IF" FC >= 0 "THEN" FB <= 0 "ELSE" FB >= 0     
 "END" ZEROINRAT; 
         "EOP"    
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 1      
        
        
        
 AUTHOR: J.C.P.BUS.         
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED: 740218.
        
        
 BRIEF DESCRIPTION:         
        
     THIS  SECTION  CONTAINS  TWO  PROCEDURES  FOR  SOLVING  SYSTEMS  OF      
     NON-LINEAR EQUATIONS, OF WHICH THE JACOBIAN IS KNOWN TO BE A BAND        
     MATRIX.      
     QUANEWBND  ASKS FOR AN APPROXIMATION OF THE JACOBIAN AT THE INITIAL      
     GUESS.       
     QUANEWBND1  CALCULATES AN APPROXIMATION  OF THE JACOBIAN AT THE
     INITIAL GUESS, USING FORWARD DIFFERENCES.  
        
        
 KEYWORDS:        
        
     NON-LINEAR EQUATIONS,  
     SYSTEMS OF EQUATIONS,  
     NO EXPLICIT JACOBIAN.  
        
        
        
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 2      
        
        
        
 SUBSECTION: QUANEWBND.     
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS :      
        
     "PROCEDURE" QUANEWBND(N, LW, RW, X, F, JAC, FUNCT, IN, OUT);   
     "VALUE" N, LW, RW; "INTEGER" N, LW, RW;    
     "ARRAY" X, F, JAC, IN, OUT; "BOOLEAN" "PROCEDURE" FUNCT;       
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE  NUMBER  OF   INDEPENDENT  VARIABLES;   THE  NUMBER  OF      
             EQUATIONS SHOULD ALSO BE EQUAL TO N;         
     LW:     <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF CODIAGONALS  TO THE LEFT OF THE MAIN DIAGONAL      
             OF THE JACOBIAN;         
     RW:     <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF CODIAGONALS TO THE RIGHT OF THE MAIN DIAGONAL      
             OF THE JACOBIAN;         
     X:      <ARRAY IDENTIFIER>;      
             "ARRAY" X[1:N];
             ENTRY:  AN INITIAL ESTIMATE  OF THE SOLUTION  OF THE SYSTEM      
                     THAT HAS TO BE SOLVED;     
             EXIT:   THE CALCULATED SOLUTION OF THE SYSTEM;         
     F:      <ARRAY IDENTIFIER>;      
             "ARRAY" F[1:N];
             ENTRY:  THE  VALUES  OF  THE  FUNCTION  COMPONENTS  AT  THE      
                     INITIAL GUESS;   
             EXIT:   THE  VALUES  OF  THE  FUNCTION  COMPONENTS  AT  THE      
                     CALCULATED SOLUTION;       
     JAC:    <ARRAY IDENTIFIER>;      
             "ARRAY" JAC[1:(LW + RW) * (N - 1) + N];      
             ENTRY:  AN APPROXIMATION  OF  THE JACOBIAN  AT  THE INITIAL      
                     ESTIMATE OF THE SOLUTION;  
             EXIT:   AN APPROXIMATION OF THE JACOBIAN  AT THE CALCULATED      
                     SOLUTION;        
             AN APPROXIMATION OF THE (I, J)-TH  ELEMENT  OF THE JACOBIAN      
             IS GIVEN IN JAC[(LW + RW) * (I - 1) + J], FOR I = 1, ..., N      
             AND J = MAX(1, I - LW), ..., MIN(N, I + RW); 
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 3      
        
        
        
     FUNCT:  <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS :        
             "BOOLEAN" "PROCEDURE" FUNCT(N, L, U, X, F);  
             "VALUE" N, L, U; "INTEGER" N, L, U; "ARRAY" X, F;      
             THE MEANING OF THE FORMAL PARAMETERS IS :    
             N:  <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF INDEPENDENT VARIABLES OF THE FUNCTION  F;      
             L, U: <ARITHMETIC EXPRESSION>;     
                 LOWER  AND  UPPER  BOUND   OF  THE  FUNCTION  COMPONENT      
                 SUBSCRIPT; 
             X:  <ARRAY IDENTIFIER>;  
                 THE INDEPENDENT VARIABLES ARE GIVEN IN X[1:N];     
             F:  <ARRAY IDENTIFIER>;  
                 AFTER  A CALL OF  FUNCT  THE FUNCTION COMPONENTS  F[I],      
                 I = L, ..., U, SHOULD BE GIVEN IN F[L:U];
             IF  THE VALUE  OF  THE PROCEDURE IDENTIFIER  EQUALS  FALSE,      
             THEN THE EXECUTION OF  QUANEWBND  WILL BE TERMINATED, WHILE      
             THE VALUE OF OUT[5] IS SET EQUAL TO 2;       
     IN:     <ARRAY IDENTIFIER>;      
             "ARRAY" IN[0:4];         
             ENTRY :        
             IN  THIS  AUXILIARY  ARRAY  ONE  SHOULD  GIVE THE FOLLOWING      
             VALUES FOR CONTROLLING THE PROCESS:
             IN[0]:  THE MACHINE PRECISION;     
             IN[1]:  THE RELATIVE PRECISION ASKED FOR;    
             IN[2]:  THE ABSOLUTE PRECISION  ASKED FOR;  IF  THE  VALUE,      
                     DELIVERED IN  OUT[5]  EQUALS ZERO,  THEN  THE  LAST      
                     CORRECTION VECTOR  D, SAY,  WHICH  IS A MEASURE FOR      
                     THE ERROR IN THE SOLUTION,  SATIFIES THE INEQUALITY      
                     NORM(D) <= NORM(X) * IN[1] + IN[2],  
                     WHEREBY X DENOTES THE CALCULATED SOLUTION, GIVEN IN      
                     ARRAY  X  AND  NORM(.)  DENOTES THE EUCLIDIAN NORM;      
                     HOWEVER,WE CAN NOT GUARANTEE THAT THE TRUE ERROR IN      
                     THE SOLUTION SATISFIES THIS INEQUALITY,  ESPECIALLY      
                     IF  THE  JACOBIAN  IS  (NEARLY)  SINGULAR   AT  THE      
                     SOLUTION;        
             IN[3]:  THE  MAXIMUM  VALUE  OF  THE NORM  OF  THE RESIDUAL      
                     VECTOR ALLOWED; IF  OUT[5] = 0,  THEN THIS RESIDUAL      
                     VECTOR F, SAY, SATIFIES: NORM(F) <= IN[3];     
             IN[4]:  THE   MAXIMUM   NUMBER    OF   FUNCTION   COMPONENT      
                     EVALUATIONS ALLOWED;  L - U + 1  FUNCTION COMPONENT      
                     EVALUATIONS    ARE    COUNTED    EACH    CALL    OF      
                     FUNCT(N, L, U, X, F); IF OUT[5]=1, THEN THE PROCESS      
                     IS TERMINATED,  BECAUSE  THE NUMBER  OF EVALUATIONS      
                     EXCEEDED THE VALUE GIVEN IN IN[4];   
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 4      
        
        
        
     OUT:    <ARRAY IDENTIFIER>;      
             "ARRAY" OUT[1:5];        
             EXIT :         
             OUT[1]: THE EUCLIDIAN NORM OF THE LAST STEP ACCEPTED;  
             OUT[2]: THE EUCLIDIAN NORM  OF  THE RESIDUAL VECTOR  AT THE      
                     CALCULATED SOLUTION;       
             OUT[3]: THE  NUMBER   OF   FUNCTION  COMPONENT  EVALUATIONS      
                     PERFORMED;       
             OUT[4]: THE NUMBER OF ITERATIONS CARRIED OUT;
             OUT[5]: THE INTEGER VALUE DELIVERED IN  OUT[5]  GIVES  SOME      
                     INFORMATION  ABOUT THE TERMINATION  OF THE PROCESS;      
                     OUT[5] = 0: THE PROCESS  IS TERMINATED  IN A NORMAL      
                                 WAY; THE LAST STEP AND  THE NORM OF THE      
                                 RESIDUAL VECTOR SATISFY  THE CONDITIONS      
                                 (SEE IN[2], IN[3]);      
                     IF  OUT[5] ^= 0,  THEN  THE PROCESS  IS  TERMINATED      
                     PREMATURALY;     
                     OUT[5] = 1: THE   NUMBER   OF   FUNCTION  COMPONENT      
                                 EVALUATIONS  EXCEEDS THE VALUE GIVEN IN      
                                 IN[4];         
                     OUT[5] = 2: A CALL OF  FUNCT  DELIVERED  THE  VALUE      
                                 FALSE;         
                     OUT[5] = 3: THE  APPROXIMATION   OF   THE  JACOBIAN      
                                 MATRIX TURNS OUT TO BE SINGULAR.   
        
        
 PROCEDURES USED: 
        
     MULVEC      = CP31020, 
     DUPVEC      = CP31030, 
     VECVEC      = CP34010, 
     ELMVEC      = CP34020, 
     DECSOLBND   = CP34322. 
        
        
 EXECUTION FIELD LENGTH:    
        
     THE MAXIMUM NUMBER OF WORDS, NECESSARY  FOR THE ARRAYS  DECLARED IN      
     QUANEWBND EQUALS  MAX(N * 3 + (N - 1) * (LW + RW), 4N).        
        
        
 RUNNING TIME: PROPORTIONAL TO  N * LW * ( LW + RW + 1).  
        
        
 LANGUAGE:   ALGOL 60.      
        
        
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 5      
        
        
        
 METHOD AND PERFORMANCE:    
        
     THE METHOD USED IN QUANEWBND IS THE SAME AS GIVEN IN [1];  THE SAME      
     PROBLEMS HAVE BEEN TESTED AND THE RESULTS ARE  THE SAME  OR  BETTER      
     THAN THOSE REPORTED IN  [1];  CITING  [1],  WE CAN SAY  THAT  "THIS      
     METHOD OFFERS A USEFUL IF MODEST IMPROVEMENT UPON  NEWTON'S METHOD,      
     BUT THIS IMPROVEMENT TENDS TO VANISH  AS THE NONLINEARITIES  BECOME      
     MORE PRONOUNCED".      
        
        
 EXAMPLE OF USE: SEE QUANEWBND1 (THIS SECTION). 
        
        
 REFERENCES:      
        
     [1] BROYDEN C.G.       
         THE CONVERGENCE  OF  AN ALGORITHM  FOR SOLVING SPARSE NONLINEAR      
         SYSTEMS. 
         MATH. COMP., VOL.25 (1971).  
        
        
        
 SUBSECTION: QUANEWBND1.    
        
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE READS :      
        
     "PROCEDURE" QUANEWBND1(N, LW, RW, X, F, FUNCT, IN, OUT);       
     "VALUE" N, LW, RW; "INTEGER" N, LW, RW;    
     "ARRAY" X, F, IN, OUT; "BOOLEAN" "PROCEDURE" FUNCT;  
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE  NUMBER  OF   INDEPENDENT  VARIABLES;   THE  NUMBER  OF      
             EQUATIONS SHOULD ALSO BE EQUAL TO N;         
     LW:     <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF CODIAGONALS  TO THE LEFT OF THE MAIN DIAGONAL      
             OF THE JACOBIAN;         
     RW:     <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF CODIAGONALS TO THE RIGHT OF THE MAIN DIAGONAL      
             OF THE JACOBIAN;         
     X:      <ARRAY IDENTIFIER>;      
             "ARRAY" X[1:N];
             ENTRY:  AN INITIAL ESTIMATE  OF THE SOLUTION  OF THE SYSTEM      
                     THAT HAS TO BE SOLVED;     
             EXIT:   THE CALCULATED SOLUTION OF THE SYSTEM;         
     F:      <ARRAY IDENTIFIER>;      
             "ARRAY" F[1:N];
             EXIT:   THE  VALUES  OF  THE  FUNCTION  COMPONENTS  AT  THE      
                     CALCULATED SOLUTION;       
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 6      
        
        
        
     FUNCT:  <PROCEDURE IDENTIFIER>;  
             THE HEADING OF THIS PROCEDURE READS :        
             "BOOLEAN" "PROCEDURE" FUNCT(N, L, U, X, F);  
             "VALUE" N, L, U; "INTEGER" N, L, U; "ARRAY" X, F;      
             THE MEANING OF THE FORMAL PARAMETERS IS :    
             N:  <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF INDEPENDENT VARIABLES OF THE FUNCTION  F;      
             L, U: <ARITHMETIC EXPRESSION>;     
                 LOWER  AND  UPPER  BOUND   OF  THE  FUNCTION  COMPONENT      
                 SUBSCRIPT; 
             X:  <ARRAY IDENTIFIER>;  
                 THE INDEPENDENT VARIABLES ARE GIVEN IN X[1:N];     
             F:  <ARRAY IDENTIFIER>;  
                 AFTER  A CALL OF  FUNCT  THE FUNCTION COMPONENTS  F[I],      
                 I = L, ..., U, SHOULD BE GIVEN IN F[L:U];
             IF  THE VALUE  OF  THE PROCEDURE IDENTIFIER  EQUALS  FALSE,      
             THEN THE EXECUTION OF  QUANEWBND  WILL BE TERMINATED, WHILE      
             THE VALUE OF OUT[5] IS SET EQUAL TO 2;       
     IN:     <ARRAY IDENTIFIER>;      
             "ARRAY" IN[0:4];         
             ENTRY :        
             IN  THIS  AUXILIARY  ARRAY  ONE  SHOULD  GIVE THE FOLLOWING      
             VALUES FOR CONTROLLING THE PROCESS:
             IN[0]:  THE MACHINE PRECISION;     
             IN[1]:  THE RELATIVE PRECISION ASKED FOR;    
             IN[2]:  THE ABSOLUTE PRECISION  ASKED FOR;  IF  THE  VALUE,      
                     DELIVERED IN  OUT[5]  EQUALS ZERO,  THEN  THE  LAST      
                     CORRECTION VECTOR  D, SAY,  WHICH  IS A MEASURE FOR      
                     THE ERROR IN THE SOLUTION,  SATIFIES THE INEQUALITY      
                     NORM(D) <= NORM(X) * IN[1] + IN[2],  
                     WHEREBY X DENOTES THE CALCULATED SOLUTION, GIVEN IN      
                     ARRAY  X  AND  NORM(.)  DENOTES THE EUCLIDIAN NORM;      
                     HOWEVER,WE CAN NOT GUARANTEE THAT THE TRUE ERROR IN      
                     THE SOLUTION SATISFIES THIS INEQUALITY,  ESPECIALLY      
                     IF  THE  JACOBIAN  IS  (NEARLY)  SINGULAR   AT  THE      
                     SOLUTION;        
             IN[3]:  THE  MAXIMUM  VALUE  OF  THE NORM  OF  THE RESIDUAL      
                     VECTOR ALLOWED; IF  OUT[5] = 0,  THEN THIS RESIDUAL      
                     VECTOR F, SAY, SATIFIES: NORM(F) <= IN[3];     
             IN[4]:  THE   MAXIMUM   NUMBER    OF   FUNCTION   COMPONENT      
                     EVALUATIONS ALLOWED;  L - U + 1  FUNCTION COMPONENT      
                     EVALUATIONS    ARE    COUNTED    EACH    CALL    OF      
                     FUNCT(N, L, U, X, F); IF OUT[5]=1, THEN THE PROCESS      
                     IS TERMINATED,  BECAUSE  THE NUMBER  OF EVALUATIONS      
                     EXCEEDED THE VALUE GIVEN IN IN[4];   
             IN[5]:  THE  JACOBIAN  MATRIX  AT  THE  INITIAL  GUESS   IS      
                     APPROXIMATED  USING  FORWARD DIFFERENCES,  WITH  AN      
                     FIXED INCREMENT TO EACH VARIABLE  THAT  EQUALS  THE      
                     VALUE GIVEN IN IN[5];      
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 7      
        
        
        
     OUT:    <ARRAY IDENTIFIER>;      
             "ARRAY" OUT[1:5];        
             EXIT :         
             OUT[1]: THE EUCLIDIAN NORM OF THE LAST STEP ACCEPTED;  
             OUT[2]: THE EUCLIDIAN NORM  OF  THE RESIDUAL VECTOR  AT THE      
                     CALCULATED SOLUTION;       
             OUT[3]: THE  NUMBER   OF   FUNCTION  COMPONENT  EVALUATIONS      
                     PERFORMED;       
             OUT[4]: THE NUMBER OF ITERATIONS CARRIED OUT;
             OUT[5]: THE INTEGER VALUE DELIVERED IN  OUT[5]  GIVES  SOME      
                     INFORMATION  ABOUT THE TERMINATION  OF THE PROCESS;      
                     OUT[5] = 0: THE PROCESS  IS TERMINATED  IN A NORMAL      
                                 WAY; THE LAST STEP AND  THE NORM OF THE      
                                 RESIDUAL VECTOR SATISFY  THE CONDITIONS      
                                 (SEE IN[2], IN[3]);      
                     IF  OUT[5] ^= 0,  THEN  THE PROCESS  IS  TERMINATED      
                     PREMATURALY;     
                     OUT[5] = 1: THE   NUMBER   OF   FUNCTION  COMPONENT      
                                 EVALUATIONS  EXCEEDS THE VALUE GIVEN IN      
                                 IN[4];         
                     OUT[5] = 2: A CALL OF  FUNCT  DELIVERED  THE  VALUE      
                                 FALSE;         
                     OUT[5] = 3: THE  APPROXIMATION   OF   THE  JACOBIAN      
                                 MATRIX TURNS OUT TO BE SINGULAR.   
        
        
 PROCEDURES USED: 
        
     JACOBNBNDF  = CP34439, 
     QUANEWBND   = CP34430. 
        
        
 EXECUTION FIELD LENGTH:    
        
     QUANEWBND1  DECLARES AN AUXILIARY ARRAY OF DIMENSION ONE  AND ORDER      
     N + (N - 1) * (LW + RW).         
        
        
 RUNNING TIME: PROPORTIONAL TO  N * LW * ( LW + RW + 1).  
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     QUANEWBND1  USES  JACOBNBNDF  (SECTION  4.3.2.1)  TO  CALCULATE  AN      
     INITIAL APPROXIMATION OF THE JACOBIAN MATRIX   AT THE INITIAL GUESS      
     GIVEN  IN  X[1:N]  AND  SOLVES  THE  NONLINEAR  SYSTEM  BY  CALLING      
     QUANEWBND (THIS SECTION).        
        
        
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 8      
        
        
        
 EXAMPLE OF USE:  
        
     LET THE FUNCTION F BE DEFINED BY (SEE [1]):
     F[1] = (3 - 2 * X[1]) * X[1] + 1 - 2 * X[2],         
     F[I] = (3 - 2 * X[I]) * X[I] + 1 - X[I - 1] - 2 * X[I + 1],  I = 2,      
     ..., N - 1,  
     F[N] = (3 - 2 * X[N]) * X[N] + 1 - X[N - 1];         
     LET AN INITIAL ESTIMATE OF THE SOLUTION OF THE SYSTEM  F(X) = 0  BE      
     GIVEN BY X[I] = -1, I = 1, ..., N;  THEN THE FOLLOWING PROGRAM  MAY      
     SOLVE THIS SYSTEM FOR N = 600 AND PRINTS SOME RESULTS.         
        
     "BEGIN"      
         "PROCEDURE" QUANEWBND1(N, L, R, X, F, FU, I, O); "CODE" 34431;       
         "BOOLEAN" "PROCEDURE" FUN(N, L, U, X, F); "VALUE" N, L, U; 
         "INTEGER" N, L, U; "ARRAY" X, F;       
         "BEGIN" "INTEGER" I; "REAL" X1, X2, X3;
             X1:= "IF" L = 1 "THEN" 0 "ELSE" X[L - 1]; X2:= X[L];   
             X3:= "IF" L = N "THEN" 0 "ELSE" X[L + 1];    
             "FOR" I:= L "STEP" 1 "UNTIL" U "DO"
             "BEGIN" F[I]:= (3 - 2 * X2) * X2 + 1 - X1 - X3 * 2;    
                 X1:= X2; X2:= X3;    
                 X3:= "IF" I <= N - 2 "THEN" X[I + 2] "ELSE" 0      
             "END"; FUN:= "TRUE"      
         "END" FUN;         
        
         "INTEGER" I; "ARRAY" X, F[1:600], IN[0:5], OUT[1:5];       
         "FOR" I := 1 "STEP" 1 "UNTIL" 600 "DO" X[I]:= -1;
         IN[0]:= "-14; IN[1]:= IN[2]:= IN[3]:= "-6; IN[4]:= 20000;  
         IN[5]:= 0.001;     
         QUANEWBND1(600, 1, 1, X, F, FUN, IN, OUT);       
         OUTPUT(71, "("//,"(" NORM RESIDUALVECTOR: ")"+.15D"+3D,/,  
         "(" LENGTH OF LAST STEP: ")"+.15D"+3D,/,         
         "(" NUMBER OF FUNCTION COMPONENT EVALUATIONS: ")"5ZD,/,    
         "(" NUMBER OF ITERATIONS: ")"4ZD/,"("REPORT: ")"D/")",     
         OUT[2], OUT[1], OUT[3], OUT[4], OUT[5])
     "END"        
        
     RESULTS:     
        
     NORM RESIDUALVECTOR: +.221010684482660"-006
     LENGTH OF LAST STEP: +.302712457332660"-006
     NUMBER OF FUNCTION COMPONENT EVALUATIONS:   6598     
     NUMBER OF ITERATIONS:     7      
     REPORT: 0    
        
        
 REFERENCES:      
        
     [1] BROYDEN C.G.       
         THE CONVERGENCE  OF  AN ALGORITHM  FOR SOLVING SPARSE NONLINEAR      
         SYSTEMS. 
         MATH. COMP., VOL.25 (1971).  
        
        
        
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                      PAGE 9      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE" 34430;    
     "PROCEDURE" QUANEWBND(N, LW, RW, X, F, JAC, FUNCT, IN, OUT);   
     "VALUE" N, LW, RW; "INTEGER" N, LW, RW;    
     "ARRAY" X, F, JAC, IN, OUT; "BOOLEAN" "PROCEDURE" FUNCT;       
     "BEGIN" "INTEGER" L, IT, FCNT, FMAX, ERR, B;         
         "REAL" MACHEPS, RELTOL, ABSTOL, TOLRES, ND, MZ, RES;       
         "ARRAY" DELTA[1:N];
        
         "REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
         "REAL" "PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X); "CODE" 34020;       
         "PROCEDURE" MULVEC(L, U, SHIFT, A, B, X); "CODE" 31020;    
         "PROCEDURE" DUPVEC(L, U, SHIFT, A, B); "CODE" 31030;       
         "REAL" "PROCEDURE" DECSOLBND(A, N, LW, RW, AUX, B);        
         "CODE" 34322;      
        
        
         "REAL" "PROCEDURE" EVALUATE(N, X, F); "VALUE" N; 
         "INTEGER" N; "ARRAY" X, F;   
         "BEGIN" FCNT:= FCNT + N; "IF" ^ FUNCT(N, 1, N, X, F) "THEN"
             "BEGIN" ERR:= 2; "GOTO" EXIT "END";
             "IF" FCNT > FMAX "THEN" ERR:= 1;   
             EVALUATE:= SQRT(VECVEC(1, N, 0, F, F))       
         "END" EVAL;        
        
         "BOOLEAN" "PROCEDURE" DIRECTION;       
         "BEGIN" "ARRAY" LU[1:L], AUX[1:5]; AUX[2]:= MACHEPS;       
             MULVEC(1, N, 0, DELTA, F, -1); DUPVEC(1, L, 0, LU, JAC);         
             DECSOLBND(LU, N, LW, RW, AUX, DELTA);        
             DIRECTION:= AUX[3] = N   
         "END" SOLLINSYS;   
        
         "BOOLEAN" "PROCEDURE" TEST(ND, TOLD, NRES, TOLRES, ERR);   
         "VALUE" ND, TOLD; "INTEGER" ERR; "REAL" ND, TOLD, NRES, TOLRES;      
         TEST:= ERR ^= 0 "OR" (NRES < TOLRES "AND" ND < TOLD);      
        
         "PROCEDURE" UPDATE JAC;      
         "BEGIN" "INTEGER"  I, J, K, R, M; "REAL" MUL, CRIT;        
             "ARRAY" PP, S[1:N];      
             CRIT:= ND * MZ;
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" PP[I]:= DELTA[I] ** 2;       
             R:= 1; K:= 1; M:= RW + 1;
             "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
             "BEGIN" MUL:= 0; "FOR" J:= R "STEP" 1 "UNTIL" M "DO"   
                 MUL:= MUL + PP[J]; J:= R - K;  
                 "IF" ABS(MUL) > CRIT "THEN"    
                 ELMVEC(K, M - J, J, JAC, DELTA, F[I] / MUL); K:= K + B;      
                 "IF" I > LW "THEN" R:= R + 1 "ELSE" K:= K - 1;     
                 "IF" M < N "THEN" M:= M + 1    
             "END"
         "END" UPDATEJAC;   
1SECTION : 5.1.1.2.2          (OCTOBER 1974)                     PAGE 10      
        
        
        
        
         MACHEPS:= IN[0]; RELTOL:= IN[1]; ABSTOL:= IN[2]; 
         TOLRES:= IN[3]; FMAX:= IN[4]; MZ:= MACHEPS ** 2; 
         IT:= FCNT:= 0; B:= LW + RW; L:= (N - 1) * B + N; B:= B + 1;
         RES:= SQRT(VECVEC(1, N, 0, F, F)); ERR:= 0;      
     ITERATE: "IF" ^ TEST(SQRT(ND), SQRT(VECVEC(1, N, 0, X, X)) * RELTOL      
         + ABSTOL, RES, TOLRES, ERR) "THEN"     
         "BEGIN" IT:= IT + 1; "IF" IT ^= 1 "THEN" UPDATEJAC;        
             "IF" ^ DIRECTION "THEN" ERR:= 3 "ELSE"       
             "BEGIN" ELMVEC(1, N, 0, X, DELTA, 1);        
                 ND:= VECVEC(1, N, 0, DELTA, DELTA);      
                 RES:= EVALUATE(N, X, F); "GOTO" ITERATE  
              "END"         
         "END";   
     EXIT: OUT[1]:= SQRT(ND); OUT[2]:=RES; OUT[3]:= FCNT; 
         OUT[4]:= IT; OUT[5]:= ERR    
     "END" QUANEWBND;       
         "EOP"    
        
 "CODE" 34431;    
     "PROCEDURE" QUANEWBND1(N, LW, RW, X, F, FUNCT, IN, OUT);       
     "VALUE" N, LW, RW; "INTEGER" N, LW, RW; "ARRAY" X, F, IN, OUT; 
     "BOOLEAN" "PROCEDURE" FUNCT;     
     "BEGIN" "INTEGER" I, K; "REAL" S;
         "PROCEDURE" QUANEWBND(N, L, R, X, F, J, G, I, O); "CODE" 34430;      
         "ARRAY" JAC[1:(LW + RW) * (N - 1) + N];
         "PROCEDURE" JACOBNBNDF(N,L, R, X, F, J, I, D, F); "CODE" 34439;      
         FUNCT(N, 1, N, X, F); S:= IN[5];       
         K:= (LW + RW)*(N - 1) + N*2 - ((LW - 1)*LW + (RW - 1)*RW) // 2;      
         IN[4]:= IN[4] - K; 
         JACOBNBNDF(N, LW, RW, X, F, JAC, I, S, FUNCT);   
         QUANEWBND(N, LW, RW, X, F, JAC, FUNCT, IN, OUT); 
         IN[4]:= IN[4] + K; OUT[3]:= OUT[3] + K 
     "END" QUANEWBND1;      
         "EOP"    
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 1      
        
        
        
 AUTHORS: J.C.P. BUS, B. VAN DOMSELAAR, J. KOK. 
        
        
 CONTRIBUTORS:  B. VAN DOMSELAAR, J. KOK.       
        
        
 INSTITUTE: MATHEMATICAL CENTRE.      
        
        
 RECEIVED:  741101.         
        
        
 BRIEF DESCRIPTION:         
        
     THIS SECTION CONTAINS TWO PROCEDURES;      
     MARQUARDT CALCULATES THE LEAST SQUARES SOLUTION OF AN
     OVERDETERMINED SYSTEM OF NONLINEAR EQUATIONS WITH MARQUARDT'S  
     METHOD;      
     GSSNEWTON CALCULATES THE LEAST SQUARES SOLUTION OF AN
     OVERDETERMINED SYSTEM OF NONLINEAR EQUATIONS WITH THE GAUSS-NEWTON       
     METHOD;      
     THE USER HAS TO PROGRAM THE EVALUATION OF THE RESIDUAL VECTOR OF         
     THE SYSTEM AND THE JACOBIAN MATRIX OF ITS PARTIAL DERIVATIVES. 
        
        
 KEYWORDS:        
        
     NONLINEAR EQUATIONS,   
     LEAST SQUARES SOLUTION,
     OVERDETERMINED SYSTEM, 
     MARQUARDT'S METHOD,    
     GAUSS-NEWTON METHOD,   
     CURVE FITTING.         
        
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 2      
        
        
        
 SUBSECTION: MARQUARDT.     
        
 CALLING SEQUENCE:
        
     THE HEADING OF THIS PROCEDURE IS:
     "PROCEDURE" MARQUARDT(M, N, PAR, RV, JJINV, FUNCT, JACOBIAN, IN,         
        OUT); "VALUE" M, N; "INTEGER" M, N;     
     "ARRAY" PAR, RV, JJINV, IN, OUT; "BOOLEAN" "PROCEDURE" FUNCT;  
     "PROCEDURE" JACOBIAN;  
     "CODE" 34440;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     M:        <ARITHMETIC EXPRESSION>;         
               THE NUMBER OF EQUATIONS;         
     N:        <ARITHMETIC EXPRESSION>;         
               THE NUMBER OF UNKNOWN VARIABLES; N SHOULD SATISFY N<=M;        
     PAR:      <ARRAY IDENTIFIER>;    
               "ARRAY" PAR[1 : N];    
               THE UNKNOWN VARIABLES OF THE SYSTEM;       
               ENTRY: AN APPROXIMATION TO A LEAST SQUARES SOLUTION  
                      OF THE SYSTEM;  
               EXIT:  THE CALCULATED LEAST SQUARES SOLUTION;        
     RV:       <ARRAY IDENTIFIER>;    
               "ARRAY" RV[1 : M];     
               EXIT:  THE RESIDUAL VECTOR AT  THE CALCULATED SOLUTION;        
     JJINV:    <ARRAY IDENTIFIER>;    
               "ARRAY" JJINV[1 : N, 1 : N];     
               EXIT:  THE INVERSE OF THE MATRIX  J' * J  WHERE J DENOTES      
                      THE MATRIX OF PARTIAL DERIVATIVES DRV[I] / DPAR[J]      
                      (I=1,...,M; J=1,...,N) AND  J'  DENOTES THE   
                      TRANSPOSE OF J. 
     FUNCT:    <PROCEDURE IDENTIFIER>;
               THE HEADING OF THIS PROCEDURE SHOULD BE:   
               "BOOLEAN" "PROCEDURE" FUNCT(M, N, PAR, RV); "VALUE" M, N;      
               "INTEGER" M, N; "ARRAY" PAR, RV; 
               ENTRY: M, N, PAR;      
                      M, N  HAVE THE SAME MEANING AS IN THE PROCEDURE         
                      MARQUARDT;      
                      "ARRAY" PAR[1:N] CONTAINS THE CURRENT VALUES OF         
                      THE UNKNOWNS AND SHOULD NOT BE ALTERED;       
               EXIT:  "ARRAY" RV[1 : M];        
                      UPON COMPLETION OF A CALL OF FUNCT, THIS ARRAY  RV      
                      SHOULD CONTAIN THE RESIDUAL VECTOR, OBTAINED WITH       
                      THE CURRENT VALUES OF THE UNKNOWNS; 
                      E.G. IN CURVE FITTING PROBLEMS:     
                      RV[I] := THEORETICAL VALUE  F(X[I], PAR) -    
                              OBSERVED VALUE  Y[I];       
               AFTER A SUCCESSFUL CALL OF FUNCT, THE BOOLEAN PROCEDURE        
               SHOULD DELIVER THE VALUE TRUE;   
               HOWEVER, IF FUNCT DELIVERS THE VALUE FALSE, THEN IT IS         
               ASSUMED THAT THE CURRENT ESTIMATES OF THE UNKNOWNS LIE         
               OUTSIDE A FEASIBLE REGION AND THE PROCESS IS TERMINATED        
               (SEE OUT[1]);
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 3      
        
        
        
               HENCE, PROPER PROGRAMMING OF FUNCT MAKES IT POSSIBLE TO        
               AVOID CALCULATION OF A RESIDUAL VECTOR WITH VALUES OF THE      
               UNKNOWN VARIABLES WHICH MAKE NO SENSE OR WHICH EVEN MAY        
               CAUSE OVERFLOW IN THE COMPUTATION;         
     JACOBIAN: <PROCEDURE IDENTIFIER>;
               THE HEADING OF THIS PROCEDURE SHOULD BE:   
               "PROCEDURE" JACOBIAN(M, N, PAR, RV, JAC, LOCFUNCT);  
               "VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV, JAC;  
               "PROCEDURE" LOCFUNCT;  
               ENTRY: M, N, PAR, RV, LOCFUNCT;  
                      FOR M,N,PAR SEE: FUNCT;   
                      RV CONTAINS THE RESIDUAL VECTOR OBTAINED WITH THE       
                      CURRENT VALUES OF THE UNKNOWNS AND SHOULD NOT BE        
                      ALTERED;        
                      A CALL OF LOCFUNCT(M,N,PAR,RV) IS EQUIVALENT WITH       
                      A CALL OF THE USER-DEFINED PROCEDURE
                      FUNCT(M,N,PAR,RV), BUT, IN ADDITION, THIS CALL IS       
                      COUNTED TO THE TOTAL NUMBER OF CALLS OF FUNCT 
                      (SEE OUT[4]) AND, MOREOVER, IF FUNCT DELIVERS THE       
                      VALUE FALSE THEN THE PROCESS IS TERMINATED;   
               EXIT:  "ARRAY" JAC[1 : M, 1 : N];
                      UPON COMPLETION OF A CALL OF JACOBIAN, JAC SHOULD       
                      CONTAIN THE PARTIAL DERIVATIVES  DRV[I] / DPAR[J],      
                      OBTAINED WITH THE CURRENT VALUES OF THE UNKNOWN         
                      VARIABLES GIVEN IN PAR[1:N];        
               IT IS A PREREQUISITE FOR THE PROPER OPERATION OF THE 
               PROCEDURE MARQUARDT THAT THE PRECISION OF THE ELEMENTS OF      
               THE MATRIX JAC  IS AT LEAST THE PRECISION DEFINED BY 
               IN[3] AND IN[4];       
     IN:       <ARRAY IDENTIFIER>;    
               "ARRAY" IN[0 : 6];     
               ENTRY:  IN THIS ARRAY THE USER SHOULD GIVE SOME DATA TO        
               CONTROL THE PROCESS;   
               IN[0]:  THE MACHINE PRECISION;   
                       FOR THE CYBER 73 A SUITABLE VALUE IS "-14;   
               IN[1], IN[2] ARE NOT USED BY THE PROCEDURE  MARQUARDT;         
               IN[3], IN[4]:
                       THE RELATIVE AND ABSOLUTE TOLERANCE FOR THE  
                       DIFFERENCE BETWEEN THE EUCLIDEAN NORM OF THE 
                       ULTIMATE AND PENULTIMATE RESIDUAL VECTOR;    
                       THE PROCESS IS TERMINATED IF THE IMPROVEMENT OF        
                       THE SUM OF SQUARES IS LESS THAN    
                       IN[3] * (SUM OF SQUARES) + IN[4] * IN[4];    
                       THESE TOLERANCES SHOULD BE CHOSEN GREATER THAN         
                       THE CORRESPONDING ERRORS OF THE CALCULATED   
                       RESIDUAL VECTOR;         
                       NOTE THAT THE EUCLIDEAN NORM OF THE RESIDUAL 
                       VECTOR IS DEFINED AS THE SQUARE ROOT OF THE SUM        
                       OF SQUARES;    
               IN[5]:  THE MAXIMUM NUMBER OF CALLS OF FUNCT ALLOWED;
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 4      
        
        
        
               IN[6]:  A STARTING VALUE USED FOR THE RELATION BETWEEN         
                       THE GRADIENT AND THE GAUSS-NEWTON DIRECTION (SEE       
                       [2]);  IF THE PROBLEM IS WELL CONDITIONED THEN A       
                       SUITABLE VALUE FOR IN[6] WILL BE  0.01; IF THE         
                       PROBLEM IS ILL CONDITIONED THEN IN[6] SHOULD BE        
                       GREATER, BUT THE VALUE OF IN[6] SHOULD SATISFY:        
                       IN[0] < IN[6] <= 1/IN[0];
     OUT:      <ARRAY IDENTIFIER>; "ARRAY" OUT[1 : 7];    
               EXIT :  IN ARRAY  OUT  SOME BY-PRODUCTS ARE DELIVERED;         
               OUT[1]: THIS VALUE GIVES INFORMATION ABOUT THE       
                       TERMINATION OF THE PROCESS;        
                       OUT[1]=0:  NORMAL TERMINATION;     
                       OUT[1]=1:  THE PROCESS HAS BEEN BROKEN OFF,  
                                  BECAUSE THE NUMBER OF CALLS OF FUNCT        
                                  EXCEEDED THE NUMBER GIVEN IN IN[5];         
                       OUT[1]=2:  THE PROCESS HAS BEEN BROKEN OFF,  
                                  BECAUSE A CALL OF FUNCT DELIVERED THE       
                                  VALUE FALSE;  
                       OUT[1]=3:  FUNCT BECAME FALSE WHEN CALLED WITH         
                                  THE INITIAL ESTIMATES OF PAR[1:N];
                                  THE ITERATION PROCESS WAS NOT STARTED       
                                  AND SO JJINV[1:N,1:N] CAN NOT BE USED;      
                       OUT[1]=4:  THE PROCESS HAS BEEN BROKEN OFF,  
                                  BECAUSE THE PRECISION ASKED FOR CAN         
                                  NOT BE ATTAINED; THIS PRECISION IS
                                  POSSIBLY CHOSEN TOO HIGH, RELATIVE TO       
                                  THE PRECISION IN WHICH THE RESIDUAL         
                                  VECTOR IS CALCULATED (SEE IN[3]); 
               OUT[2]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR    
                       CALCULATED WITH VALUES OF THE UNKNOWNS DELIVERED;      
               OUT[3]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR    
                       CALCULATED WITH THE INITIAL VALUES OF THE    
                       UNKNOWN VARIABLES;       
               OUT[4]: THE NUMBER OF CALLS OF FUNCT NECESSARY TO OBTAIN       
                       THE CALCULATED RESULT;   
               OUT[5]: THE TOTAL NUMBER OF ITERATIONS PERFORMED; NOTE         
                       THAT IN EACH ITERATION ONE EVALUATION OF THE 
                       JACOBIAN MATRIX HAD TO BE MADE;    
               OUT[6]: THE  IMPROVEMENT OF THE EUCLIDEAN NORM OF THE
                       RESIDUAL VECTOR IN THE LAST ITERATION STEP;  
               OUT[7]: THE CONDITION NUMBER OF  J' * J , I.E. THE RATIO       
                       OF ITS LARGEST TO SMALLEST EIGENVALUES;      
        
 DATA AND RESULTS:
        
     IF THIS PROCEDURE IS USED FOR CURVE FITTING THEN THE RELATIVE  
     ACCURACY IN THE CALCULATION OF THE RESIDUAL VECTOR DEPENDS STRONGLY      
     ON THE ERRORS IN THE EXPERIMENTAL DATA AND THIS SHOULD BE REFLECTED      
     IN THE PARAMETERS IN[3] AND IN[4];         
     THE MATRIX  JJINV  CAN BE USED IF SOME STATISTICAL INFORMATION 
     ABOUT THE FITTED PARAMETERS IS REQUIRED; THE STANDARD DEVIATION,         
     COVARIANCE MATRIX AND CORRELATION MATRIX MAY BE CALCULATED EASILY        
     FROM  JJINV ;
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 5      
        
        
        
 PROCEDURES USED: 
        
     MULCOL = CP31022,      
     DUPVEC = CP31030,      
     VECVEC = CP34010,      
     MATVEC = CP34011,      
     TAMVEC = CP34012,      
     MATTAM = CP34015,      
     QRISNGVALDEC = CP34273.
        
        
 REQUIRED CENTRAL MEMORY:   
        
     EXECUTION FIELD LENGTH: ONE ARRAY OF LENGTH  M * N  AND FOUR   
                             ARRAYS OF LENGTH  N  ARE DECLARED;     
        
        
 RUNNING TIME:    
        
     THE NUMBER OF ITERATION STEPS DEPENDS STRONGLY ON THE PROBLEM TO BE      
     SOLVED; HOWEVER, THE WORK DONE EACH ITERATION STEP IS PROPORTIONAL       
     TO N CUBED;  
        
        
 LANGUAGE:   ALGOL 60.      
        
        
 METHOD AND PERFORMANCE:    
        
     MARQUARDT USES MARQUARDT'S METHOD; FOR DETAILS SEE [2];        
        
        
        
 REFERENCES:      
        
     [1] D. W. MARQUARDT,   
         AN ALGORITHM FOR LEAST-SQUARES ESTIMATION OF NONLINEAR     
         PARAMETERS,        
         J. SIAM 11 (1963), PP.431-441.         
        
     [2] J. C. P. BUS, B. VAN DOMSELAAR, J. KOK,
         NONLINEAR LEAST SQUARES ESTIMATION,    
         MATHEMATICAL CENTRE, AMSTERDAM. (TO APPEAR)      
        
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 6      
        
        
        
 EXAMPLE OF USE:  
        
     THE PARAMETERS PAR[1 : 3] IN THE CURVE FITTING PROBLEM:        
     RV[I] = PAR[1] + PAR[2] * EXP(PAR[3] * X[I]) - Y[I]  
     WITH  (X[I], Y[I]), I=1,...,6  AS THE EXPERIMENTAL DATA, MAY BE
     OBTAINED BY THE FOLLOWING PROGRAM:         
        
     "BEGIN" "PROCEDURE" MARQUARDT(M,N,P,RV,JJINV,F,J,I,O); "CODE"34440;      
        "INTEGER" I;        
        "ARRAY" IN[0:6],OUT[0:7],X,Y,RV[1:6],JJINV[1:3,1:3],PAR[1:3];         
        
        "BOOLEAN" "PROCEDURE" EXPFUNCT(M,N,PAR,RV);       
            "VALUE" M,N; "INTEGER" M,N; "ARRAY" PAR,RV;   
            "BEGIN" "INTEGER" I;      
                    "FOR" I:=1 "STEP" 1 "UNTIL" M "DO"    
                    "BEGIN" "IF" PAR[3]*X[I]>680 "THEN"   
                            "BEGIN" EXPFUNCT:="FALSE";    
                                    "GOTO" STOP 
                            "END";    
                            RV[I]:=PAR[1]+PAR[2]*EXP(PAR[3]*X[I])-Y[I]        
                    "END"; EXPFUNCT:="TRUE";    
                STOP:       
           "END" EXPFUNCT;  
        
        "PROCEDURE" JACOBIAN(M,N,PAR,RV,JAC,LOCFUNCT);    
           "VALUE" M,N; "INTEGER" M,N; "ARRAY" PAR,RV,JAC;
           "PROCEDURE" LOCFUNCT;      
           "BEGIN" "INTEGER" I; "REAL" EX;      
                   "FOR" I:=1 "STEP" 1 "UNTIL" M "DO"     
                   "BEGIN" JAC[I,1]:=1;         
                           JAC[I,2]:=EX:=EXP(PAR[3]*X[I]);
                           JAC[I,3]:=X[I]*PAR[2]*EX       
                   "END"    
           "END" JACOBIAN;  
        
        IN[0]:="-14; IN[3]:="-4; IN[4]:="-1; IN[5]:=75; IN[6]:="-2; 
        "FOR" I:=1 "STEP" 1 "UNTIL" 6 "DO"      
        INPUT(60,"("2(N),/")",X[I],Y[I]);       
        PAR[1]:=580; PAR[2]:=-180; PAR[3]:=-0.160;        
        MARQUARDT(6,3,PAR,RV,JJINV,EXPFUNCT,JACOBIAN,IN,OUT);       
        OUTPUT(61,"("3/,"("X[I]     Y[I]")",/,6(B,+D,5B,3D.D,/),2/, 
        "("PARAMETERS")",/,3(+D.3D"+ZD,/),2/,"("OUT:")",/,7(5B+D.6D"+ZD,      
        /),2/,"("LAST RESIDUAL VECTOR")",/,6(6B+3Z.D,/)")",X[1],Y[1],         
        X[2],Y[2],X[3],Y[3],X[4],Y[4],X[5],Y[5],X[6],Y[6],PAR[1],PAR[2],      
        PAR[3],OUT[7],OUT[2],OUT[6],OUT[3],OUT[4],OUT[5],OUT[1],RV[1],        
        RV[2],RV[3],RV[4],RV[5],RV[6])
     "END"        
         "EOP"    
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 7      
        
        
        
     WITH THE DATA GIVEN IN THE  X - Y - TABLE THIS PROGRAM DELIVERS:         
        
     X[I]     Y[I]
      -5     127.0
      -3     151.0
      -1     379.0
      +1     421.0
      +3     460.0
      +5     426.0
        
        
     PARAMETERS   
     +5.232" +2   
     -1.568" +2   
     -1.998" -1   
        
        
     OUT:         
          +7.220828" +7     
          +1.157156" +2     
          +1.728008" -3     
          +1.654588" +2     
          +2.300000" +1     
          +2.200000" +1     
          +0.000000" +0     
        
        
     LAST RESIDUAL VECTOR   
            -29.6 
            +86.6 
            -47.3 
            -26.2 
            -22.9 
            +39.5 
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 8      
        
        
        
 SUBSECTION  :  GSSNEWTON.  
        
        
 CALLING SEQUENCE    :      
        
     THE HEADING OF THE PROCEDURE READS :       
        
     "PROCEDURE" GSSNEWTON(M, N, PAR, RV, JJINV, FUNCT, JACOBIAN, IN,         
         OUT);    
     "VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV, JJINV, IN, OUT; 
     "BOOLEAN""PROCEDURE" FUNCT; "PROCEDURE" JACOBIAN;    
     "CODE" 34441;
        
     THE MEANING OF THE FORMAL PARAMETERS IS  : 
     M       :   <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF EQUATIONS;       
     N       :   <ARITHMETIC EXPRESSION>;       
                 THE NUMBER OF UNKNOWNS IN THE M EQUATIONS  (N <= M);         
     PAR     :   <ARRAY IDENTIFIER>; "ARRAY" PAR[1 : N];  
                 THE UNKNOWNS OF THE EQUATIONS. 
                 ENTRY :  AN APPROXIMATION TO A LEAST SQUARES SOLUTION        
                     OF THE SYSTEM.   
                 EXIT  :  THE CALULATED LEAST SQUARES SOLUTION;     
     RV      :   <ARRAY IDENTIFIER>; "ARRAY" RV[1 : M];   
                 EXIT  :  THE RESIDUAL VECTOR OF THE SYSTEM AT THE  
                     CALCULATED SOLUTION;       
     JJINV   :   <ARRAY IDENTIFIER>; "ARRAY" JJINV[1 : N,1 : N];    
                 EXIT  :  THE INVERSE OF THE MATRIX  J' * J, WHERE  J         
                     IS THE JACOBIAN MATRIX AT THE SOLUTION AND J' IS         
                     J TRANSPOSED;    
        
     FUNCT :     <PROCEDURE IDENTIFIER>;        
                 THE HEADING OF THIS PROCEDURE SHOULD BE :
        
                 "BOOLEAN""PROCEDURE" FUNCT(M, N, PAR, RV); "VALUE" M,        
                 N; "INTEGER" M, N; "ARRAY" PAR, RV;      
        
                 ENTRY: M, N, PAR;    
                      M, N  HAVE THE SAME MEANING AS IN THE PROCEDURE         
                      GSSNEWTON;      
                      "ARRAY" PAR[1:N] CONTAINS THE CURRENT VALUES OF         
                      THE UNKNOWNS AND SHOULD NOT BE ALTERED.       
                 EXIT:  "ARRAY" RV[1 : M];      
                      UPON COMPLETION OF A CALL OF FUNCT, THIS ARRAY  RV      
                      SHOULD CONTAIN THE RESIDUAL VECTOR, OBTAINED WITH       
                      THE CURRENT VALUES OF THE UNKNOWNS. 
                 THE PROGRAMMER OF FUNCT MAY DECIDE THAT SOME CURRENT         
                 ESTIMATES OF THE UNKNOWNS LIE OUTSIDE A FEASIBLE   
                 REGION; IN THIS CASE  FUNCT  SHOULD DELIVER THE VALUE        
                 FALSE AND THE PROCESS IS TERMINATED (SEE OUT[1]).  
                 OTHERWISE  FUNCT  SHOULD DELIVER THE VALUE TRUE;   
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                     PAGE 9      
        
        
        
     JACOBIAN  : <PROCEDURE IDENTIFIER>;        
                 THE HEADING OF THIS PROCEDURE SHOULD BE :
        
                 "PROCEDURE" JACOBIAN(M, N, PAR, RV, JAC, LOCFUNCT);
                 "VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV, JAC;
                 "PROCEDURE" LOCFUNCT;
        
                 THE MEANING OF THE PARAMETERS OF JACOBIAN IS :     
        
                 M, N : SEE GSSNEWTON.
                 PAR : <ARRAY IDENTIFIER>; "ARRAY" PAR[1 : N];      
                     ENTRY :  CURRENT ESTIMATES OF THE UNKNOWNS.    
                     THESE VALUES SHOULD NOT BE CHANGED.  
                 RV  : <ARRAY IDENTIFIER>; "ARRAY" RV[1 : M];       
                     ENTRY :  THE RESIDUAL VECTOR OF THE SYSTEM OF  
                     EQUATIONS CORRESPONDING TO THE VECTOR OF UNKNOWNS        
                     AS GIVEN IN PAR. 
                     EXIT  :  THE ENTRY VALUES. 
                 JAC : <ARRAY IDENTIFIER>; "ARRAY" JAC[1 : M, 1 : N];         
                     EXIT  :  THE JACOBIAN MATRIX AT THE CURRENT    
                     ESTIMATES GIVEN IN  PAR, I.E. THE MATRIX OF PARTIAL      
                     DERIVATIVES      
                     D(RV)[I] / DPAR[J], I = 1(1)M, J = 1(1)N.      
                 LOCFUNCT : <PROCEDURE IDENTIFIER>; THE HEADING OF THIS       
                     PROCEDURE IS THE SAME AS THE HEADING OF FUNCT. 
        
                 A CALL OF THE PROCEDURE JACOBIAN SHOULD DELIVER THE
                 JACOBIAN MATRIX EVALUATED WITH THE CURRENT ESTIMATES         
                 OF THE UNKNOWN VARIABLES GIVEN IN  PAR   
                 IN SUCH A WAY, THAT THE PARTIAL DERIVATIVE         
                 D(RV)[I] / DPAR[J] IS DELIVERED IN JAC[I,J], I = 1(1)M,      
                 J = 1(1)N. 
                 FOR THE CALCULATION OF THE DERIVATIVES ONE CAN USE THE       
                 VALUES OF THE CURRENT ESTIMATES OF THE   
                 UNKNOWNS AS GIVEN IN PAR AND THE RESIDUAL VECTOR AS
                 GIVEN IN  RV.        
                 ONE CAN ALSO USE THE PROCEDURE  FUNCT    
                 (PARAMETER OF GSSNEWTON) THROUGH CALLS OF THE PROCEDURE      
                 LOCFUNCT (PARAMETER OF JACOBIAN). THIS PARAMETER OF
                 JACOBIAN MAY BE USED WHEN THE JACOBIAN MATRIX IS   
                 APPROXIMATED USING (FORWARD) DIFFERENCES.
                 AN APPROPRIATE PROCEDURE TO THIS PURPOSE IS  JACOBNMF        
                 (SECTION 4.3.2.1). SUCH A PROCEDURE MAY BE USED ONLY IF      
                 THE MATRIX ELEMENTS ARE COMPUTED SUFFICIENTLY ACCURATE;      
        
     IN      :   <ARRAY IDENTIFIER>; "ARRAY" IN[0 : 7];   
                 IN THIS ARRAY TOLERANCES AND CONTROL PARAMETERS SHOULD       
                 BE GIVEN.  
                 ENTRY  :   
                 IN[0] : THE MACHINE PRECISION. FOR CALCULATION ON THE        
                     CYBER 73 A SUITABLE VALUE IS "-14.   
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 10      
        
        
        
                 IN[1], IN[2] :       
                     RELATIVE AND ABSOLUTE TOLERANCE FOR THE STEP VECTOR      
                     (RELATIVE TO THE VECTOR OF CURRENT ESTIMATES IN
                     PAR).  
                     THE PROCESS IS TERMINATED IF IN SOME ITERATION (BUT      
                     NOT THE FIRST) THE EUCLIDEAN NORM OF THE CALCULATED      
                     NEWTON STEP IS LESS THAN IN[1] * NORM(PAR) + IN[2].      
                     IN[1] SHOULD NOT BE CHOSEN SMALLER THAN IN[0]. 
                 IN[3] IS NOT USED BY THE PROCEDURE  GSSNEWTON;     
                 IN[4] : ABSOLUTE TOLERANCE FOR THE EUCLIDEAN NORM OF         
                     THE RESIDUAL VECTOR. THE PROCESS IS TERMINATED WHEN      
                     THIS NORM IS LESS THAN IN[4].        
                 IN[5] : THE MAXIMUM ALLOWED NUMBER OF FUNCTION     
                     EVALUATIONS (I.E. CALLS OF FUNCT).   
                 IN[6] : THE MAXIMUM ALLOWED NUMBER OF HALVINGS OF A
                     CALCULATED NEWTON STEP VECTOR (SEE METHOD AND  
                     PERFORMANCE). A SUITABLE VALUE IS 15.
                 IN[7] : THE MAXIMUM ALLOWED NUMBER OF SUCCESSIVE IN[6]       
                     TIMES HALVED STEP VECTORS. SUITABLE VALUES ARE 1         
                     AND 2; 
        
     OUT     :   <ARRAY IDENTIFIER>; "ARRAY" OUT[1 : 9];  
                 IN ARRAY OUT INFORMATION ABOUT THE TERMINATION OF THE        
                 PROCESS IS DELIVERED.
                 EXIT  :    
                 OUT[1] :   
                     THE PROCESS WAS TERMINATED BECAUSE (OUT[1] = ) 
                   1.THE NORM OF THE RESIDUAL VECTOR IS SMALL WITH  
                     RESPECT TO IN[4],
                   2.THE CALCULATED NEWTON STEP IS SUFFICIENTLY SMALL         
                     (SEE IN[1], IN[2]),        
                   3.THE CALCULATED STEP WAS COMPLETELY DAMPED (HALVED)       
                     IN  IN[7]  SUCCESSIVE ITERATIONS,    
                   4.OUT[4] EXCEEDS IN[5], THE MAXIMUM ALLOWED NUMBER OF      
                     CALLS OF FUNCT,  
                   5.THE JACOBIAN WAS NOT FULL-RANK (SEE OUT[8]),   
                   6.FUNCT DELIVERED FALSE AT A NEW VECTOR OF       
                     ESTIMATES OF THE UNKNOWNS, 
                   7.FUNCT DELIVERED FALSE IN A CALL FROM JACOBIAN. 
                 OUT[2] : THE EUCLIDEAN NORM OF THE LAST RESIDUAL   
                     VECTOR.
                 OUT[3] : THE EUCLIDEAN NORM OF THE INITIAL RESIDUAL
                     VECTOR.
                 OUT[4] : THE TOTAL NUMBER OF CALLS OF FUNCT.       
                     OUT[4] WILL BE LESS THAN IN[5] + IN[6].        
                 OUT[5] : THE TOTAL NUMBER OF ITERATIONS. 
                 OUT[6] : THE EUCLIDEAN NORM OF THE LAST STEP VECTOR.         
                 OUT[7] : ITERATION NUMBER OF THE LAST ITERATION IN 
                     WHICH THE NEWTON STEP WAS HALVED.    
                 OUT[8], OUT[9] :     
                     RANK AND MAXIMUM COLUMN NORM OF THE JACOBIAN MATRIX      
                     IN THE LAST ITERATION, AS DELIVERED BY  LSQORTDEC        
                     (SECTION 3.1.1.2.1.1) IN AUX[3] AND AUX[5].    
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 11      
        
        
        
 DATA AND RESULTS    :      
        
     THE PROCEDURE GSSNEWTON CAN BE USED FOR APPROXIMATING AN EXACT OR A      
     LEAST SQUARES SOLUTION OF A SYSTEM OF NONLINEAR EQUATIONS. WHEN AN       
     EXACT SOLUTION IS REQUIRED, THE PROCEDURE MAY TERMINATE ONLY WITH        
     OUT[1] = 1, AND VERY SMALL VALUES SHOULD BE ASSIGNED TO IN[1] AND        
     IN[2]. WHEN A LEAST SQUARES SOLUTION IS REQUIRED, POSITIVE RESULTS       
     OF THE PROCEDURE ARE SIGNALED BY OUT[1] = 1 OR 2. WHENEVER THE 
     PROCEDURE TERMINATES WITH OUT[1] < 5, THEN THE INVERSE OF  J' * J        
     (SEE MEANING OF THE PARAMETER  JJINV) IS DELIVERED IN  JJINV. IN         
     THAT CASE THE COVARIANCE MATRIX AND THE STANDARD DEVIATIONS OF THE       
     SOLUTION CAN BE CALCULATED.      
        
     FOR A CURVE FITTING PROBLEM, SAY :         
        "ESTIMATE PARAMETERS  PAR[1], ... , PAR[N] OF A FUNCTION    
        "Y = F(X; PAR[1], ... , PAR[N]), WHEN A SET OF DATA (X[I],Y[I]),      
        "I = 1(1)M, HAS TO BE FITTED,"
     THE FOLLOWING SYSTEM OF  M  EQUATIONS IN THE  N  UNKNOWN PARAMETERS      
     PAR[1], ... , PAR[N] CAN BE DERIVED :      
        
     F(X[I]; PAR[1], ... , PAR[N]) - Y[I] = 0, I = 1(1)M. 
        
 PROCEDURES USED :
        
     VECVEC    = CP34010,   
     DUPVEC    = CP31030,   
     ELMVEC    = CP34020,   
     LSQORTDEC = CP34134,   
     LSQSOL    = CP34131,   
     LSQINV    = CP34136.   
        
 REQUIRED CENTRAL MEMORY  : 
     EXECUTION FIELD LENGTH :  AN ARRAY OF  (M + 1) * N ELEMENTS, FOUR        
     ARRAYS OF  N ELEMENTS AND ONE ARRAY OF  M ELEMENTS ARE DECLARED.         
        
 RUNNING TIME    :
     THE RUNNING TIME IS PROPORTIONAL TO THE TOTAL NUMBER OF CALLS OF         
     FUNCT. BESIDES, IN EACH ITERATION A LINEAR LEAST SQUARES SYSTEM
     IS SOLVED (NUMBER OF OPERATIONS PROPORTIONAL TO  M * (N SQUARED)).       
        
 LANGUAGE    :  ALGOL 60.   
        
        
 METHOD AND PERFORMANCE :   
     THE PROCEDURE GSSNEWTON IS BASED UPON THE GAUSS-NEWTON METHOD FOR        
     CALCULATING A LEAST SQUARES SOLUTION OF A SYSTEM OF NONLINEAR  
     EQUATIONS (SEE E.G. [1], [2]). IN SEVERAL ITERATIONS A STEP VECTOR       
     (FOR THE VECTOR OF UNKNOWNS) IS OBTAINED BY SOLVING A LINEARIZED         
     SYSTEM OF EQUATIONS. THIS STEP IS PERFORMED ONLY IF THE NORM OF THE      
     RESIDUAL VECTOR DECREASES. OTHERWISE THE NEWTON STEP VECTOR IS 
     HALVED A NUMBER OF TIMES UNTIL THE NORM OF THE RESIDUAL VECTOR IS        
     SMALLER THAN THE NORMS OF THE LAST AND SUBSEQUENT RESIDUAL VECTORS       
     (THIS PROCESS IS CALLED STEP SIZE CONTROL).
     FOR FURTHER DETAILS SEE [3] (SEE ALSO [2]).
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 12      
        
        
        
 REFERENCES  :    
        
     [1] H.O. HARTLEY :     
         THE MODIFIED GAUSS-NEWTON METHOD.      
         TECHNOMETRICS, V.3 (1961), PP. 269 - 280.        
     [2] H. SPAETH :        
         THE DAMPED TAYLOR'S SERIES METHOD FOR MINIMIZING A SUM OF  
         SQUARES AND FOR SOLVING SYSTEMS OF NONLINEAR EQUATIONS.    
         ALGORITHM 315, COLLECTED ALGORITHMS FROM  CACM,  
         COMMUNICATIONS OF THE ACM, VOL. 10 (NOV. 1967), PP. 726 - 728.       
     [3] J.C.P. BUS, B. VAN DOMSELAAR, J. KOK : 
         NONLINEAR LEAST SQUARES ESTIMATION.    
         MATHEMATICAL CENTRE (TO APPEAR).       
        
        
 EXAMPLE OF USE  :
        
     THE PARAMETERS PAR[1 : 3] IN THE CURVE FITTING PROBLEM:        
     G[I] = PAR[1] + PAR[2] * EXP(PAR[3] * X[I]) - Y[I]   
     WITH  (X[I], Y[I]), I=1,...,6  AS THE EXPERIMENTAL DATA, MAY BE
     OBTAINED BY THE FOLLOWING PROGRAM:         
        
     "BEGIN""PROCEDURE" GSSNEWTON(M, N, P, F, C, TF, JAC, I, O);    
            "CODE" 34441;   
        "INTEGER" I;        
        "ARRAY" IN[0:7], OUT[1:9], X, Y, G[1:6], V[1:3, 1:3], PAR[1:3];       
        
        "BOOLEAN""PROCEDURE" EXPFUNCT(M, N, PAR, G);      
        "VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, G;     
        "BEGIN""INTEGER" I; 
            "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" 
            "BEGIN""IF" PAR[3] * X[I] > 680 "THEN"        
                "BEGIN" EXPFUNCT:= "FALSE"; "GO TO" STOP "END";     
                G[I]:= PAR[1] + PAR[2] * EXP(PAR[3] * X[I]) - Y[I]  
            "END"; EXPFUNCT:="TRUE";  
        STOP:     
        "END" EXPFUNCT;     
        
        "PROCEDURE" JACOBIAN(M, N, PAR, G, JAC, LOCFUNCT);
        "VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, G, JAC;
        "PROCEDURE" LOCFUNCT;         
        "BEGIN" "INTEGER" I; "REAL" EX;         
            "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO" 
            "BEGIN" JAC[I,1]:=1; JAC[I,2]:= EX:= EXP(PAR[3] * X[I]);
                JAC[I,3]:= X[I] * PAR[2] * EX   
            "END" 
        "END" JACOBIAN;     
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 13      
        
        
        
        IN[0]:= "-14; IN[1]:= IN[2]:= "-6; IN[5]:= 75; IN[4]:="-10; 
        IN[6]:= 14; IN[7]:= 1;        
        "FOR" I:= 1 "STEP" 1 "UNTIL" 6 "DO"     
        INPUT(1, "("2(N),/")",X[I],Y[I]);       
        PAR[1]:= 580; PAR[2]:= - 180; PAR[3]:= - 0.160;   
        GSSNEWTON(6, 3, PAR, G, V, EXPFUNCT, JACOBIAN, IN, OUT);    
        OUTPUT(61,"("3/4B,"("X[I]     Y[I]")",/, 6(5B+D, 5B3D.D, /), 2/,      
        4B"("PARAMETERS")",/,3(4B+D.3D"+ZD,/),2/,4B"("OUT:")",/,    
        3(9B+D.6D"+ZD,/), 5(14B3ZD,/), 9B+D.6D"+ZD,2/4B,  
        "("LAST RESIDUAL VECTOR")",/,6(10B+3Z.D,/)")", X[1], Y[1],  
        X[2],Y[2],X[3],Y[3],X[4],Y[4],X[5],Y[5],X[6],Y[6],PAR[1],PAR[2],      
        PAR[3],OUT[6],OUT[2],OUT[3],OUT[4],OUT[5],OUT[1],OUT[7],OUT[8],       
        OUT[9], G[1],G[2],G[3],G[4],G[5],G[6])  
     "END"        
         "EOP"    
        
      WITH THE DATA GIVEN IN THE  X - Y - TABLE THIS PROGRAM DELIVERS:        
        
      X[I]     Y[I]         
       -5     127.0         
       -3     151.0         
       -1     379.0         
       +1     421.0         
       +3     460.0         
       +5     426.0         
        
        
      PARAMETERS  
      +5.233" +2  
      -1.569" +2  
      -1.997" -1  
        
        
      OUT:        
           +5.260478" -4    
           +1.157156" +2    
           +1.654588" +2    
                  16        
                  16        
                   2        
                   0        
                   3        
           +2.339529" +3    
        
      LAST RESIDUAL VECTOR  
             -29.6
             +86.6
             -47.3
             -26.2
             -22.9
             +39.5
        
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 14      
        
        
        
 SOURCE TEXTS    :
0"CODE" 34440;    
  "PROCEDURE" MARQUARDT(M,N,PAR,G,V,FUNCT,JACOBIAN,IN,OUT);         
  "VALUE" M,N; "INTEGER" M,N; "ARRAY" PAR,G,V,IN,OUT;     
  "BOOLEAN" "PROCEDURE" FUNCT; "PROCEDURE" JACOBIAN;      
  "BEGIN" "INTEGER" MAXFE,FE,IT,I,J,ERR;        
          "REAL" VV,WW,W,MU,RES,FPAR,FPARPRES,LAMBDA,LAMBDAMIN,     
                 P,PW,RELTOLRES,ABSTOLRES;      
          "ARRAY" EM[0:7],VAL,B,BB,PARPRES[1:N],JAC[1:M,1:N];       
        
          "PROCEDURE" MULCOL(L,U,S,T,A,B,X); "CODE" 31022;
          "PROCEDURE" DUPVEC(L,U,S,A,B); "CODE" 31030;    
          "REAL" "PROCEDURE" VECVEC(L,U,S,A,B); "CODE" 34010;       
          "REAL" "PROCEDURE" MATVEC(L,U,S,A,B); "CODE" 34011;       
          "REAL" "PROCEDURE" TAMVEC(L,U,S,A,B); "CODE" 34012;       
          "REAL" "PROCEDURE" MATTAM(L,U,S,T,A,B); "CODE" 34015;     
          "INTEGER" "PROCEDURE" QRISNGVALDEC(A,M,N,VAL,V,EM);       
          "CODE" 34273;     
        
         "PROCEDURE" LOCFUNCT(M,N,PAR,G);       
         "INTEGER" M,N; "ARRAY" PAR,G;
         "BEGIN" FE:= FE+1; "IF" FE >= MAXFE "THEN" ERR:= 1 "ELSE"  
                 "IF" "NOT" FUNCT(M,N,PAR,G) "THEN" ERR:= 2;        
                 "IF" ERR^=0 "THEN" "GOTO" EXIT 
         "END" LOCFUNCT;    
        
         VV:=10; W:=0.5; MU:= 0.01;   
         WW:=("IF" IN[6]<"-7 "THEN" "-8 "ELSE" "-1*IN[6]);
         EM[0]:=EM[2]:=EM[6]:=IN[0]; EM[4]:=10*N;         
         RELTOLRES:=IN[3]; ABSTOLRES:=IN[4]**2; MAXFE:=IN[5];       
         ERR:= 0; FE:= IT:= 1; P:=FPAR:= RES:= 0;         
         PW:=-LN(WW*IN[0])/2.30;      
        
         "IF" "NOT" FUNCT(M,N,PAR,G) "THEN"     
         "BEGIN" ERR:= 3; "GOTO" ESCAPE "END";  
         FPAR:= VECVEC(1,M,0,G,G); OUT[3]:=SQRT(FPAR);    
        
         "FOR" IT:= 1, IT+1 "WHILE" FPAR > ABSTOLRES "AND"
                       RES > RELTOLRES*FPAR+ABSTOLRES "DO"
         "BEGIN" JACOBIAN(M,N,PAR,G,JAC,LOCFUNCT);        
                 I:=QRISNGVALDEC(JAC,M,N,VAL,V,EM);       
                 "IF" IT=1 "THEN"     
                       LAMBDA:= IN[6] * VECVEC(1,N,0,VAL,VAL) "ELSE"
                 "IF" P =0 "THEN" LAMBDA:= LAMBDA*W "ELSE" P:= 0;   
        
                 "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"       
                 B[I]:=VAL[I]*TAMVEC(1,M,I,JAC,G);        
                                                               "COMMENT"      
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 15      
                                                                 ;  
        
        
            L:   "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"       
                 BB[I]:=B[I]/(VAL[I]*VAL[I]+LAMBDA);      
                 "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"       
                 PARPRES[I]:= PAR[I] - MATVEC(1,N,I,V,BB);
                 LOCFUNCT(M,N,PARPRES,G);       
                 FPARPRES:= VECVEC(1,M,0,G,G);  
                 RES:=FPAR-FPARPRES;  
                 "IF" RES < MU * VECVEC(1,N,0,B,BB) "THEN"
                 "BEGIN" P:= P+1; LAMBDA:= VV * LAMBDA;   
                     "IF" P=1 "THEN"  
                     "BEGIN" LAMBDAMIN:= WW * VECVEC(1,N,0,VAL,VAL);
                         "IF" LAMBDA<LAMBDAMIN "THEN" LAMBDA:= LAMBDAMIN      
                     "END"; 
                     "IF" P<PW "THEN" "GOTO" L "ELSE"     
                     "BEGIN" ERR:= 4; 
                             "GOTO" EXIT        
                     "END"; 
                 "END";     
        
                 DUPVEC(1,N,0,PAR,PARPRES);     
                 FPAR:=FPARPRES       
           "END" ITERATION; 
        
      EXIT:       
           "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"   
           MULCOL(1,N,I,I,JAC,V,1/(VAL[I]+IN[0]));        
           "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"   
           "FOR" J:=1 "STEP" 1 "UNTIL" I "DO"   
           V[I,J]:= V[J,I]:= MATTAM(1,N,I,J,JAC,JAC);     
        
           LAMBDA:= LAMBDAMIN:= VAL[1];         
           "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO"  
           "IF" VAL[I]>LAMBDA    "THEN" LAMBDA   := VAL[I] "ELSE"   
           "IF" VAL[I]<LAMBDAMIN "THEN" LAMBDAMIN:= VAL[I];         
        
           OUT[7]:=(LAMBDA/(LAMBDAMIN+IN[0]))**2;         
           OUT[2]:=SQRT(FPAR);        
           OUT[6]:=SQRT(RES+FPAR)-OUT[2];       
      ESCAPE:     
           OUT[4]:=FE;      
           OUT[5]:=IT-1;    
           OUT[1]:=ERR      
      "END" MARQUARDT;      
         "EOP"    
        
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 16      
        
        
        
     "CODE" 34441;
     "PROCEDURE" GSSNEWTON(M, N, PAR, RV, JJINV, FUNCT, JACOBIAN,   
       IN, OUT);  
     "VALUE" M, N; "INTEGER" M, N;    
     "ARRAY" PAR, RV, JJINV, IN, OUT; 
     "BOOLEAN" "PROCEDURE" FUNCT;     
     "PROCEDURE" JACOBIAN;  
        
     "BEGIN" "INTEGER" I, J, INR, MIT, TEXT,    
         IT, ITMAX, INRMAX, TIM, FEVAL, FEVALMAX;         
         "REAL" RHO, RES1, RES2, RN, RELTOLPAR, ABSTOLPAR, ABSTOLRES,         
            STAP, NORMX;    
         "BOOLEAN" CONV, TESTTHF, DAMPING ON;   
         "ARRAY" JAC[1:M + 1,1:N], PR, AID, SOL[1 : N], FU2[1 : M], 
            AUX[2 : 5];     
         "INTEGER""ARRAY" CI[1:N];    
        
         "REAL""PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010; 
         "PROCEDURE" DUPVEC(L, U, S, A, B); "CODE" 31030; 
         "PROCEDURE" ELMVEC(L, U, S, A, B, X); "CODE" 34020;        
         "PROCEDURE" LSQORTDEC(A, M, N, AUX, AID, CI); "CODE" 34134;
         "PROCEDURE" LSQSOL(A, M, N, AID, CI, B); "CODE" 34131;     
         "PROCEDURE" LSQINV(A, N, AID, CI); "CODE" 34136; 
        
         "BOOLEAN""PROCEDURE" LOC FUNCT(M, N, PAR, RV);   
         "VALUE" M, N; "INTEGER" M, N; "ARRAY" PAR, RV;   
         "BEGIN" LOC FUNCT:= TEST THF:= FUNCT(M, N, PAR, RV)        
             "AND" TEST THF; FEVAL:= FEVAL + 1  
         "END" LOC FUNCT;   
        
        
         ITMAX:= FEVALMAX:= IN[5]; AUX[2]:= N * IN[0]; TIM:= IN[7]; 
         RELTOLPAR:= IN[1] ** 2; ABSTOLPAR:= IN[2] ** 2;  
         ABSTOLRES:= IN[4] ** 2; INRMAX:= IN[6];
         DUPVEC(1, N, 0, PR, PAR);    
         "IF" M < N "THEN"  
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" JAC[M + 1, I]:= 0;     
         TEXT:= 4; MIT:= 0; TEST THF:= "TRUE";  
         RES2:= STAP:= OUT[5]:= OUT[6]:= OUT[7]:= 0;      
         FUNCT(M, N, PAR, FU2); RN:= VECVEC(1, M, 0, FU2, FU2);     
         OUT[3]:= SQRT(RN); FEVAL:= 1; DAMPING ON:= "FALSE";        
         "FOR" IT:= 1, IT + 1 "WHILE" IT <= ITMAX "AND"   
                FEVAL < FEVALMAX "DO" 
         "BEGIN" OUT[5]:= IT; JACOBIAN(M, N, PAR, FU2, JAC, LOCFUNCT);        
            "IF" "NOT" TEST THF "THEN"
            "BEGIN" TEXT:= 7; "GO TO" FAIL "END";         
            LSQORTDEC(JAC, M, N, AUX, AID, CI); 
            "IF" AUX[3] ^= N "THEN"   
            "BEGIN" TEXT:= 5; "GO TO" FAIL "END";         
            LSQSOL(JAC, M, N, AID, CI, FU2); DUPVEC(1, N, 0, SOL, FU2);       
            STAP:= VECVEC(1, N, 0, SOL, SOL);   
            RHO:= 2; NORMX:= VECVEC(1, N, 0, PAR, PAR);   
                                                               "COMMENT"      
1SECTION : 5.1.3.1.3         (DECEMBER 1975)                    PAGE 17      
                                                                 ;  
        
        
            "IF" STAP > RELTOLPAR * NORMX + ABSTOLPAR     
               "OR" IT = 1 "AND" STAP > 0 "THEN"
            "BEGIN""FOR" INR:= 0, INR + 1       
               "WHILE""IF" INR = 1 "THEN" DAMPING ON "OR" RES2 >= RN
               "ELSE""NOT" CONV "AND" (RN <= RES1 "OR" RES2 < RES1) "DO"      
               "BEGIN""COMMENT" DAMPING STOPS WHEN        
                  R0 > R1 "AND" R1 <= R2 (BEST RESULT IS X1, R1)    
                  WITH  X1 = X0 + I * DX, I:= 1, .5, .25, .125, ETC.  ;       
                  RHO:= RHO / 2; "IF" INR > 0 "THEN"      
                  "BEGIN" RES1:= RES2; DUPVEC(1, M, 0, RV, FU2);    
                     DAMPING ON:= INR > 1       
                  "END";    
                  "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"     
                  PR[I]:= PAR[I] - SOL[I] * RHO;
                  FEVAL:= FEVAL + 1;  
                  "IF" "NOT" FUNCT(M, N, PR, FU2) "THEN"  
                  "BEGIN" TEXT:= 6; "GO TO" FAIL "END";   
                  RES2:= VECVEC(1, M, 0, FU2, FU2); CONV:= INR >= INRMAX      
               "END" DAMPING OF STEP VECTOR;    
               "IF" CONV "THEN"       
               "BEGIN""COMMENT" RESIDUE CONSTANT; MIT:= MIT + 1;    
                  "IF" MIT < TIM "THEN" CONV:= "FALSE"    
               "END" "ELSE" MIT:= 0;  
               "IF" INR > 1 "THEN"    
               "BEGIN" RHO:= RHO * 2; ELMVEC(1, N, 0, PAR, SOL, - RHO);       
                  RN:= RES1; "IF" INR > 2 "THEN" OUT[7]:= IT        
               "END""ELSE"  
               "BEGIN" DUPVEC(1, N, 0, PAR, PR); RN:= RES2;         
                  DUPVEC(1, M, 0, RV, FU2)      
               "END";       
        
               "IF" RN <= ABSTOLRES "THEN"      
               "BEGIN" TEXT:= 1; ITMAX:= IT "END""ELSE"   
               "IF" CONV "AND" INRMAX > 0 "THEN"
               "BEGIN" TEXT:= 3; ITMAX:= IT "END"         
               "ELSE" DUPVEC(1, M, 0, FU2, RV)  
            "END" ITERATION WITH DAMPING AND TESTS "ELSE" 
            "BEGIN" TEXT:= 2; RHO:= 1; ITMAX:= IT "END"   
         "END" OF ITERATIONS;         
        
         LSQINV(JAC, N, AID, CI);     
         "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"    
         "BEGIN" JJINV[I,I]:= JAC[I,I];         
             "FOR" J:= I + 1 "STEP" 1 "UNTIL" N "DO"      
             JJINV[I,J]:= JJINV[J,I]:= JAC[I,J] 
         "END" CALCULATION OF INVERSE MATRIX OF NORMAL EQUATIONS;   
     FAIL :       
         OUT[6]:= SQRT(STAP) * RHO; OUT[2]:= SQRT(RN); OUT[4]:= FEVAL;        
         OUT[1]:= TEXT; OUT[8]:= AUX[3]; OUT[9]:= AUX[5]  
     "END" GSSNEWTON;       
         "EOP"