|
8 | 8 | * |
9 | 9 | * |
10 | 10 | * IDENTIFICATION |
11 | | - * $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.108 2010/02/26 02:01:08 momjian Exp $ |
| 11 | + * $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.109 2010/08/03 21:21:03 tgl Exp $ |
12 | 12 | * |
13 | 13 | *------------------------------------------------------------------------- |
14 | 14 | */ |
@@ -5410,3 +5410,63 @@ plist_same(int npts, Point *p1, Point *p2) |
5410 | 5410 |
|
5411 | 5411 | return FALSE; |
5412 | 5412 | } |
| 5413 | + |
| 5414 | + |
| 5415 | +/*------------------------------------------------------------------------- |
| 5416 | + * Determine the hypotenuse. |
| 5417 | + * |
| 5418 | + * If required, x and y are swapped to make x the larger number. The |
| 5419 | + * traditional formula of x^2+y^2 is rearranged to factor x outside the |
| 5420 | + * sqrt. This allows computation of the hypotenuse for significantly |
| 5421 | + * larger values, and with a higher precision than when using the naive |
| 5422 | + * formula. In particular, this cannot overflow unless the final result |
| 5423 | + * would be out-of-range. |
| 5424 | + * |
| 5425 | + * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) ) |
| 5426 | + * = x * sqrt( 1 + y^2/x^2 ) |
| 5427 | + * = x * sqrt( 1 + y/x * y/x ) |
| 5428 | + * |
| 5429 | + * It is expected that this routine will eventually be replaced with the |
| 5430 | + * C99 hypot() function. |
| 5431 | + * |
| 5432 | + * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the |
| 5433 | + * case of hypot(inf,nan) results in INF, and not NAN. |
| 5434 | + *----------------------------------------------------------------------- |
| 5435 | + */ |
| 5436 | +double |
| 5437 | +pg_hypot(doublex,doubley) |
| 5438 | +{ |
| 5439 | +doubleyx; |
| 5440 | + |
| 5441 | +/* Handle INF and NaN properly */ |
| 5442 | +if (isinf(x)||isinf(y)) |
| 5443 | +returnget_float8_infinity(); |
| 5444 | + |
| 5445 | +if (isnan(x)||isnan(y)) |
| 5446 | +returnget_float8_nan(); |
| 5447 | + |
| 5448 | +/* Else, drop any minus signs */ |
| 5449 | +x=fabs(x); |
| 5450 | +y=fabs(y); |
| 5451 | + |
| 5452 | +/* Swap x and y if needed to make x the larger one */ |
| 5453 | +if (x<y) |
| 5454 | +{ |
| 5455 | +doubletemp=x; |
| 5456 | + |
| 5457 | +x=y; |
| 5458 | +y=temp; |
| 5459 | +} |
| 5460 | + |
| 5461 | +/* |
| 5462 | + * If y is zero, the hypotenuse is x. This test saves a few cycles in |
| 5463 | + * such cases, but more importantly it also protects against |
| 5464 | + * divide-by-zero errors, since now x >= y. |
| 5465 | + */ |
| 5466 | +if (y==0.0) |
| 5467 | +returnx; |
| 5468 | + |
| 5469 | +/* Determine the hypotenuse */ |
| 5470 | +yx=y /x; |
| 5471 | +returnx*sqrt(1.0+ (yx*yx)); |
| 5472 | +} |