'BIGNUM.BAS v0.n
'Sep-Dec 1996 by Marc Kummel aka Treebeard.
'Contact mkummel@rain.org, http://www.rain.org/~mkummel/
' https://web.archive.org/web/20200220020034/http://www.rain.org/~mkummel/tbvault.html
'Big number arithmetic routines, done with a crypto project in mind. It's
'an engaging task to work out these algorithms in BASIC, since we know them
'so well on paper. My aim is to minimize string space; then speed. Not
'quite "arbitrary precision" since numbers and answer must fit in BASIC
'string space, 64K with PDS far strings, much less for QBasic. Don't use
'the same string twice in one call! bMul(s$,s$,s$) should figure s=s*s, but
'it won't work because passed strings are modified (and restored) during the
'call, so use temporary strings or pass by value. See bModPower() for an
'example. Can be much optimized, but it works and makes sense.
'
'---------------------------------------------------------------------------
'useful globals
Const false
= 0, true
= Not false
, nada
= -1 Const zero$
= "0", one$
= "1", two$
= "2", three$
= "3" Const four$
= "4", five$
= "5" Const dp$
= ".", neg$
= "-", asc0
= 48 Const basechr$
= "B", basesep$
= "," 'number bases Const negative
= -1, positive
= 1 'returned by bComp() Const maxlongdig
= 8 'max digits in long var& Const maxsegment
= 65536 'memory segment size Const memget
= 0, memput
= 1, memclr
= 2 'memory operations Const maxstack
= 10 'stack depth Const maxmem
= 14 'memory slots (0-9 + constants) Const pimem
= 10, pi2mem
= 11, emem
= 12 'pi, 2pi, e Const ln10mem
= 13, ln2mem
= 14 'natural logs of 10 & 2
'useful shared stuff
Dim Shared topline%
, botline%
, topview%
, botview%
'Trick to give bmem$() its own 64k segment in PDS (but not QB).
'Watch out for overflows or conflicts with multiple modules.
Common bmem$
(), zmem$
(), cname$
() Dim bmem$
(maxmem
), zmem$
(1 To maxstack
), cname$
(10 To maxmem
)
esc$
= Chr$(27) 'useful stringsnull$ = ""
error$
= "error! " 'unhelpful error messages abort$ = "<abort>"
op$ = "+" 'default action
digits% = 100 'digits for division etc
bignumpath$ = null$ 'path for files (or current dir if null)
prmcntfile$ = "BNPRMCNT.DAT" 'prime count table
logfile$ = "BNLOG.LOG" 'log file
ReDim prmcnt&
(0) 'stub for prime count table
'useful constants
cname$(pimem) = "ã": bmem$(pimem) = "3.14159265358979323846264338327"
cname$(pi2mem) = "2ã": bmem$(pi2mem) = "6.28318530717958647692528676654"
cname$(emem) = "e": bmem$(emem) = "2.71828182845904523536028747135"
cname$(ln10mem) = "ln(10)": bmem$(ln10mem) = "2.30258509299404568401799145468"
cname$(ln2mem) = "ln(2)": bmem$(ln2mem) = ".693147180559945309417232121458"
'screen size and window
topline% = 1
botline% = 25
topview% = topline%
botview% = botline% - 2
'screen buffer for ScreenSave() and ScreenRestore()
screenbufsize% = 80 * botline% * 2 'screen size, 2 bytes per character
Dim screenbuf%
(screenbufsize%
- 1) 'buffer screenrow% = 1 'cursor row
Def Seg = 0 'video segment for mono or color If Peek(&H463) = &HB4 Then screenseg&
= &HB000 Else screenseg&
= &HB800 'color
'trap ctrl-esc for emergency exit. Speed and size hog, but useful.
'On Key(15) GoSub EventTrap
'start up
' LoadPrimeTable 'prime count table
' ShowHelp 'something to look at
'========================================================================
n$ = "1"
m$ = "99999999999999999999999999999999999999999999999999"
r$ = ""
m$ = r$ ' you can't use the destination variable in the source
'========================================================================
's = |s|
'
'return true if s is negative
'
bIsNeg%
= (Left$(s$
, 1) = neg$
)
'return sign of number (-1 or +1)
'
If bIsNeg
(s$
) Then bSign%
= negative
Else bSign%
= positive
'Strip leading 0s and final "." (but leave something)
'
n% = 1
n% = n% + 1
'Strip trailing 0s to "." (but leave something)
'
n% = n% - 1
'Strip s$ to whole number and base 10 integer logarithm and sign. Decimal
'point is implied after the first digit, and slog% counts places left or
'right. bLogPut() reverses the process, and bLogDp() gives info on the
'decimals. Tricky, but it works and simplifies dividing and multipling.
'eg s$ -> s$ , slog%
' 660 -> 66 , 2 (6.6 * 10^ 2) (or 660,2 if zeroflag%=false)
' 6.6 -> 66 , 0 (6.6 * 10^ 0)
' .066 -> 66 , -2 (6.6 * 10^-2)
'bDiv(), bMul(), and bSqr() use this to trim unnecessary zeros and to locate
'decimal point. These set zeroflag% to trim trailing zeros, but bDivIntMod()
'must set it false in order to figure remainder of division. A kludge.
'
Sub bLogGet
(s$
, slog%
, sign%
, zeroflag%
) bStripZero s$
n% = dpt% + 1
n% = n% + 1
slog% = dpt% - n%
s$
= Left$(s$
, dpt%
- 1) + Mid$(s$
, dpt%
+ 1) slog% = dpt% - 2
'remove trailing 0's if zeroflag%
If zeroflag%
Then bStripTail s$
'Strip a number to "standard form" with no leading or trailing 0s and no
'final "." All routines should return all arguments in this form.
'
bStripZero s$
'Restore a number from the integer and log figured in bLogGet(). s$ is taken
'as a number with the decimal after first digit, and decimal is moved slog%
'places left or right, adding 0s as required. Called by bDiv() and bMul().
'
Sub bLogPut
(s$
, slog%
, sign%
) s$ = zero$
s$
= dp$
+ String$(-slog%
- 1, zero$
) + s$
s$
= s$
+ String$(slog%
- last%
+ 1, zero$
) + dp$
s$
= Left$(s$
, slog%
+ 1) + dp$
+ Mid$(s$
, slog%
+ 2) bClean s$
If sign%
= negative
Then s$
= neg$
+ s$
'shift decimal n% digits (minus=left), i.e multiply/divide by 10.
'
bLogGet s$, slog%, sign%, false
bLogPut s$, slog% + n%, sign%
's = -s
'
'Take whole number and log from bLogGet() and return number of decimal
'places in the expanded number; OR take string and number of decimal points
'desired and return the log. It works both ways.
'
bLogDp%
= Len(s$
) - 1 - logdp%
'out = s1 / s2 using fast long-integer algorithm. s2$ must be <= 8 digits.
's1$ and s2$ must be stripped first, no decimals.
'
Sub bDivLong
(s1$
, s2$
, quotient$
, remainder$
) quotient$ = null$
remainder& = 0
dividend&
= remainder&
* 10&
+ Val(Mid$(s1$
, i%
, 1)) dig% = dividend& \ divisor&
quotient$
= quotient$
+ Chr$(asc0
+ dig%
) remainder& = dividend& - dig% * divisor&
'out = s1 / s2 using character algorithm, digit by digit, slow but honest.
's1$ and s2$ must be stripped first, no decimals.
'
Sub bDivChar
(s1$
, s2$
, quotient$
, remainder$
) last1%
= Len(s1$
) 'length of the dividend last2%
= Len(s2$
) 'length of the divisor quotient$ = null$
remainder$ = null$
'get next digit of dividend or zero$ if past end
dvd$
= remainder$
+ Mid$(s1$
, i%
, 1) dvd$ = remainder$ + zero$
'if dividend < divisor then digit%=0 else have to calculate it.
'do fast compare using string operations. see bComp%()
bStripZero dvd$
If (ldvd%
< last2%
) Or ((ldvd%
= last2%
) And (dvd$
< s2$
)) Then 'divisor is bigger, so digit is 0, easy!
dig% = 0
remainder$ = dvd$
'dividend is bigger, but no more than 9 times bigger.
'subtract divisor until we get remainder less than divisor.
'time hog, average is 5 tries through j% loop. There's a better way.
remainder$ = null$
borrow% = 0
n% = last2% - j%
n%
= Val(Mid$(dvd$
, ldvd%
- j%
, 1)) - n%
- borrow%
If n%
>= 0 Then borrow%
= 0 Else borrow%
= 1: n%
= n%
+ 10 remainder$
= Chr$(asc0
+ n%
) + remainder$
'if remainder < divisor then exit
bStripZero remainder$
dvd$ = remainder$
quotient$
= quotient$
+ Chr$(asc0
+ dig%
)
'out = s1 / s2
'
'strip divisor
t$ = s2$
bLogGet t$, slog2%, sign2%, true
'divide by zero?
'do powers of 10 with shifts
'the hard way
'strip all
s2$ = t$: t$ = null$
bLogGet s1$, slog1%, sign1%, true
'figure decimal point and sign of answer
outlog% = slog1% + bLogDp(s2$, slog2%)
If sign1%
<> sign2%
Then outsign%
= negative
Else outsign%
= positive
'bump digits past leading zeros and always show whole quotient
olddigits% = digits%
digits%
= digits%
+ Len(s2$
) If digits%
< outlog%
+ 1 Then digits%
= outlog%
+ 1
'do it, ignore remainder
'clean up
bLogPut
out$
, outlog%
, outsign%
bLogPut s1$, slog1%, sign1%
bLogPut s2$, slog2%, sign2%
digits% = olddigits%
