This one is the version to be assembled an put into LIB or RAM.
Code:
// Double Exponential Integration for the wp34s calculator
//
// v1.0l-386 (20170323) by M. César Rodríguez GGPL
//
// USAGE:
// Inputs:
// - integration limits in X (upper limit) and Y (lower limit)
// - name of integrand program in the alpha register
// Outputs:
// - approximate value of integral in X
// - updates L
// - corrupts stack (works in both SSIZE8 and SSIZE4)
//
// Remarks:
// - accepts +infinite and -infinite as integration limits. If none
// of the limits is infinite, the method applied is the tanh-sinh
// one. If both are infinite then the sinh-sinh method is selected.
// Otherwise the exp-sinh method is used
// - the program tries to get all displayed digits correct, so
// increasing the number of displayed digits will (hopefully) give
// closer approximations (needing more time). But, in DBLON mode, the
// program ignores the number of displayed digits and does its best
// to get as much correct digits as possible
// - during program execution the succesive aproximations are
// displayed (as with the built-in integration program). Pressing
// [<-] will interrupt the program & the last approximation
// will be returned
// - all internal computations are carried out with flag D set (no
// matter the value the user set to it), thus many errors when
// evaluating the integrand (say 0/0 or 1/0, overflows or ...) will
// not raise an error, but return infinite or NaN instead. Such
// conditions will be detected by this program that will continue
// the computation simply ignoring the offending point
// - keep in mind that, when both integration limits are infinite (or
// one is infinite and the other 0) the integrand may be evaluated
// at really bigs/small points (up to 1e±201 in DBLOFF mode and
// 1e±3088 in DBLON mode). Ensure that your integrand program
// behaves as expected for such arguments (remember that this program
// works with flag D set, as stated previously)
// - the double exponential method relies on the function to be
// integrated being analytic over the integration interval, except,
// perhaps, at the interval ends. Thus, if it is known that the
// integrand, _or any of its derivatives_, has a discontinuity at some
// point inside the integration interval, it is advised to break the
// interval at such point and compute two separate integrals, one
// at each side of the problematic point. Better results will be
// get this way. Thus, beware of abs, fp, ip and such non analytic
// functions inside the integrand program. Other methods (Romberg)
// may behave better in this respect, but, on the other side, the
// double exponential method manages discontinuities (of the
// integrand or its derivatives) at the interval edges better and
// is usually way faster that the Romberg method
// - as many other quadrature methods (Romberg included), the
// double exponential algorithm does not manage well highly
// oscillatory integrands. Here, highly oscillatory means that the
// integrand changes sign many (or infinite) times along the
// integration interval
//
// Method used:
// The algorith used here is adapted from that described in:
// H. Takahasi and M. Mori. "Double Exponential Formulas for Numerical
// Integration." Publications of RIMS, Kyoto University 9 (1974),
// 721–741.
// Another useful reference may be:
// David H. Bailey, Xiaoye S. Li and Karthik Jeyabalan, "A comparison
// of three high-precision quadrature schemes," Experimental
// Mathematics, vol. 14 (2005), no. 3, pg 317-329.
//
// 252 steps
// 14 local registers
LBL'DEI' // Double-Exponential Integration
INTM? // check for invalid mode
ERR 13 // ERR_BAD_MODE
LocR 14 // Local registers and flags:
// .00 X2L to return it in L
// .01 bma2 (b - a)/2, a & b are the
// integration limits
// .02 bpa2 (b + 2)/2
// .03 eps epsilon
// .04 thr convergence threshold
// .05 lvl level counter
// .06 tm maximum abscissa in the t domain
// .07 h space between abscissas in
// the t domain
// .08 ssp partial sum at each level
// .09 j abscissa counter
// .10 ch cosh(t)
// .11 w weight
// .12 rp abscissa (r) or fplus + fminus (p)
// .13 ss value of integral
// F.00 DF backup of D flag
// F.01 rev b is < a
// F.02 ES expsinh mode
// F.03 bpa2z bpa2 is zero in sinhsinh
// or expsinh modes
// F.04 left expsinh case with
// -infinite limit
// F.05 TS tanhsinh mode
// F.06 lg0 true if level > 0
// check arguments ************************************
STO .00 // save X (for L at exit) (X2L)
NaN? // check for invalid limits
JMP DEI_bad_rng // b is NaN, exit
x<>y
NaN?
JMP DEI_bad_rng // a is NaN, exit
x=? Y // check for equal limits
JMP DEI_0 // a == b, return 0
// set mode flags, bma2 and bpa2 **********************
FS?S D // backup & set flag D, specials are not errors here
SF .00 // DF
x>? Y // a > b?
SF .01 // yes, flag it (rev)
INF? // a == ±inf?
SF .02 // yes, flag it (temporarily in ES)
x<>y
INF? // b == ±inf?
SF .03 // yes, flag it (temporarily in bpa2z)
// set the ES, bpa2z, left & TS flags according to the
// a & b values ++++++++++++++++++++++++++++++++++++++++
FS? .02 // was a == ±inf? (ES)
FC? .03 // was b != ±inf? (bpa2z)
SKIP 004 // skip if a != ±inf or b != ±inf
# 000 // here in the sinhsinh case --------------
# 002 // 2*bpa2 = 0, 2*bma2 = 2
CF .02 // as this flag was dirty (ES)
JMP DEI_ba2_rdy // done with sinhsinh
FC? .02 // was a != ±inf? (ES)
FS? .03 // was b == ±inf? (bpa2z)
SKIP 007 // skip if a == ±inf or b == ±inf
cENTER // here in the tanhsinh case --------------
+ // 2*bpa2 ready
x<> Z
-
ABS // 2*bma2 ready
SF .05 // flag the tanhsinh case (TS)
JMP DEI_ba2_rdy // done with tanhsinh
INF? // here in the expsinh case ----------------
x<>y // now X is the finite limit
x>? Y // finite limit > infinite one?
SF .04 // yes, left case (left)
RCL+ X // 2*bpa2 ready
CF .03 // as this flag was dirty (bpa2z)
x=0? // finite limit == 0?
SF .03 // yes, flag it (bpa2z)
SF .02 // flag the expsinh case (ES)
# 002 // 2*bma2 ready
DEI_ba2_rdy:: # 002 // here with 2*bpa & 2*bma2 on stack
/ // and mode flags ready
STO .01 // save bma2
x<>y
RCL/ L
STO .02 // save bpa2
// done with mode flags, bpa2 and bma2 ****************
// compute precision parameters ************************
# 034 // 34 digits in DBLON mode
+/-
DBL? // in DBL mode use all digits
SKIP 007
# 009 // otherwise use displayed digits + 3,
1/x // let's get the number of displayed digits
ROUND
RCL- L
EXPT // now X is -(number of displayed digits + 1)
# 002
- // -(the number of displayed digits + 3)
ENTER // here X = -(number of displayed + 3) (or -34)
10^x
STO .03 // save epsilon (eps)
SQRT // compute the convergence threshold, thr,
SDL 001 // = 10*sqrt(epsilon)
STO .04 // save the convergence threshold (thr)
DROP
ABS // compute the maximum level
LB
ROUNDI
# 002
+ // maxlevel = round(log2(digits)) + 2
STO .05 // save level (level)
FC? .03 // the sinhsinh case and the expsinh case
SKIP 004 // when the finite limit == 0 can use a
# 001 // smaller epsilon
# 000
NEIGHB
STO .03 // save epsilon for such cases (eps)
// precision parameters (epsilon, thr, level) ready ****
// compute maximum t value (tm) ************************
// tm is computed in such a way that all abscissas are
// representable at the given precision and that the
// integrand is never evaluated at the limits
FC? .05 // tanhsinh mode? (TS)
SKIP 004 // no, check other modes
# 001 // start tm computation for tanh-sinh mode
RCLMIN .01 // (bma2)
RCL+ X // X = 2*MIN(1, bma2)
SKIP 008 // goto remaining tm calculation
FS? .03 // neither SS nor ES with 0 limit? (bpa2z)
SKIP 004 // no, must be SS or ES with 0 limit
# 1/2 // start tm computation for ES with != 0 limit
RCL/ .02 // (bpa2)
ABS // X = ABS(1/(2*bpa2))
SKIP 002 // goto remaining tm calculation
RCL .03 // start tm computation for SS & ES0 cases (eps)
SQRT // X = SQRT(eps)
RCL/ .03 // continue tm computation for all cases (eps)
LN
FC? .05 // not in TS mode? (TS)
RCL+ X // 2*ln(...) for all cases except TS
RCL+ X // 4*ln(...) for all cases except TS (2*LN(...))
PI
/
LN
STO .06 // tm done
// maximum t (tm) ready ********************************
// level loop ******************************************
# 002
STO .07 // h = 2 (h)
DEI_lvl_loop:: # 000
STO .08 // ssp = 0 (ssp)
# 001
STO .09 // j = 1 (j)
# 002
STO/ .07 // h /= 2
RCL .07 // X = t
// j loop ++++++++++++++++++++++++++++++++++++++++++++++
// compute abscissas and weights ----------------------
DEI_j_loop:: COSH // cosh(t) (cosh is much faster than sinh/tanh)
STO .10 // save for later (ch)
x^2
DEC X
SQRT
PI
*
# 002
/ // pi/2*sqrt(cosh(t)^2 - 1) = pi/2*sinh(t)
FS? .02 // ES mode? (ES)
SKIP 002 // yes, want EXP
COSH // no, want COSH
SKIP 001
EXP
FILL
FC? .05 // TS mode? (TS)
SKIP 002 // no, done with weight (for ES, SS cases)
x^2 // end weight computation for TS mode
1/x
STO .11 // save weight (w)
DROP // stack filled with exp/cosh of pi/2*sinh(t)
FS? .02 // ES mode? (ES)
SKIP 003 // yes, skip this part
x^2 // no ES mode, get r (the abscissa in the
DEC X // normalized domain)
SQRT
FS? .05 // TS mode? (TS)
RCL/ Y // yes, adjust r
FS? .04 // ES mode -infinity? (left)
+/- // yes, adjust r
STO .12 // save normalized abscissa (rp)
// done with abscissas and weights --------------------
// evaluate integrand ----------------------------------
RCL* .01 // r*(b - a)/2 (bma2)
RCL+ .02 // (b + a)/2 + r*(b - a)/2 (bpa2)
FILL
XEQa // f(bpa2 + bma2*r)
SPEC? // do not stop in error (flag D is set)
# 000
RCL* .11 // fplus*w (w)
x<> .12 // p = fplus*w stored, r in X (rp)
RCL .02 // (b + a)/2 (bpa2)
RCL .01 // (b - a)/2 (bma2)
FS? .02 // ES mode? (ES)
SKIP 003
RCL* Z // no ES mode, "normal" abscissa
-
SKIP 002
RCL/ Z // ES mode, "inverse" abscissa
+ // X = (b + a)/2 ± (b - a)/2*(r or 1/r)
FILL
XEQa // f(bpa2 ± bma2*(r or 1/r))
SPEC? // do not stop in error (flag D is set)
# 000
FS? .02 // ES mode? (ES)
SKIP 002
RCL* .11 // no ES mode, normal weight, fminus*w (w)
SKIP 001
RCL/ .11 // ES mode, inverse weight, fminus/w (w)
RCL+ .12 // p = fplus*w + fminus*(w or 1/w) (rp)
RCL* .10 // X = ch*p (ch)
STO+ .08 // ssp += p
// done with integrand --------------------------------
// level 0 is special ----------------------------------
FC? .06 // is level == 0? (lg0)
JMP DEI_updte_j // yes, go to next level
INC .09 // no, sweep odd abscissas, j += 1 (j)
ABS
x<? .03 // is abs(p) < epsilon?
JMP DEI_j_exit // done with level
// and check user interrupt
KEY? Y // any key?
SKIP 005 // no, continue
# 035 // yes, check keycode for <-
x!=? Z // key is not [<-]?
SKIP 002 // it is not [<-], continue
RCL .13 // it is [<-], recall ss and exit (ss)
JMP DEI_converged
// if level > 0 done ----------------------------------
DEI_updte_j:: INC .09 // j += 1 (j)
RCL .09 // X = j (j)
RCL* .07 // X = t = j*h (h)
x<=? .06 // t <= tm? (tm)
JMP DEI_j_loop // yes, continue j loop
// done with j loop ++++++++++++++++++++++++++++++++++++
DEI_j_exit:: RCL .08 // no, update ss. X = ssp (ssp)
STO+ .13 // ss += ssp (ss)
FS? .06 // is level > 0 (lg0) ----------------------
JMP DEI_chk_conv // yes, progress to next level
# 001 // no, add 1st series term
FS? .04 // left? (left)
+/-
RCL .02 // (b + a)/2 (bpa2)
FS? .02 // ES mode? (ES)
+
FILL
XEQa // f(bpa2 (+/-1 in ES mode))
SPEC? // do not stop in error (flag D is set)
# 000
STO+ .13 // ss += f(bpa2) (ss)
// done with level 0 ----------------------------------
// show progress --------------------------------------
DEI_chk_conv:: GSB DEI_factors // apply factor to sum
TOP? // show value of approximation
MSG 25 // MSG_INTEGRATE
// done with display ----------------------------------
// check convergence ----------------------------------
RCL .08 // ssp (ssp)
RCL+ X // 2*ssp
RCL- .13 // 2*ssp - ss (ss)
ABS // ABS(2*ssp - ss)
RCL .04 // thr (thr)
RCL* .13 // thr*ss (ss)
ABS // ABS(thr*ss)
x>? Y // ABS(2*ssp - ss) < ABS(thr*ss)?
SKIP 003 // yes, computation done
SF .06 // no, mark level 0 done (lg0)
DSL .05 // update level (lvl)
JMP DEI_lvl_loop // loop
// done with level loop ********************************
// normal exit ****************************************
DEI_converged:: GSB DEI_factors // apply factors to sum
FC? .00 // restore user's D flag (DF)
CF D
RCL .00 // RCL initial X (X2L)
STO L
x<>y // X = ss
RTN // bye
// exit cause a == b **********************************
DEI_0:: DROP // here if limits are equal
# 000 // result is 0
SKIP 002
// exit cause a or b is NaN ****************************
DEI_bad_rng:: DROP // here if limit is NaN
# NaN // result is NaN
y<> .00 // restore argument (X2L)
y<> L // save in L
x<> Y // restore stack
DROP
RTN // bye
// routine to apply constants to series ****************
DEI_factors:: RCL .13 // get ss (ss)
RCL* .01 // ss*bma2 (bma2)
RCL* .07 // ss*bma2*h (h)
PI
* // ss*bma2*h*pi
# 002
/ // ss*bma2*h*pi/2
FS? .01 // reverse? (rev)
+/-
RTN // done with constants
END
Search
Member List
Calendar
Help
/
(34S) Integration (using the double-exponential method)
