Reply by Keith E. Larson ●November 12, 20032003-11-12
Hello Dan and Martin

The GasLight library by Gary Sitton is a great library. Writing math
functions is what I cut my teeth on!

Keep in mind that over time these functions have evolved a bit and are in
use with todays compiler. Basically to invoke 40b math you need to define
your variables as 'long double' instead of double. The only trouble I can
think of is wether or not the 'long double' library contains
transcendentals. The 40b basics like mpy, add, sub and divide are
definitely there, but Im not sure if the transcendentals are just calls 32b
versions.

If you want to try your hand writing a math function here are some helpful hints

If you define and link a math routine of your own making *before* it is
linked from the RTS library and the compiler will igonore subsequent finds.
That is... calls your function instead.

If interested, I started but did not finish a 96 bit (36 exp, 64 mant)
library for the C3x. It is not well tested, but the list of functions is to
the point that other functions can be derived using C code. The list
includes TADD, TSUB, TCOMP, TMPY, TRCPF,TSQRT, TEXP, TLOG, TLOG2, TSIN. I
would like to add this to the DSK package, but I have not had the time to
put the final package together. If anyone is interested I might just have
to reopen this issue!

Depending on the nature of the math you are doing, you may find that
precision may not be needed. An example would be Log Differential
Compression (LDC) that works with audio and video. Basically if you express
the first derivitive of a signal in floating point it is almost like ADPCM
(constrained to powers of 2) so you dont need to carry much precision. At
the backend, you simply re-integrate the signal. But what is way cool is
that by being a linear process you can perform DSP math on the compressed
data stream (hence the TI patent)! Milage will obviously vary, but it might
be usefull.

DSK3 Code on TI-Web
-------------------
The TI-web page is back up so I will just be happy!

The eventual download mechanism will be a hyperlink from within a PDF file
(maybe the HL also needs to be inside DSK3DW?), but at the moment the
following should work for external download. You will also note that it
appears as a 'document' rather than a DSK3Vxxx.ZIP file. This is because
within TI the mechanism for updating 'documents' is a lot easier and a lot
more timely than 'software'.

Quite a bit has been done in the last few months so read the help file
version information. In particular, just this last week a new STDOUT
feature and set of applications have been added for string printing and
character displays. Basically rather than carrying the rather hefty STDIO
library, a set of functions were created that will string print float,
decimal and hex values and copy them in the properly packed character form
for display on the PC.

http://dspvillage.ti.com/docs/catalog/resources/appnoteabstract.jhtml?family
Id2&abstractName=spra976

Best regards,
Keith Larson
+-----------+
|Keith Larson |
|Member Group Technical Staff |
|Texas Instruments Incorporated |
| |
| 281-274-3288 |
| |
| www.micro.ti.com/~klarson |
|-----------+
| TMS320C3x/C4x/VC33 Applications |
| |
| TMS320VC33 |
| The lowest cost and lowest power 500 w/Mflop |
| floating point DSP on the planet! |
+-----------+

Reply by Martin Lopez Vela ●November 12, 20032003-11-12


Hello again Dan,

Yes it is, TI and Prentice Hall, is from a "Digital Signal Processing
Series"

Regards

Mart Lez
CIATEC -----Mensaje original-----
De: Dan Eyer [mailto:]
Enviado el: Martes, 11 de Noviembre de 2003 06:27 p.m.
Para: Martin Lopez Vela
Asunto: RE: [c3x] long double math on VC33? Thanks, this will help a lot!

Is this book a TI publication, or was it from
someone else? --- Martin Lopez Vela <> wrote:
> Hello Dan,
>
> In the book "Digital Signal processing Applications
> with the TMS320
> Family, vol. 3 " appear the some extended-precision,
> floating-point (40
> bit) math function:
>
> SINX
> COSX
> EXPX
> LNX
> ATANX
> SQRTX
> FPINVX
> FDIVX
> FMULTX
>
> And this is the code:
>


__________________________________


Reply by Martin Lopez Vela ●November 12, 20032003-11-12
Hello Dan,

In the book "Digital Signal processing Applications with the TMS320
Family, vol. 3 " appear the some extended-precision, floating-point (40
bit) math function:

SINX
COSX
EXPX
LNX
ATANX
SQRTX
FPINVX
FDIVX
FMULTX

And this is the code:

*********************************************************
* PROGRAM: SINX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PRECISION SINE FUNCTION: R0 <= SIN(R0). *
* *
* APPROXIMATE ACCURACY: 9 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: NONE. *
* REGISTERS FOR INPUT: R0 (ARGUMENT IN RADIANS). *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: AR0, IR0, AND R0-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: FMULTX. *
* EXECUTION CYCLES (MIN, MAX): 160, 160. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL SINX
.GLOBL ECOSX
.GLOBL FMULTX

; INTERNAL CONSTANTS

.DATA

; SCALING COEFFS. FOR SIN(X)

NRM2 .WORD 00000006FH ; BOTTOM OF 2/PI
NRM1 .WORD 0FF22F983H ; TOP OF 2/PI

; POLYNOMIAL COEFFS. FOR SIN(X)

SHF2 .WORD 0000000A3H ; BOTTOM OF C1 (PI/2)
SHF1 .WORD 000490FDAH ; TOP OF C1 (PI/2)
.WORD 0000000D1H ; BOTTOM OF C3
.WORD 0FFDAA218H ; TOP OF C3
.WORD 0000000E3H ; BOTTOM OF C5
.WORD 0FC2335E0H ; TOP OF C5
.WORD 0F8E69754H ; TOP OF C7
.WORD 0F3280B28H ; TOP OF C9
COF .WORD 0ED9997B4H ; TOP OF C11

ACOF .WORD COF ; ADDRESS OF COEFFS.

CON .FLOAT -1.0, 0.0, 1.0, 0.0 ; MAPPING CONSTS.

ACON .WORD CON ; ADDRESS OF CONSTS.

.TEXT

; START OF SINX PROGRAM

SINX:

PUSH DP ; SAVE DP
LDP @NRM1 ; LOAD DATA PAGE POINTER

; COSX ENTRY POINT

ECOSX:

; SCALE AND MAP VARIABLE X

PUSHF R0 ; SAVE ORIGINAL X
ABSF R0 ; R0 <= |X|
LDF @NRM1,R1 ; R1 <= TOP OF 2/PI
OR @NRM2,R1 ; OR IN BOTTOM OF 2/PI
CALL FMULTX ; R0 <= |X|*2/PI
FIX R0,IR0 ; IR0 <= INTEGER QUADRANT Q
FLOAT IR0,R1 ; R1 <= FLOATING QUADRANT Q
SUBF R1,R0 ; R0 <= X, -1 < X < 1
NEGF R0,R3 ; R3 <= -X
ADDI 1,IR0 ; R2 <= Q + 1
AND 3,IR0 ; IR0 <= TABLE INDEX
TSTB 2,IR0 ; LOOK AT 2ND LSB
LDFNZ R3,R0 ; IF 1 THEN R0 <= -X
LDP @ACON ; LOAD DATA PAGE POINTER
LDI @ACON,AR0 ; AR0 -> CONST. TABLE
ADDF *+AR0(IR0),R0 ; FINAL MAPPING, R0 <= X + C
NEGF R0,R3 ; R3 <= -X
LDI @ACOF,AR0 ; AR0 -> COEFF. TABLE

; EVALUATE TRUNCATED SERIES

LDF R0,R1 ; R1 <= X
CALL FMULTX ; R0 <= X**2
LDF R0,R1 ; R1 <= X**2

MPYF *AR0--,R1,R0 ; R0 <= X**2*C11
ADDF *AR0--,R0 ; R0 <= C9 + R0

MPYF R1,R0 ; R0 <= X**2*(C9 + R0)
ADDF *AR0--,R0 ; R0 <= C7 + R0

MPYF R1,R0 ; R0 <= X**2*(C7 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C5
OR *AR0--,R2 ; OR IN BOTTOM OF C5
ADDF R2,R0 ; R0 <= C5 + R0

CALL FMULTX ; R0 <= X**2*(C5 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C3
OR *AR0--,R2 ; OR IN BOTTOM OF C3
ADDF R2,R0 ; R0 <= C3 + R0

CALL FMULTX ; R0 <= X**2*(C3 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C1
OR *AR0,R2 ; OR IN BOTTOM OF C1
ADDF R2,R0,R1 ; R1 <= C1 + R0

; TEST FOR X < 0 AND RETURN

NEGF R3,R0 ; R0 <= X
BRD FMULTX ; R0 <= X*R1 = SIN(X), (DELAYED)
POPF R5 ; TEST ORIGINAL X
LDFN R3,R0 ; IF X < 0 THEN R0 <= -X
POP DP ; UNSAVE DP

; RETURN OCCURS FROM FMULTX !

*********************************************************
* PROGRAM: COSX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PRECISION COSINE FUNCTION: R0 <= COS(R0). *
* *
* APPROXIMATE ACCURACY: 9 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: NONE. *
* REGISTERS FOR INPUT: R0 (ARGUMENT IN RADIANS). *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: AR0, IR0, AND R0-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: ECOSX (SINX). *
* EXECUTION CYCLES (MIN, MAX): 165, 165. *
* *
* NOTE: USES SHF1 AND SHF2 FROM SINX PROGRAM! *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL COSX
.GLOBL ECOSX

.TEXT

; START OF COSX PROGRAM

COSX:

PUSH DP ; SAVE DP
LDP @NRM1 ; LOAD DATA PAGE POINTER

BRD ECOSX ; R0 <= COS(X) = SIN(X'), (DELAYED)
LDF @SHF1,R1 ; R1 <= TOP OF PI/2
OR @SHF2,R1 ; OR IN BOTTOM OF PI/2
ADDF R1,R0 ; R0 <= X' = X + PI/2

; RETURN OCCURS FROM SINX (ALIAS FMULTX) !

*********************************************************
* PROGRAM: EXPX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PREC. EXPONENTIAL: R0 <= EXP(R0). *
* *
* APPROXIMATE ACCURACY: 9 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: |R0| <= 88.0. *
* REGISTERS FOR INPUT: R0. *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: AR0 AND R0-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: FMULTX AND FPINVX. *
* EXECUTION CYCLES (MIN, MAX): 115 (R0 <=0 ), 160. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL EXPX
.GLOBL FMULTX
.GLOBL FPINVX

; INTERNAL CONSTANTS

.DATA

; SCALING COEFFS. FOR 2-**X

ENRM2 .WORD 000000029H ; BOTTOM OF 1/LN(2)
ENRM1 .WORD 00038AA3BH ; TOP OF 1/LN(2)

; POLYNOMIAL COEFFS. FOR 2**-X, 0 <= X < 1.

.WORD 000000000H ; C0 (1.0)
.WORD 00000000AH ; BOTTOM OF C1
.WORD 0FFCE8DE8H ; TOP OF C1
.WORD 00000006EH ; BOTTOM OF C2
.WORD 0FD75FDEDH ; TOP OF C2
.WORD 000000046H ; BOTTOM OF C3
.WORD 0FB9CA833H ; TOP OF C3
.WORD 0F91D8C56H ; TOP OF C4
.WORD 0F6D1E7A9H ; TOP OF C5
.WORD 0F31AA7D7H ; TOP OF C6
C7 .WORD 0EFC9BD9CH ; TOP OF C7

AC7 .WORD C7

.TEXT

; START OF EXPX PROGRAM

EXPX:

; SCALE VARIABLE X

PUSH DP ; SAVE DP
LDP @AC7 ; LOAD DATA PAGE POINTER
NEGF R0,R2 ; R2 <= -X
LDF R0,R1 ; R1 <= X
LDFN R2,R0 ; IF X < 0 THEN R1 <= |X|
LDF @ENRM1,R1 ; R1 <= TOP OF 1/LN(2)
OR @ENRM2,R1 ; OR IN BOTTOM OF 1/LN(2)
CALL FMULTX ; R0 <= X = |X|/LN(2)
FIX R0,R3 ; R3 <= I = INTEGER OF X
FLOAT R3,R1 ; R1 <= FLT. PT. I
SUBF R1,R0,R1 ; R1 <= FRACTION OF |X|, 0 <= X < 1
NEGI R3 ; R3 <= -I
LSH 24,R3 ; MOVE -I TO EXP.
PUSH R3 ; SAVE AS INT.
POPF R3 ; R3 <= FLT. PT. 2**-I
LDI @AC7,AR0 ; AR0 -> COEFF. TABLE
POP DP ; UNSAVE DP

; EVALUATE TRUNCATED SERIES

MPYF *AR0--,R1,R0 ; R0 <= X*C7
ADDF *AR0--,R0 ; R0 <= C6 + R0

MPYF R1,R0 ; R0 <= X*(C6 + R0)
ADDF *AR0--,R0 ; R0 <= C5 + R0

MPYF R1,R0 ; R0 <= X*(C5 + R0)
ADDF *AR0--,R0 ; R0 <= C4 + R0

MPYF R1,R0 ; R0 <= X*(C4 + R0)
LDF *AR0--,R4 ; R4 <= TOP OF C3
OR *AR0--,R4 ; OR IN BOTTOM OF C3
ADDF R4,R0 ; R0 <= C3 + R0

MPYF R1,R0 ; R0 <= X*(C3 + R0)
LDF *AR0--,R4 ; R4 <= TOP OF C2
OR *AR0--,R4 ; OR IN BOTTOM OF C2
ADDF R4,R0 ; R0 <= C2 + R0

CALL FMULTX ; R0 <= X*(C2 + R0)
LDF *AR0--,R4 ; R4 <= TOP OF C1
OR *AR0--,R4 ; OR IN BOTTOM OF C1
ADDF R4,R0 ; R0 <= C1 + R0

CALL FMULTX ; R0 <= X*(C1 + R0)

; TEST FOR X < 0 AND RETURN

LDF R2,R2 ; TEST ORIGINAL -X
BND FPINVX ; IF -X < 0 THEN R0 <= 1/X, (DELAYED)
ADDF *AR0,R0,R1 ; R1 <= 2**-X = C0 + R0
MPYF R3,R1,R0 ; R0 <= 2**-(I + X) TRUNC.
LDM R1,R0 ; R0 <= FULL MANTISSA

RETS ; RETURN (IF NO FPINVX BRANCH)

*********************************************************
* PROGRAM: LNX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PREC. LOGARITHM BASE E: R0 <= LN(R0). *
* *
* APPROXIMATE ACCURACY: 8 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: R0 > 0.0. *
* REGISTERS FOR INPUT: R0. *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: AR0 AND R0-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: FMULTX. *
* EXECUTION CYCLES (MIN, MAX): 193, 193. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL LNX
.GLOBL FMULTX

; INTERNAL CONSTANTS

.DATA

; SCALING COEFFS. FOR LN(1+X)

LNRM2 .WORD 0000000F7H ; BOTTOM OF LN(2)
LNRM1 .WORD 0FF317217H ; TOP OF LN(2)

; POLYNOMIAL COEFFS. FOR LN(1+X), 0 <= X < 1.

C0 .FLOAT 1.0 ; C0 (1.0)

.WORD 0000000FFH ; BOTTOM OF C1
.WORD 0FF7FFFC3H ; TOP OF C1
.WORD 0000000B4H ; BOTTOM OF C2
.WORD 0FE80107FH ; TOP OF C2
.WORD 0000000D9H ; BOTTOM OF C3
.WORD 0FE29E18FH ; TOP OF C3
.WORD 000000097H ; BOTTOM OF C4
.WORD 0FD897D13H ; TOP OF C4
.WORD 000000041H ; BOTTOM OF C5
.WORD 0FD2BAD82H ; TOP OF C5
.WORD 0000000E7H ; BOTTOM OF C6
.WORD 0FCBCC3F1H ; TOP OF C6
.WORD 000000043H ; BOTTOM OF C7
.WORD 0FB13D187H ; TOP OF C7
C8 .WORD 0F8AC87BFH ; TOP OF C8

AC8 .WORD C8

.TEXT

; START OF LNX PROGRAM

LNX:

LDF R0,R0 ; TEST X
RETSLE ; RETURN NOW IF X <= 0

; SCALE VARIABLE X

PUSH DP ; SAVE DP
LDP @AC8 ; LOAD DATA PAGE POINTER
PUSHF R0 ; SAVE AS FLT. PT.
POP R3 ; R3 <= INTEGER FORMAT
ASH -24,R3 ; R3 <= E = SIGNED EXP.
FLOAT R3,R1 ; R1 <= FLT. PT. E VALUE
LDF @C0,R2 ; R2 <= 1.0
LDE R2,R0 ; EXP. R0 <= 0 (1 <= X < 2)
SUBRF R0,R2 ; R2 <= X - 1 (0 <= X < 1)
LDF @LNRM1,R0 ; R0 <= TOP OF LN(2)
OR @LNRM2,R0 ; OR IN BOTTOM OF LN(2)
CALL FMULTX ; R0 <= E*LN(2)
LDF R0,R3 ; R3 <= E*LN(2)
LDI @AC8,AR0 ; AR0 -> COEFF. TABLE
POP DP ; UNSAVE DP

; EVALUATE TRUNCATED SERIES

LDF R2,R1 ; R1 <= X
MPYF *AR0--,R1,R0 ; R0 <= X*C8
LDF *AR0--,R2 ; R2 <= TOP OF C7
OR *AR0--,R2 ; OR IN BOTTOM OF C7
ADDF R2,R0 ; R0 <= C7 + R0

MPYF R1,R0 ; R0 <= X*(C7 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C6
OR *AR0--,R2 ; OR IN BOTTOM OF C6
ADDF R2,R0 ; R0 <= C6 + R0

MPYF R1,R0 ; R0 <= X*(C6 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C5
OR *AR0--,R2 ; OR IN BOTTOM OF C5
ADDF R2,R0 ; R0 <= C5 + R0

CALL FMULTX ; R0 <= X*(C5 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C4
OR *AR0--,R2 ; OR IN BOTTOM OF C4
ADDF R2,R0 ; R0 <= C4 + R0

CALL FMULTX ; R0 <= X*(C4 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C3
OR *AR0--,R2 ; OR IN BOTTOM OF C3
ADDF R2,R0 ; R0 <= C3 + R0

CALL FMULTX ; R0 <= X*(C3 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C2
OR *AR0--,R2 ; OR IN BOTTOM OF C2
ADDF R2,R0 ; R0 <= C2 + R0

CALL FMULTX ; R0 <= X*(C2 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C1
OR *AR0--,R2 ; OR IN BOTTOM OF C1
ADDF R2,R0 ; R0 <= C1 + R0

CALL FMULTX ; R0 <= X*(C1 + R0)

; ADD IN SCALED EXPONENT.

ADDF R3,R0 ; R0 <= LN(X) + E*LN(2)

RETS ; RETURN

*********************************************************
* PROGRAM: ATANX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PRECISION ARC TANGENT: R0 <= ATAN(R0). *
* *
* APPROXIMATE ACCURACY: 8 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: NONE. *
* REGISTERS FOR INPUT: R0. *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: AR0, IR0, AND R0-7. *
* REGISTERS FOR OUTPUT: R0 (IN RADIANS). *
* ROUTINES NEEDED: FMULTX, AND FDIVX. *
* EXECUTION CYCLES (MIN, MAX): 210 (|ATANX|<=1), 332. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL ATANX
.GLOBL FMULTX
.GLOBL FDIVX

; INTERNAL CONSTANTS

.DATA

; SCALING COEFFS. FOR ATAN(X)

.WORD 00000005DH ; BOTTOM OF -PI/4
.WORD 0FFB6F025H ; TOP OF -PI/4
.WORD 0000000A2H ; BOTTOM OF PI/4
.WORD 0FF490FDAH ; TOP OF PI/4
.WORD 000000000H ; BOTTOM OF ZERO
.WORD 080000000H ; TOP OF ZERO

; POLYNOMIAL COEFFS. FOR ATAN(X), -1 <= X <= 1.

C1 .WORD 000000000H ; TOP OF C1 (1.0)
.WORD 00000006EH ; BOTTOM OF C3
.WORD 0FED55594H ; TOP OF C3
.WORD 0000000D9H ; BOTTOM OF C5
.WORD 0FD4CBBE4H ; TOP OF C5
.WORD 0000000FFH ; BOTTOM OF C7
.WORD 0FDEE8038H ; TOP OF C7
.WORD 000000056H ; BOTTOM OF C9
.WORD 0FC5A3D83H ; TOP OF C9
.WORD 000000093H ; BOTTOM OF C11
.WORD 0FCE5CE8BH ; TOP OF C11
.WORD 0000000BFH ; BOTTOM OF C13
.WORD 0FB2FC1FDH ; TOP OF C13
.WORD 0FAFB91FEH ; TOP OF C15
C17 .WORD 0F73BD74AH ; TOP OF C17

AC17 .WORD C17

.TEXT

; START OF ATANX PROGRAM

ATANX:

; SCALE VARIABLE X

PUSH DP ; SAVE DP
LDP @AC17 ; LOAD DATA PAGE POINTER
ABSF R0,R2 ; R2 <= |X|
SUBF @C1,R2 ; R2 <= |X| - 1
BLED SKIP ; IF |X| > 1 THEN SCALE (DELAYED)
LDF R0,R3 ; R3 <= X
LDF R0,R1 ; R1 <= X
LDI 0,IR0 ; IR0 <= 0, POST SCALE INDEX

; SCALE FOR |X| > 1

PUSHF R0 ; SAVE X
ABSF R0,R1 ; R1 <= |X|
ADDF @C1,R1 ; R1 <= |X| + 1
LDF R2,R0 ; R0 <= |X| - 1
CALL FDIVX ; R0 <= (|X| - 1)/(|X| + 1)

; TEST FOR X' < 0

POPF R4 ; GET ORIGINAL X
BGED SKIP ; IF X < 0 THEN R0 <= -X' (DELAYED)
LDF R0,R3 ; R3 <= X'
LDF R0,R1 ; R1 <= X'
SUBI 2,IR0 ; IR0 <= -2, (PI/4)

NEGF R0,R3 ; R3 <= -X'
SUBI 2,IR0 ; IR0 <= -4, (-PI/4)

SKIP: CALL FMULTX ; R0 <= X**2
LDI @AC17,AR0 ; AR0 -> COEFF. TABLE
POP DP ; UNSAVE DP

; EVALUATE TRUNCATED (ODD) SERIES

LDF R0,R1 ; R1 <= X**2
MPYF *AR0--,R1,R0 ; R0 <= X**2*C17
ADDF *AR0--,R0 ; R0 <= C15 + R0

MPYF R1,R0 ; R0 <= X**2*(C15 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C13
OR *AR0--,R2 ; OR IN BOTTOM OF C13
ADDF R2,R0 ; R0 <= C13 + R0

MPYF R1,R0 ; R0 <= X**2*(C13 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C11
OR *AR0--,R2 ; OR IN BOTTOM OF C11
ADDF R2,R0 ; R0 <= C11 + R0

CALL FMULTX ; R0 <= X**2*(C11 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C9
OR *AR0--,R2 ; OR IN BOTTOM OF C9
ADDF R2,R0 ; R0 <= C9 + R0

CALL FMULTX ; R0 <= X**2*(C9 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C7
OR *AR0--,R2 ; OR IN BOTTOM OF C7
ADDF R2,R0 ; R0 <= C7 + R0

CALL FMULTX ; R0 <= X**2*(C7 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C5
OR *AR0--,R2 ; OR IN BOTTOM OF C5
ADDF R2,R0 ; R0 <= C5 + R0

CALL FMULTX ; R0 <= X**2*(C5 + R0)
LDF *AR0--,R2 ; R2 <= TOP OF C3
OR *AR0--,R2 ; OR IN BOTTOM OF C3
ADDF R2,R0 ; R0 <= C3 + R0

CALL FMULTX ; R0 <= X**2*(C3 + R0)

; FINISH UP

ADDF *AR0--,R0,R1 ; R1 <= C1 + R0
LDF R3,R0 ; R0 <= X (SIGNED)
CALL FMULTX ; R0 <= ATANX(X) = X*(1 + R0)
NOP *AR0++(IR0) ; AR0 -> C (0.0, PI/4 OR -PI/4)

; ADD IN POST SCALE VALUE C AND RETURN

POP R4 ; R4 <= RETURN ADDRESS
BUD R4 ; RETURN (DELAYED)
LDF *AR0--,R1 ; R1 <= TOP OF C
OR *AR0,R1 ; OR IN BOTTOM OF C
ADDF R1,R0 ; R0 <= ATAN(X) + C

*********************************************************
* PROGRAM: SQRTX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* APPROXIMATE ACCURACY: 10 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: R0 >= 0.0. *
* REGISTERS FOR INPUT: R0. *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: R0-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: FMULTX. *
* EXECUTION CYCLES (MIN, MAX): 138, 138. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL SQRTX
.GLOBL FMULTX

; INTERNAL CONSTANTS

.DATA

CNST1 .SET 0.5
CNST2 .SET 1.5
CNST3 .FLOAT 1.103553391 ; ADJUSTED 1.0
CNST4 .FLOAT 0.780330086 ; ADJUSTED SQRT(1/2)

SMSK .WORD 0FF7FFFFFH

.TEXT

; START OF SQRTX PROGRAM.

SQRTX:

LDF R0,R3 ; TEST AND SAVE V
RETSLE ; RETURN NOW IF V <= 0

; GET APPROXIMATION TO 1/V. FOR V = (1+M)*2**E
; AND 0 <= M < 1, FOR E EVEN: X[0] = (1-M/2)*2**-E/2
; AND FOR E ODD: X[0] = SQRT(1/2)*(1-M/2)*2**-E/2

PUSH DP ; SAVE DP
LDP @SMSK ; LOAD DATA PAGE POINTER
PUSHF R0 ; SAVE V AS FLT. PT. V = (1+M)*2**E
POP R4 ; R4 <= V AS INTEGER
XOR @SMSK,R4 ; R4 <= COMPLEMENT ALL BUT SIGN
LDI R4,R1 ; R1 <= (1-M/2)*2**-E
LDI R4,R5 ; R5 <= R1
LSH 8,R1 ; R1 <= R1 EXP. REMOVED
ASH -1,R4 ; R4 <= R4 WITH -E/2 EXP.
PUSH R4 ; SAVE R4 AS INTEGER
POPF R4 ; R4 <= FLT. PT.
LDE R4,R1 ; R1 <= (1-M/2)*2**-E/2
LDF @CNST3,R2 ; R2 <= 1.1... FOR ODD E
LSH 7,R5 ; TEST LSB OF E (AS SIGN)
LDFNN @CNST4,R2 ; IF E EVEN R2 <= 0.78...
MPYF R2,R1 ; R1 <= CORRECTED ESTIMATE

; GENERATE V/2 (USES MPYF).

MPYF CNST1,R0 ; R0 <= V/2 TRUNC.
LDM R3,R0 ; R0 <= V/2 FULL PREC.

; NEWTON ITERATION FOR Y(X) = X - V**-2 = 0 ...

MPYF R1,R1,R2 ; R2 <= X[0]**2
MPYF R0,R2 ; R2 <= (V/2) * X[0]**2
SUBRF CNST2,R2 ; R2 <= 1.5 - (V/2) * X[0]**2
MPYF R2,R1 ; R1 <= X[1] = X[0] * (1.5 -
(V/2)*X[0]**2)
;
MPYF R1,R1,R2 ; R2 <= X[1]**2
MPYF R0,R2 ; R2 <= (V/2) * X[1]**2
SUBRF CNST2,R2 ; R2 <= 1.5 - (V/2) * X[1]**2
MPYF R2,R1 ; R1 <= X[2] = X[1] * (1.5 -
(V/2)*X[1]**2)
;
MPYF R1,R1,R2 ; R2 <= X[2]**2
MPYF R0,R2 ; R2 <= (V/2) * X[2]**2
SUBRF CNST2,R2 ; R2 <= 1.5 - (V/2) * X[2]**2
MPYF R2,R1 ; R1 <= X[3] = X[2] * (1.5 -
(V/2)*X[2]**2)
;
LDF R0,R2 ; R2 <= V/2
LDF R1,R0 ; R0 <= X[3]
CALL FMULTX ; R0 <= X[3]**2
LDF R1,R4 ; R4 <= X[3]
LDF R2,R1 ; R1 <= V/2
LDF R4,R2 ; R2 <= X[3]
CALL FMULTX ; R0 <= (V/2) * X[3]**2
SUBRF CNST2,R0 ; R0 <= 1.5 - (V/2) * X[3]**2
LDF R2,R1 ; R1 <= X[3]
CALL FMULTX ; R0 <= X[4] = X[3] * (1.5 -
(V/2)*X[3]**2)

; INVERT FINAL RESULT AND RETURN

BRD FMULTX ; R0 = SQRT(V) = V*SQRT(1/V) (DELAYED)
LDF R3,R1 ; R1 <= V
POP DP ; UNSAVE DP
NOP ; DEAD CYCLE

; RETURN OCCURS FROM FMULTX !

*********************************************************
* PROGRAM: FDIVX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PRECISION DIVIDE: R0 <= R0/R1. *
* *
* APPROXIMATE ACCURACY: 10 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: R1 != 0.0. *
* REGISTERS FOR INPUT: R0 (DIVIDEND) AND R1 (DIVISOR).*
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: R0-7. *
* REGISTERS FOR OUTPUT: R0 (QUOTIENT). *
* ROUTINES NEEDED: FMULTX AND FPINVX. *
* EXECUTION CYCLES (MIN, MAX): 107, 107. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL FDIVX
.GLOBL FPINVX
.GLOBL FMULTX

.TEXT

; START OF FDIVX PROGRAM

FDIVX:

LDF R0,R3 ; R3 <= X
LDF R1,R0 ; R1 <= Y
CALL FPINVX ; R0 <= 1/Y
LDF R3,R1 ; R1 <= X
BR FMULTX ; R0 <= X/Y

; RETURN OCCURS FROM FMULTX !

*********************************************************
* PROGRAM: FMULTX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PRECISION MULTIPLY: R0 <= R0*R1. *
* *
* APPROXIMATE ACCURACY: 10 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: NONE. *
* REGISTERS FOR INPUT: R0. *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: R0 AND R4-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: NONE. *
* EXECUTION CYCLES (MIN, MAX): 20, 20. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL FMULTX

.TEXT

; START OF FMULTX PROGRAM

FMULTX:

ABSF R0,R4 ; R4 <= |XA|
XOR R1,R0 ; R0 <= SIGN INFO.
ABSF R1,R7 ; R7 <= |XB|
MPYF R4,R7,R6 ; R6 <= A*B
LDF R4,R5 ; R5 <= |XA|
ANDN 0FFH,R5 ; R5 <= A = XA - EA*2**-24
SUBRF R4,R5 ; R5 <= EA*2**-24
MPYF R7,R5 ; R5 <= B*EA*2**-24
ADDF R6,R5 ; R5 <= A*B + B*EA*2**-24
LDF R7,R6 ; R6 <= |XB|
ANDN 0FFH,R6 ; R6 <= B = XB - EB*2**-24
SUBRF R7,R6 ; R6 <= EB*2**-24
MPYF R4,R6 ; R6 <= A*EB*2**-24
ADDF R6,R5 ; R5 <= |XA*XB| = A*B +
(B*EA+A*EB)*2**-24
NEGF R5,R6 ; R6 <= - |XA*XB|

; TEST FOR XA*XB < 0 AND RETURN

POP R4 ; R4 <= RETURN ADDRESS
BUD R4 ; RETURN (DELAYED)
LDF R0,R0 ; TEST ORIGINAL (XA ^ XB)
LDFN R6,R5 ; IF XA*XB < 0 THEN R5 <= -|XA*XB|
LDF R5,R0 ; R0 <= XA*XB

*********************************************************
* PROGRAM: FPINVX *
* *
* WRITTEN BY: GARY A. SITTON *
* GAS LIGHT SOFTWARE *
* HOUSTON, TEXAS *
* MARCH 1989. *
* *
* EXTENDED PREC. FLT. PT. INVERSE: R0 <= 1/R0. *
* *
* APPROXIMATE ACCURACY: 10 DECIMAL DIGITS. *
* INPUT RESTRICTIONS: R0 != 0.0. *
* REGISTERS FOR INPUT: R0. *
* REGISTERS USED AND RESTORED: DP AND SP. *
* REGISTERS ALTERED: R0-1 AND R4-7. *
* REGISTERS FOR OUTPUT: R0. *
* ROUTINES NEEDED: FMULTX. *
* EXECUTION CYCLES (MIN, MAX): 76, 76. *
*********************************************************

; EXTERNAL PROGRAM NAMES

.GLOBL FPINVX
.GLOBL FMULTX

; INTERNAL CONSTANTS

.DATA

ONE .SET 1.0
TWO .SET 2.0

MSK .WORD 0FF7FFFFFH

.TEXT

; START OF FPINVX PROGRAM

FPINVX:

LDF R0,R0 ; TEST F
RETSZ ; RETURN NOW IF F = 0

; GET APPROXIMATION TO 1/F. FOR F = (1+M) * 2**E
; AND 0 <= M < 1, USE: X[0] = (1-M/2) * 2**-E

PUSH DP ; SAVE DP
LDP @MSK ; LOAD DATA PAGE POINTER
PUSHF R0 ; SAVE AS FLT. PT. F = (1+M) * 2**E
POP R1 ; FETCH BACK AS INTEGER
XOR @MSK,R1 ; COMPLEMENT E & M BUT NOT SIGN BIT
PUSH R1 ; SAVE AS INTEGER, AND BY MAGIC...
POPF R1 ; R1 <= X[0] = (1-M/2) * 2**-E.
POP DP ; UNSAVE DP

; NEWTON ITERATION FOR: Y(X) = X - 1/F = 0 ...

MPYF R1,R0,R4 ; R4 <= F * X[0]
SUBRF TWO,R4 ; R4 <= 2 - F * X[0]
MPYF R4,R1 ; R1 <= X[1] = X[0] * (2 - F * X[0])

MPYF R1,R0,R4 ; R4 <= F * X[1]
SUBRF TWO,R4 ; R4 <= 2 - F * X[1]
MPYF R4,R1 ; R1 <= X[2] = X[1] * (2 - F * X[1])

MPYF R1,R0,R4 ; R4 <= F * X[2]
SUBRF TWO,R4 ; R4 <= 2 - F * X[2]
MPYF R4,R1 ; R1 <= X[3] = X[2] * (2 - F * X[2])

; FOR THE LAST ITERATION: X[4] = (X[3] * (1 - (F * X[3]))) + X[3]

CALL FMULTX ; R0 <= F * X[3] = 1 + EPS
SUBRF ONE,R0 ; R0 <= 1 - F * X[3] = EPS
CALL FMULTX ; R0 <= X[3] * EPS
ADDF R1,R0 ; R0 <= X[4] = (X[3]*(1 - (F*X[3]))) +
X[3]

RETS ; RETURN

.END

Maybe it can help you.

Mart Lez
CIATEC.
Omega 201, Fracc. Delta
Le, Gto., C.P. 37545
Mico
www.ciatec.mx

-----Mensaje original-----
De: dan_eyer [mailto:]
Enviado el: Martes, 11 de Noviembre de 2003 12:41 p.m.
Para:
Asunto: [c3x] long double math on VC33? Hi,

We're trying to port some code from a processor
with 64-bit floating-point variables to a VC33,
which only has 32-bit and 40-bit floats, and we
need to maintain as much precision as possible.

Does anyone know if the math library functions
like sqrt support the "long double" (i.e. the
40-bit) floats or not? Or if there are functions
meant just for the long doubles? I know about
MPY_LD and DIV_LD, but we need sqrt and other
functions with the full 40 bit precision.

The prototypes for sqrt etc only seem to list
"double" (32-bit on VC33) parameters, but it
seems strange to have a higher-precision format
and not support doing math in it.

Maybe there are third-party libraries for doing
64-bit floating-point math on the VC33?

Thanks,
Dan Eyer _____________________________________
Note: If you do a simple "reply" with your email client, only the author
of this message will receive your answer. You need to do a "reply all"
if you want your answer to be distributed to the entire group.

_____________________________________
About this discussion group:

To Join: Send an email to

To Post: Send an email to

To Leave: Send an email to

Archives: http://groups.yahoo.com/group/c3x

More Groups: http://www.dsprelated.com ">http://docs.yahoo.com/info/terms/



Reply by dan_eyer ●November 11, 20032003-11-11

Hi,

We're trying to port some code from a processor
with 64-bit floating-point variables to a VC33,
which only has 32-bit and 40-bit floats, and we
need to maintain as much precision as possible.

Does anyone know if the math library functions
like sqrt support the "long double" (i.e. the
40-bit) floats or not? Or if there are functions
meant just for the long doubles? I know about
MPY_LD and DIV_LD, but we need sqrt and other
functions with the full 40 bit precision.

The prototypes for sqrt etc only seem to list
"double" (32-bit on VC33) parameters, but it
seems strange to have a higher-precision format
and not support doing math in it.

Maybe there are third-party libraries for doing
64-bit floating-point math on the VC33?

Thanks,
Dan Eyer