Support Forum / Форум поддержки

Code: Select all

/* Vulgar Arithmetic Library */

' Vulgar fractions - also know as finite precision rational
'                    or slash numbers.

' This library implements a fixed slash number system with 
' 26 bit Numerator and Denominator and a Sign bit. It does not 
' support an Exact bit.

' Numerator and denominator are stored as real and imaginary 
' numbers respectively in complex number variables.

' Mediant rounding is used to approximate vulgar numbers with
' numerators or denominators too large to represent exactly.

' The sign of a vulgar number is stored in the numerator.

' Two vulgar numbers are "almost equal" if the difference
' between them is too small to be repesented exactly.

' Vulgar numbers too large to be represented are rounded to 
' positive or negative infinity, i.e. numerical overflow. 

' For an introduction to slash number systems see:
' Finite Precision Rational Arithmetic: Slash Number Systems
' by D.W. Matula & P. Kornerup 1983 ISSN 0105-8517

' Gordon Charlton 2015


/* Glossary */

' maxint - Largest integer exactly represented in Smart BASIC.


' gcd(n1,n2) - Returns the greatest common denominator of two 
'              integers.


' v.maxsize - When the numerator or denominator of a vulgar 
'             fraction exceeds this, the fraction is rounded.

' v.zero    - Zero, expressed as a vulgar fraction. Vulgar 
'             numbers too small to be represented are rounded
'             to zero.

' v.one     - One, expressed as a vulgar fraction.
' v.inf     - Positive number too large to be represented.
' v.neginf  - Negative number too large to be represented.


' numerator(v)   - Returns the numerator of a vulgar number.
' denominator(v) - Returns the denominator of a vulgar number.


' vround(v,size) - Returns the best approximation to vulgar
'                  number v with values less than size.


' vulgar(n,d)    - Returns n/d as a vulgar number, reduced to
'                  its simplest terms. (i.e. 48/96 => 1/2)
'                  If it is not representable exactly without
'                  exceeding v.maxsize, it is rounded.
'                  n and d should not exceed maxint
'                  (= 9,007,199,254,740,991 ≈ 9 quadrillion).


' vabs(v)        - Returns absolute value of vulgar number v.
' vnegate(v)     - Returns v*-1.
' vreciprocal(v) - Returns 1/v.


' vsign(v) - Returns vulgar 0 if v is 0, vulgar -1 if v is 
'            negative and vulgar 1 if v is positive.


' vinteg(v) - Returns the integer part of v as a vulgar number.
' vfrac(v)  - Returns the fractional part of v as a vulgar
'             number.


' vadd(v1,v2)  - Returns v1+v2, with rounding.
' vsub(v1,v2)  - Returns v1-v2, with rounding.
' vmult(v1,v2) - Returns v1*v2, with rounding.
' vdiv(v1,v2)  - Returns v1/v2, with rounding.


' veeq(v1,v2) - Returns true (1) if v1 is exactly equal to v2.
'               Otherwise returns false (0).

' veq(v1,v2)  - Returns true if v1 is equal or almost equal to
'               v2. Otherwise returns false.

' vgt(v1,v2)  - Returns true if v1 is greater than v2. 
'               Otherwise returns false. v1 is not greater
'               than v2 if v1 and v2 are almost equal.

' vlt(v1,v2)  - Returns true if v1 is less than v2. 
'               Otherwise returns false. v1 is not less
'               than v2 if v1 and v2 are almost equal.

' vge(v1,v2)  - Returns true if v1 is greater than, equal or
'               almost equal to v2. Otherwise returns false.

' vle(v1,v2)  - Returns true if v1 is less than, equal or
'               almost equal to v2. Otherwise returns false.

' vmax(v1,v2) - Returns v1 if it is larger that v2, otherwise 
'               returns v2.

' vmin(v1,v2) - Returns v1 if it is less thsn v2, otherwise 
'               returns v2.


' vsqrt(v)     - Returns the square root of v, with rounding.
' vpower(v1,n) - Returns v1^n, with rounding. n is an integer.

' vrndfareyinc(size) - Returns a vulgar number v, 0<=v<=1 with 
'                      denominator <= size. All valid vulgar 
'                      numbers in the range are returned with 
'                      equal probability.This is equivalent
'                      to randomly selecting an element from a
'                      Farey sequence of order size. The values
'                      returned are unevenly distributed in the
'                      interval 0 to 1, especially for small 
'                      values of size.

' vrndfareyex(size)  - Like vrdfareyinc, except 0<v<1.

' vrndscaleinc(size) - Generates a large random integer n 
'                      scaled to 0<=n<=1 and returns it as a 
'                      vulgar approximation v with denominator 
'                      <= size. The values returned are evenly 
'                      distributed in the interval 0 to 1. All 
'                      valid vulgar numbers in the range are 
'                      not returned with equal probability, 
'                      especially for small values of size.

' vrndscaleex(size)  - Like vrndrealinc, except 0<v<1.

' v$(v) - Returns v expressed as a string. 


/* Library Code */

maxint = 2^53-1

def gcd(a,b)
    a = abs(a)
    b = abs(b)
    while b<>0 
          temp = a%b ! a = b ! b = temp 
    end while
    return a
end def

v.maxsize = integ(sqrt(maxint))
                  
  v.zero = 0+1i 
   v.one = 1+1i
   v.inf = 1+0i
v.neginf = -1+0i 

def numerator(vnum)
    return real(vnum)
end def

def denominator(vnum)
    return imag(vnum)
end def

def vround(vnum,size)       
    n = abs(numerator(vnum)) 
    d = denominator(vnum)    
    if (n<size) and (d<size) then
       lcf = gcd(n,d)
       return (numerator(vnum)/lcf)+(d/lcf*1i)
    else
       thatnum  = 1 ! othernum = 0
       thatden  = 0 ! otherden = 1
       while d<>0
             confrac = integ(n/d)
             thisnum = confrac * thatnum + othernum 
             if thisnum>size then 
       break
             end if
             thisden = confrac * thatden + otherden
             if thisden>size then 
       break
             end if
             temp = n%d ! n = d ! d = temp 
             othernum = thatnum ! thatnum  = thisnum
             otherden = thatden ! thatden  = thisden
       end while
       thatnum *= sign(numerator(vnum))
       return thatnum+thatden*1i
    end if
end def

def vulgar(num,den)
    return vround(num+den*1i,v.maxsize)
end def

def vabs(vnum)
    return abs(numerator(vnum))+1i*denominator(vnum)
end def

def vnegate(vnum)
    return 0-numerator(vnum)+1i*denominator(vnum)
end def

def vreciprocal(vnum)
    if numerator(vnum)<0 then
       sign = -1
    else 
       sign = 1
    end if
    return sign*denominator(vnum)+1i*abs(numerator(vnum))
end def

def vsign(vnum)
    return sign(vnum)+1i
end def

def vinteg(vnum)
    return integ(numerator(vnum)/denominator(vnum))+1i
end def

def vfrac(vnum)
    return vsub(vnum,vinteg(vnum))
end def

def vadd(vnum1,vnum2)
    n1 = numerator(vnum1)
    d1 = denominator(vnum1)
    n2 = numerator(vnum2)
    d2 = denominator(vnum2)
    return vulgar((n1*d2)+(n2*d1),(d1*d2))
end def

def vsub(vnum1,vnum2)
    return vadd(vnum1,vnegate(vnum2))
end def

def vmult(vnum1,vnum2)
    n1 = numerator(vnum1)
    d1 = denominator(vnum1)
    n2 = numerator(vnum2)
    d2 = denominator(vnum2)
    return vulgar(n1*n2,d1*d2)
end def

def vdiv(vnum1,vnum2)
    return vmult(vnum1,vreciprocal(vnum2))
end def

def veeq(vnum1,vnum2)
    if vnum1=vnum2 then
       return 1
    else 
       return 0
    endif
end def

def veq(vnum1,vnum2)
    return veeq(vsub(vnum1,vnum2),v.zero)
end def

def vgt(vnum1,vnum2)
    if sign(vsub(vnum1,vnum2))=1 then
       return 1
    else
       return 0
    endif
end def

def vlt(vnum1,vnum2)
    return vgt(vnum2,vnum1)
    end def

def vge(vnum1,vnum2)
    return 1-vlt(vnum1,vnum2)
end def

def vle(vnum1,vnum2)
    return 1-vgt(vnum1,vnum2)
end def 

def vmax(vnum1,vnum2)
    if vgt(vnum1,vnum2) then 
       return vnum1
    else
       return vnum2
    end if
end def

def vmin(vnum1,vnum2)
    if vlt(vnum1,vnum2)then
       return vnum1
    else
       return vnum2
    end if
end def

def vsqrt(vnum)
    num = numerator(vnum)*denominator(vnum)
    sq = sqrt(num)
    estimate = integ(sq)
    numer = 0
    denom = 1
    if fract(sq)<>0 then
       remainder = num-estimate^2
       if remainder>estimate then 
          estimate += 1 
          remainder = num-estimate^2
       end if 
       est2 = estimate*2
       do
          prevnumer = numer
          prevdenom = denom
          numer = prevdenom*remainder
          denom = prevdenom*est2+prevnumer
          lcf = gcd(numer,denom)
          numer /= lcf
          denom /= lcf
       until (numer>v.maxsize)or(denom>v.maxsize)
       numer = prevnumer
       denom = prevdenom
    end if
    sqrtnumden = vadd(estimate+1i,numer+denom*1i)
    return vdiv(sqrtnumden,denominator(vnum)+1i)
end def

def vpower(vnum,n)
    answer = v.one
    exponent = abs(integ(n))
    while exponent>0
          if odd(exponent) then
             answer = vmult(answer,vnum)
          end if
          exponent = integ(exponent/2)
          vnum = vmult(vnum,vnum)
    end while
    if n>=0 then 
       return answer
    else 
       return vreciprocal(answer)
    end if 
end def

def vrndfareyinc(size)
    do
      denom = rnd(size)+1
      numer = rnd(denom+1)
    until gcd(numer,denom) = 1
    return numer+denom*1i
end def

def vrndfareyex(size)
    do
      answer = vrndfareyinc(size)
    until not(veeq(answer,v.zero) or veeq(answer,v.one))
    return answer
end def

def vrndscaleinc(size)
    largernd = rnd(v.maxsize+1)*(v.maxsize+1)
    largernd += rnd(v.maxsize+1)
    return vround(largernd+maxint*1i,size)
end def

def vrndscaleex(size)
    do
      answer = vrndscaleinc(size)
    until not(veeq(answer,v.zero) or veeq(answer,v.one))
    return answer
end def

def v$(v)
    return str$(numerator(v),"#")&"/"&str$(denominator(v),"#")
end def