home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
World of Shareware - Software Farm 2
/
wosw_2.zip
/
wosw_2
/
QBAS
/
UBAS830.ZIP
/
DETP.ASM
< prev
next >
Wrap
Assembly Source File
|
1990-01-28
|
6KB
|
382 lines
;detp.asm
;determinant of matrix det P
;1990 by Yuji KIDA
comment &
usage :
dim Detp%(500):bload "detp",Detp%()
N=(P-1)\2-2:NN=(N+7)\8*8
dim A%(N-1,NN-1),W%(N-1,NN-1) ;column size = multiple of 8
call Detp%(W%(0,0),N,PP,DP)
original BASIC program is :
640 *DETP(P):' assume every elements >=0 and <p
650 local I%,J%,K%,DP,A
660 DP=1:' determinant mod p
670 for J%=0 to N-2
680 for I%=J% to N-1:' search non_zero
690 if W(I%,J%) then cancel for:goto 720
700 next
710 cancel for:DP=0:goto 840
720 if I%<>J% then neg DP
730 :swap block W(I%,J%..N-1),W(J%,J%..N-1)
740 A=W(J%,J%)
750 DP=DP*A@P
760 Inv=modinv(W(J%,J%),P)
770 for I%=J%+1 to N-1
780 C=W(I%,J%)
790 if C then C=C*Inv@P
800 :for K%=J%+1 to N-1:W(I%,K%)=(W(I%,K%)-C*W(J%,K%))@P:next
810 next
820 next
830 DP=DP*W(N-1,N-1)@P
840 return(DP)
&
INCLUDE UB.MAC
jmp start0
even
base_segment dw ?
matrixsize dw ?
modulus dw ?
segstep dw ?
var_J dw ?
invword dw ?
multiplier dw ?
counter dw ?
detoffset dw ?
detsegment dw ?
det dw ?
modinvword:
; inp : ax (GCD with 'modulus' must be 1)
; out : ax = 1/ax mod 'modulus'
; destroy : nothing
; cmp ax,1
; je modinvwordret
push bx
push dx
push si
push di
cmp ax,bx ;bx=modulus
jb modinvw10
xor dx,dx
div bx
mov ax,dx
cmp ax,1
je modinvwout
modinvw10:
xchg ax,bx ;##
xor si,si ;coef1
mov di,1 ;coef2
modinvw20:
xor dx,dx
div bx
push dx ;remainder
push di ;coef2
mul di
div cs:[modulus]
mov ax,si
sub ax,dx
jae modinvw30
add ax,cs:[modulus]
modinvw30:
mov di,ax ;new coef2=coef1-quotient*coef2
pop si ;new coef1=old coef2
mov ax,bx
pop bx
cmp bx,1
jne modinvw20 ;GCD must 1 otherwise endlessloop
mov ax,di
modinvwout:
pop di
pop si
pop dx
pop bx
modinvwordret:
ret
getvalue:
lodsw
cmp ax,1
jne getvalueilg
lodsw
clc
ret
getvalueilg:
stc
ret
detpilg:
mov bx,AR1
mov word ptr cs:[bx],8001h ;error
jmp detpret
start0:
mov ax,cs
mov ds,ax
mov es,ax
lds_si v2
call getvalue
jc detpilg
mov cs:[matrixsize],ax
lds_si v3
call getvalue
jc detpilg
mov cs:[modulus],ax
lds_si v4
mov cs:[detoffset],si
mov cs:[detsegment],ds
mov ax,cs
mov ds,ax
mov es,ax
mov [det],1 ;set initial det
mov bx,V1+2
mov ax,[bx] ;segment of W()
mov ds,ax
xor si,si
mov ax,[si] ;ax=dimension
cmp ax,2
jne detpilg
mov si,8
mov ax,[si] ;size of 2 dim -1
inc ax
xor dx,dx
mov cx,8
div cx
or dx,dx
jnz detpilg
mov cs:[segstep],ax
mov ax,ds
add ax,arrayheadseg
mov ds,ax
mov es,ax
mov cs:[base_segment],ax
mov bx,cs:[modulus]
mov cx,cs:[matrixsize]
detp10:
push cx
xor si,si
mov cx,cs:[matrixsize]
detp20:
mov ax,[si]
test ah,80h
jz detp30
and ah,7fh
neg ax
detp30:
add ax,7fffh
xor dx,dx
div bx ;modulus
mov [si],dx
add si,2
loop detp20
mov ax,ds
add ax,cs:[segstep]
mov ds,ax
mov es,ax
pop cx
loop detp10
; 670 for J%=0 to N-2
mov ax,cs:[base_segment]
mov ds,ax
xor si,si ;ds:si = adr of W(0,J)
mov ax,cs:[matrixsize]
dec ax
jnz detp105
jmp detp200
detp105:
mov cs:[counter],ax
mov cs:[var_J],0
;main loop
detp110:
mov ax,ds
mov es,ax
; 680 for I%=J% to N-1:' search non_zero
; 690 if W(I%,J%) then cancel for:goto 720
; 700 next
; 710 cancel for:DP=0:goto 840
mov cx,cs:[matrixsize]
sub cx,cs:[var_J]
xor ax,ax
detp120:
cmp es:[si],ax
jne detp130
mov dx,es
add dx,cs:[segstep]
mov es,dx
loop detp120
jmp detpequ0 ;det=0
; 720 if I%<>J% then neg DP
detp130:
mov ax,ds
mov dx,es
cmp ax,dx
je detp150
mov ax,bx ;neg det
sub ax,cs:[det] ;
mov cs:[det],ax ;
; 730 :swap block W(I%,J%..N-1),W(J%,J%..N-1)
mov cx,cs:[matrixsize] ;exchange rows
sub cx,cs:[var_J]
mov di,si
detp140:
mov ax,es:[di]
xchg [di],ax
mov es:[di],ax
add di,2
loop detp140
mov ax,ds
mov es,ax
; 740 A=W(J%,J%)
; 750 DP=DP*A@P
detp150:
mov ax,[si] ;calc new det
mul cs:[det]
div bx ;##
mov cs:[det],dx
; 770 for I%=J%+1 to N-1
mov cx,cs:[matrixsize]
sub cx,cs:[var_J]
dec cx
jcxz detp190
; 760 Inv=modinv(W(J%,J%),P)
mov ax,[si]
cmp ax,1
je detp155
call modinvword
detp155:
mov cs:[invword],ax
; 780 C=W(I%,J%)
; 790 if C then C=C*Inv@P
detp160:
mov ax,es
add ax,cs:[segstep]
mov es,ax
;740
mov ax,es:[si]
or ax,ax
jz detp185
; 800 :for K%=J%+1 to N-1:W(I%,K%)=(W(I%,K%)-C*W(J%,K%))@P
; :next
; 810 next
push cx
mul cs:[invword]
div bx ;##
mov ax,dx
or ax,ax
jz detp170
neg ax
add ax,bx ;##
detp170:
mov cs:[multiplier],ax
mov cx,cs:[matrixsize]
sub cx,cs:[var_J]
mov di,si
detp180:
mov ax,[di]
mul cs:[multiplier]
div bx ;##
mov ax,dx
xor dx,dx
add ax,es:[di]
adc dx,0
div bx ;##
mov es:[di],dx
add di,2
loop detp180
pop cx
detp185:
loop detp160
; 820 next
detp190:
inc cs:[var_J]
mov ax,ds
add ax,cs:[segstep]
mov ds,ax
add si,2
dec cs:[counter]
jz detp200
jmp detp110
detpequ0:
mov cs:[det],0
; 830 DP=DP*W(N-1,N-1)@P
detp200:
mov ax,[si] ;final element
mul cs:[det]
div bx
mov ax,dx
lds si,dword ptr cs:[detoffset]
or ax,ax
jz detp210
mov word ptr [si],1
add si,2
detp210:
mov [si],ax
; 840 return(DP)
detpret:
return