# HP Forums

You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2
This is my implementation of the double exponential integration method for the wp34s.

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

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

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

Pauli
My thumb is up. It could add to the list of WP34s 'uniquenesses'.
BTW, the license and copyright header reads: "M. César". That "M" must stand for marvellous... or may be mathematician... or both in one single letter, by some sort of arcane&syncopated algebraic expression ;O)
(03-25-2017 12:33 AM)Paul Dale Wrote: [ -> ]What do people think about including the last one in the 34S firmware?

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

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

It does the very same thing, in this respect, that the built in integrator: corrupt all the stack. I was able to write a version (not to be used in XROM) preserving the stack, but I found some trouble when porting it to XROM. Perhaps I need some advice here. My approach was:
• Use 2 local register to backup Z & T
• When returning, build a stack with the backed up T & Z, the result in Y and a dummy value in X
• Then do a DROP and return
It worked if the user was working in SSIZE4, but not in SSIZE8.

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

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

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

Jaja, Luigi, la M solo es de Manuel.
(03-25-2017 11:52 PM)emece67 Wrote: [ -> ]It does the very same thing, in this respect, that the built in integrator: corrupt all the stack.

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

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

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

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

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

- Pauli
(03-24-2017 02:23 PM)emece67 Wrote: [ -> ]cut

I wish more programs posted here would have exhaustive comments like yours instead of just the program (that in RPN / RPL is unlikely to be "self descriptive")
(03-24-2017 02:21 PM)emece67 Wrote: [ -> ]This is my implementation of the double exponential integration method for the wp34s.

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

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

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

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

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

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

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

Thus, I've modified my program to behave more or less in the same way, so it now (v1.1r and v1.1x):
1. Returns the computed result in X (as before)
2. Returns in Y the previous approximation to the result
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L (as before)
5. Keeps ABCD untouched (no matter if the user is in SSIZE4 or SSIZE8)

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

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

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

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

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

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

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

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

In my local builds that's not true, cause I have a modified CONST menu, so I try to avoid such constructs.
(03-26-2017 06:34 PM)emece67 Wrote: [ -> ]Surely that's the behavior of the 34C, but the 34s:
1. Returns the computed result in X
2. Returns in Y what seems to be some kind of a previous approximation to the result (it is not much smaller than the result, as expected for an error estimation, but approximately equal to it)
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L

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

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

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

(03-26-2017 06:34 PM)emece67 Wrote: [ -> ]Thus, I've modified my program to behave more or less in the same way, so it now (v1.1r and v1.1x):
1. Returns the computed result in X (as before)
2. Returns in Y the previous approximation to the result
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L (as before)
5. Keeps ABCD untouched (no matter if the user is in SSIZE4 or SSIZE8)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Regards.

Edited: removing the comment about ABCD remaining untouched.
(03-27-2017 05:42 PM)emece67 Wrote: [ -> ]Keeps ABCD untouched

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

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

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

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

- Pauli
Program description modified removing the statement about ABCD remaining untouched.
It would be prudent to also remove the ->A..D and A..D-> commands as well.

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

Regards,
Pages: 1 2
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :