(34S) Integration (using the double-exponential method)

[emece67]

This is my implementation of the double exponential integration method for the wp34s.

As stated in other threads, this method copes with discontinuities of the integrand or its derivatives at the integration interval ends, accepts infinities as interval ends and is, usually, much faster than the Romberg method.

There are here 3 versions of this program (the 3 versions are not identical in behaviour):

1. The 1st one (v1.2r) is suited to be entered by hand (uses local labels for long jumps)
2. The 2nd one (v1.0l) must be assembled (contains JMPs and GSBs) and is suited to be stored in RAM or LIB
3. The 3rd one (v1.3x) is aimed to replace the built in integrator, must be build with the whole firmware
   

This is the version suited to be keyed in by hand.

Code:


// Double Exponential Integration for the wp34s calculator
//
// v1.2r-393 (20170327) 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
//    - estimated error in Y
//    - upper integration limit in Z
//    - lower integration limit in T
//    - updates L (with upper limit)
//    - may corrupt ABCD
//
//  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. The estimated error will be 0 in this case
//    - if the computed error estimation is not much smaller than the computed
//      result, it is assumed that all digits of the result are corrupted by
//      roundoff. In such cases, the reported result is 0 and the reported
//      error estimation equals the sum of the absolute values of the computed
//      result and error. This usually happens when the integral evaluates to 0
//    - 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.
//
// 275 steps
// 16 local registers

        LBL'DEI'          // Double-Exponential Integration
                INTM?       // check for invalid mode
                  ERR 13    // ERR_BAD_MODE
                LocR 16     // Local registers and flags:

                // local registers
                // .00    XB      backup of X (b)
                // .01    YB      backup of Y (a)
                // .02    bma2    (b - a)/2, a & b are the integration limits
                // .03    bpa2    (b + 2)/2
                // .04    eps     epsilon
                // .05    thr     convergence threshold
                // .06    lvl     level counter
                // .07    tm      maximum abscissa in the t domain
                // .08    h       space between abscissas in the t domain
                // .09    ssp     partial sum at each level
                // .10    j       abscissa counter
                // .11    ch      cosh(t)
                // .12    w       weight
                // .13    rp      abscissa (r) or fplus + fminus (p)
                // .14    ss      value of integral
                // .15    ss1     previous ss

                // local flags
                // 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

                // save stack **************************************************
                FS?S D      // backup & set flag D, specials are not errors here
                  SF .00
                cSTO .00    // save XY (ba)
                // A..D[->]    // save ABCD
                // check arguments  ************************************
                NaN?        // check for invalid limits
                  GTO 95    // b is NaN, exit
                x<>y
                NaN?
                  GTO 95    // a is NaN, exit
                x=? Y       // check for equal limits
                  GTO 90    // a == b, return 0
                // set mode flags, bma2 and bpa2  **********************
                x>? Y       // a > b?
                  SF .01    // yes, flag it
                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?
                  FC? .03   // was b != ±inf?
                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
                GTO 10      // done with sinhsinh
                FC? .02     // was a != ±inf?
                  FS? .03   // was b == ±inf?
                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
                GTO 10      // 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
                RCL+ X      // 2*bpa2 ready
                CF .03      // as this flag was dirty
                x=0?        // finite limit == 0?
                  SF .03    // yes, flag it
                SF .02      // flag the expsinh case
                # 002       // 2*bma2  ready
        LBL 10
                # 002       // here with 2*bpa & 2*bma2 on stack
                /           // and mode flags ready
                STO .02     // save bma2
                x<>y
                RCL/ L
                STO .03     // 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 .04     // save epsilon
                SQRT        // compute the convergence threshold, .05,
                SDL 001     // = 10*sqrt(epsilon)
                STO .05     // save the convergence threshold
                DROP
                ABS         // compute the maximum level
                LB
                ROUNDI
                # 002
                +           // maxlevel = round(log2(digits)) + 2
                STO .06     // save 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 .04     // save epsilon for such cases
                // 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?
                  SKIP 004  // no, check other modes
                # 001       // start tm computation for tanh-sinh mode
                RCLMIN .02
                RCL+ X      // X = 2*MIN(1, bma2)
                SKIP 008    // goto remaining tm calculation
                FS? .03     // neither SS nor ES with 0 limit?
                  SKIP 004  // no, must be SS or ES with 0 limit
                # 1/2       // start tm computation for ES with != 0 limit
                RCL/ .03    // RCP/ bpa2
                ABS         // X = ABS(1/(2*bpa2))
                SKIP 002    // goto remaining tm calculation
                RCL .04     // start tm computation for SS & ES0 cases
                SQRT        // X = SQRT(eps)
                RCL/ .04    // continue tm computation for all cases
                LN
                FC? .05     // not in TS mode?
                  RCL+ X    // 2*ln(...) for all cases except TS
                RCL+ X      // 4*ln(...) for all cases except TS (2*LN(...))
                PI
                /
                LN
                STO .07     // tm done
                // maximum t (tm) ready ********************************
                // level loop ******************************************
                # 002
                STO .08     // h = 2
        LBL 20
                # 000
                STO .09     // ssp = 0
                # 001
                STO .10     // j = 1
                # 002
                STO/ .08    // h /= 2
                RCL .08     // X = t
                // j loop ++++++++++++++++++++++++++++++++++++++++++++++
                // compute abscissas and weights  ----------------------
        LBL 30
                COSH        // cosh(t) (cosh is much faster than sinh/tanh)
                STO .11     // save for later
                x^2
                DEC X
                SQRT
                PI
                *
                # 002
                /           // pi/2*sqrt(cosh(t)^2 - 1) = pi/2*sinh(t)
                FS? .02     // ES mode?
                  SKIP 002  // yes, want EXP
                COSH        // no, want COSH
                SKIP 001
                EXP
                FILL
                FC? .05     // TS mode?
                  SKIP 002  // no, done with weight (for ES, SS cases)
                x^2         // end weight computation for TS mode
                1/x
                STO .12     // save weight
                DROP        // stack filled with exp/cosh of pi/2*sinh(t)
                FS? .02     // ES mode?
                  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?
                  RCL/ Y    // yes, adjust r
                FS? .04     // ES mode -infinity?
                  +/-       // yes, adjust r
                STO .13     // save normalized abscissa
                // done with abscissas and weights  --------------------
                // evaluate integrand ----------------------------------
                RCL* .02    // r*(b - a)/2
                RCL+ .03    // (b + a)/2 + r*(b - a)/2
                FILL
                XEQa        // f(bpa2 + bma2*r)
                SPEC?       // do not stop in error (flag D is set)
                  # 000
                RCL* .12    // fplus*w
                x<> .13     // p = fplus*w stored, r in X
                RCL .03     // (b + a)/2
                RCL .02     // (b - a)/2
                FS? .02     // ES mode?
                  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?
                  SKIP 002
                RCL* .12    // no ES mode, normal weight, fminus*w
                SKIP 001
                RCL/ .12    // ES mode, inverse weight, fminus/w
                RCL+ .13    // p = fplus*w + fminus*(w or 1/w)
                RCL* .11    // X = ch*p
                STO+ .09    // ssp += p
                // done with integrand  --------------------------------
                // level 0 is special ----------------------------------
                FC? .06     // is level == 0?
                  GTO 40    // yes, go to next level
                INC .10     // no, sweep odd abscissas, j += 1
                ABS         // X = ABS(p)
                RCL .09     // ssp
                RCL* .04    // eps
                ABS         // X = ABS(ssp*eps), Y = ABS(p)
                x>=? Y      // ABS(p) <= ABS(eps*eps)?
                  GTO 50    // then 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 .15     // it is, recall the last computation
                GTO 70      // and exit
                // if level > 0 done  ----------------------------------
        LBL 40
                INC .10     // j += 1
                RCL .10     // X = j
                RCL* .08    // X = t = j*h
                x<=? .07    // t <= tm?
                  GTO 30    // yes, continue j loop
                // done with j loop ++++++++++++++++++++++++++++++++++++
        LBL 50
                RCL .09     // no, update ss. X = ssp
                STO+ .14    // ss += ssp
                FS? .06     // is level > 0 ----------------------------
                  GTO 60    // yes, progress to next level
                # 001       // no, add 1st series term
                FS? .04     // left?
                  +/-
                RCL .03     // (b + a)/2
                FS? .02     // ES mode?
                  +
                FILL
                XEQa        // f(bpa2 (+/-1 in ES mode))
                SPEC?       // do not stop in error (flag D is set)
                  # 000
                STO+ .14    // ss += f(bpa2)
                // done with level 0  ----------------------------------
                // apply constant coeffs to sum ------------------------
        LBL 60
                RCL .14
                RCL* .02    // ss*bma2
                RCL* .08    // ss*bma2*h
                PI
                *           // ss*bma2*h*pi
                # 002
                /           // ss*bma2*h*pi/2
                FS? .01     // reverse?
                  +/-       // yes,so change sign
                // done with constant coeffs  --------------------------
                // show progress  --------------------------------------
                TOP?        // show value of approximation
                  MSG 25    // MSG_INTEGRATE
                // check convergence  ----------------------------------
                RCL .09
                RCL+ X      // 2*ssp
                RCL- .14    // 2*ssp - ss
                ABS         // ABS(2*ssp - ss)
                RCL .05
                RCL* .14    // thr*ss
                ABS         // ABS(thr*ss)
                x>? Y       // ABS(2*ssp - ss) < ABS(thr*ss)?
                  SKIP 004  // yes, computation done
                z<> .15     // no, save ss (with all coeffs applied)
                SF .06      // mark level 0 done,
                DSL .06     // update level &...
                  GTO 20    // loop.
                RCL Z       // recall result
                // done with level loop ********************************
                // compute error estimation and check for bad results **
                // here with result in X
        LBL 70
                RCL X       // stack: ss-ss-?-?
                ABS         // stack: |ss|-ss-?-?
                RCL Y       // stack: ss-|ss|-ss-?
                RCL- .15    // stack: (ss-ss1)-|ss|-ss-?
                ABS         // stack: err-|ss|-ss-?
                RCL X       // stack: err-err-|ss|-ss
                SDL 001     // stack: 10*err-err-|ss|-ss
                x<? Z       // 10*err < |ss|?
                  SKIP 004  // yes, assume result is OK
                x<> Z       // no, bad result. stack: |ss|-err-10*err-ss
                +           // stack: new_err-10*err-ss-ss
                # 000       // stack: 0-new_err-10*err-ss
                STO T       // stack: 0-new_err-10*err-0
                x<> T       // stack: 0-err-... or ss-err-...
                // done with bad results  ******************************
                //  exit  **********************************************
        LBL 80
                TOP?        // show new computed value (0)
                  MSG 25    // MSG_INTEGRATE
                cRCL .00    // stack: b-a-s-ss1
                // [->]A..D    // try to restore ABCD
                STO L       // b into L
                cx<> Z      // stack: s-ss1-b-a
                FC? .00     // restore user's D flag (DF)
                  CF D
              RTN           // bye
                // exit cause a == b  **********************************
        LBL 90
                # 000       // error & result both 0
                SKIP 002
                // exit cause a or b is NaN ****************************
        LBL 95
                # NaN       // error & result both NaN
                RCL X
                GTO 80      // exit
              END

