Ошибка в программе вычисления интеграла, метод трапеций Опции

[Hey]

Приветствую,

пишу программу вычисления интеграла функции 1/ln x методом трапеций, для начала взял n=20 (затем буду увеличивать). Файл .exe создан, но выдает ошибку и вылетает. При просмотре в дебаггере вижу, что еще на этапе вычисления шага в регистре ST0 оказывается отрицательная величина, еще через пару шагов вылетает. В чем тут дело?

Заранее спасибо.

.486p
.model small 
.stack 100h
.data 
a dw 2
b dw 5
x dq 2
n dw 20
h dq ?
y dq ?

.code 
main proc
mov ax, @data
mov ds, ax
finit
fild b
fisub a
fidiv n
fst h  ;вычисляем шаг h 


fld1
fild b
fyl2x
fldln2
fmul
fld1
fdiv
fADD y
fild b
fdivr st(1), st ;вычисляем (ln a + ln b)/2

cycl:
fadd y
fst y 

fild x 
fadd h
fst x ;производим приращение х и каждый раз сохраняем
ficom b
je mult_  ;пока не равно b(=5),продолжаем 

fld x
fyl2x
fldln2
fmul
fld1 
fdiv ;
loop cycl ;

mult_: 
fld y ;
fld h ; 
fmul   ;когда х достиг значения 5, перемножаем на h

exit:   
mov ax, 4c00h ;
int 21h ;
main endp
end main

[IUnknown]

Hey, я опять не понял, зачем тебе 2 цикла, в каждом из которых ты вычисляешь значения интегралов? Давай запишем алгоритм который тебе нужен, в псевдокоде. Итак, чтобы определить, при каком N значение интеграла найдено с точностью Eps, тебе достаточно:

   n <- 0
   y <- 0.0;
   repeat
      prev <- y;
      n <- n + DeltaN;
      y <- Integral(a, b, n);
   until abs(y - z) < eps;
   print "finally Int = ", y, " when n = ", n

Всё, никаких двойных циклов, все прекрасно делается одним. А теперь смотри, как это записывается на ассемблере (я заменил переменную z на prev):

main proc
	mov ax, @data
	mov ds, ax
	
outer_:
	fld y
	fstp prev
	; fldz
	; fstp y

	; вычисляем шаг для текущего N
	Сalc_h
	
	iCalcF b	; 1 / ln b
	fstp y
	iCalcF a	; 1 / ln a
	fadd y	; Складываем
	fdiv x	; Делим сумму на 2
	fstp y	; Сохр. начальное приближение
	
	; внутренний цикл - вычисляем интеграл
	fld x		; X в стек сопроцессора
inner_loop:
	fadd h	; X <- X + h
	CalcF	; F(X)
	fadd y
	fstp y	; Сохраняем измененное значение интеграла
	loop inner_loop
	
	fld y
	fmul h	; нашли окончательное значение при текущем N
	fst y		; сохранили его назад в Y для последующего сравнения с prev
	
	; дублируем значение с вершины стека
	; fld st
	
	; выводим значение интеграла и N
	call outfloat
	call OutInt ; выводим актуальное n
	
	fld y
	fsub prev
	fabs
	fcomp eps
	fstsw ax
	sahf
	jp uncomp
	jz m1
	jc exit
m1:
	jmp outer_
uncomp:
	mov dx, offset mes
	mov ah, 09
	int 21h

exit:
	mov ax, 4c00h ;
	int 21h ;
main endp

Проще, правда, чем делать дважды одно и то же?

Теперь насчет

Вывод n также не работает.

Процедуру OutInt проверять не пробовал?

Вообще-то, она должна выглядеть так (вообще-то, жестко задавать в программе, что процедура выводит именно N - это бред, лучше перед ее вызовом занести ax <- N):

outint  proc    near
        push    cx
        push    dx
        push    bx
        push    ax
	
	mov ax, n  ; <--- !!! А этого лучше не делать !!!

; Проверяем число на знак.
        test    ax, ax
        jns     short @oi1
; Если оно отрицательное, выведем минус и оставим его модуль.
        mov     ah, 02h
        mov     dl, '-'
        int     21h
        pop     ax
        push    ax
        neg     ax
; Количество цифр будем держать в CX.
@oi1:   xor     cx, cx
        mov     bx, 10
@oi2:   xor     dx, dx
        div     bx
; Делим число на десять. В остатке получается последняя цифра.
; Сразу выводить её нельзя, поэтому сохраним её в стэке.
        push    dx
        inc     cx
; А с частным повторяем то же самое, отделяя от него очередную
; цифру справа, пока не останется ноль, что значит, что дальше
; слева только нули.
        test    ax, ax
        jnz     short @oi2
; Теперь приступим к выводу.
        mov     ah, 02h
@oi3:   pop     dx
; Извлекаем очередную цифру, переводим её в символ и выводим.
        add     dl, 30h
        int     21h
; Повторим ровно столько раз, сколько цифр насчитали.
        loop    @oi3
	
	; переводим строку после выведенного числа
	mov ah, 9
	mov dx, offset NL
	int 21h

        pop     ax
        pop     bx
        pop     dx
        pop     cx
        ret
outint  endp

С этими изменениями программа прекрасно отрабатывает:

F:\Asm30>trapece
2.5912
20
2.5898
40
2.5896
60
2.5895
80

F:\Asm30>