Re: Symbolic expand function doesn't always do anything

[Paul Kienzle]

One more point of concern: sometimes floating point
processors keep guard bits during computations.  It
might be that these extra bits are converted when the
double is converted to int64, but I wouldn't bet on it.

You can force the value to be written to and read from
memory which will strip any guard bits:

  volatile double dminusone=fabs(d)-1.;
  if (floor(d)==d && fabs(d) != dminusone) ...

I don't know  if this is necessary in practice on the
architectures that  octave runs on.  Given that
accuracy is not critical in this operation, I would
be happy with the floor(d)==d test by itself.  Maybe you
want a function which always calls
        GiNaC::numeric(double)
even if the value looks like an integer?

Hopefully the conversion between double and int64
will be efficient even on intel architectures.  We had
to do silly things with rand to get 53 bit mantissas.

- Paul

On Sep 14, 2004, at 9:47 PM, Paul Kienzle wrote:

> Better make that:
>
> if (floor(d) == d && fabs(d)-1 != fabs(d))
>         number = GiNaC::numeric(octave_int64_t(d));
> else
>         number = GiNaC::numeric(d);
>
>
> On Sep 14, 2004, at 8:32 PM, Paul Kienzle wrote:
>
> > How about:
> >
> > if (floor(d) == d && d+1 != d)
> >         number = GiNaC::numeric(int(d));
> > else
> >         number = GiNaC::numeric(d);
> >
> > - Paul
> >
> > On Sep 14, 2004, at 11:22 AM, Benjamin Sapp wrote:
> >
> > > Paul,
> > >
> > > What the best way of doing that?  Here's an idea,
> > >
> > > double d;
> > > int id;
> > > if (ceil(d) == floor(d)) {
> > >         id = (int)d;
> > >         number = GiNaC::numeric(id);
> > >         // printf("it's an integer: %d\n", id);
> > > }else{
> > >         // printf(" not an integer: %f\n", d);
> > >         number = GiNaC::numeric(d);
> > > }
> > >
> > > The problem with this is what happens when we get very large numbers?
> > > If the number is large enough then the number will always fall on an
> > > integer value...  right?  Should we also check how large the value is
> > > and only do it for small numbers?  where is that cutoff? (The 2^53 
> > > you
> > > referred to maybe?)  Maybe I need to go look at the standard.  Is 
> > > that
> > > the only other gotcha?  I don't feel confident enough in my 
> > > knowledge of
> > > floating point representation to say for sure.
> > >
> > > Thanks,
> > > Ben.
> > >
> > > On Mon, 2004-09-13 at 19:36, Paul Kienzle wrote:
> > > > Ben,
> > > >
> > > > My understanding of IEEE 754 is that integers are exactly 
> > > > representable
> > > > up
> > > > to about 2^53.  Couldn't the symbolic package treat integer values 
> > > > as
> > > > exact vpa values, and only force the user to mark the cases where 
> > > > the
> > > > value is meant to be an approximation?  I think this would capture 
> > > > more
> > > > of the common cases correctly.
> > > >
> > > > - Paul
> > > >
> > > >
> > > > On Sep 13, 2004, at 12:01 PM, Benjamin Sapp wrote:
> > > >
> > > > > Hi,
> > > > >
> > > > > Indeed, you are correct it is an octave integration peculiarity.  
> > > > > You
> > > > > see, octave interprets any thing that can be a number as a double
> > > > > first.  Then, when a double is used in a symbolic expression it is
> > > > > converted to the exact form like in ginsh.  By that time it's too 
> > > > > late
> > > > > because a 1 is already double approximation of a 1 rather than 
> > > > > exactly
> > > > > a
> > > > > 1.   You can perform your example from ginsh in octave by 
> > > > > surrounding
> > > > > each exact number with quotes and vpa().  For example on my 
> > > > > computer:
> > > > >
> > > > > octave:47> expand( (x-vpa("1"))^vpa("3"))
> > > > > ans =
> > > > >
> > > > > -1+x^3-3*x^2+3*x
> > > > >
> > > > > It's a bit cumbersome but I don't think there's a better solution.
> > > > >
> > > > > Good luck,
> > > > > Ben.
> > > > >
> > > > > On Mon, 2004-09-13 at 08:28, edA-qa mort-ora-y wrote:
> > > > > > I'm using the OctaveForge integration of the GiNaC symbolic 
> > > > > > library
> > > > > > and
> > > > > > having some trouble using the "expand" function.
> > > > > >
> > > > > > I want is to perform the polynomial expansion:
> > > > > >
> > > > > >         >x = sym( "x" )
> > > > > >         >expand( (x-1)^2 )
> > > > > >         ans = (x-1)^2
> > > > > >
> > > > > >         I want to get
> > > > > >         ans = 1 + x^2 - 2*x
> > > > > >
> > > > > > If you do this form, it works:
> > > > > >         >expand( (x-1) * (x-1) )
> > > > > >         ans = 1.0+x^2-(2.0)*x
> > > > > >
> > > > > > Using ginsh with expand produces the desired results, so I'm 
> > > > > > assuming
> > > > > > it
> > > > > > is an Octave integration peculariaty:
> > > > > >
> > > > > >         ginsh>  expand( (x-1)^3 );
> > > > > >         -1-3*x^2+x^3+3*x
> > > > >
> > > > >
> > > > >
> > > > > -------------------------------------------------------------
> > > > > Octave is freely available under the terms of the GNU GPL.
> > > > >
> > > > > Octave's home on the web:  http://www.octave.org
> > > > > How to fund new projects:  http://www.octave.org/funding.html
> > > > > Subscription information:  http://www.octave.org/archive.html
> > > > > -------------------------------------------------------------
> >
> >
> >
> > -------------------------------------------------------------
> > Octave is freely available under the terms of the GNU GPL.
> >
> > Octave's home on the web:  http://www.octave.org
> > How to fund new projects:  http://www.octave.org/funding.html
> > Subscription information:  http://www.octave.org/archive.html
> > -------------------------------------------------------------



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------