QB64.org Forum
Active Forums => Programs => Topic started by: luke on June 05, 2021, 11:21:04 am
-
Supports integer operations of practically unlimited size. Addition, subtraction multiplication and division of positive and negative numbers are available, as is an efficient exponentiation. Lightly tested, as always.
Full list of commands and their usage:
Dim As bn_int a, b, q, r
Dim As String s
Dim As _Unsigned Long n
+----------------+------------+----------------------------+---------------------------------+-------------------------------------------------+
| Function | Arguments | Full Name | Operation | Notes |
+----------------+------------+----------------------------+---------------------------------+-------------------------------------------------+
| bni_new | a | New | Initialise memory | Must be called before use |
| bni_del | a | Delete | Free memory | Must be called when finished |
| bni_fromstr | a, s | From String | a = Val(s) | Accepts optional sign |
| bni_tostr$ | a | To String | Return Str$(a) | Always adds sign unless 0 |
| bni_fromnative | a, n | From Native | a = n | |
| bni_add | a, b, r | Add | r = a + b | |
| bni_addu | a, b, r | Add Unsigned | r = a + b, treat a & b unsigned | |
| bni_sub | a, b, r | Subtract | r = a - b | |
| bni_subu | a, b, r | Subtract Unsigned | r = a - b, treat a & b unsigned | If a < b, r is the base complement and negative |
| bni_neg | a | Negate | a = -a | |
| bni_mul | a, b, r | Multiply | r = a * b | a, b, r must not be aliased |
| bni_mulp | a, b | Multiply in-place | a = a * b | |
| bni_shl | a, n | Shift Left | a = a * 2^n | Requires n <= FOOT_SIZE * 8 |
| bni_shr | a, n | Shift Right | a = a / 2^n | Requires n <= FOOT_SIZE * 8 |
| bni_basemul | a, n | Multiply by Base | a = a * FOOT_BASE ^ n | |
| bni_div | a, b, q, r | Divide | q = a / b, remainder r | a, b, q, r must not be aliased. r >= 0 |
| bni_divp | a, b, r | Divide in-place | a = a / b, remainder r | r >= 0 |
| bni_pow | a, n, r | Power | r = a^n | a, r must not be aliased |
| bni_powp | a, n | Power in-place | a = a^n | |
| bni_tz | a | Test Zero | a = 0? | |
| bni_tn | a | Test Negative | a < 0? | |
| bni_tp | a | Test Positive | a > 0? | |
| bni_teq | a, b | Test Equal | a = b? | |
| bni_tlt | a, b | Test Less Than | a < b? | |
| bni_tltu | a, b | Test Less Than Unsigned | Abs(a) < Abs(b)? | |
| bni_tgt | a, b | Test Greater Than | a > b? | |
| bni_tgtu | a, b | Test Greater Than Unsigned | Abs(a) > Abs(b)? | |
+----------------+------------+----------------------------+---------------------------------+-------------------------------------------------+
Functions and demo program:
Option _Explicit
_Define A-Z As _UNSIGNED LONG
Const FOOT_SIZE = 4
Const FOOT_MAX = 2 ^ 32 - 1
Const FOOT_BASE = 2 ^ 32
Type bn_int
length As Long
feet As _MEM
sign As Integer '-1 = negative, 1 = positive, 0 = zero
End Type
Dim As bn_int a, b, c, d
bni_new a
bni_new b
bni_new c
bni_new d
bni_fromnative a, 10
bni_powp a, 2784
bni_fromstr b, "999999999999999999999998999999999999999999999999"
Dim t#
t# = Timer(0.001)
bni_div a, b, c, d
'The 115th Fibonacci number
Print Mid$(bni_tostr$(c), 2715)
Print "Time: "; Timer(0.001) - t#
'Initialise a new biginteger with value 0
Sub bni_new (a As bn_int)
a.feet = _MemNew(FOOT_SIZE * 4) 'Start with 4 feet, somewhat arbitrarily
If a.feet.SIZE <> FOOT_SIZE * 4 Then Error 257
a.length = 1
set_foot a, 0, 0
End Sub
'Delete a biginteger
Sub bni_del (a As bn_int)
_MemFree a.feet
End Sub
'Parse string of digits into biginteger
Sub bni_fromstr (a As bn_int, original_in$)
'We will read 7 digits at a time, because that is guaranteed to be at most a 32-bit value
Const DIGITS_PER_ROUND = 7
Dim i, sign, multiplier As bn_int, temp As bn_int, s$, in$
bni_new temp
bni_new multiplier
bni_fromnative a, 0
bni_fromnative multiplier, 10 ^ DIGITS_PER_ROUND
'Remove any sign prefix if present
If Left$(original_in$, 1) = "-" Then
sign = -1
in$ = Mid$(original_in$, 2)
ElseIf Left$(original_in$, 1) = "+" Then
sign = 1
in$ = Mid$(original_in$, 2)
Else
in$ = original_in$
End If
For i = 1 To Len(in$) Step DIGITS_PER_ROUND
s$ = Mid$(in$, i, DIGITS_PER_ROUND)
bni_fromnative temp, Val(s$)
If Len(s$) < DIGITS_PER_ROUND Then
'This only occurs on the final loop if the number of digits is not a multiple of DIGITS_PER_ROUND
bni_fromnative multiplier, 10 ^ Len(s$)
End If
bni_mulp a, multiplier
bni_add a, temp, a
Next i
If bni_tz(a) Then
a.sign = 0
ElseIf sign = -1 Then
a.sign = -1
Else
a.sign = 1
End If
bni_del temp
bni_del multiplier
End Sub
'b = a
Sub bni_copy (a As bn_int, b As bn_int)
If a.feet.OFFSET = b.feet.OFFSET Then Exit Sub
If b.length Then bni_del b
b.length = a.length
b.sign = a.sign
b.feet = _MemNew(a.feet.SIZE)
_MemCopy a.feet, a.feet.OFFSET, a.length * FOOT_SIZE To b.feet, b.feet.OFFSET
End Sub
'Store an _unsigned long in a biginteger
Sub bni_fromnative (a As bn_int, b)
'No need to grow because all bn_int are at least 4 feet
a.length = 1
a.sign = Sgn(b) 'Note: b is _unsigned so this is a zero test
set_foot a, 0, b
End Sub
'Store an _unsigned _integer64 in a biginteger
Sub bni_fromnative_long (a As bn_int, b)
'No need to grow because all bn_int are at least 4 feet
set_foot a, 0, b And FOOT_MAX
Dim t
t = _SHR(b, FOOT_SIZE * 8)
If t Then
set_foot a, 1, t
a.length = 2
Else
a.length = 1
End If
a.sign = Sgn(b)
End Sub
Function bni_tostr$ (a As bn_int)
Dim human_base As bn_int, q As bn_int, r As bn_int, s$, d$
bni_new human_base
'Limited by the size of an _unsigned long
bni_fromnative human_base, 10000000
bni_new q
bni_new r
bni_copy a, q
Do
bni_divp q, human_base, r
d$ = LTrim$(Str$(get_foot(r, 0)))
If Not bni_tz(q) Then d$ = String$(7 - Len(d$), "0") + d$
s$ = d$ + s$
Loop While Not bni_tz(q)
If a.sign = -1 Then
s$ = "-" + s$
ElseIf a.sign = 1 Then
s$ = "+" + s$
Else
s$ = " " + s$
End If
bni_del r
bni_del q
bni_del human_base
bni_tostr$ = s$
End Function
Sub bni_add (a As bn_int, b As bn_int, r As bn_int)
If a.sign = 0 Then
'quick return on 0
bni_copy b, r
ElseIf b.sign = 0 Then
bni_copy a, r
ElseIf a.sign = b.sign Then
'pos + pos = pos; neg - neg = -(pos + pos)
bni_addu a, b, r
r.sign = a.sign
ElseIf a.sign = -1 Then
'neg + pos = pos - pos
If bni_tgtu(a, b) Then
bni_subu a, b, r
bni_neg r
Else
bni_subu b, a, r
End If
ElseIf b.sign = -1 Then
'pos + neg = pos - pos
If bni_tgtu(b, a) Then
bni_subu b, a, r
bni_neg r
Else
bni_subu a, b, r
End If
End If
End Sub
'r = a + b for bigintegers
Sub bni_addu (a As bn_int, b As bn_int, r As bn_int)
Dim i, carry, a_foot, b_foot, r_foot
'The final result can be up to one foot longer than the longest addend.
bni_grow r, max(a.length, b.length) + 1
For i = 0 To max(a.length, b.length) - 1
If i < a.length Then
a_foot = get_foot(a, i)
Else
a_foot = 0
End If
If i < b.length Then
b_foot = get_foot(b, i)
Else
b_foot = 0
End If
If carry Then
r_foot = a_foot + b_foot + 1
If r_foot <= a_foot Then carry = 1 Else carry = 0
Else
r_foot = a_foot + b_foot
If r_foot < a_foot Then carry = 1 Else carry = 0
End If
set_foot r, i, r_foot
Next i
set_foot r, i, carry
r.length = max(a.length, b.length) + carry
r.sign = 1
End Sub
'r = a - b for bigintegers
Sub bni_sub (a As bn_int, b As bn_int, r As bn_int)
If bni_tz(a) Then
'0 - b = -b
bni_copy b, r
bni_neg r
ElseIf bni_tz(b) Then
'a - 0 = a
bni_copy a, r
ElseIf bni_tn(a) And bni_tn(b) Then
'neg(a) - neg(b) = neg(a) + pos(b) = pos(b) - pos(a)
If bni_tlt(a, b) Then
bni_subu a, b, r
bni_neg r
Else
bni_subu b, a, r
End If
ElseIf bni_tn(a) Then
'neg(a) - pos(b) = -(pos(a) + pos(b))
bni_addu a, b, r
bni_neg r
ElseIf bni_tn(b) Then
'pos(a) - neg(b) = pos(a) + pos(b)
bni_addu a, b, r
ElseIf bni_tlt(a, b) Then
'a - b < 0, so a - b = -(b - a)
bni_subu b, a, r
bni_neg r
Else
bni_subu a, b, r
End If
End Sub
'r = a - b
'In general assumes a >= b. If a < b then r is negative and
'its value is the base complement of the real answer.
Sub bni_subu (a As bn_int, b As bn_int, r As bn_int)
Dim i, borrow, a_foot, b_foot, r_foot, r_len
bni_grow r, max(a.length, b.length)
r_len = 1 'Don't write directly to r.length in case a = r or b = r
For i = 0 To max(a.length, b.length) - 1
If i < a.length Then
a_foot = get_foot(a, i)
Else
a_foot = 0
End If
If i < b.length Then
b_foot = get_foot(b, i)
Else
b_foot = 0
End If
r_foot = a_foot - b_foot + borrow
If borrow Then
r_foot = a_foot - b_foot - 1
borrow = r_foot >= a_foot
Else
r_foot = a_foot - b_foot
borrow = r_foot > a_foot
End If
set_foot r, i, r_foot
If r_foot <> 0 Then r_len = i + 1
Next i
If r_len = 1 And get_foot(r, 0) = 0 Then r.sign = 0 Else r.sign = 1
r.length = r_len
If borrow Then r.sign = -1 'Turns out a < b
End Sub
'a = -a
Sub bni_neg (a As bn_int)
a.sign = -Sgn(a.sign)
End Sub
'r = a * b for bigintegers
'Requires r <> b, see bni_mulp for alternative
'Based on Knuth 4.3.1 Algorithm M
Sub bni_mul (a As bn_int, b As bn_int, r As bn_int)
Dim i, j, r_foot, carry
'a_foot and b_foot are an int64 so a 64-bit multiplication is performed.
'Otherwise a 32-bit multiplication would drop the higher bits.
Dim a_foot As _Unsigned _Integer64, b_foot As _Unsigned _Integer64
Dim temp As _Unsigned _Integer64
'The final result can be up to the sum of the length of the multiplicands
bni_grow r, a.length + b.length
'Clear the destination; if you don't, you get a multiply & add instead.
bni_zero r
For i = 0 To a.length - 1
carry = 0
a_foot = get_foot(a, i)
For j = 0 To b.length - 1
b_foot = get_foot(b, j)
r_foot = get_foot(r, i + j)
temp = a_foot * b_foot + r_foot + carry
set_foot r, i + j, temp And FOOT_MAX
carry = _SHR(temp, FOOT_SIZE * 8)
Next j
set_foot r, i + j, carry
Next i
r.length = a.length + b.length
If carry = 0 Then r.length = r.length - 1
r.sign = a.sign * b.sign
End Sub
'a = a * b for bigintegers
Sub bni_mulp (a As bn_int, b As bn_int)
Dim result As bn_int
bni_new result
bni_mul a, b, result
bni_del a
a = result
End Sub
'a = a * 2^b, b <= FOOT_SIZE * 8
Sub bni_shl (a As bn_int, b)
Dim i, foot, mask, mask_shift, carry, carry2
If b > FOOT_SIZE * 8 Then Error 5
bni_grow a, a.length + 1
'bitmask for bits that will be carried between feet
mask_shift = FOOT_SIZE * 8 - b
mask = _SHL(2 ^ b - 1, mask_shift)
For i = 0 To a.length - 1
foot = get_foot(a, i)
carry = foot And mask
foot = _SHL(foot, b) + carry2
carry2 = _SHR(carry, mask_shift)
set_foot a, i, foot
Next i
If carry2 Then
'add a new foot
a.length = a.length + 1
set_foot a, a.length - 1, carry2
End If
End Sub
'a = a / 2^b, b <= FOOT_SIZE * 8
Sub bni_shr (a As bn_int, b)
Dim i, foot, mask, mask_shift, carry, carry2
If b > FOOT_SIZE * 8 Then Error 5
'bitmask for bits that will be carried between feet
mask_shift = FOOT_SIZE * 8 - b
mask = 2 ^ b - 1
For i = a.length - 1 To 0 Step -1
foot = get_foot(a, i)
carry = foot And mask
foot = _SHR(foot, b) + carry2
carry2 = _SHL(carry, mask_shift)
set_foot a, i, foot
Next i
End Sub
'a = a * b^n
'Assumes n >= 0
Sub bni_basemul (a As bn_int, n)
bni_grow a, a.length + n
_MemCopy a.feet, a.feet.OFFSET, a.length * FOOT_SIZE To a.feet, a.feet.OFFSET + n * FOOT_SIZE
_MemFill a.feet, a.feet.OFFSET, n * FOOT_SIZE, 0 As _UNSIGNED LONG
a.length = a.length + n
End Sub
Sub bni_div (a As bn_int, b As bn_int, q As bn_int, r As bn_int)
If bni_tz(b) Then
Error 5
ElseIf bni_tlt(a, b) Then
bni_fromnative q, 0
bni_copy a, r
ElseIf b.length = 1 Then
bni_div_short a, b, q, r
Else
bni_div_big a, b, q, r
End If
End Sub
'a = a / b, remainder r for bigintegers
Sub bni_divp (a As bn_int, b As bn_int, r As bn_int)
Dim result As bn_int
bni_new result
bni_div a, b, result, r
bni_del a
a = result
End Sub
'q = a / b, remainder r
'Requires b <= FOOT_MAX, gives r >= 0
Sub bni_div_short (a As bn_int, b As bn_int, q As bn_int, r As bn_int)
Dim i, a_foot, b_foot, d
Dim temp As _Unsigned _Integer64
If b.length > 1 Then Error 5 'Use bni_bigdiv for bigger divisors
bni_grow q, a.length 'This might be one foot too large, we'll deal with it at the end
'remainder may be as big as the divisor
bni_grow r, b.length
bni_fromnative r, 0
b_foot = get_foot(b, 0)
For i = a.length - 1 To 0 Step -1
a_foot = get_foot(a, i)
temp = (temp * FOOT_BASE + a_foot)
d = temp \ b_foot
set_foot q, i, d
temp = temp Mod b_foot
Next i
q.length = a.length
'Trim leading 0
If get_foot(q, q.length - 1) = 0 And q.length > 1 Then q.length = q.length - 1
set_foot r, 0, temp
r.length = 1
If temp Then r.sign = 1 Else r.sign = 0
If q.length = 1 And get_foot(q, 0) = 0 Then q.sign = 0 Else q.sign = a.sign * b.sign
End Sub
'q = a / b, remainder r
'Requires a >= b, b > FOOT_MAX
'Based on Knuth 4.3.1 Algorithm D
Sub bni_div_big (a As bn_int, b As bn_int, q As bn_int, r As bn_int)
Dim d, j, n, m
Dim As bn_int u, v, temp
Dim As _Unsigned _Integer64 temp2, qhat, rhat
bni_new temp
bni_new u
bni_new v
bni_copy a, u
bni_copy b, v
n = v.length
m = u.length - v.length
bni_grow q, m + 1 'TODO: Shrink if too big
d = bni_div_big_normalise(u, v)
For j = m To 0 Step -1
'Do not combine these next few lines, it is separate to force 64-bit calculations
temp2 = get_foot(u, j + n)
temp2 = temp2 * FOOT_BASE
temp2 = temp2 + get_foot(u, j + n - 1)
qhat = temp2 \ get_foot(v, n - 1)
rhat = temp2 Mod get_foot(v, n - 1)
'If qhat < foot_base then qhat * v_{n-2} will not overflow
While qhat >= FOOT_BASE Or qhat * get_foot(v, n - 2) > FOOT_BASE * rhat + get_foot(u, j + n - 2)
qhat = qhat - 1
rhat = rhat + get_foot(v, n - 1)
If rhat >= FOOT_BASE Then Exit While
Wend
bni_fromnative_long temp, qhat
bni_mulp temp, v
bni_basemul temp, j
bni_subu u, temp, u
If bni_tn(u) Then
'qhat was too big, add back one divisor
qhat = qhat - 1
bni_copy v, temp
bni_basemul temp, j
bni_addu u, temp, u
'This will cause a carry into the next digit that we don't care about
u.length = u.length - 1
End If
set_foot q, j, qhat
Next j
q.length = m + 1
If get_foot(q, 0) = 0 Then q.length = q.length - 1
q.sign = a.sign * b.sign
bni_shr u, d
bni_copy u, r
End Sub
'Normalise a & b for division and return the normalisation constant.
Function bni_div_big_normalise (a As bn_int, b As bn_int)
'Need high foot of b to be >= FOOT_BASE/2, which can be done by
'bitshifting until the high bit is 1.
Dim foot, mask, d, old_length
foot = get_foot(b, b.length - 1)
mask = FOOT_BASE / 2
'Find highest 1 bit
Do While (foot And mask) = 0
d = d + 1
mask = _SHR(mask, 1)
Loop
old_length = a.length
bni_shl a, d
If a.length = old_length Then
'Algorithm requires a have a leading foot, even if it's zero
bni_grow a, a.length + 1
a.length = a.length + 1
set_foot a, a.length - 1, 0
End If
bni_shl b, d
bni_div_big_normalise = d
End Function
'r = a^b
'Based on Knuth 4.6.3 Algorithm A
Sub bni_pow (a As bn_int, b, r As bn_int)
Dim n, parity
Dim As bn_int z, t
n = b
bni_new t
bni_new z
bni_fromnative r, 1
bni_copy a, z
Do
parity = n And 1
n = n \ 2
If parity Then
bni_mulp r, z
If n = 0 Then Exit Sub
End If
bni_mulp z, z
Loop
End Sub
'a = a^b
Sub bni_powp (a As bn_int, b)
Dim result As bn_int
bni_new result
bni_pow a, b, result
bni_copy result, a
bni_del result
End Sub
'a = 0?
Function bni_tz (a As bn_int)
bni_tz = a.sign = 0
End Function
'a < 0?
Function bni_tn (a As bn_int)
bni_tn = a.sign = -1
End Function
'a > 0?
Function bni_tp (a As bn_int)
bni_tp = a.sign = 1
End Function
'a = b?
Function bni_teq (a As bn_int, b As bn_int)
If a.length <> b.length Or a.sign <> b.sign Then Exit Function
Dim i
For i = 0 To a.length - 1
If get_foot(a, i) <> get_foot(b, i) Then Exit Function
Next i
bni_teq = -1
End Function
'a < b?
Function bni_tlt (a As bn_int, b As bn_int)
'a.sign = b.sign so just a.sign throughout
If a.sign > b.sign Or a.sign * a.length > a.sign * b.length Then Exit Function
If a.sign < b.sign Or a.sign * a.length < a.sign * b.length Then
bni_tlt = -1
Exit Function
End If
Dim i, cmp
For i = a.length - 1 To 0 Step -1
cmp = a.sign * get_foot(a, i) - a.sign * get_foot(b, i)
If cmp < 0 Then
bni_tlt = -1
Exit Function
ElseIf cmp > 0 Then
Exit Function
End If
Next i
End Function
'abs(a) < abs(b)?
Function bni_tltu (a As bn_int, b As bn_int)
If a.length > b.length Then Exit Function
If a.length < b.length Then
bni_tltu = -1
Exit Function
End If
Dim i, cmp
For i = a.length - 1 To 0 Step -1
cmp = get_foot(a, i) - get_foot(b, i)
If cmp < 0 Then
bni_tltu = -1
Exit Function
ElseIf cmp > 0 Then
Exit Function
End If
Next i
End Function
'a > b?
Function bni_tgt (a As bn_int, b As bn_int)
'a.sign = b.sign so just a.sign throughout
If a.sign < b.sign Or a.sign * a.length < a.sign * b.length Then Exit Function
If a.sign > b.sign Or a.sign * a.length > a.sign * b.length Then
bni_tgt = -1
Exit Function
End If
Dim i, cmp
For i = a.length - 1 To 0 Step -1
cmp = a.sign * get_foot(a, i) - a.sign * get_foot(b, i)
If cmp > 0 Then
bni_tgt = -1
Exit Function
ElseIf cmp < 0 Then
Exit Function
End If
Next i
End Function
'abs(a) > abs(b)?
Function bni_tgtu (a As bn_int, b As bn_int)
If a.length < b.length Then Exit Function
If a.length > b.length Then
bni_tgtu = -1
Exit Function
End If
Dim i, cmp
For i = a.length - 1 To 0 Step -1
cmp = get_foot(a, i) - get_foot(b, i)
If cmp > 0 Then
bni_tgtu = -1
Exit Function
ElseIf cmp < 0 Then
Exit Function
End If
Next i
End Function
'Set all feet in a biginteger to 0, regardless of length
Sub bni_zero (a As bn_int)
_MemFill a.feet, a.feet.OFFSET, a.feet.SIZE, 0 As _UNSIGNED LONG
End Sub
'Ensure a biginteger has at least min_capacity feet
Sub bni_grow (a As bn_int, min_capacity)
Dim current_capacity As _Offset, new_capacity As _Offset, new_feet As _MEM
current_capacity = a.feet.SIZE / FOOT_SIZE
If current_capacity >= min_capacity Then Exit Sub
'Get new capacity with headroom by doublings (ensuring it's always an exact number of feet)
new_capacity = current_capacity
While new_capacity < min_capacity
new_capacity = new_capacity * 2
Wend
'Create new mem block
new_feet = _MemNew(FOOT_SIZE * new_capacity)
If new_feet.SIZE <> FOOT_SIZE * new_capacity Then Error 257
'Copy old feet across
_MemCopy a.feet, a.feet.OFFSET, a.length * FOOT_SIZE To new_feet, new_feet.OFFSET
'And put it all in-place
_MemFree a.feet
a.feet = new_feet
End Sub
Function get_foot (a As bn_int, p)
get_foot = _MemGet(a.feet, a.feet.OFFSET + p * FOOT_SIZE, _Unsigned Long)
End Function
Sub set_foot (a As bn_int, p, v)
_MemPut a.feet, a.feet.OFFSET + p * FOOT_SIZE, v
End Sub
Function min (a, b)
If a < b Then min = a Else min = b
End Function
Function max (a, b)
If a < b Then max = b Else max = a
End Function
-
you have the right idea, string math is woefully slow and cumbersome
-
It's not really as slow as you make it (unless you are running on a Commodore 64!) We are talking nano-seconds with today's CPU's.
you have the right idea, string math is woefully slow and cumbersome
-
It's not really as slow as you make it (unless you are running on a Commodore 64!) We are talking nano-seconds with today's CPU's.
I really wish I could agree with George but after fixing subtr$() Function, it takes 30 secs now to calculate inverse of STx number for 1st 115 terms in Fibonacci. ref: https://www.qb64.org/forum/index.php?topic=2921.msg133102#msg133102
@luke I am looking forward to seeing more.
-
Out of interest I compared these routines against the code at https://www.qb64.org/forum/index.php?topic=2921.msg121797#msg121797 by calculating some factorials.
My driver programs:one$ = "1"
factorial$ = "1"
limit$ = "10000"
v$ = "1"
t# = Timer(0.001)
Do
v$ = add$(v$, one$)
factorial$ = mult$(factorial$, v$)
Loop Until v$ = limit$
t# = Timer(0.001) - t#
Print factorial$
Print t#
Dim As bn_int factorial, v, limit, one
bni_new factorial
bni_new v
bni_new one
bni_new limit
bni_fromnative one, 1
bni_fromnative factorial, 1
bni_fromnative limit, 10000
bni_fromnative v, 1
Dim t#, result$
t# = Timer(0.001)
Do
bni_add v, one, v
bni_mulp factorial, v
Loop Until bni_teq(v, limit)
result$ = bni_tostr$(factorial)
t# = Timer(0.001) - t#
Print result$
Print t#It's written in such a stilted fashion because I don't have subtraction routines yet.
For smaller numbers (~1000!) the stringless approach was faster by about an order of magnitude: 0.03s vs 0.37s. For 10000!, the stringless approach took about 3.7 seconds*. I killed the string method after a minute with no result. So I think this approach has some merit!
* The actual factorial calculation takes about 3.0 seconds. The rest of the time is spent converting the number into a string for display, which is a non-trivial operation involving lots of divisions.
EDIT: 100000! took 6.5 minutes5.5 minutes with $checking:off. I don't know the answer's correct but it had the right number of digits, so I'll call that a success.
-
Impressive but not the print out here:
35660 digits confirmed on internet.
-
There we go, finally got subtraction and full division working, and powers too. As a demo the program calculates 10^2784 / 999999999999999999999998999999999999999999999999 to find the 115th Fibonacci number (after https://www.qb64.org/forum/index.php?topic=2921.msg133102#msg133102 ). On my computer this operation (including converting the result to base 10 and PRINTing it) takes 0.01 seconds, which is somewhat quicker than the string math programs :)
-
"somewhat"
String math is doing 10000! in 3.25 mins * 60 / 3.7 secs = 53 x's faster
Inverse STx now taking about 2 secs, say 2.1 on outside (on a good run you can get it under 2) / .01 = 210 x's faster
Why mult by 10^2784? You can find Fibonacci terms from decimal to right. Does that put the 115th term just before the decimal?
-
Why mult by 10^2784? You can find Fibonacci terms from decimal to right. Does that put the 115th term just before the decimal?
Because this is all strictly integers, just doing 1/STx would give 0. I worked out the 115th number appears after 2784 decimal places, so I instead calculate 10^2784/STx which puts the 115th term just before the decimal.
To keep things neat I extract it with MID$, but the full number result is indeed all the Fibonacci numbers up to 115th.
-
Because this is all strictly integers, just doing 1/STx would give 0. I worked out the 115th number appears after 2784 decimal places, so I instead calculate 10^2784/STx which puts the 115th term just before the decimal.
To keep things neat I extract it with MID$, but the full number result is indeed all the Fibonacci numbers up to 115th.
Oh you've not gotten to the Mr$() part. Mr$ stands for math regulator that handles all the stuff for signs and decimals.
myAnswer (as ?) = Mr(float1 as ?, opCode, float2 as ?) ' opCode for +-*/
but how would you get back to a number? You pretty much have to go to a string form sometime in process, so you could pass and print numbers.
-
I currently have no plans to extend this to floating point, because that's a can of worms I'm not quite ready to open yet.
-
Oh you've not gotten to the Mr$() part. Mr$ stands for math regulator that handles all the stuff for signs and decimals.
myAnswer (as ?) = Mr(float1 as ?, opCode, float2 as ?) ' opCode for +-*/
but how would you get back to a number? You pretty much have to go to a string form sometime in process, so you could pass and print numbers.
I fiddled around with a non-string method last year, for a bit, but came to the same conclusion. For me, it was easier and more dependable to stick with strings.
Pete
-
but how would you get back to a number? You pretty much have to go to a string form sometime in process, so you could pass and print numbers.
This is handled by the bni_tostr$ and bni_fromstr functions.