[emece67]

This one is the version to be assembled an put into LIB or RAM. Beware that this version does not exactly reproduce the behavior of the other two.

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

[emece67]

And this last one is that aimed to replace the built-in integration program.

Code:


// Double Exponential Integration for the wp34s calculator
//
// v1.3x-398 (20170104) by M. César Rodríguez GGPL
//
// USAGE:
//  Inputs:
//    - integration limits in X (upper limit) and Y (lower limit)
//  Outputs:
//    - approximate value of integral in X
//    - estimated error in Y
//    - upper integration limit in Z
//    - lower integration limit in T
//    - updates L (with upper limit)
//    - in SSIZE8 mode, all ABCD will be corrupted
//
//  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. The estimated error will be 0 in this case
//    - if the computed error estimation is not much smaller than the computed
//      result, it is assumed that all digits of the result are corrupted by
//      roundoff. In such cases, the reported result is 0 and the reported
//      error estimation equals the sum of the absolute values of the computed
//      result and error. This usually happens when the integral evaluates to 0
//    - if the user has set the D flag, 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±199 in DBLOFF mode and
//      1e±3088 in DBLON mode). Ensure that your integrand program
//      behaves as expected for such arguments (you may consider to set
//      flag D in some circumstances)
//    - 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.
//
// 263 steps
// 16 local registers

              XLBL"INTEGRATE" // Double-Exponential Integration
                INTM?       // check for invalid mode
                  ERR ERR_BAD_MODE
                LocR 16     // Local registers and flags:

