\ Cauchy principal value integration using
\ 6-point Gauss-Legendre quadrature
\


\ ---------------------------------------------------
\     (c) Copyright 1999  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------

\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependencies:
\   Assumes independent floating point stack
\   Uses a FORmula TRANslator for clarity
\
\ Description of algorithm:
\   integrate f(x)/(x - x0) from x0-delta to x0+delta
\   using even-order Gauss-Legendre quadrature formula

\ Usage example:
\   15 set-precision  ok
\   use( fexp 1e0 0e0 )pval fs. 2.11450175075134E0  ok

MARKER  -pval

: undefined    BL WORD  FIND  NIP  0=  ;

undefined  f"  [IF]  include  ftran111.f  [THEN]

FVARIABLE wt1    FVARIABLE xi1      FVARIABLE y1    FVARIABLE y1m
FVARIABLE wt2    FVARIABLE xi2      FVARIABLE y2    FVARIABLE y2m
FVARIABLE wt3    FVARIABLE xi3      FVARIABLE y3    FVARIABLE y3m

FVARIABLE x0     FVARIABLE delta

0.238619186083197e0  xi1 F!
0.661209386466265e0  xi2 F!
0.932469514203152e0  xi3 F!

0.467913934572691e0  wt1 F!
0.360761573048139e0  wt2 F!
0.171324492379170e0  wt3 F!

v: fdummy

: )pval    ( xt -- )   ( f: delta x0 -- integral)
        defines fdummy          \ vector function name
        x0 F!   delta F!
        f" y1 = x0 + delta*xi1" \ scale interval
        f" y2 = x0 + delta*xi2"
        f" y3 = x0 + delta*xi3"
        f" y1m = x0 - delta*xi1"
        f" y2m = x0 - delta*xi2"
        f" y3m = x0 - delta*xi3"

        f" wt1 * ( fdummy(y1) - fdummy(y1m) ) / xi1 +
           wt2 * ( fdummy(y2) - fdummy(y2m) ) / xi2 +
           wt3 * ( fdummy(y3) - fdummy(y3m) ) / xi3 "
;