Details | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 4973 | right-hear | 1 | /* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */ |
| 2 | /* |
||
| 3 | * hypot() function for DJGPP. |
||
| 4 | * |
||
| 5 | * hypot() computes sqrt(x^2 + y^2). The problem with the obvious |
||
| 6 | * naive implementation is that it might fail for very large or |
||
| 7 | * very small arguments. For instance, for large x or y the result |
||
| 8 | * might overflow even if the value of the function should not, |
||
| 9 | * because squaring a large number might trigger an overflow. For |
||
| 10 | * very small numbers, their square might underflow and will be |
||
| 11 | * silently replaced by zero; this won't cause an exception, but might |
||
| 12 | * have an adverse effect on the accuracy of the result. |
||
| 13 | * |
||
| 14 | * This implementation tries to avoid the above pitfals, without |
||
| 15 | * inflicting too much of a performance hit. |
||
| 16 | * |
||
| 17 | */ |
||
| 18 | |||
| 19 | #include |
||
| 20 | #include |
||
| 21 | #include |
||
| 22 | |||
| 23 | /* Approximate square roots of DBL_MAX and DBL_MIN. Numbers |
||
| 24 | between these two shouldn't neither overflow nor underflow |
||
| 25 | when squared. */ |
||
| 26 | #define __SQRT_DBL_MAX 1.3e+154 |
||
| 27 | #define __SQRT_DBL_MIN 2.3e-162 |
||
| 28 | |||
| 29 | double |
||
| 30 | hypot(double x, double y) |
||
| 31 | { |
||
| 32 | double abig = fabs(x), asmall = fabs(y); |
||
| 33 | double ratio; |
||
| 34 | |||
| 35 | /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|). */ |
||
| 36 | if (abig < asmall) |
||
| 37 | { |
||
| 38 | double temp = abig; |
||
| 39 | |||
| 40 | abig = asmall; |
||
| 41 | asmall = temp; |
||
| 42 | } |
||
| 43 | |||
| 44 | /* Trivial case. */ |
||
| 45 | if (asmall == 0.) |
||
| 46 | return abig; |
||
| 47 | |||
| 48 | /* Scale the numbers as much as possible by using its ratio. |
||
| 49 | For example, if both ABIG and ASMALL are VERY small, then |
||
| 50 | X^2 + Y^2 might be VERY inaccurate due to loss of |
||
| 51 | significant digits. Dividing ASMALL by ABIG scales them |
||
| 52 | to a certain degree, so that accuracy is better. */ |
||
| 53 | |||
| 54 | if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX) |
||
| 55 | return abig * sqrt(1.0 + ratio*ratio); |
||
| 56 | else |
||
| 57 | { |
||
| 58 | /* Slower but safer algorithm due to Moler and Morrison. Never |
||
| 59 | produces any intermediate result greater than roughly the |
||
| 60 | larger of X and Y. Should converge to machine-precision |
||
| 61 | accuracy in 3 iterations. */ |
||
| 62 | |||
| 63 | double r = ratio*ratio, t, s, p = abig, q = asmall; |
||
| 64 | |||
| 65 | do { |
||
| 66 | t = 4. + r; |
||
| 67 | if (t == 4.) |
||
| 68 | break; |
||
| 69 | s = r / t; |
||
| 70 | p += 2. * s * p; |
||
| 71 | q *= s; |
||
| 72 | r = (q / p) * (q / p); |
||
| 73 | } while (1); |
||
| 74 | |||
| 75 | return p; |
||
| 76 | } |
||
| 77 | } |
||
| 78 | |||
| 79 | #ifdef TEST |
||
| 80 | |||
| 81 | #include |
||
| 82 | |||
| 83 | int |
||
| 84 | main(void) |
||
| 85 | { |
||
| 86 | printf("hypot(3, 4) =\t\t\t %25.17e\n", hypot(3., 4.)); |
||
| 87 | printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", hypot(3.e+150, 4.e+150)); |
||
| 88 | printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", hypot(3.e+306, 4.e+306)); |
||
| 89 | printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n", |
||
| 90 | hypot(3.e-320, 4.e-320)); |
||
| 91 | printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n", |
||
| 92 | hypot(0.7*DBL_MAX, 0.7*DBL_MAX)); |
||
| 93 | printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", hypot(DBL_MAX, 1.0)); |
||
| 94 | printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", hypot(1.0, DBL_MAX)); |
||
| 95 | printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", hypot(0.0, DBL_MAX)); |
||
| 96 | |||
| 97 | return 0; |
||
| 98 | } |
||
| 99 | |||
| 100 | #endif |