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

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
GaussKronrod integration program. Even with no floatingpoint
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 backwardcompatible with the
builtin 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 builtin 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 1e7
STD > relative tolerance of 1e12
FIX 6 > absolute tolerance of 1e7
Examples of tolerance from ITOL or stack:
0.001 > relative tolerance of 0.001
6 > relative tolerance of 1e7
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 APPhysics 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 runtimes, it might be better to just default to 1e12.
 Making FIX mode specify absolute tolerance makes it slightly
incompatible with the 49/50 builtin method. Is this worth
resurrecting from the past?
 Any other suggestions or thoughts?
wes

Re: hpgcc integration program
>  Are the tolerance options too confusing or unnecessary? With the
> very short runtimes, it might be better to just default to 1e12.
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

Re: hpgcc integration program
On Sun, 3 Feb 2008 07:15:42 0800 (PST), Wes
wrote:
> With the
>very short runtimes, it might be better to just default to 1e12.
Yes, it is the best.
Damir

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 runtimes, it might be better to just default to 1e12.
>
> 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(ba) / 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))

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 GaussKronrod 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 1e12
S(0,2,exp(x^2),x)
Adaptive Simpson's Rule: 4065 evaluations
my GaussKronrod: 21 (minimum possible)
S(0,sqrt(4*pi),sin(x^2),x)
Adaptive Simpson's Rule: 19797 evaluations
my GaussKronrod: 147
Maxima's GaussKronrod: 63
Another important things is that with GaussKronrod, the endpoints
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, GaussKronrod is the way to go.
wes

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 GaussKronrod: 147
> Maxima's GaussKronrod: 63
>
> Another important things is that with GaussKronrod, the endpoints
> 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, GaussKronrod is the way to go.
>
> wes
Just out of curiosity, would you try S(0,1,(sqrt(x)/(x1)  1/
ln(x)),x)?
For comparison, the actual result to 15 figures is
3.64899739785761e2.
Thanks,
Gerson.

Re: hpgcc integration program
> Just out of curiosity, would you try S(0,1,(sqrt(x)/(x1)  1/
> ln(x)),x)?
> For comparison, the actual result to 15 figures is
> 3.64899739785761e2.
>
> Thanks,
>
> Gerson.
3.64899739786e2 + 2.59101085497e13 in 0.3365 seconds (TEVAL)
for comparison with builtin:
3.64899739831E2 + 3.64899739883E13 in 2261.9767 seconds.
Interestingly, the builtin's interval does not contain the actual
value.
wes

Re: hpgcc integration program
BEGIN PGP SIGNED MESSAGE
Hash: SHA1
Wes wrote:
>> Just out of curiosity, would you try S(0,1,(sqrt(x)/(x1)  1/
>> ln(x)),x)?
>> For comparison, the actual result to 15 figures is
>> 3.64899739785761e2.
>>
>> Thanks,
>>
>> Gerson.
>
> 3.64899739786e2 + 2.59101085497e13 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 builtin:
> 3.64899739831E2 + 3.64899739883E13 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 ReplyTo address.
BEGIN PGP SIGNATURE
Version: GnuPG v2.0.4svn0 (GNU/Linux)
Comment: Using GnuPG with SUSE  http://enigmail.mozdev.org
iD8DBQFHq1K+6ricFBSet+gRAkBWAKCYv3eOIq5h0JoGOzxDm8 TxK4LgigCff5P1
/XOGisRf2LapYuI9HFmtFTo=
=mfUB
END PGP SIGNATURE

Re: hpgcc integration program
>
> 3.64899739786e2 + 2.59101085497e13 in 0.3365 seconds (TEVAL)
>
> for comparison with builtin:
> 3.64899739831E2 + 3.64899739883E13 in 2261.9767 seconds.
>
> Interestingly, the builtin'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/cgisys/cgiw...gi?read=109056
Gerson.

Re: hpgcc integration program
On Feb 4, 7:48 pm, TW wrote:
> >  Are the tolerance options too confusing or unnecessary? With the
> > very short runtimes, it might be better to just default to 1e12.
>
> 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 1e12 with optional arguments to override the
default sounds like the best option.
wes

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

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  rungekutta or relative.

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  rungekutta 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.
> rungekutta or relative.
I'm familiar with rungekutta, but what is "relative"?
wes

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  rungekutta 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.
>
> > rungekutta or relative.
>
> I'm familiar with rungekutta, but what is "relative"?
>
> wes
"Relative of the RungeKutta method"  merson or uncle or motherin
law

Re: hpgcc integration program
mjc wrote in
news:95c281c95fe3491d85cbc4804f0b89e0@c4g2000hsg.googlegroups.com:
>>
>> > rungekutta or relative.
>>
>> I'm familiar with rungekutta, but what is "relative"?
>>
>> wes
>
> "Relative of the RungeKutta method"  merson or uncle or motherin
> law
>
I'm familar with uncle or motherinlaw, but what is "merson"?
;)

Re: hpgcc integration program
On Feb 13, 3:56 pm, Joe Veazey wrote:
> mjc wrote innews:95c281c95fe3491d85cbc4804f0b89e0@c4g2000hsg.googlegroups.com:
>
>
>
> >> > rungekutta or relative.
>
> >> I'm familiar with rungekutta, but what is "relative"?
>
> >> wes
>
> > "Relative of the RungeKutta method"  merson or uncle or motherin
> > law
>
> I'm familar with uncle or motherinlaw, but what is "merson"?
> ;)
GIYF
A variant of the RungeKutta method.
See http://eom.springer.de/K/k056060.htm

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 buildin 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:0cad709a1643423c989635d9bbdf3ea8@v46g2000hsv.googlegroups.com...
> On Feb 4, 7:48 pm, TW wrote:
>> >  Are the tolerance options too confusing or unnecessary? With the
>> > very short runtimes, it might be better to just default to 1e12.
>>
>> 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 1e12 with optional arguments to override the
> default sounds like the best option.
>
> wes

Re: hpgcc integration program
On Aug 23, 9:17 am, "VeliPekka 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
HP34C, 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 buildin 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 builtin
integrator. Since the whole point of my program was to implement a
GaussKronrod 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 1e12
and a default max depth of 200 and the result is then rounded to the
50g's max relative tolerance of 1e11.
> 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 hpgcc3beta 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, "VeliPekka 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 buildin 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:0cad709a1643423c989635d9bbdf3ea8@v46g2000hsv.googlegroups.com...
>
> > On Feb 4, 7:48 pm, TW wrote:
> >> >  Are the tolerance options too confusing or unnecessary? *With the
> >> > very short runtimes, it might be better to just default to 1e12.
>
> >> 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 1e12 with optional arguments to override the
> > default sounds like the best option.
>
> > wes