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 |
|
long double math on VC33?
Started by ●November 11, 2003
Reply by ●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 ●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 ●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! | +-----------+ |