"General Purpose" Interpolation subroutines

AllyCat

Senior Member
#1
Hi,

PICaxe Basic struggles to implement the more complex mathematical calculations (such as trigonometric or logarithmic functions) because the interpreter is slow and doesn't support signed or fractional (non-integer) numbers. One solution is to use a "lookup table" of results, but normally a table cannot contain every possible value that might be required. Therefore, it is necessary to "interpolate" values between the sample-points which have been stored in the table.

The mathematical function may be considered as a graph of Y versus X coordinates, with a lookup table containing sample points along the curve. The simplest form of interpolation assumes that adjacent points are joined by a straight line, i.e. that the smooth curve is replaced by a "piecewise linear" representation. In this case, the "target" result can be obtained from the coordinates of the start and end points of the appropriate line, then calculating the gradient (slope) of the line and the proportional distance of the target point along that line.

The two subroutines to be described in this thread are part of a family of routines which have a consistent structure as described in the thread here. To minimise the size of the algorithm and the number of (RAM) variables required, a few constraints are placed on the range of the numerical values. The primary limitations here are that the table may contain a maximum of 128 sample points (normally it will be much less) which must be equally-spaced along the X-axis, with 255 intermediate interpolated values available (i.e. defined by a single byte) along each segment between the reference points. Each sample point may have a maximum value in the Y-axis of 16-bits (i.e. two bytes = 65535).

Another constraint is that the interpolating lines must always have a positive (or zero) gradient, i.e. each subsequent Y-value in the table may not be smaller than the previous point. In practice this is not a severe restriction because cyclic trigonometric functions can be considered as separate segments "mirrored" in the X or Y axes, with any negative gradients handled by simply reversing the numerical values along the X-axis. Such a restriction also ensures that the "reverse" lookup (i.e. from y-axis to X-axis) using the same data produces a unique solution (or at worst a contiguous range).

Thus, each interpolation line fits (as the diagonal) across a rectangular "cell" which has a width of 256 units and a maximum height of either 4095 or 65535 units (depending on the algorithm employed). Typically there will be 16 cells (i.e. rectangles with their NE and SW corners touching), to give a total horizontal (x-axis) range of 0 - 4095 and vertical (Y-axis) range of 0 - 65535. These integer values can of course be scaled and/or offset, to suit the numerical range and resolution of the dataset.

These requirements dictate that most cells (and the graph of the whole dataset) will be principally "tall and thin" (e.g. 65536 high by 4096 wide). So, if the dataset converts from a "large" range of numerical values to a smaller range, such as a logarithmic function, then it is necessary to perform a "reverse lookup" function, i.e.. from the Y-axis to the X-axis. This is the reason that two different subroutines are required, one to convert from a point on the X-axis to the corresponding value on the Y-axis and a second to accept a point on the Y-axis and calculate the value along the X-axis. The latter routine is slightly more complex because it needs to "search" through some or all of the (Y) entries in the table, to find the cell which contains the target datapoint.

The lookup tables are "indexed" so that up to 8 different datasets may be used, for example (arc)sine/cos, (arc)tan{2}, exponential/log and NTC resistance/temperature, etc. This is why the number of datapoints per function will be typically limited to around 16, i.e. 16 words of 2 bytes x 8 functions = 256 bytes, or one "page" in the EEPROM (table). The table index number has a value from 0 to 7, but is stored in the top 3 bits of b1 (to avoid conflicts with the flags available in the division subroutine described in the thread linked above). This number is converted by a preliminary lookup table to give a "pointer" address of the first entry in the appropriate table. Typically, the number of cells defined by each table will be a binary power (e.g. 8, 16 or 32), but any value may be used, within the constraints of available memory and scale factors.

So here is the first subroutine to calculate the value on the Y-axis which corresponds to a point on the curve/graph directly above a supplied value along the X-axis. This first example includes a Sine function table for one quadrant, divided into 16 piecewise-linear sections.

Code:
picaxe 08m2						; Or most others
#no_data

symbol TOM = 32					; Multiplier for Lookup Table number (i.e. flags * 32)
symbol WIDTH = 256				; Piecewise-Linear horizontal segment length
symbol tempb = b1					; Temporary (local) and parameter-passing byte 
symbol tempw = w1					; Temporay (loacl) word and parameter-passing word 
symbol wx 	= w2					; Interpolation X-word and Numerator High word in ddiv
symbol wy  	= w3					; Interpolation Y-word and Divisor in ddiv

symbol TABBASE = 0				; Table of table base addressess (max 8 bytes)
symbol SINBASE = 8				; Start address of Sine table (length 34 bytes)
symbol SININDEX = 0 * TOM			; Sine/Cos Lookup table pointer
data TABBASE,(SINBASE)
data SINBASE,(0,0,141,12,252,24,45,37,2,49,94,60,39,71,66,81,144,90)
data (1,99,125,106,243,112,83,118,144,122,157,125,117,127,255,127)   ; Sine table ($8000 = +1)

