|
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