|
9 | 9 | * ====================================================
|
10 | 10 | */
|
11 | 11 |
|
12 |
| -/// pow(x,y) return x**y |
13 |
| -/// |
14 |
| -/// n |
15 |
| -/// Method: Let x = 2 * (1+f) |
16 |
| -/// 1. Compute and return log2(x) in two pieces: |
17 |
| -/// log2(x) = w1 + w2, |
18 |
| -/// where w1 has 53-24 = 29 bit trailing zeros. |
19 |
| -/// 2. Perform y*log2(x) = n+y' by simulating muti-precision |
20 |
| -/// arithmetic, where |y'|<=0.5. |
21 |
| -/// 3. Return x**y = 2**n*exp(y'*log2) |
22 |
| -/// |
23 |
| -/// Special cases: |
24 |
| -/// 1. (anything) ** 0 is 1 |
25 |
| -/// 2. 1 ** (anything) is 1 |
26 |
| -/// 3. (anything except 1) ** NAN is NAN |
27 |
| -/// 4. NAN ** (anything except 0) is NAN |
28 |
| -/// 5. +-(|x| > 1) ** +INF is +INF |
29 |
| -/// 6. +-(|x| > 1) ** -INF is +0 |
30 |
| -/// 7. +-(|x| < 1) ** +INF is +0 |
31 |
| -/// 8. +-(|x| < 1) ** -INF is +INF |
32 |
| -/// 9. -1 ** +-INF is 1 |
33 |
| -/// 10. +0 ** (+anything except 0, NAN) is +0 |
34 |
| -/// 11. -0 ** (+anything except 0, NAN, odd integer) is +0 |
35 |
| -/// 12. +0 ** (-anything except 0, NAN) is +INF, raise divbyzero |
36 |
| -/// 13. -0 ** (-anything except 0, NAN, odd integer) is +INF, raise divbyzero |
37 |
| -/// 14. -0 ** (+odd integer) is -0 |
38 |
| -/// 15. -0 ** (-odd integer) is -INF, raise divbyzero |
39 |
| -/// 16. +INF ** (+anything except 0,NAN) is +INF |
40 |
| -/// 17. +INF ** (-anything except 0,NAN) is +0 |
41 |
| -/// 18. -INF ** (+odd integer) is -INF |
42 |
| -/// 19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer) |
43 |
| -/// 20. (anything) ** 1 is (anything) |
44 |
| -/// 21. (anything) ** -1 is 1/(anything) |
45 |
| -/// 22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) |
46 |
| -/// 23. (-anything except 0 and inf) ** (non-integer) is NAN |
47 |
| -/// |
48 |
| -/// Accuracy: |
49 |
| -/// pow(x,y) returns x**y nearly rounded. In particular |
50 |
| -/// pow(integer,integer) |
51 |
| -/// always returns the correct integer provided it is |
52 |
| -/// representable. |
53 |
| -/// |
54 |
| -/// Constants : |
55 |
| -/// The hexadecimal values are the intended ones for the following |
56 |
| -/// constants. The decimal values may be used, provided that the |
57 |
| -/// compiler will convert from decimal to binary accurately enough |
58 |
| -/// to produce the hexadecimal values shown. |
59 |
| -/// |
60 |
| -
|
| 12 | +// pow(x,y) return x**y |
| 13 | +// |
| 14 | +// n |
| 15 | +// Method: Let x = 2 * (1+f) |
| 16 | +// 1. Compute and return log2(x) in two pieces: |
| 17 | +// log2(x) = w1 + w2, |
| 18 | +// where w1 has 53-24 = 29 bit trailing zeros. |
| 19 | +// 2. Perform y*log2(x) = n+y' by simulating muti-precision |
| 20 | +// arithmetic, where |y'|<=0.5. |
| 21 | +// 3. Return x**y = 2**n*exp(y'*log2) |
| 22 | +// |
| 23 | +// Special cases: |
| 24 | +// 1. (anything) ** 0 is 1 |
| 25 | +// 2. 1 ** (anything) is 1 |
| 26 | +// 3. (anything except 1) ** NAN is NAN |
| 27 | +// 4. NAN ** (anything except 0) is NAN |
| 28 | +// 5. +-(|x| > 1) ** +INF is +INF |
| 29 | +// 6. +-(|x| > 1) ** -INF is +0 |
| 30 | +// 7. +-(|x| < 1) ** +INF is +0 |
| 31 | +// 8. +-(|x| < 1) ** -INF is +INF |
| 32 | +// 9. -1 ** +-INF is 1 |
| 33 | +// 10. +0 ** (+anything except 0, NAN) is +0 |
| 34 | +// 11. -0 ** (+anything except 0, NAN, odd integer) is +0 |
| 35 | +// 12. +0 ** (-anything except 0, NAN) is +INF, raise divbyzero |
| 36 | +// 13. -0 ** (-anything except 0, NAN, odd integer) is +INF, raise divbyzero |
| 37 | +// 14. -0 ** (+odd integer) is -0 |
| 38 | +// 15. -0 ** (-odd integer) is -INF, raise divbyzero |
| 39 | +// 16. +INF ** (+anything except 0,NAN) is +INF |
| 40 | +// 17. +INF ** (-anything except 0,NAN) is +0 |
| 41 | +// 18. -INF ** (+odd integer) is -INF |
| 42 | +// 19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer) |
| 43 | +// 20. (anything) ** 1 is (anything) |
| 44 | +// 21. (anything) ** -1 is 1/(anything) |
| 45 | +// 22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) |
| 46 | +// 23. (-anything except 0 and inf) ** (non-integer) is NAN |
| 47 | +// |
| 48 | +// Accuracy: |
| 49 | +// pow(x,y) returns x**y nearly rounded. In particular |
| 50 | +// pow(integer,integer) |
| 51 | +// always returns the correct integer provided it is |
| 52 | +// representable. |
| 53 | +// |
| 54 | +// Constants : |
| 55 | +// The hexadecimal values are the intended ones for the following |
| 56 | +// constants. The decimal values may be used, provided that the |
| 57 | +// compiler will convert from decimal to binary accurately enough |
| 58 | +// to produce the hexadecimal values shown. |
| 59 | +// |
61 | 60 | use super::{fabs, get_high_word, scalbn, sqrt, with_set_high_word, with_set_low_word};
|
62 | 61 |
|
63 | 62 | const BP: [f64; 2] = [1.0, 1.5];
|
|
0 commit comments