cosine:							; Enter with angle in w1 (2^16 = 360 degrees)
	w1 = w1 + $4000				; Add one quadrant and fall through into sine 
sine:								; Enter with angle in w1 (2^16 = 360 degrees)
	if w1 < $8000 then 				; Shift down if in 2nd or 3d quadrants (180 - 360 degs)	
posine:	
		if w1 > $3FFF then			; Reverse 2nd and 4th quadrants
			w1 = $8000 - w1
		endif
		w1 = w1 / 4				; Max Table X-Axis = 16 * 256 = 4095 (Top 4 bits spare)
		b1 = SININDEX				; Select Table and fall into interpXtoY	
interpXtoY:	 						; Retain flags in .0 to .3
		wx = w1					; Frees up b2 and b3
		b2 = b1 / TOM + tabbase	; Table address pointer in top 3 bits
		read b2,b2					; Address of start (word) of table
		b2 = b2 + b3 + b3			; Add 2* wx/256 (b3 still contains hi byte of w1) = base of cell
		b3 = b2 + 2				; Point to top of cell
		read b3,WORD wy
		call readw1
		wy = wy - w1 				; Calculate height of cell (wy - wx)
		wx = wx / 256 * 65280 + wx	; Horizontal distance across cell (65280 = -1 * 256)
		w1 = wx * 256 ** wy + w1	; Height inside cell (0<wx<256) * gradient + height of cell base
		w1 = w1 max $7FFF
	else
		w1 = w1 - $8000			; Invert befeore and after positive calculation
		call posine
		w1 = -w1		
	endif 		
return							; From interpolation or inversion routines

readw1:							; PE does not accept READ w1,WORD w1 directly (b2 is pointer)
	b3 = b2 + 1					; Point to high byte
	read b2,b2						; Read low byte
	read b3,b3						; Read high byte
return							; Result in w1
The "reverse" functions (arcsin/cos using interpolation from Y to X) are rather more complex because the table needs to be searched and then a "high resolution" division may be needed, so those subroutines will be held over for the next post.

Also, apart for the sine/cos function, a test harness is not included yet, but a very similar subroutine is included in the Sunrise and Sunset calculation in the thread here.

Cheers, Alan.
 

AllyCat

Senior Member
#2
Hi,

A recent topic in the Active Forum reminded me that I didn't "finish" this thread, so here are a few more details. Not yet the promised "reverse lookup" (X from Y) algorithm, but code for ARCTAN and an example of how the subroutine can use multiple tables (e.g. for SIN and ATAN).

Unlike SIN and COS, which have a maximum range from zero to 1.0 (or down to -1.0 if angles outside the "first quadrant" are considered), the TAN function spans all the way from zero up to "infinity" (and the negative equivalent in other quadrants) which makes it difficult to create a complete integer-based lookup table. However, angles between 0 and 45 degrees do convert to a range from 0 to 1.0 (i.e. the same range as for SIN and COS). Then, angles between 45 and 90 degrees can be calculated from the inverse representation (i.e. by dividing unity by the complementary angle). Thus, {A}TAN is best implemented by using 8 octants, each spanning 45 degrees, whilst SIN and COS need only 4 quadrants of 90 degrees each.

Normally, ATAN would be considered the "reverse" (or inverse) function, but ATAN (or more specifically ATAN2 which can receive two independent directional vectors) is generally more useful than TAN. The gradient of the lookup curve (table) doesn't change greatly over the range of 0 to 45 degrees (i.e. the X and Y axes can have comparable resolutions), so it's practical to create a table using ATAN for the forward (X to Y) conversion direction. TAN can then be calculated by the reverse lookup algorithm, or from the knowledge that TAN(x) = SIN(x) divided by COS(x).

Accuracy:


The thread linked at the top of this post paid particular attention to the accuracy of the ATAN function, so let's consider what compromises are necessary or worthwhile here. The available PICAxe table space is 256 bytes (in either EEPROM or TABLE memory), i.e. an average of 16 words per table if all 8 potential tables were implemented. However, the two tables below can handle all the "Trig" functions, so probably a total of 3 or 4 tables is as many as might be needed in practice, so a larger number of points per table can be considered.

Three basic factors determine the "accuracy" of the results: the "resolution" (i.e. numerical range) of the tables' "output" values (Y), the resolution of the input values (X) and the error due to the degree of "curvature" of the table values (i.e. how far the actual function deviates from a straight line joining each pair of reference points). Also the "scaling" of the input/output values (e.g. representing 1.0 by an integer value of 10,000 , or 36,000 to represent 360 degrees, etc.) may affect the resolution/accuracy. But here, using the full range of a Word (i.e. 0 - 65535) to represent one revolution (360 degrees) and the range -1.0 to +1.0, makes the best use of the available number system.

The table output values (Y) are stored in full Words (2 bytes), so are unlikely to limit the accuracy/resolution any more than any subsequent numerical processing. However, the chosen subroutine limits the number of interpolation points to 256 within each "line segment", so the use of 16 lines (17 lookup points) limits the "input" resolution to 4096 values. This is the limiting factor for the fairly linear {A}TAN function, so it may be worthwhile to increase the table to 32 segments (33 words, or 66 bytes). It's probably unnecessary to go any further because the whole table represents only one octant (45 degrees) so the effective output value also has a nominal resolution of 8192 points. For compatibility, the Y values have been scaled to a value of 32768 = Unity, but the "error" value at the centre of the "worst" line is no more than 2 in the Least Significant digit, so becomes negligible when scaled to a full revolution.

The SIN function has considerably more curvature so the worst interpolation error is about 5 in the LS digit, when using a 32 section lookup table. Therefore, it may be worthwhile to increase the table to 64 segments (130 bytes) to reduce the worst interpolation errors to around 1 in the LS digit. However, even the 5 digits error (for 32 lines) scales to less than 0.0002 (i.e. 5 / 32768) , or a little more if alternate data values are simply deleted from the table (because the remaining points no longer have a symmetrical best fit to the maximum errors).


So here are sample subroutines for SIN, COS and ARCTAN with a test harness and conversion factors for scaling to integer values of 10,000 (i.e. 4 decimal places) and 36,000 (i.e. 360 degrees with two decimal places). Also the octant-handling for ATAN2 is shown, but adding a test harness and division routine for this will exceed the 10,000 character forum limit.

Code:
; SIN,COS and ARCTAN functions,	AllyCat,	March 2018  
#picaxe 08m2		; And most others

symbol TOM = 32				; Multiplier for Lookup Table number (i.e. flags * 32)
symbol tempb = b1			; Temporary (local) and parameter-passing Byte 
symbol tempw = w1			; Temporary (local) and parameter-passing Word 
symbol wx = w2				; Interpolation X-word and Numerator High word in ddiv
symbol wy = w3				; Interpolation Y-word and Divisor in ddiv
symbol NEGX = 32768
symbol degs = w9
symbol tan = w10

symbol TABBASE = 0			; Table of table base addressess (max 8 bytes)
symbol SINBASE = 8			; Start address of Sine table (length 130 bytes)
symbol SINTAB = 0 * TOM			; Sine/Cos Lookup table pointer
symbol ATANBASE = 138
symbol ATANTAB = 1 * TOM

data TABBASE,(SINBASE,ATANBASE)

data SINBASE,(0,0,36,3,72,6,107,9,140,12,171,15,200,18,226,21,249,24)	; 64 segment SINE table
data (12,28,26,31,36,34,41,37,39,40,32,43,18,46,252,48)
data (223,51,187,54,142,57,88,60,24,63,207,65,124,68,30,71)
data (181,73,65,76,193,78,53,81,156,83,247,85,68,88,131,90)
data (181,92,216,94,237,96,243,98,234,100,209,102,168,104,111,106)
data (37,108,203,109,96,111,228,112,86,114,183,115,6,117,67,118)
data (109,119,134,120,139,121,126,122,94,123,43,124,229,124,140,125)
data (31,126,158,126,11,127,99,127,168,127,218,127,247,127,2,128)

data ATANBASE,(0,0,23,5,44,10,61,15,69,20,68,25,54,30,26,35,238,39)		; 32 segment ATAN table
data (176,44,94,49,247,53,122,58,230,62,56,67,114,71,146,75)
data (152,79,132,83,85,87,12,91,168,94,42,98,147,101,225,104)
data (24,108,53,111,59,114,41,117,255,119,192,122,107,125,1,128)