// local registers
#define XB    .00   // backup of X (b)
#define YB    .01   // backup of Y (a)
#define bma2  .02   // (b - a)/2, a & b are the integration limits
#define bpa2  .03   // (b + 2)/2
#define eps   .04   // epsilon
#define thr   .05   // convergence threshold
#define lvl   .06   // level counter
#define tm    .07   // maximum abscissa in the t domain
#define h     .08   // space between abscissas in the t domain
#define ssp   .09   // partial sum at each level
#define j     .10   // abscissa counter
#define ch    .11   // cosh(t)
#define w     .12   // weight
#define rp    .13   // abscissa (r) or fplus + fminus (p)
#define ss    .14   // value of integral
#define ss1   .15   // previous ss

// local flags
#define DF    .00   // backup of D flag
#define rev   .01   // b is < a
#define ES    .02   // expsinh mode
#define bpa2z .03   // bpa2 is zero in sinhsinh or expsinh modes
#define left  .04   // expsinh case with -Infinite limit
#define TS    .05   // tanhsinh mode
#define lg0   .06   // true if level > 0

                FS?S D      // backup & set flag D, specials are not errors here
                  SF DF
                cSTO XB     // save XY (ba)
                // check arguments  ************************************
                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  **********************
                x>? Y       // a > b?
                  SF rev    // yes, flag it
                INF?        // a == ±inf?
                  SF ES     // yes, flag it (temporarily in ES)
                x<>y
                INF?        // b == ±inf?
                  SF bpa2z  // yes, flag it (temporarily in bpa2z)
                // set the ES, bpa2z, left & TS flags according to the
                // a & b values ++++++++++++++++++++++++++++++++++++++++
                FS? ES      // was a == ±inf?
                  FC? bpa2z // was b != ±inf?
                    JMP DEI_inf_1   // skip if a != ±inf or b != ±inf
                _INT 000    // here in the sinhsinh case  --------------
                _INT 002    // 2*bpa2 = 0, 2*bma2 = 2
                CF ES       // as this flag was dirty
                JMP DEI_ba2_rdy     // done with sinhsinh
