Page 1 of 1

Matlib (modified)

Posted: Sat Sep 14, 2013 12:23 pm
by Henko
Modified due to:
- extensions of the standard math functions of Smart Basic -> removed from matlib
- couple of bugs in prior version
- some additional functions
- short explanation about each function

Code: Select all

' IMPORTANT: all functions are based upon OPTION BASE 1
'
' vec_in (n,v,x,y)
' vec_out (n$,n,v,x,y,l,d)
' vec_zero (n,v)
' vec_one (n,v)
' vec_rnd (n,v,mini,maxi)
' vec_copy (n,from,into)
' vec_scal (n,v,s)
' vec_len (n,v)
' vec_norm (n,v)
' vec_plus (n,v1,v2)
' vec_min (n,v1,v2)
' inprod (n,v1,v2)
' mat_in (n,m,mat,x,y)
' mat_out (n$,n,m,mat,x,y,l,d)
' mat_zero (n,m,mat)
' mat_unit (n,mat)
' mat_rnd (n,m,mat,mini,maxi)
' mat_scal (n,m,mat,s)
' mat_copy (n,m,from,into)
' mat_trans (n,m,mat,matt)
' mat_plus (n,m,mata,matb)
' mat_min (n,m,mata,matb)
' mat_mul (n,m,mata,matb,matc)
' mat_inv (n,mata,ainv)
' mat_vec (n,m,mat,vin,vout)
' def ls (n,m,x,y,c)
' def rot2(deg,vin(),vout())
' def rot_x(deg,vin(),vout())
' def rot_y(deg,vin(),vout())
' def rot_z(deg,vin(),vout())
' def rot_xyz(dgx,dgy,dgz,vin(),vout())
' def poly (n,coef,x)
' def surf_n(n,c(,))
' def surf_3(a(),b())
' def mod (a,m)
' def rest (a,b)
' def ln(x)
' def log10(x)
' def pi
' def rad

' input a n-sized vector v() at screen position x,y
'
def vec_in (n,v(),x,y)
dim b$(20)
for i=1 to n
b$(i)="a" & i ! field b$(i) text b$(i) at x,y+30*(i-1) size 80,25
next i
som=0
loop:
for i=1 to n
if field_changed(b$(i)) then
  v(i)=field_text$(b$(i)) ! field b$(i) delete
  draw text n2a$(v(i),6,2) at x,y+30*(i-1)
  som=som+1
end if
next i
if som<n then goto loop
end def

' print a vector at position x,y with number length l and precision d
'
def vec_out (n$,n,v(),x,y,l,d)
ys=y-30
if n$<>"" then
  draw text n$ at x,y
  ys=ys+30
end if
for i=1 to n
  a$=n2a$(i,3,0) & " " & n2a$(v(i),l,d)
  draw text a$ at x-12,ys+25*i
next i
end def

' produce a vector filled with zero's
'
def vec_zero (n,v())
for i=1 to n ! v(i)=0 ! next i
end def

' produce a vector filled with one's
'
def vec_unit (n,v())
for i=1 to n ! v(i)=1 ! next i
end def

' produce a vector filled with random numbers between mini and maxi
'
def vec_rnd (n,v(),mini,maxi)
for i=1 to n ! v(i)=mini+(maxi-mini)*rnd(1) ! next i
end def

' copy vector from() into vector into()
'
def vec_copy (n,from(),into())
for i=1 to n ! into(i)=from(i) ! next i
end def

' multiply a vector with scalar s
'
def vec_scal (n,v(),s)
for i=1 to n ! v(i)=s*v(i) ! next i
end def

' produce the length of a vector
'
def vec_len (n,v()) = sqrt(inprod(n,v,v))

' normalize a vector to unit length
'
def vec_norm (n,v())
vec_scal(n,v,1/vec_len(n,v))
end def

' add vector v2() to vector v1()
'
def vec_plus (n,v1(),v2())
for i=1 to n ! v1(i)=v1(i)+v2(i) ! next i
end def

' subtract vector v2() from vector v1()
'
def vec_min (n,v1(),v2())
for i=1 to n ! v1(i)=v1(i)-v2(i) ! next i
end def

' produce the scalar product of two vectors
'
def inprod (n,v1(),v2())
som=0
for i=1 to n ! som=som+v1(i)*v2(i) ! next i
inprod=som
end def

' input a nxm matrix at position x,y (left upper corner)
'   n is the number of rows, m is the number of columns
'
def mat_in (n,m,mat(,),x,y)
dim b$(100)
b$(1)=0
for i=1 to n
  for j=1 to m
    k=m*(i-1)+j ! b$(k)="     a" & i & j ! tt$=b$(k)
    field tt$ text tt$ at x+100*(j-1),y+30*(i-1) size 80,25
  next j
next i
som=0
loop1:
for i=1 to n
  for j=1 to m
    k=m*(i-1)+j
    if field_changed(b$(k)) then
      mat(i,j)=field_text$(b$(k))
      field b$(k) delete
      draw text n2a$(mat(i,j),6,2) at x+100*(j-1),y+30*(i-1)
      som=som+1
    end if
  next j
next i
if som<n*m then goto loop1
end def

' produce a nxm matrix with zero elements
'
def mat_nul (n,m,mat(,))
for i=1 to n
  for j=1 to m ! mat(i,j)=0 ! next j
next i
end def

' produce a unit matrix with one's in the diagonal and zeros elsewhere
' 
def mat_unit (n,mat(,))
for i=1 to n
  for j=1 to n ! if i=j then mat(i,j)=1 else mat(i,j)=0 ! next j
next i
end def

' produce a matrix with random elements between mini and maxi
'
def mat_rnd (n,m,mat(,),mini,maxi)
for i=1 to n
  for j=1 to m ! mat(i,j)=mini+(maxi-mini)*rnd(1) ! next j
next i
end def

' multiply all elements of a matrix with scalair s
'
def mat_scal (n,m,mat(,),s)
for i=1 to n
  for j=1 to m ! mat(i,j)=s*mat(i,j) ! next j
next i
end def

' copy matrix from() into matrix into()
'
def mat_copy (n,m,from(,),into(,))
for i=1 to n
  for j=1 to m ! into(i,j)=from(i,j) ! next j
next i
end def

' produce the transpose of matrix mat() into matrix matt()
'
def mat_trans (n,m,mat(,),matt(,))
for i=1 to n
  for j=1 to m ! matt(j,i)=mat(i,j) ! next j
next i
end def

' add matb() to mata()
'
def mat_plus (n,m,mata(,),matb(,))
for i=1 to n
  for j=1 to m ! mata( i,j)=mata(i,j)+matb(i,j) ! next j
next i
end def

' subtract matb() from mata()
'
def mat_min (n,m,mata(,),matb(,))
for i=1 to n
  for j=1 to m ! mata( i,j)=mata(i,j)-matb(i,j) ! next j
next i
end def

' produce product of mata() and matb() giving matc()
def mat_mul (n,m,mata(,),matb(,),matc(,))
for i=1 to n
  for j=1 to n
    tot=0
    for k=1 to m ! tot=tot+mata(i,k)*matb(k,j) ! next k
    matc(i,j)=tot
    next j
  next i
end def

' produce the inverse of square matrix a() giving matrix ainv()
'
def mat_inv (nvar,a(,),ainv(,))
dim w(25,50)                      
for i=1 to nvar                 
  for j=1 to nvar ! w(i,j)=a(i,j) ! w(i,j+nvar)=0  ! next j
  w(i,i+nvar)=1
  next i
for piv=1 to nvar
  fac=w(piv,piv)
  for j=piv to piv+nvar ! w(piv,j)=w(piv,j)/fac ! next j
  for i=1 to nvar
    if i<>piv then
      fac=w(i,piv)
      for j=piv to piv+nvar ! w(i,j)=w(i,j)-fac*w(piv,j) ! next j
      endif
    next i
  next piv
for i=1 to nvar
  for j=1 to nvar ! ainv(i,j)=w(i,j+nvar) ! next j
  next i
end def

' print a nxm matrix at position x,y left upper corner
' each element having total length l and d decimals
' n$ is an text, printed to identify the matrix (may be empty)
'
def mat_out (n$,n,m,mat(,),x,y,l,d)
ys=y-30
if n$<>"" then
  draw text n$ at x,y
  ys=ys+30
end if
for i=1 to n
  for j=1 to m
    a$=n2a$(mat(i,j),l,d) ! draw text a$ at x+12*l*(j-1),ys+25*i
  next j
next i
end def

' produce product of a matrix and a vector vin(), giving vout()
' the matrix has size nxm, the vin() has size m, vout() has size n
'
def mat_vec (n,m,mat(,),vin(),vout())
dim v(20)
for i=1 to n
  tot=0
  for j=1 to m ! tot=tot+mat( i,j)*vin(j) ! next j
  v(i)=tot
next i
for i=1 to m ! vout(i)=v(i) ! next i
end def

