(49G) Gauss-Lobatto-Kronrod 4-7-13 Integration

[Gerald H]

For input

Programme taking one argument from & returning one argument to the stack
Lower bound of integration
Upper bound of integration

the programme GK returns a numerical integral accurate to screen display setting.

eg For STD display

<< LN INV >>
2
98765432167

the programme returns

4069312104.01

in 33.12 sec

& for 4 FIX

406959338.70

in 5.74 sec.

For sub-programmes see below GK.

Code:


GK

::
  CK3&Dispatch
  # 811
  ::
    ID x0AE
    4UNROLL
    CK3
    3DROP
  ;
;

x0AE

::
  2DUP
  2DUPSWAP
  2DUP
  %+
  %2
  %/
  3UNROLL
  %-
  %2
  %/
  4PICK
  8PICK
  EVAL
  DUPTYPEREAL?
  NcaseTYPEERR
  4PICK
  9PICK_
  EVAL
  %0
  DUPDUP
  '
  NULLLAM
  BINT9
  NDUPN
  DOBIND
  % .942882415695
  % .816496580928
  % .641853342346
  % .4472135955
  % .236383199662
  %0
  BINT6
  ONE_DO
  INDEX@
  #2*
  PICK
  %CHS
  LOOP
  BINT11
  ZERO_DO
  6GETLAM
  %*
  7GETLAM
  %+
  BINT14
  PICK
  EVAL
  BINT11
  UNROLL
  LOOP
  8PICK
  5PICK
  %+
  %5
  %*
  5GETLAM
  4GETLAM
  %+
  DUP
  7PUTLAM
  %+
  6GETLAM
  %*
  %6
  %/
  3PUTLAM
  10PICK_
  3PICK
  %+
  % 432.
  %*
  9PICK_
  6PICK
  %+
  % 625.
  %*
  %+
  7PICK
  % 672.
  %*
  %+
  7GETLAM
  % 77.
  %*
  %+
  6GETLAM
  %*
  % 1470.
  %/
  2PUTLAM
  BINT11
  ROLL
  %+
  % 9.42738402189E-2
  %*
  SWAP
  BINT10
  ROLL
  %+
  % .155071987337
  %*
  %+SWAP
  8ROLL
  %+
  % .18882157396
  %*
  %+SWAP
  6ROLL
  %+
  % .199773405227
  %*
  %+
  SWAP4ROLL
  %+
  % .224926465333
  %*
  %+SWAP
  % .242611071901
  %*
  %+
  7GETLAM
  % 1.58271919735E-2
  %*
  %+
  6GETLAM
  %*
  DUP1PUTLAM_
  2GETLAM
  %-
  %ABS
  3GETLAM
  1GETLAM
  %-
  %ABS
  DUP%0=
  ITE
  ::
    2DROP
    %1
  ;
  %/
  ID x0B1
  % .00000000001
  %MAX
  SWAPDUP
  %0>
  OVER
  %1
  %<
  ANDITE
  %/
  DROP
  9GETLAM
  8GETLAM
  5GETLAM
  4GETLAM
  1GETLAM
  DUP%0=
  IT
  ::
    DROP
    8GETLAM
    9GETLAM
    %-
  ;
  %ABS
  ABND
  9PICK_
  7ROLL
  TRUE
  {
    LAM 3
    LAM 2
    LAM 1
  }
  BIND
  ID x0AF
  1GETABND
  ?SEMI
  ID x0B0
  5UNROLL
;

x0B1

::
  %.1
  BINT0
  BINT4
  ZERO_DO
  #2*
  BINT48
  INDEX@#-
  TestSysFlag
  IT
  #1+
  LOOP
  ::
    DUP#0=csedrp
    BINT12
  ;
  FPTR2 ^RP#
;

x0AF

::
  4PICK
  6PICK
  2DUP
  %+
  %2
  %/
  3UNROLL
  %-
  %2
  %/
  2DUP
  2DUP
  % .816496580928
  %*
  2DUP
  %-
  5UNROLL
  %+
  3UNROLL
  % .4472135955
  %*
  2DUP
  %-
  4UNROLL
  %+
  4PICK
  4PICK
  8PICK
  5PICK
  5PICK
  BINT5
  ZERO_DO
  LAM 3
  EVAL
  5UNROLL
  LOOP
  '
  NULLLAM
  BINT16
  NDUPN
  DOBIND
  14GETLAM
  13GETLAM
  %+
  4GETLAM
  1GETLAM
  %+
  2DUP
  %5
  %*
  %+
  10GETLAM
  %*
  %6
  %/
  SWAP
  % 625.
  %*
  ROT
  % 77.
  %*
  %+
  5GETLAM
  2GETLAM
  %+
  % 432.
  %*
  %+
  3GETLAM
  % 672.
  %*
  %+
  10GETLAM
  %*
  % 1470.
  %/
  DUPROT
  %-
  %ABS
  12GETLAM
  LAM 2
  %*
  %<=
  9GETLAM
  16GETLAM
  %<=
  OR
  15GETLAM
  7GETLAM
  %<=
  DUPUNROT
  ORcase
  ::
    ABND
    LAM 1
    AND
    NOT?SEMI
    ID x0B0
    DISPROW5
    FALSE'
    LAM 1
    STO
  ;
  2DROP
  16GETLAM
  9GETLAM
  14GETLAM
  5GETLAM
  9GETLAM
  8GETLAM
  5GETLAM
  4GETLAM
  8GETLAM
  11GETLAM
  4GETLAM
  3GETLAM
  11GETLAM
  6GETLAM
  3GETLAM
  1GETLAM
  6GETLAM
  7GETLAM
  1GETLAM
  2GETLAM
  7GETLAM
  15GETLAM
  2GETLAM
  13GETLAM
  %0
  BINT6
  ZERO_DO
  5UNROLL
  12GETLAM
  ID x0AF
  %+
  LOOP
  ABND
;

x0B0

"Out of Tolerance"