Calculate Logarithms with PICAXE math

westaust55

Moderator
There have been a number of requests on the PICAXE forum for a routine to calculate logarithms using PICAXE math functions.

While there is the MicroMega uM-FPU chips which come at an extra cost but can perform these functions relatively quickly and easily
if one is only after a log function, then a PICAXE based routine may suffice.

The below code will accept an integer number from 1 to 65535 and calculate the Log to base 10.

For this code, I have gone back to first principles and used the information from the book:
The Art of Computer Programming Volume 1 - by Donald E. Knuth

The concept is extremely easy and the major hurdle is handling very large numbers often far greater than 65536.
In fact numbers up to around 100 Million are involved in an attempt to maintain accuracy.
That said, the accuracy I have observed by entering quite a few numbers and comparison with MS Excel
indicates that results are usually within ~0.0002, but occasionally may be out by ~0.002.


Code:
; Calculate logarithms (base 10) for an integer value from 1 to 65535
;
; The routine is based on the text in Chapter 1.2.2 in the book
; The Art of Computer  Programming by Donald E. Knuth   
;
; Log10(x) = n + b1/2 + b2/4 + b3/8 + b4/16 . . . + b16/65536
;
; Picaxe routine created from the above formula by
; Westaust55 (Perth, Western Australia)
;
; Accuracy by observation comaped with MS Excel  is that results
; are within 0.0002 but occasionally may be out by ~0.002
;
; Date: 28 Oct 2011
; Rev: 0
;
;
#PICAXE 40X1
#NO_DATA
#NO_TABLE

; define variable alias names
SYMBOL Reserved = b0 ; reserved for bit values
SYMBOL  LogMbit = bit0

SYMBOL       x = w4  ; w4 = b9:b8
SYMBOL MantSrc = w5  ; w5 = b11:b10
SYMBOL HiMword = w6  ; w6 = b13:b12
SYMBOL LoMWord = w7  ; w7 = b15:b14
SYMBOL LogMant = w8  ; w8 = b17:b16
SYMBOL Part1a  = w9  ; w9 = b19:b18 
SYMBOL Part1b  = w10 ; w10 = b21:b20
SYMBOL Part2   = w11 ; w11 = b23:b22

SYMBOL Charact = b26
SYMBOL Counter = b27  

; start of the program example
Main:
; Entry of a value
DataEntry:
	SERTXD  ( "Input an value:")
	SERRXD   [65000], #x		; early 40X1 parts Must have the timeout - PE tells you so!
	SERTXD (cr,lf,"Entered Value: ", #x, cr,lf)


; Calculate the Integer Part (aka Characteristic) of the logarithm
MantSrc = x 				; original value forms the source for the mantissa	

Charact = 0					; clear the Characteristic/Integer part in case second call
	DO WHILE MantSrc >= 10		; for each decade
		INC Charact   		; increment the Characterstic
		MantSrc = MantSrc / 10	; reduce the input value by one decade
	LOOP
	
; Calculate Fractional part (aka Mantissa) of the logarithm	
; here we scale the value by 1000 to have an integer value of reasonable size to work with
	IF x > 9999 THEN
		MantSrc = x / 10
	ELSEIF x > 999 THEN
		MantSrc = x * 1
	ELSEIF x >99 THEN
		MantSrc = x * 10
	ELSEIF x > 9 THEN
		MantSrc = x * 100
	ELSE
		MantSrc = x * 1000
	ENDIF

	LogMant = 0							; clear the Log10 mantissa in case second routine call

; determine "bx" for each of b1 to b16	
	FOR Counter = 1 to 16					; for each bit for a full word value
		HiMword = MantSrc ** MantSrc
		Part1a  = HiMword * 6				; divide 1000 and multiply by 65536 over next 3 lines
		Part1b  = HiMword * 5
		Part2 = HiMword * 53 /100 			; in reality this is mult by 65.536 to 2 decimal places
		LoMWord = HiMword * 6 / 1000
		
;		HiMword = 
		
		LoMWord = MantSrc *  MantSrc / 1000		; calc low word and divide low word by 1000

; determine each "bx" bit for the log mantissa		
		LogMbit = 0						; clear the "bx" bit by default - then set below as required
		IF Part1A > 999 THEN				; overflow would occur so in effect divide by 10 here
			MantSrc = Part1b + Part2 + LoMWord / 10 + Part1a + 2
			LogMbit = 1	
		ELSE
			MantSrc = Part1a * 10 + Part1b + Part2 + LoMWord + 2 ; calc new Mantissa value - add 2 for integer rounding error
			IF MantSrc >= 10000 THEN		; if  over 10,000 then reduce by a decade and set log "bx" bit
				MantSrc = MantSrc/10
				LogMbit = 1
			ENDIF
		ENDIF

		LogMant = LogMant * 2 + LogMbit		; shift left one bit and include the latest "bx" bit
	NEXT Counter						; contine until 16 bits are calculated

; For the mantissa we now have a value where bit15 = 0.50, bit14 = 0.250, bit13 = 0.125, etc
; Move the value 4 decimal places to the left to create an integer value that we can understand
; so for example, 0.3010 appears as %0101100101110111 (=dec 119) but after "shift" will become 3010 
	LogMant = LogMant ** 10000				;
		
; now print out the result for Log10(x) for the entered value for x	
	SERTXD ("Log(",#x,") = ",#Charact,".")
	IF LogMant <10 THEN
		SERTXD ("000")   
	ELSEIF LogMant <100 THEN
		SERTXD ("00")   
	ELSEIF LogMant <1000 THEN
		SERTXD ("0")   
	ENDIF
	SERTXD (#LogMant, cr,lf,cr,lf)		

	GOTO Main
Note that from Log10, values to base 2 or natural logs, if required, can be calculated using a fixed scaling factor.

The number of BASIC variables in the above routine can be reduced by carefully reworking the code, but presented here, so folks may better understand how it works.
 
Last edited:

elektriker

New Member
Hello westaust55.
you wrote:

The routine is based on the text in Chapter 1.2.2 in the book
; The Art of Computer Programming by Donald E. Knuth
;
; Log10(x) = n + b1/2 + b2/4 + b3/8 + b4/16 . . . + b16/65536

can you tell me what is in this formula n , b1 ....b16 ?

Thank you
 

westaust55

Moderator
Looking through the code I posted:
n = characteristic/integer part of the answer.
b1 through b16 are the bit values (1 or 0) of the 16-bit word being the fractional part (= mantissa)

for &#8220;n&#8221;, this is simply incremented by 1 for each decade, so for the input value in the range 1 to 9.99 = = > n = 0, 10 to 99.99 = = > n = 1, 100 to 999.9 = = > n = 2, etc

for the mantissa part I do some pre-scaling to get a larger value for working in the PICAXE integer only basic before starting to calculate the log of the mantissa.
Due to the 16 bit limitation of the PICAXE BASIC I have used several word variables as a multi-bit value to maximise accuracy.

In summary, none of these values, &#8220;n&#8221;, &#8220;b1&#8221;, &#8220;b2&#8221;, &#8230; &#8220;b16&#8221; are fixed values but are calculated for each entry or use of the routine.
 
Top