如何使用Assembly MIPS {formula provided}計算sin（x）？

Question

我只是堅持瞭如何使用以下公式來計算大會MIPS的sin（x）

http://i.stack.imgur.com/YmxL4.jpg

PLZ如果您有任何想法寫下來..

Answer 1

這個人是從Matthew Stanton

## Matthew Stanton 
## Finds the value of sin(x) 
## Register Use: 
## $t0 value of n 
## $f0 (Series*x^2)/(n(n-1)) 
## $f1 absolute value of (x^2)/(n(n-1)) 
## $f2 holds x^2 
## $f3 holds remainders +or-(x^2)/(n(n-1)) 
## $f4 accuracey 
## $f12 Holds sin(x) 

.text 
.globl main 

main: 
li $t0,3 # Initilize N 
li.s $f4,1.0e-6 # Set Accuracey 
li $v0,4 # syscall for Print String 
la $a0, prompt1 # load address of prompt 
syscall # print the prompt 
li $v0,6 # Reads user number 
syscall 
mul.s $f2,$f0,$f0 # x^2 
mov.s $f12,$f0 # Answer 
for: 
abs.s $f1,$f0 # compares to the non-negative value of the series 
c.lt.s $f1,$f4 # is number < 1.0e-6? 
bc1t endfor 
subu $t1,$t0,1 # (n-1) 
mul $t1,$t1,$t0 # n(n-1) 
mtc1 $t1, $f3 # move n(n-1) to a floating register 
cvt.s.w $f3, $f3 # converts n(n-1) to a float 
div.s $f3,$f2,$f3 # (x^2)/(n(n-1)) 
neg.s $f3,$f3 # -(x^2)/(n(n-1)) 
mul.s $f0,$f0,$f3 # (Series*x^2)/(n(n-1)) 

add.s $f12,$f12,$f0 # Puts answer into $f12 

addu $t0,$t0,2 # Increment n 
b for # Goes to the beggining of the loop 
endfor: 
li $v0,2 # Prints answer in $f12 
syscall 
li $v0,10 # code 10 == exit 
syscall # Halt the program. 




.data 
prompt1: .asciiz "Program will calculate sin(x). Please input a value for x!" 

This code is from

; FILE: Source:sinegen.ASM   REV: 31 --- 16-bit sinetable generator 
; History 
; 31  18th September 1998: 1st version. 
; 

    IFGT 0 

Inspiration for this document and source came from PAC/#amycoders 
who needed good&short sinetable generator. My friend ArtDent coded 
this kind of routine years ago, but unfortunately he didn't backup 
his amiga sources when he went pc. Anyways he still remembered the 
principle well and he pointed me the algorithm to use. This whole 
document and source was written by me (Piru) in 5 hours. 

sine&cosine table generation 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
Lets have a look at sine and cosine graph: 

     pi 2pi 
    _ | | 
    |/|\| | | 
--/-+-\-+-/-- 
    | | |\|/| 
    0 | T 
    | | 
    1/2pi 3/2pi 

     pi 3/2pi 
    _ | | _ 
    |\| | |/| 
--+-\-+-/-+-- 
    | |\_/| | 
    0 |  | 
    1/2pi 2pi 


We notice that sine is phase shifted 90 degrees compared to 
cosine. Also we notice that both sine and cosine are symmetrical 
to 1/2pi and pi, thus can be easily mirrored. So we need to 
calculate only 90 degrees of either sine or cosine and we can 
derive whole table from it and also the other function. 

These are the formulas to calculate sin x and cos x: 

sin x = x - x^3/3! + x^5/5! - x^7/7! + ... 

cos x = 1 - x^2/2! + x^4/4! - x^6/6! + ... 

x is real, 0 <= x <= 1/2pi 


Out of these two the latter (cos x) is easier to calculate. 

You can save space by combining sine and cosine tables. Just 
take last 90 degrees of cosine before cosine table and you 
have sinetable at table - 90 degrees. :) 

So after thinking a while I came up with this pseudocode 
routine that calculates 90 degrees of sine + 360 degrees 
cosine: 

