| OLD | NEW |
| 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_log10l.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_log10l.c */ |
| 2 /* | 2 /* |
| 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> | 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
| 4 * | 4 * |
| 5 * Permission to use, copy, modify, and distribute this software for any | 5 * Permission to use, copy, modify, and distribute this software for any |
| 6 * purpose with or without fee is hereby granted, provided that the above | 6 * purpose with or without fee is hereby granted, provided that the above |
| 7 * copyright notice and this permission notice appear in all copies. | 7 * copyright notice and this permission notice appear in all copies. |
| 8 * | 8 * |
| 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES | 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
| 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF | 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
| (...skipping 42 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 53 * | 53 * |
| 54 * ERROR MESSAGES: | 54 * ERROR MESSAGES: |
| 55 * | 55 * |
| 56 * log singularity: x = 0; returns MINLOG | 56 * log singularity: x = 0; returns MINLOG |
| 57 * log domain: x < 0; returns MINLOG | 57 * log domain: x < 0; returns MINLOG |
| 58 */ | 58 */ |
| 59 | 59 |
| 60 #include "libm.h" | 60 #include "libm.h" |
| 61 | 61 |
| 62 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 62 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 63 long double log10l(long double x) | 63 long double log10l(long double x) { |
| 64 { | 64 return log10(x); |
| 65 » return log10(x); | |
| 66 } | 65 } |
| 67 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 | 66 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
| 68 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) | 67 /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) |
| 69 * 1/sqrt(2) <= x < sqrt(2) | 68 * 1/sqrt(2) <= x < sqrt(2) |
| 70 * Theoretical peak relative error = 6.2e-22 | 69 * Theoretical peak relative error = 6.2e-22 |
| 71 */ | 70 */ |
| 72 static const long double P[] = { | 71 static const long double P[] = { |
| 73 4.9962495940332550844739E-1L, | 72 4.9962495940332550844739E-1L, 1.0767376367209449010438E1L, |
| 74 1.0767376367209449010438E1L, | 73 7.7671073698359539859595E1L, 2.5620629828144409632571E2L, |
| 75 7.7671073698359539859595E1L, | 74 4.2401812743503691187826E2L, 3.4258224542413922935104E2L, |
| 76 2.5620629828144409632571E2L, | 75 1.0747524399916215149070E2L, |
| 77 4.2401812743503691187826E2L, | |
| 78 3.4258224542413922935104E2L, | |
| 79 1.0747524399916215149070E2L, | |
| 80 }; | 76 }; |
| 81 static const long double Q[] = { | 77 static const long double Q[] = { |
| 82 /* 1.0000000000000000000000E0,*/ | 78 /* 1.0000000000000000000000E0,*/ |
| 83 2.3479774160285863271658E1L, | 79 2.3479774160285863271658E1L, 1.9444210022760132894510E2L, |
| 84 1.9444210022760132894510E2L, | 80 7.7952888181207260646090E2L, 1.6911722418503949084863E3L, |
| 85 7.7952888181207260646090E2L, | 81 2.0307734695595183428202E3L, 1.2695660352705325274404E3L, |
| 86 1.6911722418503949084863E3L, | 82 3.2242573199748645407652E2L, |
| 87 2.0307734695595183428202E3L, | |
| 88 1.2695660352705325274404E3L, | |
| 89 3.2242573199748645407652E2L, | |
| 90 }; | 83 }; |
| 91 | 84 |
| 92 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), | 85 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), |
| 93 * where z = 2(x-1)/(x+1) | 86 * where z = 2(x-1)/(x+1) |
| 94 * 1/sqrt(2) <= x < sqrt(2) | 87 * 1/sqrt(2) <= x < sqrt(2) |
| 95 * Theoretical peak relative error = 6.16e-22 | 88 * Theoretical peak relative error = 6.16e-22 |
| 96 */ | 89 */ |
| 97 static const long double R[4] = { | 90 static const long double R[4] = { |
| 98 1.9757429581415468984296E-3L, | 91 1.9757429581415468984296E-3L, -7.1990767473014147232598E-1L, |
| 99 -7.1990767473014147232598E-1L, | 92 1.0777257190312272158094E1L, -3.5717684488096787370998E1L, |
| 100 1.0777257190312272158094E1L, | |
| 101 -3.5717684488096787370998E1L, | |
| 102 }; | 93 }; |
| 103 static const long double S[4] = { | 94 static const long double S[4] = { |
| 104 /* 1.00000000000000000000E0L,*/ | 95 /* 1.00000000000000000000E0L,*/ |
| 105 -2.6201045551331104417768E1L, | 96 -2.6201045551331104417768E1L, 1.9361891836232102174846E2L, |
| 106 1.9361891836232102174846E2L, | 97 -4.2861221385716144629696E2L, |
| 107 -4.2861221385716144629696E2L, | |
| 108 }; | 98 }; |
| 109 /* log10(2) */ | 99 /* log10(2) */ |
| 110 #define L102A 0.3125L | 100 #define L102A 0.3125L |
| 111 #define L102B -1.1470004336018804786261e-2L | 101 #define L102B -1.1470004336018804786261e-2L |
| 112 /* log10(e) */ | 102 /* log10(e) */ |
| 113 #define L10EA 0.5L | 103 #define L10EA 0.5L |
| 114 #define L10EB -6.5705518096748172348871e-2L | 104 #define L10EB -6.5705518096748172348871e-2L |
| 115 | 105 |
| 116 #define SQRTH 0.70710678118654752440L | 106 #define SQRTH 0.70710678118654752440L |
| 117 | 107 |
| 118 long double log10l(long double x) | 108 long double log10l(long double x) { |
| 119 { | 109 long double y, z; |
| 120 » long double y, z; | 110 int e; |
| 121 » int e; | |
| 122 | 111 |
| 123 » if (isnan(x)) | 112 if (isnan(x)) |
| 124 » » return x; | 113 return x; |
| 125 » if(x <= 0.0) { | 114 if (x <= 0.0) { |
| 126 » » if(x == 0.0) | 115 if (x == 0.0) |
| 127 » » » return -1.0 / (x*x); | 116 return -1.0 / (x * x); |
| 128 » » return (x - x) / 0.0; | 117 return (x - x) / 0.0; |
| 129 » } | 118 } |
| 130 » if (x == INFINITY) | 119 if (x == INFINITY) |
| 131 » » return INFINITY; | 120 return INFINITY; |
| 132 » /* separate mantissa from exponent */ | 121 /* separate mantissa from exponent */ |
| 133 » /* Note, frexp is used so that denormal numbers | 122 /* Note, frexp is used so that denormal numbers |
| 134 » * will be handled properly. | 123 * will be handled properly. |
| 135 » */ | 124 */ |
| 136 » x = frexpl(x, &e); | 125 x = frexpl(x, &e); |
| 137 | 126 |
| 138 » /* logarithm using log(x) = z + z**3 P(z)/Q(z), | 127 /* logarithm using log(x) = z + z**3 P(z)/Q(z), |
| 139 » * where z = 2(x-1)/x+1) | 128 * where z = 2(x-1)/x+1) |
| 140 » */ | 129 */ |
| 141 » if (e > 2 || e < -2) { | 130 if (e > 2 || e < -2) { |
| 142 » » if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ | 131 if (x < SQRTH) { /* 2(2x-1)/(2x+1) */ |
| 143 » » » e -= 1; | 132 e -= 1; |
| 144 » » » z = x - 0.5; | 133 z = x - 0.5; |
| 145 » » » y = 0.5 * z + 0.5; | 134 y = 0.5 * z + 0.5; |
| 146 » » } else { /* 2 (x-1)/(x+1) */ | 135 } else { /* 2 (x-1)/(x+1) */ |
| 147 » » » z = x - 0.5; | 136 z = x - 0.5; |
| 148 » » » z -= 0.5; | 137 z -= 0.5; |
| 149 » » » y = 0.5 * x + 0.5; | 138 y = 0.5 * x + 0.5; |
| 150 » » } | 139 } |
| 151 » » x = z / y; | 140 x = z / y; |
| 152 » » z = x*x; | 141 z = x * x; |
| 153 » » y = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); | 142 y = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3)); |
| 154 » » goto done; | 143 goto done; |
| 155 » } | 144 } |
| 156 | 145 |
| 157 » /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ | 146 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ |
| 158 » if (x < SQRTH) { | 147 if (x < SQRTH) { |
| 159 » » e -= 1; | 148 e -= 1; |
| 160 » » x = 2.0*x - 1.0; | 149 x = 2.0 * x - 1.0; |
| 161 » } else { | 150 } else { |
| 162 » » x = x - 1.0; | 151 x = x - 1.0; |
| 163 » } | 152 } |
| 164 » z = x*x; | 153 z = x * x; |
| 165 » y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7)); | 154 y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7)); |
| 166 » y = y - 0.5*z; | 155 y = y - 0.5 * z; |
| 167 | 156 |
| 168 done: | 157 done: |
| 169 » /* Multiply log of fraction by log10(e) | 158 /* Multiply log of fraction by log10(e) |
| 170 » * and base 2 exponent by log10(2). | 159 * and base 2 exponent by log10(2). |
| 171 » * | 160 * |
| 172 » * ***CAUTION*** | 161 * ***CAUTION*** |
| 173 » * | 162 * |
| 174 » * This sequence of operations is critical and it may | 163 * This sequence of operations is critical and it may |
| 175 » * be horribly defeated by some compiler optimizers. | 164 * be horribly defeated by some compiler optimizers. |
| 176 » */ | 165 */ |
| 177 » z = y * (L10EB); | 166 z = y * (L10EB); |
| 178 » z += x * (L10EB); | 167 z += x * (L10EB); |
| 179 » z += e * (L102B); | 168 z += e * (L102B); |
| 180 » z += y * (L10EA); | 169 z += y * (L10EA); |
| 181 » z += x * (L10EA); | 170 z += x * (L10EA); |
| 182 » z += e * (L102A); | 171 z += e * (L102A); |
| 183 » return z; | 172 return z; |
| 184 } | 173 } |
| 185 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 | 174 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
| 186 // TODO: broken implementation to make things compile | 175 // TODO: broken implementation to make things compile |
| 187 long double log10l(long double x) | 176 long double log10l(long double x) { |
| 188 { | 177 return log10(x); |
| 189 » return log10(x); | |
| 190 } | 178 } |
| 191 #endif | 179 #endif |
| OLD | NEW |