DEI_inf_1::     FC? ES      // was a != ±inf?
                  FS? bpa2z // was b == ±inf?
                    JMP DEI_inf_2   // skip if a == ±inf or b == ±inf
                cENTER      // here in the tanhsinh case  --------------
                +           // 2*bpa2 ready
                x<> Z
                -
                ABS         // 2*bma2 ready
                SF TS       // flag the tanhsinh case
                JMP DEI_ba2_rdy     // done with tanhsinh
DEI_inf_2::     INF?        // here in the expsinh case ----------------
                  x<>y      // now X is the finite limit
                x>? Y       // finite limit > infinite one?
                  SF left   // yes, left case
                RCL+ X      // 2*bpa2 ready
                CF bpa2z    // as this flag was dirty
                x=0?        // finite limit == 0?
                  SF bpa2z  // yes, flag it
                SF ES       // flag the expsinh case
                _INT 002    // 2*bma2 ready
DEI_ba2_rdy::   _INT 002    // here with 2*bpa & 2*bma2 on stack
                /           // and mode flags ready
                STO bma2    // save bma2
                x<>y
                RCL/ L
                STO bpa2    // save bpa2
                // done with mode flags, bpa2 and bma2  ****************
                // compute precision parameters ************************
                _INT 034    // 34 digits in DBLON mode
                +/-
                DBL?        // in DBL mode use all digits
                  JMP DEI_dbl
                _INT 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)
                _INT 002
                -           // -(the number of displayed digits + 3)
DEI_dbl::       ENTER       // here X = -(number of displayed + 3) (or -34)
                10^x
                STO eps     // save epsilon
                SQRT        // compute the convergence threshold, thr,
                SDL 001     // = 10*sqrt(epsilon)
                STO thr     // save the convergence threshold
                DROP
                ABS         // compute the maximum level
                LB
                ROUNDI
                _INT 002
                +           // maxlevel = round(log2(digits)) + 2
                STO lvl     // save level
                FC? bpa2z   // the sinhsinh case and the expsinh case
                  JMP DEI_tanhsinh  // when the finite limit == 0 can use a
                _INT 001    // smaller epsilon
                _INT 000
                NEIGHB
                STO eps     // save epsilon for such cases
                // 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
DEI_tanhsinh::  FC? TS      // tanhsinh mode?
                  JMP DEI_ss_es     // no, check other modes
                _INT 001    // start tm computation for tanh-sinh mode
                RCLMIN bma2
                RCL+ X      // X = 2*MIN(1, bma2)
                JMP DEI_cont        // goto remaining tm calculation
DEI_ss_es::     FS? bpa2z   // neither SS nor ES with 0 limit?
                  JMP DEI_ss_es0    // no, must be SS or ES with 0 limit
                _INT 1/2    // start tm computation for ES with != 0 limit
                RCL/ bpa2
                ABS         // X = ABS(1/(2*bpa2))
                JMP DEI_cont        // goto remaining tm calculation
DEI_ss_es0::    RCL eps     // start tm computation for SS & ES0 cases
                SQRT        // X = SQRT(eps)
DEI_cont::      RCL/ eps    // continue tm computation for all cases
                LN
                FC? TS      // not in TS mode?
                  RCL+ X    // 2*ln(...) for all cases except TS
                RCL+ X      // 4*ln(...) for all cases except TS (2*LN(...))
                PI
                /
                LN
                STO tm      // tm done
                // maximum t (tm) ready ********************************
                // level loop ******************************************
                _INT 002
                STO h       // h = 2
DEI_lvl_loop::  _INT 000
                STO ssp     // ssp = 0
                _INT 001
                STO j       // j = 1
                _INT 002
                STO/ h      // h /= 2
                RCL h       // X = t
                // j loop ++++++++++++++++++++++++++++++++++++++++++++++
                // compute abscissas and weights  ----------------------
DEI_j_loop::    COSH        // cosh(t) (cosh is much faster than sinh/tanh)
                STO ch      // save for later
                x^2
                DEC X
                SQRT
                Num PIon2
                *           // pi/2*sqrt(cosh(t)^2 - 1) = pi/2*sinh(t)
                FS? ES      // ES mode?
                  EXP       // yes, want EXP
                FC? ES
                  COSH      // no, want COSH
                FILL
                FC? TS      // TS mode?
                  JMP DEI_not_ts    // no, done with w (for ES, SS cases)
                x^2         // end weight computation for TS mode
                1/x
DEI_not_ts::    STO w       // save weight
                DROP        // stack filled with exp/cosh of pi/2*sinh(t)
                FS? ES      // ES mode?
                  JMP DEI_es_skip   // yes, skip this part
                x^2         // no ES mode, get r (the abscissa in the
                DEC X       // normalized domain)
                SQRT
