это от Мэтью Стэнтон
## 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