I was bored on a snowy March day, so I tried typing this up, not sure if anyone else will find this interesting.
I am doing this in Pascal, well for clarity and for the mere fact that I like Pascal.
The first step is to calculate the square root of 10, then its square root and so on, down to at least the tenth square root. Also, we need to calculate the power of 2^n. I am going to place them in an array. Sure I could use the standard square root function, but somehow that seems like cheating. But there is a well known routine to calculate a square root and it is actually pretty efficient.
Function MySqrt(f: Double): Double;
Var
Root, k: Double;
Begin
Root := f / 2.0;
Repeat
If (Root<>0) Then
Root := 0.5 * (Root + (f / Root));
Until abs(f-(Root*Root)) < 1e-10;
MySqrt := Root;
End;
It is known as the Heron method and it can be read about
here:
https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method
So first of all, we need to create an array of Square Roots and Divisors, I am building both of them into one array structure.
Const
KkSize = 50;
Type
CellType = Record
SqrtVal: Double;
Divisor: Double;
End;
Var
i: LongInt;
Divisor: Double;
WorkArray: Array[1..KkSize] of CellType;
...
Value := 10.0;
Divisor := 2.0;
For i := 1 to KkSize Do
Begin
Value := MySqrt(Value);
WorkArray[i].SqrtVal := Value;
WorkArray[i].Divisor := Divisor;
Divisor := Divisor * 2.0;
End;
Now finally comes calculating the actual log function. First, we set LogNbr to zero and DivVal to the number we want to convert to its log. Now one has to change the n value so that it is greater than zero and less than 10. The Significant will be the value of LogNbr that is to the left of the decimal point. These are the instructions that do that:
If DivVal > 10 Then
Repeat
LogNbr := LogNbr + 1;
DivVal := DivVal / 10;
until (DivVal <= 10);
If DivVal < 1 Then
Repeat
LogNbr := LogNbr - 1;
DivVal := DivVal * 10;
until (DivVal >= 1);
Finally, we calculate the Mantissa with these instructions:
For i:=1 To KkSize Do
Begin
If DivVal >= WorkArray[i].SqrtVal Then
Begin
LogNbr := LogNbr + (1.0 / WorkArray[i].Divisor);
DivVal := DivVal / WorkArray[i].SqrtVal;
End;
End;
The whole idea is to check If DivVal is greater than the sqrt, if so then LogNbr := LogNbr + (1/Divisor).
Divisor is going to be a power of two. Then we divide the DivVal by SqrtVal.
I was able to piece this logic together from the book, “An Introduction to the History of Mathematics” by Howard Eves which was published in 1964. My parents had purchased this used book from a library and I had discovered it in the 1980s and found it interesting.
Here is the complete function:
Function MyLog(n: Double): Double;
Var
LogNbr, DivVal: Double;
i: LongInt;
Begin
LogNbr := 0.0;
DivVal := n;
If DivVal > 10 Then
Repeat
LogNbr := LogNbr + 1;
DivVal := DivVal / 10;
until (DivVal <= 10);
If DivVal < 1 Then
Repeat
LogNbr := LogNbr - 1;
DivVal := DivVal * 10;
until (DivVal >= 1);
For i:=1 To KkSize Do
Begin
If DivVal >= WorkArray[i].SqrtVal Then
Begin
LogNbr := LogNbr + (1.0 / WorkArray[i].Divisor);
DivVal := DivVal / WorkArray[i].SqrtVal;
End;
End;
MyLog := LogNbr;
End;
No comments:
Post a Comment