hpgcc integration program - Hewlett Packard

This is a discussion on hpgcc integration program - Hewlett Packard ; Hi All, A few weeks ago, there was a discussion here about numeric integration on the 50g in which I questioned the appropriateness of the Romberg method. Romberg was perfect for the tiny memory available on the 34C, but now ...

+ Reply to Thread
Results 1 to 18 of 18

Thread: hpgcc integration program

  1. hpgcc integration program

    Hi All,

    A few weeks ago, there was a discussion here about numeric integration
    on the 50g in which I questioned the appropriateness of the Romberg
    method. Romberg was perfect for the tiny memory available on the 34C,
    but now with more memory available, perhaps it's time to move to
    another method. There are libraries already available, but none that
    I'm aware of that use ARM code. With the power of HPGCC and HPStack/
    HPParser, I wanted to see what a C program could do.

    I've been playing around with it and I've now got a (mostly) working
    Gauss-Kronrod integration program. Even with no floating-point
    hardware, the speed is quite good. For instance, S(0,1,'1/
    SQRT(X))','X') which normally takes 25 minutes in STD mode, and
    S(1,0,'LN(X)^5','X') which normally takes 36 minutes in STD mode, each
    take less than 0.4 seconds. (I'm referring to numeric integration in
    approx mode, not symbolic in exact mode.)

    The program still needs some work, but I wanted to run my ideas about
    its interface past the group. I want to keep the default options as
    simple as possible and be more or less backward-compatible with the
    built-in numeric integrator, but with some useful options. I'd
    welcome comments on the following interface ideas.

    Stack diagram:
    lower upper expr var --> area error
    lower upper expr {var tolerance } --> area error
    lower upper expr {var tolerance maxdepth} --> area error

    The error can be stored in IERR and LASTARG via an RPL wrapper
    program. Like the built-in integration, if convergence is not
    achieved, then the returned error is negative.

    The display mode is used to determine the tolerance (like the built-
    in); or you may override the display by using the ITOL variable; or
    you may override ITOL by specifying the tolerance on the stack.

    If the tolerance value is positive, it will be interpreted as a
    relative tolerance, while a value <= 0 will be interpreted as an
    absolute tolerance. A specified tolerance value N >= 1 will be
    interpreted the same as SCI N. If no tolerance is specified, SCI N or
    ENG N will specify a relative tolerance of 10^-(N+1), STD will specify
    a relative tolerance of 10^-12, and FIX N will specify an absolute
    tolerance of 10^-(N+1) (a la 41C Advantage Pac, 15C, 34C).

    Examples of tolerance from Display Mode:
    SCI 6 --> relative tolerance of 1e-7
    STD --> relative tolerance of 1e-12
    FIX 6 --> absolute tolerance of 1e-7

    Examples of tolerance from ITOL or stack:
    0.001 --> relative tolerance of 0.001
    6 --> relative tolerance of 1e-7
    -0.01 --> absolute tolerance of 0.01 (to the hundredths place)
    -100 --> absolute tolerance of 100 (to the hundreds place)
    0 --> absolute tolerance of 0 (nearest integer)

    Most people probably will just use the Display Mode and not bother
    with anything else. Others may want to set a default value in ITOL in
    the HOME directory and that becomes the default tolerance. The AP-
    Calculus exam requires answers to a certain number of decimal places
    (absolute error), while the AP-Physics exam wants the proper number of
    significant figures (relative error). A student taking both courses
    may set different ITOL variables in his CALC and PHYS directories as
    defaults and override them on the stack when necessary.

    The default maximum depth of subdivisions can be overridden by an
    IDEPTH variable, which can be overridden by the specifying maxdepth on
    the stack.

    Questions:
    - Are the tolerance options too confusing or unnecessary? With the
    very short run-times, it might be better to just default to 1e-12.
    - Making FIX mode specify absolute tolerance makes it slightly
    incompatible with the 49/50 built-in method. Is this worth
    resurrecting from the past?
    - Any other suggestions or thoughts?

    -wes

  2. Re: hpgcc integration program

    > - Are the tolerance options too confusing or unnecessary? With the
    > very short run-times, it might be better to just default to 1e-12.


    For me that would probably be best. What is the difference in size of
    code if you remove it? My experience has been that it generally
    doesn't make a whole lot of difference. With the speed, will anyone
    ever specify anything else? Any instances where you wouldn't want
    full out to 12?

    TW

  3. Re: hpgcc integration program

    On Sun, 3 Feb 2008 07:15:42 -0800 (PST), Wes
    wrote:

    > With the
    >very short run-times, it might be better to just default to 1e-12.


    Yes, it is the best.

    Damir

  4. Re: hpgcc integration program

    On Feb 4, 9:54 am, Damir wrote:
    > On Sun, 3 Feb 2008 07:15:42 -0800 (PST), Wes
    > wrote:
    >
    > > With the
    > >very short run-times, it might be better to just default to 1e-12.

    >
    > Yes, it is the best.
    >
    > Damir


    It might be interesting to compare something as simple as an adaptive
    Simpson's rule.

    Here is Python code from Wikipedia:

    def simpsons_rule(f,a,b):
    c = (a+b) / 2.0
    h3 = abs(b-a) / 6.0
    return h3*(f(a) + 4.0*f(c) + f(b))

    def recursive_asr(f,a,b,eps,sum):
    "Recursive implementation of adaptive Simpson's rule."
    c = (a+b) / 2.0
    left = simpsons_rule(f,a,c)
    right = simpsons_rule(f,c,b)
    if abs(left + right - sum) <= 15*eps:
    return left + right + (left + right - sum)/15
    return recursive_asr(f,a,c,eps/2,left) + recursive_asr(f,c,b,eps/
    2,right)

    def adaptive_simpsons_rule(f,a,b,eps):
    "Calculate integral of f from a to b with max error of eps."
    return recursive_asr(f,a,b,eps,simpsons_rule(f,a,b))

  5. Re: hpgcc integration program

    On Feb 4, 9:40 pm, mjc wrote:

    > It might be interesting to compare something as simple as an adaptive
    > Simpson's rule.
    >
    > Here is Python code from Wikipedia:


    One nice thing about Gauss-Kronrod is that it takes an amazingly few
    function evaluations to come up with very good estimates. I entered
    the given Adaptive Simpson's Rule code into Maxima for a comparison
    with a tolerance of 1e-12

    S(0,2,exp(-x^2),x)
    Adaptive Simpson's Rule: 4065 evaluations
    my Gauss-Kronrod: 21 (minimum possible)


    S(0,sqrt(4*pi),sin(x^2),x)
    Adaptive Simpson's Rule: 19797 evaluations
    my Gauss-Kronrod: 147
    Maxima's Gauss-Kronrod: 63


    Another important things is that with Gauss-Kronrod, the end-points
    are never evaluated. This allows you to integrate S(a,b,f(x),x) even
    if f(a) or f(b) are undefined, such an asymptote or hole, such as
    S(0,1,1/sqrt(x),x).

    From everything I've read, Gauss-Kronrod is the way to go.

    -wes

  6. Re: hpgcc integration program

    On Feb 4, 6:42*pm, Wes wrote:

    >
    > S(0,sqrt(4*pi),sin(x^2),x)
    > Adaptive Simpson's Rule: 19797 evaluations
    > my Gauss-Kronrod: 147
    > Maxima's Gauss-Kronrod: 63
    >
    > Another important things is that with Gauss-Kronrod, the end-points
    > are never evaluated. *This allows you to integrate S(a,b,f(x),x) even
    > if f(a) or f(b) are undefined, such an asymptote or hole, such as
    > S(0,1,1/sqrt(x),x).
    >
    > From everything I've read, Gauss-Kronrod is the way to go.
    >
    > -wes


    Just out of curiosity, would you try S(0,1,(sqrt(x)/(x-1) - 1/
    ln(x)),x)?
    For comparison, the actual result to 15 figures is
    3.64899739785761e-2.

    Thanks,

    Gerson.

  7. Re: hpgcc integration program

    > Just out of curiosity, would you try S(0,1,(sqrt(x)/(x-1) - 1/
    > ln(x)),x)?
    > For comparison, the actual result to 15 figures is
    > 3.64899739785761e-2.
    >
    > Thanks,
    >
    > Gerson.


    3.64899739786e-2 +- 2.59101085497e-13 in 0.3365 seconds (TEVAL)

    for comparison with built-in:
    3.64899739831E-2 +- 3.64899739883E-13 in 2261.9767 seconds.

    Interestingly, the built-in's interval does not contain the actual
    value.

    -wes

  8. Re: hpgcc integration program

    -----BEGIN PGP SIGNED MESSAGE-----
    Hash: SHA1

    Wes wrote:
    >> Just out of curiosity, would you try S(0,1,(sqrt(x)/(x-1) - 1/
    >> ln(x)),x)?
    >> For comparison, the actual result to 15 figures is
    >> 3.64899739785761e-2.
    >>
    >> Thanks,
    >>
    >> Gerson.

    >
    > 3.64899739786e-2 +- 2.59101085497e-13 in 0.3365 seconds (TEVAL)


    I assume you are calling sys_slowOff() (i.e. running @75 MHz) ?

    You will then enjoy the 192 MHz overclocking feature in HPGCC3, which
    will boost the performance by x ~2.5 ...

    Our beta testers just reported this factor on average with their
    computing intensive applications like matrix inversion or factorization.

    We are pleased to see another great application made with HPGCC!


    >
    > for comparison with built-in:
    > 3.64899739831E-2 +- 3.64899739883E-13 in 2261.9767 seconds.
    >


    > -wes



    - --
    Ingo Blank

    http://hpgcc.org
    http://blog.hpgcc.org

    To reply by email, please replace the word "spam" with my first name in
    the Reply-To address.
    -----BEGIN PGP SIGNATURE-----
    Version: GnuPG v2.0.4-svn0 (GNU/Linux)
    Comment: Using GnuPG with SUSE - http://enigmail.mozdev.org

    iD8DBQFHq1K+6ricFBSet+gRAkBWAKCYv3eOIq5h0JoGOzxDm8 TxK4LgigCff5P1
    /XOGisRf2LapYuI9HFmtFTo=
    =mfUB
    -----END PGP SIGNATURE-----

  9. Re: hpgcc integration program

    >
    > 3.64899739786e-2 +- 2.59101085497e-13 in 0.3365 seconds (TEVAL)
    >
    > for comparison with built-in:
    > 3.64899739831E-2 +- 3.64899739883E-13 in 2261.9767 seconds.
    >
    > Interestingly, the built-in's interval does not contain the actual
    > value.
    >
    > -wes


    This is really a tough integral for numeric integrators:

    http://www.cs.berkeley.edu/~wkahan/Math128/INTGTkey.pdf

    The exact result is 2 - ln(4) - y, where y is Euler's constant
    (0.577215664901533..). Details in this draft:

    http://www.hpmuseum.org/cgi-sys/cgiw...gi?read=109056

    Gerson.


  10. Re: hpgcc integration program

    On Feb 4, 7:48 pm, TW wrote:
    > > - Are the tolerance options too confusing or unnecessary? With the
    > > very short run-times, it might be better to just default to 1e-12.

    >
    > For me that would probably be best. What is the difference in size of
    > code if you remove it? My experience has been that it generally
    > doesn't make a whole lot of difference.


    Probably not a significant difference. I had compatibility in mind,
    but that seems less and less important the more I think about it.
    Perhaps the excitement of figuring out how to determine the display
    settings from within my C program blinded me to whether or not it
    _should_ be done. Honestly, I've always thought it rather annoying
    to have the display settings determine the integration tolerance.

    > With the speed, will anyone ever specify anything else?
    > Any instances where you wouldn't want full out to 12?


    Probably not too often, but sometimes top speed is more important than
    accuracy. For example, suppose you wanted to graph
    f(x)=S(0,x,g(t),t). With the low screen resolution, the tolerance
    could be quite lax so as to increase plotting speed.

    Another time you might want to control the tolerance is experimenting
    with badly behaved functions.

    I think defaulting to 1e-12 with optional arguments to override the
    default sounds like the best option.

    -wes

  11. Re: hpgcc integration program

    On Feb 7, 9:49 pm, Ingo Blank wrote:

    > I assume you are calling sys_slowOff() (i.e. running @75 MHz) ?


    I'm using a custom crt0.o in which I commented out the sys_slowOn().
    When graphing a C function, there's much less screen flicker than just
    using sys_slowOff().

    > You will then enjoy the 192 MHz overclocking feature in HPGCC3, which
    > will boost the performance by x ~2.5 ...


    Sounds intriguing. Can you still access the keyboard at that speed?
    Do you get much screen flicker switching back and forth?

    > We are pleased to see another great application made with HPGCC!


    Well, it's got a ways to go before it's a finished application. And
    as far as it being "great," that remains to be seen. But I'm very
    encouraged by the early results.

    -wes

  12. Re: hpgcc integration program

    On Feb 8, 7:02 am, Wes wrote:
    ......
    > Probably not too often, but sometimes top speed is more important than
    > accuracy. For example, suppose you wanted to graph
    > f(x)=S(0,x,g(t),t). With the low screen resolution, the tolerance
    > could be quite lax so as to increase plotting speed.

    ......
    > -wes


    In this case, it is probably better to integrate it as a differential
    equation - runge-kutta or relative.

  13. Re: hpgcc integration program

    On Feb 10, 3:36 am, mjc wrote:
    > On Feb 8, 7:02 am, Wes wrote:
    > .....
    >
    > > Probably not too often, but sometimes top speed is more important than
    > > accuracy. For example, suppose you wanted to graph
    > > f(x)=S(0,x,g(t),t). With the low screen resolution, the tolerance
    > > could be quite lax so as to increase plotting speed.

    > .....
    > > -wes

    >
    > In this case, it is probably better to integrate it as a differential
    > equation - runge-kutta or relative.


    Yes, you're quite right. That example would be better graphed as a
    Diff.Eq.

    Last month I was working on a physics problem which involved something
    like:

    f(x) = integral from 0 to h(x) of g(x,t) dt

    I don't think this one could be done in this way (?), or at least not
    easily.

    > runge-kutta or relative.


    I'm familiar with runge-kutta, but what is "relative"?

    -wes

  14. Re: hpgcc integration program

    On Feb 12, 11:25 am, Wes wrote:
    > On Feb 10, 3:36 am, mjc wrote:
    >
    > > On Feb 8, 7:02 am, Wes wrote:
    > > .....

    >
    > > > Probably not too often, but sometimes top speed is more important than
    > > > accuracy. For example, suppose you wanted to graph
    > > > f(x)=S(0,x,g(t),t). With the low screen resolution, the tolerance
    > > > could be quite lax so as to increase plotting speed.

    > > .....
    > > > -wes

    >
    > > In this case, it is probably better to integrate it as a differential
    > > equation - runge-kutta or relative.

    >
    > Yes, you're quite right. That example would be better graphed as a
    > Diff.Eq.
    >
    > Last month I was working on a physics problem which involved something
    > like:
    >
    > f(x) = integral from 0 to h(x) of g(x,t) dt
    >
    > I don't think this one could be done in this way (?), or at least not
    > easily.
    >
    > > runge-kutta or relative.

    >
    > I'm familiar with runge-kutta, but what is "relative"?
    >
    > -wes


    "Relative of the Runge-Kutta method" - merson or uncle or mother-in-
    law

  15. Re: hpgcc integration program

    mjc wrote in
    news:95c281c9-5fe3-491d-85cb-c4804f0b89e0@c4g2000hsg.googlegroups.com:

    >>
    >> > runge-kutta or relative.

    >>
    >> I'm familiar with runge-kutta, but what is "relative"?
    >>
    >> -wes

    >
    > "Relative of the Runge-Kutta method" - merson or uncle or mother-in-
    > law
    >


    I'm familar with uncle or mother-in-law, but what is "merson"?
    ;-)

  16. Re: hpgcc integration program

    On Feb 13, 3:56 pm, Joe Veazey wrote:
    > mjc wrote innews:95c281c9-5fe3-491d-85cb-c4804f0b89e0@c4g2000hsg.googlegroups.com:
    >
    >
    >
    > >> > runge-kutta or relative.

    >
    > >> I'm familiar with runge-kutta, but what is "relative"?

    >
    > >> -wes

    >
    > > "Relative of the Runge-Kutta method" - merson or uncle or mother-in-
    > > law

    >
    > I'm familar with uncle or mother-in-law, but what is "merson"?
    > ;-)


    GIYF

    A variant of the Runge-Kutta method.

    See http://eom.springer.de/K/k056060.htm

  17. Re: hpgcc integration program

    While the news message is relatively old
    I'd still like to present my thoughts:

    1) it should be totally compatible with the build-in integrator
    including the FIX mode
    2) optional variables could be located in CASDIR directory

    Nice work...
    have you tested double or triple integral speeds?
    with this in mind and considering plotting
    are you now willing to think about using the old "display method"
    as The Method to specify accuracy?
    Thank you!
    VPN

    "Wes" wrote in message
    news:0cad709a-1643-423c-9896-35d9bbdf3ea8@v46g2000hsv.googlegroups.com...
    > On Feb 4, 7:48 pm, TW wrote:
    >> > - Are the tolerance options too confusing or unnecessary? With the
    >> > very short run-times, it might be better to just default to 1e-12.

    >>
    >> For me that would probably be best. What is the difference in size of
    >> code if you remove it? My experience has been that it generally
    >> doesn't make a whole lot of difference.

    >
    > Probably not a significant difference. I had compatibility in mind,
    > but that seems less and less important the more I think about it.
    > Perhaps the excitement of figuring out how to determine the display
    > settings from within my C program blinded me to whether or not it
    > _should_ be done. Honestly, I've always thought it rather annoying
    > to have the display settings determine the integration tolerance.
    >
    >> With the speed, will anyone ever specify anything else?
    >> Any instances where you wouldn't want full out to 12?

    >
    > Probably not too often, but sometimes top speed is more important than
    > accuracy. For example, suppose you wanted to graph
    > f(x)=S(0,x,g(t),t). With the low screen resolution, the tolerance
    > could be quite lax so as to increase plotting speed.
    >
    > Another time you might want to control the tolerance is experimenting
    > with badly behaved functions.
    >
    > I think defaulting to 1e-12 with optional arguments to override the
    > default sounds like the best option.
    >
    > -wes




  18. Re: hpgcc integration program

    On Aug 23, 9:17 am, "Veli-Pekka Nousiainen"
    wrote:

    > with this in mind and considering plotting are you now willing to think about
    > using the old "display method" as The Method to specify accuracy?


    The more I've thought about the issue, the more I agree with Tim's
    comment: "With the speed, will anyone ever specify anything else?"
    Also, I don't really like the idea of the display mode changing the
    answer of a calculation. We certainly wouldn't like it if SQRT or SIN
    behaved this way. I can understand the original logic with the
    HP-34C, but with today's speed, I think it's no longer needed. The
    defaults can still be overridden by optional arguments, so the trade-
    off with speed and accuracy is still available.

    > 1) it should be totally compatible with the build-in integrator
    > including the FIX mode


    To be totally compatible, I'd have to use the exact same Romberg
    implementation so as to get the same result as the built-in
    integrator. Since the whole point of my program was to implement a
    Gauss-Kronrod method, I won't be getting the exact same results
    anyway, so I might as well use the better results. With that in mind,
    I'm currently calculating with a default relative tolerance of 1e-12
    and a default max depth of 200 and the result is then rounded to the
    50g's max relative tolerance of 1e-11.

    > 2) optional variables could be located in CASDIR directory


    How about looking in the current directory, then PATH, then CASDIR.
    That way a user could have different variable values in different
    directories if he wanted. I have separate PHYSICS and CALCULUS
    directories -- physics class requires relative tolerance (significant
    figures) while the calculus AP exam requires absolute tolerance
    (decimal places).

    > have you tested double or triple integral speeds?


    No, but that's a good suggestion. The way the program is currently
    written, I'd have to have separate double integral and triple integral
    programs. I think a better way would be to have a generic program
    that could do any level of nested integrals. You never know when you
    might want that quintuple integral.

    I used the program in class quite successfully all last semester.
    There are a few bugs with the called HPStack/HPParser libraries that
    need to be fixed before I can release anything. (I also recently
    signed up for the hpgcc3-beta and am looking forward using it for this
    program.)

    I've sometimes wondered about the need for such a program. However,
    at the end of last year, my students were taking a standardized exam,
    one section of which required them to answer 17 question in 50
    minutes. One of the questions required numeric integration of a
    function containing a cusp. In standard mode, the 50g takes 33
    minutes to complete. Of course, they didn't need such high precision
    and using FIX 5 would give an acceptable answer in 30 seconds, but
    it's quite likely that many (most?) 50g users don't know about this
    feature as it is not documented except in the AUR. (FWIW, my program
    takes 0.22 seconds.)

    Thanks for taking interest. Now that school has started up again, I
    have some motivation to dust off the code and make some improvements.

    -wes


    On Aug 23, 9:17*am, "Veli-Pekka Nousiainen"
    wrote:
    > While the news message is relatively old
    > I'd still like to present my thoughts:
    >
    > 1) it should be totally compatible with the build-in integrator
    > * * including the *FIX mode
    > 2) optional variables could be located in CASDIR directory
    >
    > Nice work...
    > have you tested double or triple integral speeds?
    > with this in mind and considering plotting
    > are you now willing to think about using the old "display method"
    > as The Method to specify accuracy?
    > Thank you!
    > VPN
    >
    > "Wes" wrote in message
    >
    > news:0cad709a-1643-423c-9896-35d9bbdf3ea8@v46g2000hsv.googlegroups.com...
    >
    > > On Feb 4, 7:48 pm, TW wrote:
    > >> > - Are the tolerance options too confusing or unnecessary? *With the
    > >> > very short run-times, it might be better to just default to 1e-12.

    >
    > >> For me that would probably be best. *What is the difference in size of
    > >> code if you remove it? *My experience has been that it generally
    > >> doesn't make a whole lot of difference.

    >
    > > Probably not a significant difference. *I had compatibility in mind,
    > > but that seems less and less important the more I think about it.
    > > Perhaps the excitement of figuring out how to determine the display
    > > settings from within my C program blinded me to whether or not it
    > > _should_ be done. *Honestly, I've always thought it rather annoying
    > > to have the display settings determine the integration tolerance.

    >
    > >> With the speed, will anyone ever specify anything else?
    > >> Any instances where you wouldn't want full out to 12?

    >
    > > Probably not too often, but sometimes top speed is more important than
    > > accuracy. *For example, suppose you wanted to graph
    > > f(x)=S(0,x,g(t),t). *With the low screen resolution, the tolerance
    > > could be quite lax so as to increase plotting speed.

    >
    > > Another time you might want to control the tolerance is experimenting
    > > with badly behaved functions.

    >
    > > I think defaulting to 1e-12 with optional arguments to override the
    > > default sounds like the best option.

    >
    > > -wes



+ Reply to Thread