Post Reply 
(34S) Integration (using the double-exponential method)
03-29-2017, 12:03 AM
Post: #21
RE: (WP-34S) Integration (using the double-exponential method)
Thanks, I'll try to find some time to pull this into XROM and do some testing.

Pauli
Find all posts by this user
Quote this message in a reply
03-31-2017, 06:22 AM
Post: #22
RE: (WP-34S) Integration (using the double-exponential method)
(03-26-2017 11:36 AM)Dieter Wrote:  BTW, Pi/2 is directly available as CNST 86. ;-)

I XROM Num PIon2 works regardless of the constant number Smile

- Pauli
Find all posts by this user
Quote this message in a reply
03-31-2017, 06:24 AM
Post: #23
RE: (WP-34S) Integration (using the double-exponential method)
What is the rationale behind setting the D flag irrespective of the user's wishes?

- Pauli
Find all posts by this user
Quote this message in a reply
03-31-2017, 09:57 AM
Post: #24
RE: (WP-34S) Integration (using the double-exponential method)
(03-31-2017 06:24 AM)Paul Dale Wrote:  What is the rationale behind setting the D flag irrespective of the user's wishes?

It is, more or less, explained in the header. The goal is not to stop the computation with an error (getting a corrupted stack and with the local registers allocated) in case the evaluation of the integrand fails. This may happen frequently because the integrand may be evaluated at really big/small numbers. E.g.
\[\int_{-1}^\infty{\ln(x^3+1)\over x^2+1}dx\]
may raise an error as, when the integrand is evaluated at 1e+199, the 3rd power will overflow (although the integral is convergent and the program, in its current form, seems to compute it correctly at \(\approx2.73777317\ldots\)

In case the integration limits are finite (tanh-sinh method), the integrand will be evaluated very near the limits (care was taken for the correct computation of tm, if such computation is correct, the integrand may be evaluated at the neighbour of the limits, but not at the limits). If the integrand goes to infinity at one limit, it is possible for it to also evaluate as infinity in the neighbour, raising an error.

In both cases, and if the integral is convergent, the weight applied to such abscissa will "probably" be low enough for the product of the integrand and the weight to be negligible, so there's low error (or no error at all) discarding such point.

The user's D flag is restored at program end.

Regards.

EDIT. The above example is only explanatory, but I've not tested if evaluating such integral does indeed evaluate the integrand at 1e199. Perhaps the algorithm has decided, prior to that, that is has reached convergence and has stopped with a solution.

César - Information must flow.
Find all posts by this user
Quote this message in a reply
03-31-2017, 10:44 AM
Post: #25
RE: (WP-34S) Integration (using the double-exponential method)
Thanks for the enlightenment. We do not want the integration routine to error out if D is cleared.

I've got the integration routine going on the 34S -- I made a few minor changes. Getting rid of all of the SKIP commands in favour of the assembler generated JMP pseudo commands. I also avoided a few branches when checking flags and used the PI/2 constant directly. I hope I've got it right. There looks to be some possibility of changing some of the divisions to multiplications by choosing different constants (things like * 1/2 instead of / 2).

I've switched the D flag off when running the user's code if it was off on entry to integrate -- the "no surprises rule" in action, the user is expecting an error if the numbers go bad in their code. The calculator shouldn't hide it.

So who's willing to try this integration code and give it a bit of a stress?


Pauli


Code:
// Double Exponential Integration for the wp34s calculator
//
// v1.2x-393 (20170327) 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)
//    - 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.
//
// 264 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

                // save stack **************************************************
                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_eso  // 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_eso::    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 weight (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
DEI_es_mode::   RCL/ Z      // ES mode, "inverse" abscissa
                +           // X = (b + a)/2 ± (b - a)/2*(r or 1/r)
DEI_call_usr::  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
                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::    FC? DF      // honour the user's D flag setting
          CF DF
        XEQUSR
        POPUSR
                SPEC?       // do not stop in error (flag D is set)
          _INT 0
        SF D
        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
Find all posts by this user
Quote this message in a reply
04-01-2017, 01:35 AM
Post: #26
RE: (WP-34S) Integration (using the double-exponential method)
(03-31-2017 10:44 AM)Paul Dale Wrote:  In:
Code:

EI_xeq_user::    FC? DF      // honour the user's D flag setting
        CF DF
        XEQUSR
        POPUSR
        SPEC?       // do not stop in error (flag D is set)
        _INT 0
        SF D
        RTN
"CF DF" should be "CF D".

Also, my name suffered from some CP1252 to UTF8 conversion that changed "é" and "í" into some funny chars. Same happened to the dash between the page numbers in the first reference. The remarks may be changed to reflect the new flag D behavior.

In my test, this version behaves as expected.

Regards.

César - Information must flow.
Find all posts by this user
Quote this message in a reply
04-01-2017, 03:49 AM
Post: #27
RE: (WP-34S) Integration (using the double-exponential method)
Thanks for the correction. I don't know what happened to your name -- I will fix.

- Pauli
Find all posts by this user
Quote this message in a reply
04-01-2017, 04:36 AM
Post: #28
RE: (WP-34S) Integration (using the double-exponential method)
I should also confirm permission to include this in the WP 34S code.

Pauli
Find all posts by this user
Quote this message in a reply
04-01-2017, 11:06 AM (This post was last modified: 04-01-2017 11:25 AM by emece67.)
Post: #29
RE: (WP-34S) Integration (using the double-exponential method)
(04-01-2017 04:36 AM)Paul Dale Wrote:  I should also confirm permission to include this in the WP 34S code.

Of course you have such permission.

This is a new version with some minor changes just to be on the safe side, now that user's D flag is obeyed. The comments on the header have also been retouched.

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

Regards.

César - Information must flow.
Find all posts by this user
Quote this message in a reply
04-02-2017, 09:34 AM
Post: #30
RE: (WP-34S) Integration (using the double-exponential method)
The code is in -- it will appear in the next build.


Pauli
Find all posts by this user
Quote this message in a reply
04-02-2017, 01:07 PM (This post was last modified: 04-02-2017 01:12 PM by emece67.)
Post: #31
RE: (WP-34S) Integration (using the double-exponential method)
Thanks Pauli.

There's a little problem in the header you added. It says:

/* The exposed functions in this file don't use the normal prologue/epilogue
* and do not operate in double precision. They call back to user code and
* manage the stack, input, output and locals themselves.
*
* Changes to this will likely cause breakage.
*/

But it does work in double precision. In any case, it is not a big issue as this will not be read by common users.

Regards.

César - Information must flow.
Find all posts by this user
Quote this message in a reply
04-02-2017, 10:02 PM
Post: #32
RE: (WP-34S) Integration (using the double-exponential method)
I don't seem to be having a good time of this Sad

I'll update later today.

The code is in the 34S repository now.


Pauli
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)