Thursday, March 09, 2023

How to calculate Log10 without calling the Log or Ln Function

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: