Why not a hypot(x,y) in Math.chpl?

Hello. I have just found that (apparently) there is no hypot(x,y) function in Chapel's Math module. Given that there seems to be some reason for its existence in C, namely that hypot(x,y) works better/is more robust/accurate/ than sqrt(x*x + y*y), I wonder why it was not included in Math.

Not too important at all for my work: I am just curious.

Hi Nelson,

I'm not aware of any particularly deep reason. Note that it does seem
there is a version in the LinearAlgebraJama package module. If it was
something that would be useful to you to have, we can look into adding it.

I do see some past conversations where Damian has mentioned the C
library's behavior w.r.t. NaN handling, so I suspect there will be some
discussion along those lines.

Thanks,
Lydia

Hi Lydia,

Thanks for the quick reply! As I said, not terribly important for me to have it. What happened is that I was translating some code from C, and got curious if Chapel had hypot as well.
Regards,

Nelson

I wouldn't be surprised if this was just an oversight given that math.h defines it in C (my version at least). At least, I don't recall any discussion and intentional decision to omit it. I think it'd be very reasonable to open an issue requesting it. I will be curious to hear @damianmoz 's thoughts here.

Nelson, you are probably well-aware of this, but for others who might be reading, you could access it in the meantime using Chapel's interoperability features:

extern proc hypot(x: real, y: real): real;
extern proc hypotf(x: real(32), y: real(32)): real(32);

and/or, for the 32-bit case, you could define:

extern "hypotf" proc hypot(x: real(32), real(32)): real(32);

if you prefer to use Chapel overloading to resolve which to call. Here's an example demonstrating the calls: Run – Attempt This Online

-Brad

Am I correct that this syntax is not yet used by Chapel itself in Math.chpl?

It would greatly simpllfy the existing code.

1 Like

If all you want to do is use the system mathematics library's definitions, use what Brad has suggested

extern proc hypot(x : real(64), y : real(64)) : real(64);
extern "hypotf" proc hypot(x : real(32), y : real(32)) : real(32);

I do not use these routines because they fail to propagate a NaN in the presence of an infinity, i.e. if one argument is a Damian screwup, and the other is an infinity, it returns the infinity. I prefer that it tells me that I screwed up. My 2c.

What is your application? What data do you have (or more correctly are you certain of)?

I am assuming that eventually these definitions will find their way into Math.chpl

I believe that's correct and think it's simply because most of the current code predated the syntax being added.

-Brad

Hi Damian and Brad,

Thanks for the enlightening comments. As I mentioned in my first post, this is really just curiosity on my part: since Math.chpl mostly mirrors <math.h>, I wondered why not hypot. This all started because I wrote a simple Barne's analysis procedure, which calculates the value of a field z(x,y) at a point from surrounding known zb(xb_i,yb_i), weighing their values essentially by exp(-dist((x,y),(xb_i,yb_i)**2) where dist is the euclidean norm. For now I haven't even tested it, and will just write sqrt( (x-xb[i])**2 + (y-yb[i])**2). I do not expect hypot to make any big difference here.

Cheers

Nelson

I like this new syntax. If faciltitates negative programming.

Thanks for that info Brad.

The main advantage of hypot(x, y) is that is it most implementations will be accurate to one (1) unit in the last place. For your problem where you are computing a weight, your approach looks reasonable. If you are looking at a Givens transformation, that's a whole new ball game because you need robustness in the presence of totally unknown data.

These days, I tend to use my own routine

(r, s, t) = sumofsq(x, y); /// r*(s+t) = x^x + y*y, s >> t

The scale factor r ensures that neither overflow or underflow occur and I get to tweak the accuracy n subsequent calculations.

I then figure out whether I use sqrt or rsqrt on the result because I often want the reciprocal in my computations of things like the norm.

My 2c.