DEI_es_skip::   FS? TS      // TS mode?
                  RCL/ Y    // yes, adjust r
                FS? left    // ES mode -infinity?
                  +/-       // yes, adjust r
                STO rp      // save normalized abscissa
                // done with abscissas and weights  --------------------
                // evaluate integrand ----------------------------------
                RCL* bma2   // r*(b - a)/2
                RCL+ bpa2   // (b + a)/2 + r*(b - a)/2
                GSB DEI_xeq_user    // f(bpa2 + bma2*r)
                RCL* w      // fplus*w
                x<> rp      // p = fplus*w stored, r in X
                RCL bpa2    // (b + a)/2
                RCL bma2    // (b - a)/2
                FS? ES      // ES mode?
                  JMP DEI_es_mode
                RCL* Z      // no ES mode, "normal" abscissa
                -
                JMP DEI_call_usr_m
DEI_es_mode::   RCL/ Z      // ES mode, "inverse" abscissa
                +           // X = (b + a)/2 ± (b - a)/2*(r or 1/r)
DEI_call_usr_m:: GSB DEI_xeq_user   // f(bpa2 ± bma2*(r or 1/r))
                FS? ES      // ES mode?
                  RCL/ w    // ES mode, inverse weight, fminus/w
                FC? ES
                  RCL* w    // no ES mode, normal weight, fminus*w
                RCL+ rp     // p = fplus*w + fminus*(w or 1/w)
                RCL* ch     // X = ch*p
                SPEC?
                  _INT 000  // safety check for weight
                STO+ ssp    // ssp += p
                // done with integrand  --------------------------------
                // level 0 is special ----------------------------------
                FC? lg0     // is level == 0?
                  JMP DEI_updte_j   // yes, go to next level
                INC j       // no, sweep odd abscissas, j += 1
                ABS         // X = ABS(p)
                RCL ssp
                RCL* eps
                ABS         // X = ABS(ssp*eps), Y = ABS(p)
                x>=? Y      // ABS(p) <= ABS(eps*eps)?
                  JMP DEI_j_exit    // then done with level
                // and check user interrupt ............................
                KEY? Y      // any key?
                  JMP DEI_updte_j   // no, continue
                _INT 035    // yes, check keycode for [<-]
                x!=? Z      // key is not [<-]?
                  JMP DEI_updte_j   // it is not [<-], continue
                RCL ss1     // it is, recall the last computation
                JMP DEI_chk_bad     // and exit (checking big error results)
                // call user's function  -------------------------------
DEI_xeq_user::    SPEC?     // abscissa is good?
                    JMP DEI_bad_absc  // no, skip point
                  FC? DF    // honour the user's D flag setting
                    CF D
                  XEQUSR
                  POPUSR
                  SPEC?     // do not stop in error (if flag D was set)
                    _INT 0
                  SF D      // set flag D for internal calculations
                RTN
DEI_bad_absc::    _INT 000
                RTN
                // if level > 0 done  ----------------------------------
DEI_updte_j::   INC j       // j += 1
                RCL j       // X = j
                RCL* h      // X = t = j*h
                x<=? tm     // t <= tm?
                  JMP DEI_j_loop    // yes, continue j loop
                // done with j loop ++++++++++++++++++++++++++++++++++++
DEI_j_exit::    RCL ssp     // no, update ss. X = ssp
                STO+ ss     // ss += ssp
                FS? lg0     // is level > 0 ----------------------------
                  JMP DEI_chk_conv  // yes, progress to next level
                _INT 001    // no, add 1st series term
                FS? left    // left?
                  +/-
                RCL bpa2    // (b + a)/2
                FS? ES      // ES mode?
                  +
                GSB DEI_xeq_user    // f(bpa2 (+/-1 in ES mode))
                STO+ ss     // ss += f(bpa2)
                // done with level 0  ----------------------------------
                // apply constant coeffs to sum ------------------------
DEI_chk_conv::  RCL ss
                RCL* bma2 // ss*bma2
                RCL* h    // ss*bma2*h
                Num PIon2
                *         // ss*bma2*h*pi/2
                FS? rev   // reverse?
                  +/-     // yes,so change sign
                // done with constant coeffs  --------------------------
                // show progress  --------------------------------------
                TOP?        // show value of approximation
                  MSG MSG_INTEGRATE
                // check convergence  ----------------------------------
                RCL ssp
                RCL+ X      // 2*ssp
                RCL- ss     // 2*ssp - ss
                ABS         // ABS(2*ssp - ss)
                RCL thr
                RCL* ss     // thr*ss
                ABS         // ABS(thr*ss)
                x>? Y       // ABS(2*ssp - ss) < ABS(thr*ss)?
                  JMP DEI_fin       // yes, computation done
                z<> ss1     // no, save ss (with all coeffs applied)
                SF lg0      // mark level 0 done,
                DSL lvl     // update level &...
                  JMP DEI_lvl_loop  // loop.
DEI_fin::       RCL Z       // recall result
                // done with level loop ********************************
                // compute error estimation and check for bad results **
                // here with result in X
DEI_chk_bad::   RCL X       // stack: ss-ss-?-?
                ABS         // stack: |ss|-ss-?-?
                RCL Y       // stack: ss-|ss|-ss-?
                RCL- ss1    // stack: (ss-ss1)-|ss|-ss-?
                ABS         // stack: err-|ss|-ss-?
                RCL X       // stack: err-err-|ss|-ss
                SDL 001     // stack: 10*err-err-|ss|-ss
                x<? Z       // 10*err < |ss|?
                  JMP DEI_res_okay  // yes, assume result is OK
                x<> Z       // no, bad result. stack: |ss|-err-10*err-ss
                +           // stack: new_err-10*err-ss-ss
                _INT 000    // stack: 0-new_err-10*err-ss
                STO T       // stack: 0-new_err-10*err-0
DEI_res_okay::  x<> T       // stack: 0-err-... or ss-err-...
                // done with bad results  ******************************
                //  exit  **********************************************
DEI_exit::      TOP?        // show new computed value (0)
                  MSG MSG_INTEGRATE
                cRCL XB     // stack: b-a-s-ss1
                STO L       // b into L
                cx<> Z      // stack: s-ss1-b-a
                FC? DF      // restore user's D flag (DF)
                  CF D
              RTN           // bye
                // exit cause a == b  **********************************
DEI_0::         _INT 000    // error & result both 0
                JMP DEI_dup_exit
                // exit cause a or b is NaN ****************************
DEI_bad_rng::   Num NaN     // error & result both NaN
DEI_dup_exit::  RCL X
                JMP DEI_exit        // exit

#undef XB
#undef YB
#undef bma2
#undef bpa2
#undef eps
#undef thr
#undef lvl
#undef tm
#undef h
#undef ssp
#undef j
#undef ch
#undef w
#undef rp
#undef ss
#undef ss1

#undef DF
#undef rev
#undef ES
#undef bpa2z
#undef left
#undef TS
#undef lg0

[Paul Dale]

What do people think about including the last one in the 34S firmware?
Does it do the right thing with the final stack?

Pauli

[Luigi Vampa]

My thumb is up. It could add to the list of WP34s 'uniquenesses'.
BTW, the license and copyright header reads: "M. César". That "M" must stand for marvellous... or may be mathematician... or both in one single letter, by some sort of arcane&syncopated algebraic expression ;O)

[emece67]

(03-25-2017 12:33 AM)Paul Dale Wrote:  What do people think about including the last one in the 34S firmware?

Perhaps a little bit more time for people other than me to test the program is advisable.

(03-25-2017 12:33 AM)Paul Dale Wrote:  Does it do the right thing with the final stack?

It does the very same thing, in this respect, that the built in integrator: corrupt all the stack. I was able to write a version (not to be used in XROM) preserving the stack, but I found some trouble when porting it to XROM. Perhaps I need some advice here. My approach was:

• Use 2 local register to backup Z & T
• When returning, build a stack with the backed up T & Z, the result in Y and a dummy value in X
• Then do a DROP and return
  

It worked if the user was working in SSIZE4, but not in SSIZE8.

I do not want to use xIN/xOUT because it forces DBLON mode and because I'm not sure if XEQUSR/POPUSR can be used inside xIN/xOUT.

In any case, I'm still trying to fix this (looking into ->A..D->).

(03-25-2017 10:21 PM)Luigi Vampa Wrote:  BTW, the license and copyright header reads: "M. César". That "M" must stand for marvellous... or may be mathematician... or both in one single letter, by some sort of arcane&syncopated algebraic expression ;O)

Jaja, Luigi, la M solo es de Manuel.

[Paul Dale]

(03-25-2017 11:52 PM)emece67 Wrote:  It does the very same thing, in this respect, that the built in integrator: corrupt all the stack.

The built in integrator doesn't corrupt the stack. It leaves the integral in X, the error estimate in Y and the original limits of integration in Z & T and also preserves the original X value in L. Refer to the 34C manual page 211. I'm don't remember what it does to A, B, C & D with stack size 8 -- XROM switches to four level stack mode on entry so they won't be directly altered, although the fills before calling user code might.

Quote:I do not want to use xIN/xOUT because it forces DBLON mode and because I'm not sure if XEQUSR/POPUSR can be used inside xIN/xOUT.

xIN/xOUT cannot nest so they cannot be used in integrate because the user's program might use a function that uses xIN/xOUT.

Quote:In any case, I'm still trying to fix this (looking into ->ABCD->).

[->]A..D and A..D[->] stash these four registers away in temporary volatile memory. They provided a way to gain some extra working registers before we had locals. It must be assumed that the temporary registers are lost when XEQUSR is called.


- Pauli

[pier4r]

(03-24-2017 02:23 PM)emece67 Wrote:  cut

I wish more programs posted here would have exhaustive comments like yours instead of just the program (that in RPN / RPL is unlikely to be "self descriptive")

[Dieter]

(03-24-2017 02:21 PM)emece67 Wrote:  This is my implementation of the double exponential integration method for the wp34s.

Looks like a nice and carefully coded program. Let me ask a few questions:

There is a logic that calculates an epsilon value depending on the display setting. This seems to work fine in FIX mode. But how does this work in SCI/ENG or ALL mode? Note that (unlike RDP) ROUND rounds 1/9 to at most 11 places. Or a certain number of significant digits in SCI/ENG mode.

At LBL 40 the program checks if the integration variable is within the integration limits. This is done by calculating the current x and comparing it to tm in R.06. In such cases I always wonder if roundoff errors might not lead to an error here. Suppose tm is 2 and h*j rounds to 2,000000000000001 so that the test returns false. Can this possibly happen? A simple counter would be a safe solution.

Very small values of the integral, e.g. 1E–16, may be calculated because the actual value is zero, but just as well the correct value of the true integral. Is the program able to distinguish these two cases?

BTW, Pi/2 is directly available as CNST 86. ;-)

Dieter

[emece67]

(03-26-2017 12:25 AM)Paul Dale Wrote:  The built in integrator doesn't corrupt the stack. It leaves the integral in X, the error estimate in Y and the original limits of integration in Z & T and also preserves the original X value in L. Refer to the 34C manual page 211. I'm don't remember what it does to A, B, C & D with stack size 8 -- XROM switches to four level stack mode on entry so they won't be directly altered, although the fills before calling user code might.

Surely that's the behavior of the 34C, but the 34s:

1. Returns the computed result in X
2. Returns in Y what seems to be some kind of a previous approximation to the result (it is not much smaller than the result, as expected for an error estimation, but approximately equal to it)
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L
5. Corrupts ABCD with some unknown data (no matter if the user is in SSIZE4 or SSIZE8, ABCD are always corrupted)
   

Thus, I've modified my program to behave more or less in the same way, so it now (v1.1r and v1.1x):

1. Returns the computed result in X (as before)
2. Returns in Y the previous approximation to the result
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L (as before)
5. Keeps ABCD untouched (no matter if the user is in SSIZE4 or SSIZE8)
   

(03-26-2017 11:36 AM)Dieter Wrote:  There is a logic that calculates an epsilon value depending on the display setting. This seems to work fine in FIX mode. But how does this work in SCI/ENG or ALL mode? Note that (unlike RDP) ROUND rounds 1/9 to at most 11 places. Or a certain number of significant digits in SCI/ENG mode.

Later, epsilon (eps) is mainly used to compute some other values (thr & tm). Such values depend on the number of significant digits displayed thus, epsilon is not an absolute threshold for anything, but a number carrying information about the number of significant digits displayed. This way, eps is more meaningful here in ALL/SCI/ENG modes than in FIX mode.

There was also an absolute comparison with eps that allows the program not to iterate (mainly when both integration ends are infinite) through a big number of abscissas whose weights are so small than they do not contribute to the final result. This was only an heuristic to speed-up computation, not the main convergence test. In any case, thanks to your comment, I've found that such test interferes with cases where the integrand is always very small, so I have changed such test and now all tests are relative.

(03-26-2017 11:36 AM)Dieter Wrote:  At LBL 40 the program checks if the integration variable is within the integration limits. This is done by calculating the current x and comparing it to tm in R.06. In such cases I always wonder if roundoff errors might not lead to an error here. Suppose tm is 2 and h*j rounds to 2,000000000000001 so that the test returns false. Can this possibly happen? A simple counter would be a safe solution.

j is a small integer only updated with INC. h is an integer power of 2 from 2 to 1/256 updated by continued division by 2, all these powers have an exact decimal representation in the wp34s. Thus, h*j is always exactly representable in the machine, so no concerns about roundoff here.

(03-26-2017 11:36 AM)Dieter Wrote:  Very small values of the integral, e.g. 1E–16, may be calculated because the actual value is zero, but just as well the correct value of the true integral. Is the program able to distinguish these two cases?

Being now all tests relative, integrands that are always small do not pose any special problem. But the program does not distinguish if the true integral is zero and the computed, small, result is due to round-off, or the computed, small, result is a nice approximation of a small but not zero integral.

(03-26-2017 11:36 AM)Dieter Wrote:  BTW, Pi/2 is directly available as CNST 86. ;-)

In my local builds that's not true, cause I have a modified CONST menu, so I try to avoid such constructs.

[Dieter]

(03-26-2017 06:34 PM)emece67 Wrote:  Surely that's the behavior of the 34C, but the 34s:

1. Returns the computed result in X
2. Returns in Y what seems to be some kind of a previous approximation to the result (it is not much smaller than the result, as expected for an error estimation, but approximately equal to it)
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L
   