in: table, tablesize (90 degrees * 5) 

quart = tablesize/5 
x = 0; x_add = (1/2 * pi)/quart 
for q = 0 to (quart - 1) 
    fact = 1; d = 0; cosx = 1; powx = 1 
    powx_mul = - (x * x) ; rem this will magically toggle sign 
    repeat 
    powx = powx * powx_mul 
    d++; fact = fact * d 
    d++; fact = fact * d 
    cosx = cosx + powx/fact 
    until d = 12 
    table[quart - q] = cosx   ; rem /¯ 
    table[quart + q] = cosx   ; rem ¯\ 
    table[quart * 3 - q] = -cosx  ; rem  \_ 
    table[quart * 3 + q] = -cosx  ; rem  _/ 
    table[quart * 5 - q] = cosx  ; rem   /¯ 
    x = x + x_add 
endfor 

Then I just coded this in 020+ asm adding fixedpoint math 
and stuff: 

    ENDC 

TESTSINE SET 0 
    IFNE TESTSINE 

Main lea (sine,pc),a0 
    move.l #256,d0 
    bsr sinegen 
    rts 

sine ds.w 256 
cosine ds.w 256*4 
    ENDC 

; 68020+ 16:16 fixedpoint sinetable generator. 
; Coded by Harry "Piru" Sintonen. 
; Not specially optimized as usually this thing is ran only once at 
; init time. 68060 will woe on 64 bit muls & swaps - who cares ;) 

; IN: a0.l=pointer to array of word (will contain 450 degree 16-bit sinetable) 
;  d0.l=wordsper90degrees 
; OUT: d0.l=0 
sinegen 
    movem.l d0-d7/a0-a5,-(sp) 

    move.l #26353589,d1 ; pi/2*65536*256 
    divu.l d0,d1 
    move.l d1,a5 

    add.l d0,d0 
    add.l d0,a0 
    lea 0(a0,d0.l*2),a2 
    lea 0(a0,d0.l*4),a4 
    move.l a0,a1 
    move.l a2,a3 
    addq.l #2,a1  ; these two can be removed 
    addq.l #2,a2  ; really ;) 

    moveq #0,d0  ; x 

    moveq #12,d7 

.oloop move.l d0,d5 
    moveq #1,d1 
    lsr.l #8,d5 
    swap d1  ; 1<<16 = cos x 
    move.l d1,d3 

    mulu.l d5,d4:d5 
    move.w d4,d5 
    moveq #0,d2  ; d 
    swap d5 
    moveq #1,d6  ; factorial 
    neg.l d5  ; change sign of powx 

.iloop muls.l d5,d4:d3 ; calculate x^d 
    move.w d4,d3 
    swap d3 
    move.l d3,d4 

    addq.l #1,d2  ; calculate d! 
    mulu.l d2,d6 
    addq.l #1,d2 
    mulu.l d2,d6 

    divs.l d6,d4 
    add.l d4,d1  ; cos x += x^d/d! 

    cmp.l d7,d2 
    bne.b .iloop 

    lsr.l #1,d1 
    tst.w d1  ; if d1=$8000 then d1=d1-1 ;) 
    dbpl d1,.rule 
.rule 
    move.w d1,(a0)+ 
    move.w d1,-(a1) 
    move.w d1,-(a4) 
    neg.w d1 
    move.w d1,-(a2) 
    move.w d1,(a3)+ 

    add.l a5,d0 
    subq.l #1,(sp)  ; watch out - don't mess with stack:) 
    bne.b .oloop 

    movem.l (sp)+,d0-d7/a0-a5 
    rts 

Answer 2

要計算X * 3需要3次乘法。要計算X * 5，它將花費2次以上的乘法。爲了得到合理的精確度，它會加起來很多乘法。然後是等式的階乘部分 - 爲了合理的精確度，它也是很多的補充。

您無法解決查找表的性能問題;因爲查找表花費的不僅僅是查找「罪惡」表。

基本上，您需要找到適用於計算機的不同公式。

我會忍不住下手CORDIC：http://en.wikipedia.org/wiki/CORDIC

• 布倫丹