I am trying to solve a problem posed in this question which asks that a program be written that calculates logarithms of numbers without including <math.h>.
I wrote the following program:
#include <stdio.h>
#define ln2 0.693147180559945
double power (double argu, int expo)
{
double power = 1;
for (int i = 0; i < expo; i++) { power *= argu; }
return power;
}
double log_2 (double argu)
{
int i = 0;
for (; (double)(1 << i) <= argu; i++) {}
i -= 1;
double who = (double)i, fra = 0;
for (int j = 1, k = 1; k < 450; j *= -1, k += 1)
{
fra += ((double)j / (double)k) * (power (((argu - (double)(1 << i)) / (double)(1 << i)), (double)k));
}
return who + (fra / ln2);
}
double log_b (unsigned int base, double argu)
{
return ((log_2(argu)) / (log_2(base)));
}
int main (void)
{
for (;;)
{
unsigned int base;
double argu;
printf ("Give me a nonzero positive integer base, b = ");
if (scanf ("%u", &base) != 1) { printf ("Error!!"); break; }
printf ("Give me a nonzero positive rational argument, x = ");
if (scanf ("%lf", &argu) != 1) { printf ("Error!!"); break; }
printf ("log_b(x) = %.12lf\n\n", log_b (base, argu));
}
return 1;
}
I used the following method; approximated the log base 2 using left shift and Taylor expanded around that approximation, see function log_2(); then changed the base from 2 by applying the fact that ratios of logarithms are base independent, see function log_b().
My programme works for other numbers but give it an argument larger than a billion, it hangs. I suspect, as 2^32 roughly equals 109, the computer is treating the variable double argu declared in the main() as an int instead of a double. Is my suspicion correct? If yes, then why is a double being treated as an int? On the other hand, if not, then why is the program not accepting arguments larger than a billion, a double should be able to accommodate numbers smaller than or equal to 253 right?
Your interpretation for the expression (double)(1 << i) <= argu is incorrect: 1 << i is evaluated using integer arithmetics, because << requires integer operands and it uses int because the constant 1 has int type. Hence if the result does not fit in an int or if the shift amount is greater or equal to the width of type int, 32 on most current architectures, the behavior is undefined.
You can use a larger type to accommodate up to 263: 1ULL can be shifted left at least 63 positions.
Yet your approach is still flawed for numbers less or equal to 1.0 where you try and compute 1 << -1 and for numbers larger or equal to 264 that will still cause an infinite loop.
Instead of mixing integer and floating point arithmetics, you can just use double arithmetics:
double log_2(double argu) {
double who = 0;
while (argu >= 2) {
who += 1;
argu /= 2;
}
double x = argu - 1; // x > -1.0 && x < 1.0
// Compute log(1+x) using the Taylor series first discovered by Nicolaus Mercator
double frac = 0;
for (int j = 1, k = 1; k < 450; j = -j, k++) {
frac += (double)j / (double)k * power(x, k);
}
return who + frac / ln2;
}
Note that this function can be simplified, computing the power of -x incrementally, thereby removing the need for function power():
double log_2(double argu) {
double who = 0;
while (argu >= 2) {
who += 1;
argu /= 2;
}
double x = argu - 1; // x > -1.0 && x < 1.0
// Compute log(1+x) using the Taylor series first discovered by Nicolaus Mercator
double frac = 0;
double xk = x;
for (int k = 1; k < 450; k++) {
frac += xk / k;
xk *= -x;
}
return who + frac / ln2;
}
Note however that the loop can be exited after much less than 450 iterations in most cases, but if the argument is too close to the base value, it does not converge fast enough for 450 to suffice. Try 2 and 1.999 and you get 1.000302003013 which is obviously wrong!
To mitigate this problem, we can adjust the value of argu so x is in the range [-1/3 : 1/3]: this ensures quicker convergence of the loop in less than 40 iterations (less than 33 in my tests). We can also add a test to break from the loop if frac does no changes.
Here is a faster and more reliable version:
double log_2(double argu) {
double who = 0;
while (argu < 0.667) {
who -= 1;
argu *= 2;
}
while (argu > 1.334) {
who += 1;
argu /= 2;
}
double x = argu - 1; // x >= -0.33 && x <= 0.33
// Compute log(1+x) using the Taylor series
// first discovered by Nicolaus Mercator
double frac = 0;
double xk = x;
for (int k = 1; k < 40; k++) {
double frac0 = frac;
frac += xk / k;
xk *= x;
if (frac == frac0)
break;
}
return who + frac / ln2;
}