AFAIR one of the final steps of the integration routine was a STO–Y so that Y does not hold the previous approximation but the difference compared to the last one. The content of the other lower stack registers are just as they are supposed to.

(03-26-2017 06:34 PM)emece67 Wrote:  5. Corrupts ABCD with some unknown data (no matter if the user is in SSIZE4 or SSIZE8, ABCD are always corrupted)

A, B, C and D may be corrupted, yes.

(03-26-2017 06:34 PM)emece67 Wrote:  Thus, I've modified my program to behave more or less in the same way, so it now (v1.1r and v1.1x):

1. Returns the computed result in X (as before)
2. Returns in Y the previous approximation to the result
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L (as before)
5. Keeps ABCD untouched (no matter if the user is in SSIZE4 or SSIZE8)
   

I would prefer an error estimate in Y, as it is returned by the 34C and 15C. And I think this also is the way the 34s integration routine is supposed to work. Pauli?

Dieter

[Paul Dale]

(03-26-2017 06:34 PM)emece67 Wrote:  Corrupts ABCD with some unknown data (no matter if the user is in SSIZE4 or SSIZE8, ABCD are always corrupted)

I'm not seeing corruption here in SSIZE4 mode. In SSIZE8 mode, the values are whatever is left from the final function evaluation.


Pauli

[Paul Dale]

(03-26-2017 08:19 PM)Dieter Wrote:  I would prefer an error estimate in Y, as it is returned by the 34C and 15C. And I think this also is the way the 34s integration routine is supposed to work. Pauli?

That's how I thought it worked. Seems I was wrong or something got lost somewhere.


- Pauli

[Dieter]

(03-27-2017 12:10 AM)Paul Dale Wrote:  That's how I thought it worked. Seems I was wrong or something got lost somewhere.

The current integration code can be examined on sourceforge.net. The exit routine at int_done starts with the last two approximations in X and Y, then recalls the two limits and rearranges the stack so that the limits are in Z and T and the last two approximations in X and Y. And indeed there is no final STO–Y that would leave an error estimate in Y.

Edit: The switch from Gauss-Kronrod to Romberg was somewhere between v2240 and v2260. If you compare the integration code of v2240 with that of v2260 you will notice that the former indeed calculated an error estimate from the last two approximations...

Code:

            /* Set up the stack for our output */
            GSB int_restore_limits
            RCL .03
            RCL- .02                /* err, l, u, ?, G*/
            RCL .03                /* K, err, l, u, G*/
            RTN

...while already the earliest Romberg versions didn't:

Code:

/* Three matches in a row is goodness. */
int_really_done::  [cmplx]RCL cr_limits
                   STO L
                   [cmplx]x[<->] Z
                   RTN

...and this obviously has never been corrected. Since the 34s manual says the user interface is "as in the HP-15C" I think this should be done, e.g. like this:

Code:

int_done::  RCL- Y
            RCL L
            [cmplx]RCL cr_limits
            STO L
            [cmplx]x[<->] Z
            RTN

This way (last approx. – prev. approx.) is returned in Y.

Dieter

[emece67]

I'll change my program to return an error estimation in Y.

I've found that, after the change I made cause of Dieter's comment about epsilon, there's now a way to discern when a small result is due to roundoff when evaluating an integral that is zero or a small result due to an integral that is truly small, but different from 0. I'll try to add an heuristic to return 0 when the integral seems to be zero.

Regards.

[emece67]

Changes done (v1.2r and v1.2x).

Now the program:

1. Returns the computed result in X
2. Returns an estimation of the error in Y
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L
   

There's also a check around the estimated error used to distinguish between integrals that evaluate to 0 and those that evaluate to small numbers. Now, integrals suspected to be equal to 0 are reported as such, thus \(\int_{-\pi}^\pi \cos(x)dx\) will be evaluated to 0, but \(\int_{-\pi/2}^{\pi/2} 10^{-100}\cos(x)dx\) will be evaluated to \(2·10^{-100}\). Both computations will return and estimated error different from 0.

Regards.

Edited: removing the comment about ABCD remaining untouched.

[Paul Dale]

(03-27-2017 05:42 PM)emece67 Wrote:  Keeps ABCD untouched

It is not possible to use A..D-> and ->A..D to preserve the A-D registers across a call to user code. These commands stash the register values at a fixed location in volatile RAM which will be lost under some circumstances.

It will also break with nested integrals or with mixed solve and integrate.

The only way to preserve the A-D registers is to allocate four additional local registers for the purpose.


The A..D-> and ->A..D commands are left from when I implemented solve without the benefit of local registers and had to squeeze everything into a very limited number of special XROM registers. I needed more temporary space so I added these two commands. Like several other specials for XROM, they are extremely specific to my purpose at the time and come with major caveats.


- Pauli

[emece67]

Program description modified removing the statement about ABCD remaining untouched.

[Paul Dale]

It would be prudent to also remove the ->A..D and A..D-> commands as well.


Pauli

[emece67]

As you insisted, I've commented out such lines.

Regards,