[SciPy-user] how to make a sinc(x) function - "divide by zero" !

[Sebastian Haase]

Then there is a BIG (?) difference in how "where" is implemented in Numeric 
and numarray
Unless, Numeric defaults to "ignore" divide-by-zero ...

- Sebastian



On Monday 12 December 2005 11:05, Scott Ransom wrote:
> Short answer is no.  At least not for arrays:
>
> In [25]:x = asarray([-0.1, -0.05, -0.01, -0.005, -0.001, 0.0, \
> 0.001, 0.005, 0.01, 0.05, 0.1])
>
> In [26]:sinc(x)
> Out[26]:
> array([ 0.98363164,  0.99589274,  0.99983551,  0.99995888,
> 	0.99999836,  1.        ,  0.99999836,  0.99995888,
> 	0.99983551,  0.99589274,  0.98363164])
>
> Scalars, however, seem to be a different story.  I guess that I
> must have only used this function with arrays...
>
> In [27]:sinc(0.0)
> ---------------------------------------------------------------------------
> exceptions.ZeroDivisionError                         Traceback
> (most recent call last)
>
> /home/ransom/<console>
>
> /home/ransom/<console> in sinc(xs)
>
> ZeroDivisionError: float division
>
>
> Scott
>
> On Mon, Dec 12, 2005 at 10:30:08AM -0800, Sebastian Haase wrote:
> > Scott,
> > Thanks for your reply. (As Perry points out - my code actually works on
> > numarray)
> > But in your code: Don't you get a "division by zero" in Numeric here ?
> > "where" is evaluation BOTH* results FIRST (over the ENTIRE array) and
> > only AFTERWARDS "selects" based on the on the first argument
> > (* actually all "three" arrays)
> >
> > Thanks,
> > Sebastian Haase
> >
> > On Monday 12 December 2005 10:04, Scott Ransom wrote:
> > > Here is a simple version that I use.  It uses the Taylor
> > > expansion when the argument is near zero.  Note that this
> > > is using Numeric and umath right now.  This should be
> > > trivial to change to new scipy_core (which I am not using yet):
> > >
> > > def sinc(xs):
> > >     """
> > >     sinc(xs):
> > >         Return the sinc function [i.e. sin(pi * xs)/(pi * xs)]
> > >             for the values xs.
> > >     """
> > >     pxs = umath.pi*xs
> > >     return Numeric.where(umath.fabs(pxs)<1e-3,
> > >                          1.0-pxs*pxs/6.0, umath.sin(pxs)/pxs)
> > >
> > > Scott
> > >
> > > On Monday 12 December 2005 12:39 pm, Arnd Baecker wrote:
> > > > On Mon, 12 Dec 2005, Sebastian Haase wrote:
> > > > > Hi,
> > > > > I was trying to implement a "sinc" [sin(x)/x] in numarray. But the
> > > > > "where"-semantics makes it choke on the x=0: 0/0 case.
> > > > > (This should of course work for x being an array - so "if" is no
> > > > > option) The best I could come up with was:
> > > > >
> > > > > def sinc(r):
> > > > >     na.Error.pushMode(all="ignore")
> > > > >     a = na.where(r, na.divide(na.sin(r),r), 1)
> > > > >     na.Error.popMode()
> > > > >     return a
> > > > >
> > > > > but I still seem to get  a warning...
> > > > >
> > > > > >>> F.sinc(0)
> > > > >
> > > > > Warning: Encountered invalid numeric result(s)  in divide
> > > > > 1.0
> > > > >
> > > > > What is a better way of doing this ?
> > > >
> > > > Yuo could try to use `vectorize` (warning: untested code)
> > > >
> > > > def sinc(x):
> > > >     if x==0.0:                # presumably better to check for small
> > > > x return 0.0            # here ...
> > > >     else:
> > > >        return sin(x)/x
> > > >
> > > > sinc_vectorized=scipy.vectorize(sinc)
> > > >
> > > > HTH,
> > > >
> > > > Arnd
> > > >
> > > > _______________________________________________
> > > > SciPy-user mailing list
> > > > SciPy-user at scipy.net
> > > > http://www.scipy.net/mailman/listinfo/scipy-user
> >
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user at scipy.net
> > http://www.scipy.net/mailman/listinfo/scipy-user
>
> --