(34S) Integration (using the double-exponential method)
03-29-2017, 12:03 AM
Post: #21
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
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
03-31-2017, 06:22 AM
Post: #22
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
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

- Pauli
03-31-2017, 06:24 AM
Post: #23
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
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
03-31-2017, 09:57 AM
Post: #24
 emece67 Senior Member Posts: 379 Joined: Feb 2015
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.
03-31-2017, 10:44 AM
Post: #25
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
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
04-01-2017, 01:35 AM
Post: #26
 emece67 Senior Member Posts: 379 Joined: Feb 2015
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.
04-01-2017, 03:49 AM
Post: #27
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
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
04-01-2017, 04:36 AM
Post: #28
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
RE: (WP-34S) Integration (using the double-exponential method)
I should also confirm permission to include this in the WP 34S code.

Pauli
04-01-2017, 11:06 AM (This post was last modified: 04-01-2017 11:25 AM by emece67.)
Post: #29
 emece67 Senior Member Posts: 379 Joined: Feb 2015
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.
04-02-2017, 09:34 AM
Post: #30
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
RE: (WP-34S) Integration (using the double-exponential method)
The code is in -- it will appear in the next build.

Pauli
04-02-2017, 01:07 PM (This post was last modified: 04-02-2017 01:12 PM by emece67.)
Post: #31
 emece67 Senior Member Posts: 379 Joined: Feb 2015
RE: (WP-34S) Integration (using the double-exponential method)
Thanks Pauli.

/* 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.
04-02-2017, 10:02 PM
Post: #32
 Paul Dale Senior Member Posts: 1,761 Joined: Dec 2013
RE: (WP-34S) Integration (using the double-exponential method)
I don't seem to be having a good time of this

I'll update later today.

The code is in the 34S repository now.

Pauli
 « Next Oldest | Next Newest »

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