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

> - 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

On Sun, 3 Feb 2008 07:15:42 -0800 (PST), Wes <wjltemp-gg@yahoo.com>
wrote:

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

Yes, it is the best.

Damir

On Feb 4, 9:54 am, Damir <noem...@nospam.hr> wrote:
> On Sun, 3 Feb 2008 07:15:42 -0800 (PST), Wes <wjltemp...@yahoo.com>
> 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))

On Feb 4, 9:40 pm, mjc <mjco...@acm.org> 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

On Feb 4, 6:42=A0pm, Wes <wjltemp...@yahoo.com> 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. =A0This 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.

> 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

-----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+gRAkBWAKCYv3eOIq5h0JoGOzxDm8TxK4LgigCff5P1
/XOGisRf2LapYuI9HFmtFTo=
=mfUB
-----END PGP SIGNATURE-----

>
> 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/cgiwrap/hpmuseum/archv017.cgi?read=109056

Gerson.

On Feb 4, 7:48 pm, TW <timwess...@gmail.com> 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

On Feb 7, 9:49 pm, Ingo Blank <spam.and.tr...@hpgcc.org> 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

On Feb 8, 7:02 am, Wes <wjltemp...@yahoo.com> 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.

On Feb 10, 3:36 am, mjc <mjco...@acm.org> wrote:
> On Feb 8, 7:02 am, Wes <wjltemp...@yahoo.com> 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

On Feb 12, 11:25 am, Wes <wjltemp...@yahoo.com> wrote:
> On Feb 10, 3:36 am, mjc <mjco...@acm.org> wrote:
>
> > On Feb 8, 7:02 am, Wes <wjltemp...@yahoo.com> 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:)

mjc <mjcohen@acm.org> 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"?
;-)

On Feb 13, 3:56 pm, Joe Veazey <nob...@nowhere.nodomain> wrote:
> mjc <mjco...@acm.org> 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

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" <wjltemp-gg@yahoo.com> wrote in message 
news:0cad709a-1643-423c-9896-35d9bbdf3ea8@v46g2000hsv.googlegroups.com...
> On Feb 4, 7:48 pm, TW <timwess...@gmail.com> 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 

On Aug 23, 9:17 am, "Veli-Pekka Nousiainen"
<velipekka.nousiai...@saunalahti.fi> wrote:

> with this in mind and considering plotting are you now willing to think a=
bout
> 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=A0am, "Veli-Pekka Nousiainen"
<velipekka.nousiai...@saunalahti.fi> 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
> =A0 =A0 including the =A0FIX 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" <wjltemp...@yahoo.com> wrote in message
>
> news:0cad709a-1643-423c-9896-35d9bbdf3ea8@v46g2000hsv.googlegroups.com...
>
> > On Feb 4, 7:48 pm, TW <timwess...@gmail.com> wrote:
> >> > - Are the tolerance options too confusing or unnecessary? =A0With th=
e
> >> > very short run-times, it might be better to just default to 1e-12.
>
> >> For me that would probably be best. =A0What is the difference in size =
of
> >> code if you remove it? =A0My experience has been that it generally
> >> doesn't make a whole lot of difference.
>
> > Probably not a significant difference. =A0I 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. :) =A0Honestly, I've always thought it rather annoyin=
g
> > 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. =A0For example, suppose you wanted to graph
> > f(x)=3DS(0,x,g(t),t). =A0With 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