如何使用Assembly MIPS {公式提供}计算sin（x）？

## 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!\"

; 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