首先,您有一些错误:
(1) 正如 Jester 所提到的,您的“读取浮动”系统调用有一个错误。
(2) 您正在计算factorial(n)
,但泰勒余弦级数使用factorial(2n)
,因此它增加得更快。
(3) 您可能遇到的另一个问题 [将使用正确的阶乘] 是您使用整数数学作为阶乘。当计算到 10 次迭代时,阶乘将溢出一个 32 位整数。对于 10 次迭代,我们最终得到 [至少]factorial(18)
约为 6.4e15(即不适合32 整数)。对于整数数学,只要值包含 32 位,阶乘就会“振荡”,而不是稳定增加。
我采取了稍微不同的方法。
我没有像您那样预先计算所有内容,而是创建了一个使用循环的解决方案。这可能会或可能不会帮助您调试自己的实现,但如果我不建议重构,我会失职。
下面是循环实现的 asm 代码。还有一个 C 版本,我用来先调试算法。还有一个带有一些调试语句的 asm 版本。
理由:
您正在尝试cos
使用泰勒级数求和/展开 [至 10 次迭代] 进行计算。毕竟,该公式确实使用了“求和”运算符。
您正在预先计算所有x^n
项和所有阶乘项,并将它们放在不同的寄存器中。这很好,因为 MIPS 有 32 个 FP 寄存器。但是,对于双精度,我们必须使用两个寄存器,所以这意味着只有 16 个数字。在某种程度上,您所做的相当于编译器执行“循环展开”。
另一个问题是很难跟踪所有这些寄存器。什么持有什么价值。
而且,假设问题是使用 20 次迭代而不是 10 次。我们可能会用完寄存器。在实践中,这对于其他级数展开可能是必要的,因为我们可能不会很快收敛。
所以,我建议使用循环。每个幂项很容易从前一个算出。同样,对于每个阶乘。
使用循环方法的另一个优点是,我们可以监控术语值 [ ](越来越小),而不是使用固定次数的迭代,cur
以及它是否变化小于一定量 [或小于量](例如对于双精度,1e-14)我们可以停止,因为我们的值不会比浮点格式[和硬件]提供给我们的精度更好。
注意:这未在下面显示,但很容易实现。
这是asm版本。注释中使用的变量名称是指 C 代码中使用的变量名称 [随后]。
# A program to calculate cos(x) using taylor series
.data
pr1: .asciiz "Enter Float x: "
sym_fnc: .asciiz "cos(x): "
nl: .asciiz "\n"
.text
.globl main
main:
li $s7,0 # clear debug flag #+
# prompt user for x value
li $v0,4 # print string
la $a0,pr1 # prompt user for x
syscall
# read user's x value
li $v0,6 # read float
syscall
jal qcos
la $a0,sym_fnc # string
jal prtflt
li $v0,10
syscall
# qcos -- calculate cosine
#
# RETURNS:
# f0 -- cos(x)
#
# arguments:
# f0 -- x value
#
# registers:
# f2 -- x value
# f4 -- sum
# f6 -- xpow (x^n)
# f8 -- n2m1
# f10 -- factorial (nfac)
# f12 -- RESERVED (used in pflt)
# f14 -- current term
# f16 -- x^2
#
# f18 -- a one value
#
# t0 -- zero value
# t1 -- one value
#
# t6 -- negation flag
# t7 -- iteration count
qcos:
move $s0,$ra # save return address
mov.s $f2,$f0 # save x value
mul.s $f16,$f2,$f2 # get x^2
li $t0,0 # get a zero
li $t1,1 # get a one
li $t6,1 # start with positive term
# xpow = 1
mtc1 $t1,$f6 # xpow = 1
cvt.s.w $f6,$f6 # convert to float
# n2m1 = 0
mtc1 $t0,$f8 # n2m1 = 0
cvt.s.w $f8,$f8 # convert to float
# nfac = 1
mtc1 $t1,$f10 # nfac = 1
cvt.s.w $f10,$f10 # convert to float
# get a one value
mtc1 $t1,$f18 # onetmp = 1
cvt.s.w $f18,$f18 # convert to float
# zero the sum
mtc1 $t0,$f4 # sum = 0
cvt.s.w $f4,$f4 # convert to float
li $t7,10 # set number of iterations
cosloop:
div.s $f14,$f6,$f10 # cur = xpow / nfac
# apply the term to the sum
bgtz $t6,cospos # do positive? yes, fly
sub.s $f4,$f4,$f14 # subtract the term
b cosneg
cospos:
add.s $f4,$f4,$f14 # add the term
cosneg:
subi $t7,$t7,1 # bump down iteration count
blez $t7,cosdone # are we done? if yes, fly
# now calculate intermediate values for _next_ term
# get _next_ power term
mul.s $f6,$f6,$f16 # xpow *= x2
# go from factorial(2n) to factorial(2n+1)
add.s $f8,$f8,$f18 # n2m1 += 1
mul.s $f10,$f10,$f8 # nfac *= n2m1
# go from factorial(2n+1) to factorial(2n+1+1)
add.s $f8,$f8,$f18 # n2m1 += 1
mul.s $f10,$f10,$f8 # nfac *= n2m1
neg $t6,$t6 # flip sign for next time
j cosloop
cosdone:
mov.s $f0,$f4 # set return value
move $ra,$s0 # restore return address
jr $ra
# dbgflt -- debug print float number
dbgflt:
bnez $s7,prtflt
jr $ra
# dbgnum -- debug print int number
dbgnum:
beqz $s7,dbgnumdone
li $v0,1
syscall
dbgnumdone:
jr $ra
# dbgprt -- debug print float number
dbgprt:
beqz $s7,dbgprtdone
li $v0,4
syscall
dbgprtdone:
jr $ra
# prtflt -- print float number
#
# arguments:
# a0 -- prefix string (symbol name)
# f12 -- number to print
prtflt:
li $v0,4 # syscall: print string
syscall
li $v0,2 # print float
syscall
li $v0,4 # syscall: print string
la $a0,nl # print newline
syscall
jr $ra
这是C代码[它也有代码sin(x)
]:
// mipsqsin/mipstaylor -- fast sine/cosine calculation
#include <stdio.h>
#include <math.h>
#define ITERMAX 10
// qcos -- calculate cosine
double
qcos(double x)
{
int iteridx;
double x2;
double cur;
int neg;
double xpow;
double n2m1;
double nfac;
double sum;
// square of x
x2 = x * x;
// values for initial terms where n==0:
xpow = 1.0;
n2m1 = 0.0;
nfac = 1.0;
neg = 1;
sum = 0.0;
iteridx = 0;
// NOTES:
// (1) with the setup above, we can just use the loop without any special
// casing
while (1) {
// calculate current value
cur = xpow / nfac;
// apply it to sum
if (neg < 0)
sum -= cur;
else
sum += cur;
// bug out when done
if (++iteridx >= ITERMAX)
break;
// now calculate intermediate values for _next_ sum term
// get _next_ power term
xpow *= x2;
// go from factorial(2n) to factorial(2n+1)
n2m1 += 1.0;
nfac *= n2m1;
// now get factorial(2n+1+1)
n2m1 += 1.0;
nfac *= n2m1;
// flip sign
neg = -neg;
}
return sum;
}
// qsin -- calculate sine
double
qsin(double x)
{
int iteridx;
double x2;
double cur;
int neg;
double xpow;
double n2m1;
double nfac;
double sum;
// square of x
x2 = x * x;
// values for initial terms where n==0:
xpow = x;
n2m1 = 1.0;
nfac = 1.0;
neg = 1;
sum = 0.0;
iteridx = 0;
// NOTES:
// (1) with the setup above, we can just use the loop without any special
// casing
while (1) {
// calculate current value
cur = xpow / nfac;
// apply it to sum
if (neg < 0)
sum -= cur;
else
sum += cur;
// bug out when done
if (++iteridx >= ITERMAX)
break;
// now calculate intermediate values for _next_ sum term
// get _next_ power term
xpow *= x2;
// go from factorial(2n+1) to factorial(2n+1+1)
n2m1 += 1.0;
nfac *= n2m1;
// now get factorial(2n+1+1+1)
n2m1 += 1.0;
nfac *= n2m1;
// flip sign
neg = -neg;
}
return sum;
}
// testfnc -- test function
void
testfnc(int typ,const char *sym)
{
double (*efnc)(double);
double (*qfnc)(double);
double vale;
double valq;
double x;
double dif;
int iter;
switch (typ) {
case 0:
efnc = cos;
qfnc = qcos;
break;
case 1:
efnc = sin;
qfnc = qsin;
break;
default:
efnc = NULL;
qfnc = NULL;
break;
}
iter = 0;
for (x = 0.0; x <= M_PI_2; x += 0.001, ++iter) {
vale = efnc(x);
valq = qfnc(x);
dif = vale - valq;
dif = fabs(dif);
printf("%s: %d x=%.15f e=%.15f q=%.15f dif=%.15f %s\n",
sym,iter,x,vale,valq,dif,(dif < 1e-14) ? "PASS" : "FAIL");
}
}
// main -- main program
int
main(int argc,char **argv)
{
testfnc(0,"cos");
testfnc(1,"sin");
return 0;
}
这是我使用的带有调试语句的 asm 版本。这是“printf 调试”,就像在 C 中一样。
或者,模拟器mars
和spim
模拟器都具有内置的 GUI 调试功能。您可以通过单击一个按钮来单步执行代码。您可以看到所有寄存器值的实时显示。
注意:就个人而言,我更喜欢mars
. 这适用,mars
但我不知道是否spim
支持.eqv
伪操作。
# A program to calculate cos(x) using taylor series
.data
pr1: .asciiz "Enter Float x: "
sym_fnc: .asciiz "cos(x): "
nl: .asciiz "\n"
#+ddef f2,x
.eqv fpr_x $f2 #+
sym_x: .asciiz "x[f2]: " #+
#+
#+ddef f16,x2
.eqv fpr_x2 $f16 #+
sym_x2: .asciiz "x2[f16]: " #+
#+
#+ddef f4,sum
.eqv fpr_sum $f4 #+
sym_sum: .asciiz "sum[f4]: " #+
#+
#+ddef f6,xpow
.eqv fpr_xpow $f6 #+
sym_xpow: .asciiz "xpow[f6]: " #+
#+
#+ddef f8,n2m1
.eqv fpr_n2m1 $f8 #+
sym_n2m1: .asciiz "n2m1[f8]: " #+
#+
#+ddef f10,nfac
.eqv fpr_nfac $f10 #+
sym_nfac: .asciiz "nfac[f10]: " #+
#+
#+ddef f14,cur
.eqv fpr_cur $f14 #+
sym_cur: .asciiz "cur[f14]: " #+
#+
.text
.globl main
main:
#+dask
la $a0,dbgask # prompt user #+
li $v0,4 # print string #+
syscall #+
# get debug flag #+
li $v0,5 #+
syscall #+
move $s7,$v0 #+
.data #+
dbgask: .asciiz "Debug (0/1) ? " #+
.text #+
#+
# prompt user for x value
li $v0,4 # print string
la $a0,pr1 # prompt user for x
syscall
# read user's x value
li $v0,6 # read float
syscall
jal qcos
la $a0,sym_fnc # string
jal prtflt
li $v0,10
syscall
# qcos -- calculate cosine
#
# RETURNS:
# f0 -- cos(x)
#
# arguments:
# f0 -- x value
#
# registers:
# f2 -- x value
# f4 -- sum
# f6 -- xpow (x^n)
# f8 -- n2m1
# f10 -- factorial (nfac)
# f12 -- RESERVED (used in pflt)
# f14 -- current term
# f16 -- x^2
#
# f18 -- a one value
#
# t0 -- zero value
# t1 -- one value
#
# t6 -- negation flag
# t7 -- iteration count
qcos:
move $s0,$ra # save return address
mov.s $f2,$f0 # save x value
#+dflt x
la $a0,sym_x #+
mov.s $f12,fpr_x #+
jal dbgflt #+
#+
mul.s $f16,$f2,$f2 # get x^2
#+dflt x2
la $a0,sym_x2 #+
mov.s $f12,fpr_x2 #+
jal dbgflt #+
#+
li $t0,0 # get a zero
li $t1,1 # get a one
li $t6,1 # start with positive term
# xpow = 1
mtc1 $t1,$f6 # xpow = 1
cvt.s.w $f6,$f6 # convert to float
# n2m1 = 0
mtc1 $t0,$f8 # n2m1 = 0
cvt.s.w $f8,$f8 # convert to float
# nfac = 1
mtc1 $t1,$f10 # nfac = 1
cvt.s.w $f10,$f10 # convert to float
# get a one value
mtc1 $t1,$f18 # onetmp = 1
cvt.s.w $f18,$f18 # convert to float
# zero the sum
mtc1 $t0,$f4 # sum = 0
cvt.s.w $f4,$f4 # convert to float
li $t7,10 # set number of iterations
cosloop:
#+dprt "cosloop: LOOP iter="
la $a0,dprt_1 #+
jal dbgprt #+
.data #+
dprt_1: .asciiz "cosloop: LOOP iter=" #+
.text #+
#+
#+dnum $t7
move $a0,$t7 #+
jal dbgnum #+
#+
#+dprt "\n"
la $a0,dprt_2 #+
jal dbgprt #+
.data #+
dprt_2: .asciiz "\n" #+
.text #+
#+
#+dflt xpow
la $a0,sym_xpow #+
mov.s $f12,fpr_xpow #+
jal dbgflt #+
#+
#+dflt nfac
la $a0,sym_nfac #+
mov.s $f12,fpr_nfac #+
jal dbgflt #+
#+
div.s $f14,$f6,$f10 # cur = xpow / nfac
#+dflt cur
la $a0,sym_cur #+
mov.s $f12,fpr_cur #+
jal dbgflt #+
#+
# apply the term to the sum
bgtz $t6,cospos # do positive? yes, fly
#+dprt "costerm: NEG\n"
la $a0,dprt_3 #+
jal dbgprt #+
.data #+
dprt_3: .asciiz "costerm: NEG\n" #+
.text #+
#+
sub.s $f4,$f4,$f14 # subtract the term
b cosneg
cospos:
#+dprt "costerm: POS\n"
la $a0,dprt_4 #+
jal dbgprt #+
.data #+
dprt_4: .asciiz "costerm: POS\n" #+
.text #+
#+
add.s $f4,$f4,$f14 # add the term
cosneg:
#+dflt sum
la $a0,sym_sum #+
mov.s $f12,fpr_sum #+
jal dbgflt #+
#+
subi $t7,$t7,1 # bump down iteration count
blez $t7,cosdone # are we done? if yes, fly
# now calculate intermediate values for _next_ term
#+dprt "cosloop: CALC\n"
la $a0,dprt_5 #+
jal dbgprt #+
.data #+
dprt_5: .asciiz "cosloop: CALC\n" #+
.text #+
#+
# get _next_ power term
mul.s $f6,$f6,$f16 # xpow *= x2
# go from factorial(2n) to factorial(2n+1)
add.s $f8,$f8,$f18 # n2m1 += 1
#+dflt n2m1
la $a0,sym_n2m1 #+
mov.s $f12,fpr_n2m1 #+
jal dbgflt #+
#+
mul.s $f10,$f10,$f8 # nfac *= n2m1
#+dflt nfac
la $a0,sym_nfac #+
mov.s $f12,fpr_nfac #+
jal dbgflt #+
#+
# go from factorial(2n+1) to factorial(2n+1+1)
add.s $f8,$f8,$f18 # n2m1 += 1
#+dflt n2m1
la $a0,sym_n2m1 #+
mov.s $f12,fpr_n2m1 #+
jal dbgflt #+
#+
mul.s $f10,$f10,$f8 # nfac *= n2m1
#+dflt nfac
la $a0,sym_nfac #+
mov.s $f12,fpr_nfac #+
jal dbgflt #+
#+
neg $t6,$t6 # flip sign for next time
j cosloop
cosdone:
mov.s $f0,$f4 # set return value
move $ra,$s0 # restore return address
jr $ra
# dbgflt -- debug print float number
dbgflt:
bnez $s7,prtflt
jr $ra
# dbgnum -- debug print int number
dbgnum:
beqz $s7,dbgnumdone
li $v0,1
syscall
dbgnumdone:
jr $ra
# dbgprt -- debug print float number
dbgprt:
beqz $s7,dbgprtdone
li $v0,4
syscall
dbgprtdone:
jr $ra
# prtflt -- print float number
#
# arguments:
# a0 -- prefix string (symbol name)
# f12 -- number to print
prtflt:
li $v0,4 # syscall: print string
syscall
li $v0,2 # print float
syscall
li $v0,4 # syscall: print string
la $a0,nl # print newline
syscall
jr $ra
这是单个值的调试输出日志:
Debug (0/1) ? 1
Enter Float x: 0.123
x[f2]: 0.123
x2[f16]: 0.015129001
cosloop: LOOP iter=10
xpow[f6]: 1.0
nfac[f10]: 1.0
cur[f14]: 1.0
costerm: POS
sum[f4]: 1.0
cosloop: CALC
n2m1[f8]: 1.0
nfac[f10]: 1.0
n2m1[f8]: 2.0
nfac[f10]: 2.0
cosloop: LOOP iter=9
xpow[f6]: 0.015129001
nfac[f10]: 2.0
cur[f14]: 0.0075645004
costerm: NEG
sum[f4]: 0.9924355
cosloop: CALC
n2m1[f8]: 3.0
nfac[f10]: 6.0
n2m1[f8]: 4.0
nfac[f10]: 24.0
cosloop: LOOP iter=8
xpow[f6]: 2.2888667E-4
nfac[f10]: 24.0
cur[f14]: 9.536944E-6
costerm: POS
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 5.0
nfac[f10]: 120.0
n2m1[f8]: 6.0
nfac[f10]: 720.0
cosloop: LOOP iter=7
xpow[f6]: 3.4628265E-6
nfac[f10]: 720.0
cur[f14]: 4.8094813E-9
costerm: NEG
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 7.0
nfac[f10]: 5040.0
n2m1[f8]: 8.0
nfac[f10]: 40320.0
cosloop: LOOP iter=6
xpow[f6]: 5.2389105E-8
nfac[f10]: 40320.0
cur[f14]: 1.2993329E-12
costerm: POS
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 9.0
nfac[f10]: 362880.0
n2m1[f8]: 10.0
nfac[f10]: 3628800.0
cosloop: LOOP iter=5
xpow[f6]: 7.925948E-10
nfac[f10]: 3628800.0
cur[f14]: 2.1841789E-16
costerm: NEG
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 11.0
nfac[f10]: 3.99168E7
n2m1[f8]: 12.0
nfac[f10]: 4.790016E8
cosloop: LOOP iter=4
xpow[f6]: 1.1991168E-11
nfac[f10]: 4.790016E8
cur[f14]: 2.503367E-20
costerm: POS
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 13.0
nfac[f10]: 6.2270208E9
n2m1[f8]: 14.0
nfac[f10]: 8.7178289E10
cosloop: LOOP iter=3
xpow[f6]: 1.8141439E-13
nfac[f10]: 8.7178289E10
cur[f14]: 2.0809583E-24
costerm: NEG
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 15.0
nfac[f10]: 1.30767428E12
n2m1[f8]: 16.0
nfac[f10]: 2.09227885E13
cosloop: LOOP iter=2
xpow[f6]: 2.7446184E-15
nfac[f10]: 2.09227885E13
cur[f14]: 1.3117842E-28
costerm: POS
sum[f4]: 0.99244505
cosloop: CALC
n2m1[f8]: 17.0
nfac[f10]: 3.55687415E14
n2m1[f8]: 18.0
nfac[f10]: 6.4023735E15
cosloop: LOOP iter=1
xpow[f6]: 4.1523335E-17
nfac[f10]: 6.4023735E15
cur[f14]: 6.485616E-33
costerm: NEG
sum[f4]: 0.99244505
cos(x): 0.99244505
-- program is finished running --