Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use 'long double' for x,y,z arguments #93

Closed
wants to merge 1 commit into from
Closed

Conversation

dstndstn
Copy link
Collaborator

  • use cosl,sinl for radec_to_healpixlf
  • compute 1-vz in long double

Copy link
Member

@astrofrog astrofrog left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me - I don't have time to check whether this fully fixes #34, so I'll leave it up to you whether to merge or wait until we can confirm this is fixed.

@cdeil
Copy link
Member

cdeil commented Jun 26, 2018

There's a bit of a concern with this, in the sense that it probably won't fix the issue for all users: some will end up with only "double" precision still (see https://en.wikipedia.org/wiki/Long_double).

Even worse would be if there's a distribution problem where binaries (e.g. Windows from conda-forge) don't work at all for some users (e.g. some Windows machines); I don't know if that's a real concern.

I would suggest we merge tomorrow, and continue with more edits and checks in the next days.

Unless there is a more stable way to do that computation using only double.
@mreinicke - do you know if there is a trick to get this numerically stable?

@mreineck
Copy link

mreineck commented Jun 26, 2018

Hmm, as far as I understand, the problematic term is sqrt(1.-vz).

Wouldn't it be possible to rewrite this as sqrt((vx*vx+vy*vy)/(1.+vz)) ?

That should be more accurate.

@mreineck
Copy link

Going to long double is almost never a good idea, because it doesn't increase accuracy on many platforms. In most cases there is a way to rewrite the problematic expressions in a more stable fashion.

@dstndstn
Copy link
Collaborator Author

In this case it's the sin(Dec) that drops precision; Dec->90 degrees so sin(Dec)->1

@mreineck
Copy link

Yes, that's why I think it should be safer to use cos(Dec) (i.e. sqrt(x*x + y*y)) instead of sin(Dec) (i.e. z) in this situation.

I see similar code in other places in the package... the z2dec function in starutil.inc for example will fail badly near the poles. I'd recommend something like

double xyz2dec(double x, double y, double z)
  {
  return M_PI/2 - atan2(sqrt(x*x+y*y), z);
  }

It's not much slower, and it works extremely reliably.

@dstndstn
Copy link
Collaborator Author

I tried this in #94 - thanks for the tip! - seems to work in our one polar test case

@cdeil
Copy link
Member

cdeil commented Jun 26, 2018

@mreineck @dstndstn - Thank you!
Merging #94 now and closing this one without merging.

@cdeil cdeil closed this Jun 26, 2018
@cdeil cdeil deleted the long-double branch October 24, 2018 08:41
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants