# "General Purpose" Interpolation subroutines

#### AllyCat

##### Senior Member
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
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
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
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
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
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
b3 = b2 + 1				; Point to 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
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

b2 = b1 / TAM + tabbase				; b1 contains flags for table number
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
wx = -b3 * 256 + wx   				; Subtract high byte of base of cell
wy = -b3 * 256 + wy
b3 = b2 - 1					; Pointer to low byte
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
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
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
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
wx = -b3 * 256 + wx   		; Subtract high byte of base of cell
wy = -b3 * 256 + wy
b3 = b2 - 1				; Pointer to low byte
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.

#### AllyCat

##### Senior Member
Hi,

Here's an update to the above program, which may take more than an hour to fully run in the Simulator! Coming back to (use) these Subroutines after more than three years, I noticed a few omissions and details that weren't very clear. Therefore, I've revised some comments and added tables and "scaling" constants to give a range of resolutions/table sizes for all the Trigonometric functions. As before, the complete program is too large to embed fully in a single post, so this one will contain the DATA tables, the declarations and test harness. My next post will contain the fundamental Subroutines, any of which can be added into other programs as a "Library" of functions, perhaps after pre-programming the tables and then with a #NO_DATA in the header. These Subroutines can also use TABLE memory, but here they use the EEPROM, to run even in an 08M2. Optional space has been reserved at the base of the memory, for User Data, or for any other Table which might require a zero base address.

The useful ATAN2 function dates back to early computing, such as the FORTRAN language and some clarification is necessary. Conventionally, the ATAN function calculates an angle defined by Y / X , but it cannot give a unique solution (over a full cycle of 360 degrees), because -Y / -X has the same value (normally both are considered to be in the range +/- 90 degrees). To avoid this, ATAN2 keeps the X and Y values separated in their "raw" form, rather confusingly in the sequence (Y , X), with the angle specified ANTI-Clockwise from the positive X-Axis (i.e. an "Eastwards" line is the zero degrees reference). Therefore, the centre of the first Quadrant (numbered zero) and the start of the second Octant has a value of +45 degrees, or the "NE" direction. Further, the second Quadrant (numbered 1) is to the NW, with Quadrants 2 and 3 in the SW and SE directions respectively. To add to the confusion, although Q-Basic defines the parameters in the sequence (Y , X), Excel and other spreadsheets declare them as (X , Y).

Also, "Geographical" and "Meteorological" calculations generally define "North" (or "Ahead") as Zero Azimuth, with angles computed CLOCKWISE from there (i.e. Quadrants 0 , 1 , 2 and 3 are to the NE , SE , SW and NW respectively). Also, the DIVISION subroutines which I have previously described here, calculate X / Y (not Y / X), so to avoid unnecessary complication (and maybe confusion), ATAN and ATAN2 are computed here using X / Y. Fortunately, these two standards can be converted by exchanging their X and Y parameters, so the Subroutines described here can be considered as using Clockwise angles from the "North = Zero" frame of reference.
Code:
``````; General Purpose Interpolation routines with lookup tables for Trig Functions;
; Angles in Trig functions use One revolution = 65536 (16 bits) and Signed Fractional Binary notation (+/-Unity = \$8000)
; AllyCat , July 2017, updated April 2018 and December 2021 [multiple tables for SINE and (ARC)TAN]

#picaxe 20m2                ; Or X2 because M2s cannot compare the embedded X2 SIN function
; #no_data               ; Can Uncomment this line when tables have been loaded
; #define fastdiv          ; Comment out if InterpYtoX has steep gradients (> 16), i.e. X/Y => 4096/256
#define SINSEG32           ; Max 130 bytes of lookup table (EEPROM)
#define TANSEG16           ; Max 66 bytes of lookup table (EEPROM)
#define notX2              ; Don't (try to) show X2 SINs

symbol NEGX = 32768        ; Boundary of Negative numbers (Twos complement)
symbol TOM = 32        ; Lookup Table number Multiplier
symbol tempb = b1          ; Temporary (local), Parameter-passing byte and Flags
symbol tempw = w1         ; Temporary and Parameter-passing word
symbol wx     = w2          ; Interpolation X-word and Division Extension word
symbol wy      = w3            ; Interpolation Y-word and Divisor
symbol wz    = w4            ; Only required for TAN = SIN/COS in this listing
symbol wa     = w5         ; Auxiliary variable for multiplication (can overlay wx, except for TAN)
symbol carry = bit31        ; Top bit of w1
symbol TABBASE = 16      ; Table of table base pointers (max 8 bytes, 4 reserved)
symbol TANBASE = TABBASE + 4     ; Address of base of ATAN lookup table (max length 66 bytes)
symbol SINBASE = TANBASE + 54     ; Start address of Sine table (max length 130 bytes)
; Currently next spare byte = 120
symbol SINTAB = 0 * TOM        ; Sine/Cos Lookup table name (table number)
symbol TANTAB = 1 * TOM

data TABBASE,(SINBASE,TANBASE)
#IFDEF TANSEG32    ; 32 line segments (66 bytes)
symbol TANSC = 4                ; For 32 line segments
data TANBASE,(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)
#ENDIF     ; /TANSEG32
#IFDEF TANSEG16    ; 16 line segments (34 bytes)
symbol TANSC = 8                 ; For 16 line segements
data TANBASE,(0,0,45,10,70,20,56,30,241,39,98,49,126,58,61,67,151,75)        ; ATAN table (from X axis)
data (137,83,17,91,47,98,231,104,57,111,45,117,197,122,5,128)                    ; "Balanced" error curves
#ENDIF    ; /TANSEG16
#IFDEF TANSEG8        ; 8 line segments (18 bytes)
symbol TANSC = 16         ; For 8 line segments (MUST enable "fastdiv")
data TANBASE,(0,0,70,20,241,39,126,58,151,75,17,91,231,104,45,117,5,128)    ; Not optimised for 8 lines
#ENDIF    ; /TANSEG8

#IFDEF SINSEG64        ; 64 line segments (130 bytes)
symbol SINSC = 1        ; Scaling multiplier or Divisor
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)
#ENDIF    ; /SINSEG64
#IFDEF SINSEG32        ; 32 line segments (66 bytes)
symbol SINSC = 2        ; Scaling for 360 degrees = 65536 and +/-1 = 32768
data SINBASE,(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)
#ENDIF    ; /SINSEG32
#IFDEF SINSEG16        ; 16 line segments (34 bytes)
symbol SINSC = 4        ; Scale multiplier
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 (\$7FFF = 1)
#ENDIF    ; /SINSEG16

; TEST/DEMO CODE
sertxd(cr,lf,"DEG X2% .SIN   {.TAN}  .COS  {ASIN} Sin/Cos ATAN2")
for b0 = 1 to 72          ; One revolution in 5 degree increments
wx = b0 * 5           ; Degrees
sertxd(cr,lf,#wx," ")
#IFDEF notX2
sertxd("-X2 ")
#ELSE      ; It's an X2
w1 = sin wx                    ; X2 only
if w1 > 127 then      ; It's a negative value
sertxd("-")
w1 = w1 - 128     ; NOT 2's complement !
endif
sertxd(#w1,"% ")
#ENDIF    ; notX2
w1 = b0 ** 14563       ; 910.222  => 14563/65536 + 910
w1 = b0 * 910 + w1    ; Next 5 degree step
wa = w1        ; Copy to use with cosine
call SINE
wz = w1         ; Store sine to calc TAN from Sin/Cos table
call show4dpS
w1 = wa
sertxd(" {")
call TAN        ; Report the TAN value
sertxd("} ")
w1 = wa
call COSINE
wa = w1        ; Store cosine result
call show4dpS
w1 = wz
call ARCSIN
sertxd(" {")
call showdeg     ; Report the angle
sertxd("} ")
wx = wz        ; Sine of initial angle
wy = wa        ; Cosine of initial angle
call TANdiv    ; TAN = SIN / COS
wx = wz      ; Sine of initial angle
wy = wa           ; Cosine of initial angle
call ARCTAN2
call showdegs
next
end

; FORMATTING ROUTINES:
show4dpS:     ; Report (fractional) w1, signed to 4 digits (i.e. * 10,000)
; Corrupts only W1
if w1 > NEGX then
sertxd("-")
w1 = -w1
endif
w1 = w1 + 2 ** 20000         ; Convert to 4 decimal places
show4dp:               ; Unsigned
if w1 > 9999 then
sertxd("1.0000")
else
sertxd("0")
show4:
sertxd(".")
if w1 > 999 then nosupz
sertxd("0")
if w1 > 99 then nosupz
sertxd("0")
if w1 > 9 then nosupz
sertxd("0")
nosupz:
sertxd(#w1)
endif
return
showdegs:
sertxd(" ")
showdeg:
w1 = w1 ** 36000        ; Convert to hundredths of a degree (Changes W1 and B1)
b1 = w1 // 100       ; Hundredths
w1 = w1 / 100      ; Units
sertxd(#w1,".")     ; Integer part
if b1 < 10 then
sertxd("0")     ; Don't suppress leading zero in decimal
endif
sertxd(#b1)
return

; ** PASTE MAIN SUBROUTINES HERE **
; Typically :
SINE:
COSINE:
TAN:           ; Reverse Lookup (Y -> X)
TANdiv:      ; = SIN / COS
ARCTAN:
ARCTAN2:   ; X and Y (signed) inputs
ARCSIN:
FDIV:          ; Fractional Division  (i.e. +/- Unity = \$8000 = 32768)
interpXtoY:   ; Main interpolation routines
interpYtoX:
; etc.
RETURN``````
Cheers, Alan.

Last edited:

#### AllyCat

##### Senior Member
Hi,

And here should be the subroutines, but at the moment not fully tested when combined with the code in the post above. I'll revise this post later if I, or anybody else, discovers any issues.
Code:
``````; SUBROUTINES:
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 = SINTAB                 ; 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
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

TAN:        ; Uses ATAN table in reverse
b1 = w1 / 8192 + TANTAB                ; 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 fdiv    Negative        0 0 1 1 0 0 1 1 bit9
endif                                ;            Complement    0 1 1 0 0 1 1 0 bit8xorbit9
if bit9 = 1 then
sertxd("-")
endif
call interpYtoX
w1 = w1 * TANSC                ; 8 for 16 segments
if bit8 = bit9 then
sertxd("0")
else                                ; Calculate 1 / w1
wy = w1 / 2                    ; Divisor (15 steps per loop)
w1 = 32768
wx = 0
b1 = b1 & 15                ; Set division to 15 bits (counter = 000xXXXX)
call fdiv
sertxd(#w1)    ;,".")
w1 = 0                        ; Remainder in wx
b1 = b1 & 15
call fdiv                    ; Fractional part
endif
w1 = w1 ** 20000                ; Scale 10,000 = unity
call show4
return

TANdiv:            ; Uses TAN = SIN/COS with X=SIN and Y=COS already calculated
b1 = 0
if wx => NEGX then
wx = -wx
b1 = 2                        ; Sign Flag
endif
if wy => NEGX then
wy = -wy
b1 = b1 xor 2
endif
w1 = wx + wx
wx = 0
call fdiv                        ; Integer result in w1, remainder in wx
if bit9 = 1 then
sertxd("-")
endif
sertxd(#w1)                        ; Integer part
b1 = b1 & 15
w1 = 0
call fdiv                        ; Fractional part
w1 = w1 ** 20000
call show4
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 = SINTAB         ; 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 / TOM + tabbase
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
wx = -b3 * 256 + wx           ; Subtract high byte of base of cell
wy = -b3 * 256 + wy
b3 = b2 - 1                        ; Pointer to low byte
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 max cell height is restricted to 4095
wx = wx * 16                    ; Multiply numerator by 16
wy = wy / 16                    ; Divide divisor by 16
#ELSE        ; /fastdiv                ; 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    ; /fastdiv
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

FDIV:                    ; Unsigned Fractional Division 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:    ; Mathematically: TAN is Y/X (maybe (Y,X) or (X,Y)),but FDIV is X:w1/Y *!
b1 = %0000                    ; Default swaps +Y > +X (for required Y/X)
;x Clear three lowest flag bits for first octant and select 15-bit FDIV
if wx >= NEGX then        ; Value is negative    (Can use 13-bit FDIV and change TANSC)
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 9            ; Mark that wx > wy (Octant 1,2 etc) and "swap" flag
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 fdiv                    ; Result in w1=wx/wy, Can now use wx,wy as temp variables
ARCTAN:                            ; Enter with positive 15-bit value in w1 (and 3 flags set) MAX 1 octant
b1 = b1 + TANTAB            ; Select Table
w1 = w1 / TANSC            ; Correct scale: Unity in (32768) -> 4096 (16 segs)
call interpXtoY            ; Fullscale = 8192 for one octant, returns max \$8000
w1 = w1 + 2 / 4            ; Round up and scale down to Octant = 8192
if bit8 = 1 then            ; LSB of b1
w1 = 8192 - w1        ; Reverse the sector from 90 to 45 degrees
endif
w1 = b1 & 7 * 8192 + w1    ; Add the appropriate octant number (45 degrees = 4096)
return
; FORMATTING ROUTINES:     ; In previous post``````
Cheers, Alan.

Last edited: