Value of Sinc at complex infinity.

[pjgaudre]

I was wondering if we could add the values of the Sinc function at the complex infinities. In #14854, there was some talk to push this in a separate pull request. However, as @simonbyrne mentioned:

If we're going to fix it, we should also fix the cases where it returns incorrect NaNs for finite values as well (e.g. sinc(0.0 +250.0*im)).

However this is probably a bit of work: it will require a bit of complex analysis to figure out what those limits should be (unless you can find a good reference), and fix some of the other incorrect functions (e.g. sinpi(1+300*im), as well as adding tests and docs.

Although I am not quite sure how to remedy the incorrect NaNs for finite values. I have computed the complex limits for the Sinc function and found the following. Any checks or advice or collaboration would be greatly appreciated.

The Sinc function

Given that:

sinc(x+Iy) = [Fr(x,y) + I Fi(x,y)]/((x^2+y^2)*π),
where
Fr(x,y) = x*sinpi(x)*cosh(πy) + y*cospi(x)*sinh(πy),
and
Fi(x,y) = x*cospi(x)*sinh(πy) - y*sinpi(x)*cosh(πy).

We have three cases to consider

• x is identically 0, in which case, we can simplify things to:

sinc(x+Iy) = sinh(πy)/(πy).

• x is an integer,hence also a zero of sinpi(x) in which case, we can simplify things to:

sinc(x+Iy) = [Fr(x,y) + I Fi(x,y)]/((x^2+y^2)*π),
where
Fr(x,y) =  y*cospi(x)*sinh(πy),
and
Fi(x,y) = x*cospi(x)*sinh(πy).

• x+0.5 is an integer, hence also a zero of cospi(x) in which case, we can simplify things to:

sinc(x+Iy) = [Fr(x,y) + I Fi(x,y)]/((x^2+y^2)*π),
where
Fr(x,y) =  x*sinpi(x)*cosh(πy),
and
Fi(x,y) = - y*sinpi(x)*cosh(πy).

With this in mind, this gives rise to the following function definition:

function sinc(x::Complex)
xr, xi = reim(x)
  if x == 0
      return one(x)
  elseif isinf(xr) && isfinite(xi) # Real part is infinite, Imaginary part is finite.
      return zero(x)
  elseif  isfinite(xr) && isinf(xi) # Real part is finite, Imaginary part is infinite.
    Si = sign(xi)
    Sr = sign(xr)
    Sc = sign(cospi(xr))
    Ss = sign(sinpi(xr))
      if xr == 0 # Real part is 0
          return complex(Inf,0)
      elseif isinteger(xr) # Real part is an integer, hence also a zero of sinpi(xr)
          return oftype(x,Inf)*complex(Sc,Si*Sr*Sc)
      elseif isinteger(xr+0.5) # Real part+0.5 is an integer, hence also a zero of cospi(xr)
          return oftype(x,Inf)*complex(Sr*Ss,-Si*Ss)
      else
          return oftype(x,Inf)*complex(Sc,-Si*Ss)
      end
  elseif isinf(xr) && isinf(xi) # Real part and Imaginary part are infinite.
      return oftype(x,NaN)*complex(1,1) # This limit does not exist. There are infinite number of paths.
  else      # Real part and Imaginary part are finite.
      return oftype(x,sinpi(x)/(pi*x))
  end
end

Like I said, any help with this issue would be gladly appreciated! Thanks.

[simonbyrne]

cc: @stevengj, @jiahao

[simonbyrne]

2x is an integer, hence also a zero of cospi(x)

This should be x+0.5 is an integer

[pjgaudre]

@simonbyrne , That makes more sense. I wanted to represent the possibility of any real number x of the form x = n/2 where n is an odd integer. So you're right: x+0.5 is an integer is a better statement.

[pjgaudre]

Perhaps, a more compact way of writing this function is as follows:

function sinc(x::Complex)
xr, xi = reim(x)
  if x == 0
    return one(x)
  elseif  isfinite(xr) && isinf(xi)     
    Si = sign(xi)
    Sr = sign(xr)
    Sc = sign(cospi(xr))
    Ss = sign(sinpi(xr))
    return xr == 0 ? complex(Inf,0) : oftype(x,Inf)*complex( (isinteger(xr+0.5) ? Sr*Ss : Sc) , Si*(isinteger(xr) ? Sr*Sc : -Ss) )
  elseif isinf(xr)  
    return isfinite(xi) ? zero(x) : oftype(x,NaN)*complex(1,1)
  else                            
    return oftype(x,sinpi(x)/(pi*x))
  end
end

[jiahao]

For a function this complicated, writing out the if... else statements is definitely preferable than the ? : syntax. They are exact synonyms after parsing so it doesn't make the generated code more efficient one way or the other.

Here are contour plots of the complex sinc function generated by Mathematica

The plots show clearly that sinc(x+0i) is also an important special case. A proper floating point implementation should also get the correct sign of zero for the imaginary component of sinc(z) along the real axis. However, it's going to be very tricky since the sign should flip for sinc(x+0.0i) and sinc(x-0.0i), but currently doesn't. (It's also broken for every elementary complex function; see #10000)

julia> sinc(-0.5+0.0001im)
0.6366197783187171 + 0.0001273239556637434im

julia> sinc(-0.5-0.0001im)
0.6366197783187171 - 0.0001273239556637434im

julia> sinc(-0.5-0.00im)
0.6366197723675814 - 0.0im

julia> sinc(-0.5+0.00im)
0.6366197723675814 - 0.0im

[simonbyrne]

I've made some improvements to sinpi in #14922.

[jiahao]

Closed by #14922

[simonbyrne]

There are still some sinc-specific issues, e.g.

julia> sinc(0.0 +250.0*im)
Inf + NaN*im

[oscardssmith]

Closed by #37273

[musm]

There are still some sinc-specific issues, e.g.

julia> sinc(0.0 +250.0*im)
Inf + NaN*im

Looks like this is still an issue.

(should be Inf + 0 im)

[stevengj]

Arguably, the remaining problem here is in complex division, is not specific to sinc:

julia> (Inf*im) / im
Inf + NaN*im