' this is a n-degree polynomial least squares fit of a number of
' point in 2 dimensional space
' there are n points, with x- and y-coordinates in vectors x() and y()
' m is the degree of the polynomial (take 1 for straight line fit,
' m=2 for parabola, and so on)
' the coefficients of the best fit polynomial are returned in vector c()
' f.i. for m=2 : y = c(1) + c(2)*x + c(3)*x^2
'
def ls (n,m,x(),y(),c())
dim a1(20,20),a2(20,20),a3(20,20),rl(20)
m=m+1
for i=1 to n
  a1(i,1)=1
  for j=2 to m ! a1(i,j)=a1(i,j-1)*x(i) ! next j
next i
mat_trans (n,m,a1,a2)
mat_mul (m,n,a2,a1,a3)
mat_vec (m,n,a2,y,rl)
mat_inv (m,a3,a1)
mat_vec (m,n,a1,rl,c)
end def

' rotate a 2-dimensional vector vin() deg degrees, giving vector vout()
'
def rot2(deg,vin(),vout())
x=vin(1) ! vout(1)=x*cos(deg)-vin(2)*sin(deg)
vout(2)=x*sin(deg)+vin(2)*cos(deg)
end def

' rotate the 3-dim. vector vin() deg degrees about the x-axes -> vout()
' 
def rot_x(deg,vin(),vout())
vout(1)=vin(1)
y=vin(2) ! vout(2)=y*cos(deg)-vin(3)*sin(deg)
vout(3)=y*sin(deg)+vin(3)*cos(deg)
end def

' rotate the 3-dim. vector vin() deg degrees about the y-axes -> vout()
' 
def rot_y(deg,vin(),vout())
vout(2)=vin(2)
x=vin(1) ! vout(1)=x*cos(deg)-vin(3)*sin(deg)
vout(3)=x*sin(deg)+vin(3)*cos(deg)
end def

' rotate the 3-dim. vector vin() deg degrees about the z-axes -> vout()
' 
def rot_z(deg,vin(),vout())
vout(3)=vin(3)
x=vin(1) ! vout(1)=x*cos(deg)-vin(2)*sin(deg)
vout(2)=x*sin(deg)+vin(2)*cos(deg)
end def

' rotate vector vin() about all 3 axes, giving vector vout()
def rot_xyz(dgx,dgy,dgz,vin(),vout())
dim temp(3)
rot_x(dgx,vin,temp)
rot_y(dgy,temp,temp)
rot_z(dgz,temp,vout)
end def

' calculate the value of a polynomial for a give value of x
' n id the degree of the polynomial, hence n+1 coefficients must
' be passed: a0, a1, a2, ..... an, in that order
'
def poly (n,coef(),x)
res=coef(n+1)
for i=n to 1 step-1 ! res=res*x+coef(i) ! next i
poly=res
end def

' calculates the surface of a n-polygon, n>=3
' (divides the polygon in triangles and then uses surf_3 function)
' c() has size nx2 and contains the coordinates of the vertices.
'
def surf_n(n,c(,))
dim cr(20,2),u(2),v(2)
if n<3 then return
x1=c(1,1) ! y1=c(1,2) ! sum=0
for i=1 to n ! cr(i,1)=c(i,1)-x1 ! cr(i,2)=c(i,2)-y1 ! next i
for i=2 to n-1
  u(1)=cr(i,1) ! u(2)=cr(i,2) ! v(1)=cr(i+1,1) ! v(2)=cr(i+1,2)
  sum=sum+surf_3(u,v)
  next i
surf_n=sum
end def

' calculates the surface of a triangle, produced by 2 vectors
'
def surf_3(a(),b())=.5*abs(a(1)*b(2)-a(2)*b(1))

' integer remainder of a/m
'
def mod(a,m)
d=a/m
mod=m*(d-int(d))
end def

' remainder of a/m
'
def rest (a,b)
q=a/b
rest=q-int(q)
end def

' natural logarithm
'
def ln(a)=log(a)

' base 10 logarithm
'
def log10(a)=log(a)/log(10)

' value of pi
'
def pi=3.14159265

' value of 1 radian
'
def rad=180/pi

def n2a$(num,lang,dec)
b$="               "
fac=10^dec
num$=int(fac*num)/fac
tot=lang-len(num$)
if tot<1 then tot=1
a$=substr$(b$,1,tot) & num$
n2a$=a$
end def

def pre_pad$(w,a$)
sp$="               "
tot=w-len(a$)
pre_pad$=substr$(sp$,1,tot) & a$
end def