Added random numbers and tidied up code and glossary. This will be the final update until I understand Gosper's Algorithm, which, considering its complexity, may be never. When and if I do understand it, this code will be superseded by a complete rewrite.
Gosper: http://perl.plover.com/classes/cftalk/INFO/gosper.txt
update 1
The library has been extended to include more functionality present in Smart BASIC's built in maths. It currently lacks random, trigonometric and logarithmic functions and the exponent in the exponentiation function (vpower) Is an integer, not a vulgar fraction.
Test routines at the end of the library have been removed.
end updates
... or Complex Abuse.
In this collection of functions I reappropriate the Complex Number data type to implement a minimal set of Rational Number operators (also known as Vulgar Fractions or Fixed Slash Numbers.)
Rational Number representations are interesting not only because they use Primary School mathematics to represent values like 1/3 exactly rather than as an approximation, but because they are a significant part of Number Theory which is the basis for modern encryption techniques. For the mathematically inclined this is a deep and fascinating subject.
Some notes about precision and rounding. The number of rational numbers representable in a 51 bit implementation such as this (25 bits for the numerator, 25 bits for the denominator and 1 sign bit) is not 2^51-1 as rational numbers are reduced to their simplest form - the numerator and denominator are always mutually prime. For large ranges of numbers approximately 61% of number pairs are mutually prime, so approximately 1 quadrillion unique positive values and the same number of negative values can be represented.
Rational numbers are not distributed evenly along the number line. Half of the representable numbers lie in the range -1 to 1. The distance between these numbers is variable - in the worst case there is a gap of approx 3*10^-8 (three ten millionths) between two consecutive numbers, so the minumum precision in this range is seven decimal places. The other half of the representable numbers lie outside the range -1 to 0 and grow increasingly imprecise as the integer value increases. Beyond ±16,777,216 only integers can be represented.
As Knuth observes without proof in Semi-Numerical Algorithms, rational number approximations employing a mediant rounding scheme (as in this implementation) have the unusual property that rounding errors tend to cancel out rather than accumulate as is the case with floating point numbers. So, in the test code accompanying the library the square root of 4/3 is calculated with maximum precision to be 50843527/44031786, which is accurate to 15 decimal places (about one quadrillionth). If the maximum value for numerator and denominator is reduced from 67,108,864 to 7 (i.e. a 3 bit integer for the numerator and a 3 bit integer for the denominator) it calculates sqrt(4/3) as 7/6 (which is accurate to only one decimal place) but returns 7/6 squared as 4/3 -- the error is cancelled out. I did not choose 4/3 deliberately for this property, it appears to happen most if not all of the time.
Just for fun I included a fairly meaningless timing test. On my first gen iPad mini the timing loops both take about 0.8 seconds, suggesting that using the included floating point square root in a for next loop is about 640 times as fast as the vulgar square root given here in a for next loop. This is a worst case for the vulgar square root. If a square root with a rational (exact) answer is calculated it takes about half the time, and if the maximum representable number is reduced to 7 the execution time is also halved.
Further reading.
Some links that I found useful or interesting whilst writing this code.
The algorithm I used for the square root code.
http://tedmuller.us/Math/Pencil&PaperSquareRoots.htm
List of rational number functions provided in a cryptography oriented programming language. ("FLASH = Floating sLASH - similar to fixed slash.)
http://docs.certivox.com/docs/miracl/mi ... ual$chap27
This seems to be the book you need to own if you want to get serious about rational numbers.
https://books.google.co.uk/books?id=OAEgAwAAQBAJ
A lovely description of Mediant Rounding and other basic considerations from the days of the mini-computer.
http://people.csail.mit.edu/bkph/papers ... ic_OPT.pdf
Rational Numbers, Vulgar Words - a more extensive implementation I did in Forth back in the day.
http://www.forth.org/fd/FD-V15N5.pdf
Some Vulgar Functions - Extending the previous article to include trigonometric and logarithmic etc. functions.
http://www.forth.org/fd/FD-V16N2.pdf
Gordon Charlton
Copy to your dropbox from here:
https://www.dropbox.com/s/bxjhs3x9t181m ... r.txt?dl=0
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