test:  
	for tan = 0 to 10
		w1 = tan * 1000 ** 53687	; Scale to unity = 8192 (32 line segments)
		b1 = 0				; Set the ATAN2 flags to first quadrant
		call arctan
		w1 = w1 ** 36000		; Scale to hundredths of a degree
		sertxd(cr,lf,"ArcTan ")
		if tan = 10 then
			sertxd("1.0")
		else
			sertxd("0.",#tan)
		endif
		b1 = w1 // 100	
		w1 = w1 / 100
		sertxd(" = ",#w1,".")
		call showlz
		sertxd(" degrees")
	next
	for degs = 0 to 360 step 15
		w1 = degs * 100			; full scale = 36000 (hundredths of degrees)
		w1 = w1 ** 53769 + w1	; scale up to 65536 = 360 degrees
		call sine
		sertxd(cr,lf,"Sin ",#degs," = ")
		call showdec
		w1 = degs * 100
		w1 = w1 ** 53769 + w1	
		call cosine
		sertxd(" Cos = ")
		call showdec
	next
	end
	
showlz:			; Reinsert a suppressed leading zero
	if b1 < 10 then
		sertxd("0")
	endif	
	sertxd(#b1)
return
	
showdec:			; Show w1 signed with 4 decimal digits
	if w1 => NEGX then 
		sertxd("-")
		w1 = -w1
	endif	
	w1 = w1 ** 20001			; Convert 32768 to 4 decimal digits
	b1 = w1 / 10000	
	w1 = w1 // 10000
	sertxd(#b1,".")
	b1 = w1 / 100
	call showlz
	b1 = w1 // 100
	goto showlz 	; Then use its return

atan2:			; ARCTAN2 subroutine (signed vectors in wx and wy) 
	b1 = 0 				; First octant (0 - $1FFF)	
	if wx >= NEGX then 		; Value is negative
		b1 = 7
		wx = -wx			; Convert negative to positive
	endif					; WX is (now) positive
	if wy >= NEGX then 		; Value is negative
		b1 = b1 xor 3
		wy = -wy			
	endif					; WY is (now) positive
	if wx > wy then 			; WX/WY is not a proper fraction ( <1 ) so swap
		b1 = b1 xor 1			; Mark that WX > WY (Octant 1,2 etc)
		w1 = wx			; )
		wx = wy			; }Exchange WX and WY (~5 code bytes less than SWAP)
		wy = w1			; )
	endif					; Divisor larger than numerator so "fractional" result

;** Call Double-Word Division subroutine (w1 = wx:00 / wy) , then fall into arctan **

arctan:				; w1 = 8192 represents unity (for 32 segment table)
	b1 = b1 + ATANTAB		; Select Table
	call interpXtoY
	w1 = w1 + 2 / 4			; Full scale for one octant is 8192
	if bit8 = 1 then			; LSB of b1
		w1 = $2000 - w1		; Reverse the sector from 90 to 45 degrees
	endif
	w1 = b1 & 7 * $2000 + w1		; Add the other octants
return

cosine:			; Enter with angle in w1 (2^16 = 360 degrees)
	w1 = w1 + $4000				; Add one quadrant and fall through into sine 
sine:							; Enter with angle in w1 (2^16 = 360 degrees)
	if w1 < $8000 then 			; Shift down if in 2nd or 3d quadrants (180 - 360 degs)	
possine:	
		if w1 > $3FFF then		; Reverse 2nd and 4th quadrants
			w1 = $8000 - w1
		endif
		b1 = SINTAB			; Select Table 	
		call interpXtoY
	else
		w1 = w1 - $8000			; Invert before and after positive calculation
		call possine
		w1 = -w1		
	endif 		
return

interpXtoY:	 				; Retain flags in .0 to .3
	wx = w1				; Frees up b2 and b3
	b2 = b1 / TOM + TABBASE	; Table address pointer in top 3 bits
	read b2,b2				; Address of start (word) of table
	b2 = b2 + b3 + b3		; Add 2* wx/256 (b3 still contains hi byte of w1) = base of cell
	b3 = b2 + 2				; Point to top of cell
	read b3,WORD wy
	b3 = b2 + 1				; Point to high byte
	read b2,b2				; Read low byte
	read b3,b3				; Read high byte
	wy = wy - w1 				; Calculate height of cell (wy - wx)
	wx = wx / 256 * 65280 + wx		; Horizontal distance across cell (65280 = -1 * 256)
	w1 = wx * 256 ** wy + w1		; Height inside cell (0<wx<256) * gradient + height of cell base
	w1 = w1 max $7FFF
return
Cheers, Alan.
 
Last edited:

AllyCat

Senior Member
#3
Hi,

Here is the subroutine to perform interpolation in the "reverse" Lookup direction, i.e. from Y values to X values. It can use the same table(s) as above to perform "inverse" functions such as ARCSIN and ARCCOS, but also may be needed for functions with "asymmetric" scales such as Square Root or Logarithms. However, it is not the overall scales which are significant (because they can be adjusted by simple multiplication or division) but the change in the gradient of the curve (or transfer characteristic) over the range of the table.

For example with a square-law function, when the X-value changes from 0 to 1, the Y value also increases from 0 to 1 so the average gradient is 1. However, when X changes from 10 to 11, Y changes from 100 to 121 so the gradient is 21, and when X changes from 100 to 101, Y changes from 10000 to 10201 with a gradient of 201. Typically the Y-scale can be from 0 to 65535, whilst the X-scale with 16 line segments (each having 256 sub-points) would be from zero to 4095. Hence it's sensible to create a "SQUARES" table for the forward direction (even though the basic calculation doesn't need to use a lookup table) so that the larger range of values appears along the Y-axis.

The following program demonstrates the subroutine in a simple test harness. It uses only one SINE table, but with the choice of two resolutions (number of line segments), smaller than in the previous post to speed the table search and to show how the scaling needs to change. To select the lookup data from multiple tables, a subroutine is included, but this could be a Macro, or even replaced by a fixed constant assignment, if only one table is required.

The main data values are passed to and from the subroutine in W1 and may need to be scaled before and/or afterwards, depending on the number of table segments and the User-Interface parameters (degrees, radians, etc.). Also, B1 contains flags which are combined to indicate the lookup table number and the Sign, Quadrant or Octant values of trigonometric functions, etc.. The table address is recalculated after the interpolation stage to avoid the need to store it in an additional variable.

Code:
; Interpolation subroutine for "reverse" lookup through tables (X from Y) , AllyCat , April 2018
#picaxe 08m2					; Or most others
#define fastdiv
#define sinseg32

symbol NEGX = 32768				; Negative numbers boundary
symbol TAM = 32					; Lookup Table number Multiplier
symbol tempb = b1				; Temporary (local) and parameter-passing byte 
symbol tempw = w1				; Temporary (local) and parameter passing word
symbol wx = w2					; Interpolator X-word
symbol wy = w3					; Interpolator Y-word
symbol tabbase = 0				; Table of table bases (max 8 pointer bytes)
symbol sinebase = 8				; Start address of Sine table (length 66 bytes)
symbol nextbase = sinebase + 66
symbol sinindex = TAM * 0			; Sine/Cos Lookup table pointer flags
data tabbase,(sinebase,nextbase)		; Pointers to base of tables

#ifdef sinseg32
symbol SINSC = 2								; Scale multiplier for 360 degrees = 65536
data sinebase,(0,0,72,6,141,12,201,18,250,24,27,31,41,37,33,43)			; 32 Segment Sine table
data (254,48,188,54,90,60,209,65,32,71,67,76,55,81,250,85,133,90)		; With Symmetrical errors
data (219,94,246,98,212,102,114,106,206,109,232,112,187,115,71,118)
data (138,120,130,122,47,124,143,125,162,126,103,127,222,127,255,127)
#else										; 16 line segments
symbol SINSC = 4								; Scale multiplier
data sinebase,(0,0,141,12,252,24,45,37,2,49,94,60,39,71,66,81,144,90)
data (1,99,125,106,243,112,83,118,144,122,157,125,117,127,255,127)
#endif

demo:
	for b0 = 0 to 100 step 5				; Hundredths of unity
		w1 = b0 * 100 ** 53687 * 4 max 32767		; Scale Unity = 32768 
		call arcsin
		w1 = w1 ** 36000				; Scale to hundredths of a degree
		b1 = w1 // 100					; Hundredths
		w1 = w1 / 100					; Unit degrees
		sertxd(cr,lf,#b0,"% = ",#w1,".")
		w1 = b1 / 10
		b1 = b1 // 10
		sertxd(#w1,#b1," degrees")	
	next b0
end

arcsin:
	if w1 >= NEGX then 				; Value is negative
		w1 = -w1				; Convert to a positive value
		call parcsin
		w1 = -w1				; Shift result by two quadrants
	else						; Positive calculation
parcsin: 			; Calculate ArcSin of a positive number (in w1)
		b1 = sinindex 				; Select SINE table
		call interpYtoX
		w1 = w1 * SINSC				; Scale dependent on number of segments
	endif
return

arccos:
	call arcsin
	w1 = $4000 - w1					; Cos is shifted one quadrant (90 degrees)
return

getbase:						; Load b2 with base address of the required table
	b2 = b1 / TAM + tabbase				; b1 contains flags for table number
	read b2,b2
return

interpYtoX:						; Enter with table number (*32) and flags in b1 (retained)
	wx = w1						; Free up b2 and b3 to use separately
	call getbase					; Load b2 with the base address of SINE table
	do
		b2 = b2 + 2				; Move to RHS of (next) cell
		read b2,word wy				; Read Y-value at top of cell (NE corner)
	loop until wy => wx				; Repeat until target point is within the cell  	
	dec b2						; Point to High byte of LH corner (SW) of cell 
	read b2,b3					; Read High byte of base of cell
	wx = -b3 * 256 + wx   				; Subtract high byte of base of cell
	wy = -b3 * 256 + wy 
	b3 = b2 - 1					; Pointer to low byte
	read b3,b2					; Read Low byte of base of cell
	wx = wx - b2					; Subtract low byte of base of cell
	wy = wy - b2
	call getbase					; Read start of table again
	b3 = b3 - b2 				; Calculate cell number * 2 (or divide by 2 to load directly in w1)

#ifdef fastdiv						; Fast version if cell height is restricted to 4095
	wx = wx * 16					; Multiply numerator by 16
	wy = wy / 16					; Divide divisor by 16	
#else							; Scale numerator and divisor for best resolution
	for b2 = 0 to 7				 	; Multiply numerator/divisor by 256 to give result = 0 - 255
		if wx < 32768 then 			; WX can be doubled
			wx = wx + wx			; Multiply numerator by 2
		else			
			wy = wy / 2			; Divide divisor by 2	
		endif	
	next 	
#endif				
	wx = wx / wy max 255				; Horizontal distance within cell  (max = 255)
	w1 = b3 * 128 + wx				; Add the value of the LH edge of cell
return
Details of the Subroutine's Operation:

The reverse subroutine is rather more complex than the forward one for several reasons: Firstly, the algorithm needs to search through the list of Y-values to find the relevant "cell" (because they are not of equal height), whilst the constant-width of the cells in the X-direction (256 points) made calculation of the cell number particularly easy. The search process is not difficult, but it does consume more time and also it's important to ensure that the table data is constructed such that the search cannot fail (loop endlessly) because the subroutine does not set any timeout (overflow) limits.

But the main complication is the calculation of the interpolation points within each cell. In the forward direction, the "base line" (X) was fixed at 256 points, so calculating the ratio of Y/X (i.e. the line gradient) was easy because division by 256 (a byte) in binary is simple. However, in the reverse direction, the divisor (for X/Y) is not only a variable but can be very large, potentially up to most of a full Word (two bytes):

The "double-Word" division subroutine that I have previously documented can perform the required function but it has disadvantages; It is relatively slow and uses the PICaxe b1 variable which (in this family of subroutines) already contains various table and maths flags. Of course another variable (or two) could be introduced, but I have tried to keep the "RAM footprint" of these routines as small as possible. The introduction of another variable is also avoided by reading and subtracting the table's lower word one byte at a time, but this doesn't affect the executed code very much because the WORD qualifier is basically a predefined Macro anyway.

Thus, two alternative (conditionally compiled) division routines are included. The simpler assumes that the height of each cell will never exceed a value of 4095, so it is safe to multiply by 16 without causing an overflow. At the same time the divisor is divided by 16, which can't cause an overflow, but might lead to some loss of accuracy. However, these conditions cannot occur arbitrarily or "unexpectedly" because they are determined by the data in the lookup table. So it is a valid simplification provided that there are sufficient cells and/or the scale of the Y-axis is restricted. For example, a SINE curve divided into 16 cells, with unity represented as 32768 in the Y-direction, should stay just within bounds (the "tallest" cell is just less than 2048 * PI / 2).

The more "intelligent" method also employs an intrinsic multiplication by 256, but this is done by successively doubling the numerator when it can't cause a word overflow, or otherwise by halving the divisor. Then, the division can be performed by the PICaxe Basic interpreter quickly, with reasonable accuracy. This method is likely to be needed with exponential or "power" function tables, or if you're not confident about analysing the structure of the tables. ;)


Still to come: A more complete test harness to show all the Trigonometric expressions.

Cheers, Alan.
 
Last edited:

AllyCat

Senior Member
#4
Hi,

Finally, a demonstration and test harness for all the Trigonometric functions above, in 5 degree increments. It can be run in the simulator but might take half an hour or more. At 900+ bytes, it's not intended as a formal Code Snippet, but shows how the modules and tables' scale/resolution work together.

Two methods for TAN are shown, one dividing Sin / Cos (which needs one additional storage word), the other using the ATAN table in reverse. The tables have only 16 segments because the program is very close to the forum's 10,000 character posting limit.

Code:
; General Purpose Interpolation routines with lookup tables for Trig Functions;
; Angles in Trig functions use One revolution = 65536 (16 bits)
; AllyCat , July 2017, updated April 2018

#picaxe 20x2        		; Or any M2, except to compare embedded X2 SIN function
#define fastdiv			; Delete if InterpYtoX has steep gradients

symbol NEGX = 32768		; Boundary of Negative numbers (Twos complement)
symbol TAM = 32			; Lookup Table number Multiplier
symbol tempb = b1		; Temporary (local) and Parameter-passing byte 
symbol tempw = w1		; Temporary and Parameter-passing word
symbol wx 	= w2		; Interpolator X-word and Division Extension word
symbol wy  	= w3		; Interpolator Y-word and Divisor
symbol carry = bit31		; Top bit of w1
symbol tabbase = 0		; Table of table bases (max 8 bytes)
symbol atanbase = tabbase + 8	; Address of base of ATAN lookup table (length 34 bytes)
symbol sinebase = atanbase + 34	; Start address of Sine table (length 34 bytes)
symbol tanindex = TAM * 0		; ArcTan Lookup table pointer
symbol sinindex = TAM * 1		; Sine/Cos Lookup table pointer

data tabbase,(atanbase,sinebase)	; Pointer to bases of tables

; 16 line segments
symbol TANSC = 8		; Scale multiplier
data atanbase,(0,0,45,10,70,20,56,30,241,39,98,49,126,58,61,67,151,75)	; "Balanced" error curves
data (137,83,17,91,47,98,231,104,57,111,45,117,197,122,5,128)			; ATan table

symbol SINSC = 4			; Scale multiplier
data sinebase,(0,0,141,12,252,24,45,37,2,49,94,60,39,71,66,81,144,90) 
data (1,99,125,106,243,112,83,118,144,122,157,125,117,127,255,127)		; Sine table ($7FFF = 1)

; TEST CODE

	sertxd(cr,lf,"DEGs X2%_SIN  COS ASIN {Sin/Cos} ATAN2 {TAN}")

	for b0 = 0 to 72 				; One revolution in 5 degree increments
		wx = b0 * 5			; Degrees
		w1 = sin wx				; X2 only
		sertxd(cr,lf,#wx," ") 
		if w1 > 127 then			; It's a negative value
			sertxd("-")
			w1 = w1 - 128
		endif
		sertxd(#w1,"% ")
		w1 = b0 ** 14563		; 910.222  => 14563/65536 + 910
		w1 = b0 * 910 + w1		; Next 5 degree step
		w12 = w1				; Copy to use with cosine
		call sine
		w10 = w1				; Store sine result
		call showS4dp
		w1 = w12
		call cosine
		w11 = w1				; Store cosine result
		call showS4dp
		w1 = w10
		call arcsin
		call showdeg			; Report the angle
		wx = w10				; Sine of initial angle
		wy = w11				; Cosine of initial angle
		call calcTAN			; TAN = SIN / COS
		wx = w10				; Sine of initial angle
		wy = w11				; Cosine of initial angle		
		call arctan2
		call showdeg
		w1 = w12
		call TAN				; Report the TAN value
	next
end

TAN:		; Using ATAN table in reverse
	b1 = w1 / 8192 + tanindex		; Set the octant flags and table pointer
	w1 = w1 // 8192 * 4 max 32767		; Mask and scale		Octant		0 1 2 3 4 5 6 7
	if bit8 = 1 then			; Mirror about 45 degrees	Mirror		0 1 0 1 0 1 0 1 bit8
		w1 = 32767 - w1			; 32768 overflows ddiv		Negative	0 0 1 1 0 0 1 1 bit9
	endif					;				Complement	0 1 1 0 0 1 1 0 bit8xorbit9
	sertxd(" {")
	if bit9 = 1 then
		sertxd("-")
	endif	
	call interpYtoX									
	w1 = w1 * 8
	if bit8 = bit9 then
		sertxd("0.")
	else						; Calculate 1 / w1	
		wy = w1 / 2			; Divisor (15 steps per loop)
		w1 = 32768
		wx = 0
		b1 = b1 & 7			; Set division to 15 bits
		call ddiv
		sertxd(#w1,".")
		w1 = 0				; Remainder in wx
		b1 = b1 & 7
		call ddiv				; Fractional part
	endif
	w1 = w1 ** 20000			; Scale 10,000 = unity
	call show4dp
	sertxd("}")			 	
return	

calcTAN:			; Using SIN/COS
	sertxd(" {")
	b1 = 0
	if wx => NEGX then
		wx = -wx
		b1 = 2
	endif
	if wy => NEGX then
		wy = -wy
		b1 = b1 xor 2
	endif	
	w1 = wx + wx
	wx = 0
	call ddiv				; Integer result in w1, remainder in wx
showTAN:
	if bit9 = 1 then 
		sertxd("-")
	endif
	sertxd(#w1,".")			; Integer part
	b1 = b1 & 7
	w1 = 0
	call ddiv				; Fractional part
	w1 = w1 ** 20000
	call show4dp
	sertxd("}")
return		

showS4dp:					; Report (fractional) w1, signed to 4 digits (i.e. * 10,000)
	sertxd(" ") 
	if w1 > NEGX then			
		sertxd("-")
		w1 = -w1
	endif
	w1 = w1 + 2 ** 20000 		; Convert to 4 decimal places
show4dp:					; Unsigned
	if w1 > 999 then nosupz
	sertxd("0")
	if w1 > 99 then nosupz 
	sertxd("0")
	if w1 > 9 then nosupz 
	sertxd("0")	
nosupz:
	sertxd(#w1)
	return

showdeg:
	w1 = w1 ** 36000		; Convert to hundredths of a degree
	b1 = w1 // 100			; Hundredths
	w1 = w1 / 100			; Units
	sertxd(" ",#w1,".")
	if b1 < 10 then 
		sertxd("0")			; Don't suppress leading zero in decimal
	endif
	sertxd(#b1)
	return

cosine:
	w1 = w1 + $4000			; Now fall through into sine 
sine:						; Enter with angle in w1 (2^16 = 360 degrees)
	if w1 < NEGX then posine	; Shift down if 180 - 360 degrees	
		w1 = w1 - $8000
		call posine
		w1 = -w1
	return
posine:
	if w1 > $3FFF then		; Reverse 2nd and 4th quadrants
		w1 = $8000 - w1
	endif
	b1 = sinindex 			; Select Table + first quadrant (0)
	w1 = w1 / SINSC			; b1=Table number *32 (can retain flags in b1.0 to b1.4)
	goto interpXtoY			; Then use its return

interpXtoY:	 				; b1=Table number *32 (can retain flags in b1.0 to b1.4)
	wx = w1				; Free up b2 and b3. Note value in b3 used later as hi byte of w1  
	call getbase			; Get pointer to base of table into b2 
	b2 = b2 + b3 + b3		; Add 2* cellnumber (b3 is hi byte of w1) = SW corner of cell
	b3 = b2 + 2				; Point to RHS of cell
	read b3,WORD wy			; Read top (NE corner) of cell
	b3 = b2 + 1				; Point to high byte
	read b2,b2				; Read low byte of height of base of cell (for SW corner)
	read b3,b3				; Read high byte to form whole word in w1 (= height of baseline)
	wy = wy - w1 			; Calculate height of cell (range 0-65535)
	wx = wx / 256 * 65280 + wx	; Subtract Cell number * 256 from wx (65280 is -1 * 256) to give: 
	w1 = wx * 256 ** wy + w1	; Hor.Dist.inside cell (0-255) * gradient + height of cell base
	w1 = w1 max $7FFF 		; Avoid overflow into negative flag (probably unnecessary)
return 

arccos:
	call arcsin
	w1 = $4000 - w1			; Cos is shifted one quadrant (90 degrees)
return

arcsin:
	if w1 >= NEGX then 		; Value is negative
		w1 = -w1			; Convert to a positive value
		call parcsin
		w1 = -w1			; Shift result by two quadrants
	else					; Positive calculation
parcsin:
		b1 = sinindex 		; Select Sine Table
		call interpYtoX
		w1 = w1 * SINSC	
	endif
return

getbase:					; Use flags in b1 to fetch base address of required table into b2
	b2 = b1 / TAM + tabbase
	read b2,b2
return

interpYtoX:					; Enter with table number (*32) and flags in b1 (retained)
	wx = w1				; Free up b2 and b3 to use separately
	call getbase
	do
		b2 = b2 + 2			; Move to RHS of (next) cell
		read b2,word wy		; Read Y-value at top of cell (NE corner)
	loop until wy => wx		; Repeat until target point is within cell
	dec b2				; Back to high byte of LHS
	read b2,b3
	wx = -b3 * 256 + wx   		; Subtract high byte of base of cell
	wy = -b3 * 256 + wy
	b3 = b2 - 1				; Pointer to low byte
	read b3,b2
	wx = wx - b2			; Subtract low byte of base of cell
	wy = wy - b2
	call getbase
	b3 = b3 - b2 			; Cell number * 2

#ifdef fastdiv				; Fast version if cell height is restricted to 4095
	wx = wx * 16			; Multiply numerator by 16
	wy = wy / 16			; Divide divisor by 16
#else						; Scale numerator and divisor for best resolution
	for b2 = 0 to 7			; Multiply numerator/divisor by 256 to give result = 0 - 255
		if wx < 32768 then 	; WX can be doubled
			wx = wx + wx	; Multiply numerator by 2
		else
			wy = wy / 2	; Divide divisor by 2
		endif
	next 	
#endif				
	wx = wx / wy max 255		; Horizontal distance within cell  (max = 255)
	w1 = b3 * 128 + wx		; Add the value of the LH edge of cell
return

ddiv:						; Divide double-word numerator by 15-bit divisor
	for b1 = b1 to 239 step 16	; Retain starting value and repeat for each bit position
		wx = wx + wx + carry	; Rotate numerator and result left (top bit lost)
		w1 = w1 + w1		; Shift done
		if wx => wy then 		; Skip if can't subtract
			w1 = w1 + 1	; Add the flag to result (in low word)
			wx = wx - wy	; Subtract		
		endif
	next b1
	b1 = b1 & 15			; Clear the counter
return
 
arctan2:
	b1 = 48 		; Clear three lowest flag bits for first octant and select 12-bit ddiv
	if wx >= NEGX then		; Value is negative
		b1 = b1 xor 7		; Set the three flags
		wx = -wx			; Convert negative to positive
	endif					; wx is (now) positive
	if wy >= NEGX then 		; Value is negative
		b1 = b1 xor 3
		wy = -wy
	endif					; wy is (now) positive
	if wx > wy then 			; wx/wy is not a proper fraction (< 1) so swap
		b1 = b1 xor 1		; Mark that wx > wy (Octant 1,2 etc)
		w1 = wx			; )
		wx = wy			; }-Exchange wx and wy (~5 code bytes fewer than SWAP)
		wy = w1			; )
	endif					; Divisor is now larger than numerator so "fractional" result
	w1 = 0
	call ddiv				; Result in w1=wx/wy, Can now use wx,wy as temp variables
arctan:					; Enter with positive 12-bit value in w1 (and 3 flags set)
	b1 = b1 & 7 + tanindex
; Scale here if needed  	
	call interpXtoY			; Fullscale = 4096 for one octant, returns max $8000
	w1 = w1 + 2 / 4			; Round up and scale down
	if bit8 = 1 then			; LSB of b1
		w1 = $2000 - w1		; Reverse the sector from 90 to 45 degrees
	endif
	w1 = b1 & 7 * $2000 + w1	; Add the appropriate octant number (45 degrees = 4096)
return
Cheers, Alan.
 
Top