Author Topic: Stringless Math  (Read 5341 times)

0 Members and 1 Guest are viewing this topic.

Offline luke

  • Administrator
  • Seasoned Forum Regular
  • Posts: 324
    • View Profile
Stringless Math
« 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:
Code: [Select]
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:
Code: [Select]
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
« Last Edit: June 19, 2021, 03:38:09 am by luke »

Offline jack

  • Seasoned Forum Regular
  • Posts: 408
    • View Profile
Re: Stringless Math
« Reply #1 on: June 05, 2021, 11:57:04 am »
you have the right idea, string math is woefully slow and cumbersome

Offline George McGinn

  • Global Moderator
  • Forum Regular
  • Posts: 210
    • View Profile
    • Resume
Re: Stringless Math
« Reply #2 on: June 05, 2021, 01:56:00 pm »
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
____________________________________________________________________
George McGinn
Theoretical/Applied Computer Scientist
Member: IEEE, IEEE Computer Society
Technical Council on Software Engineering
IEEE Standards Association
American Association for the Advancement of Science (AAAS)

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Stringless Math
« Reply #3 on: June 05, 2021, 11:11:59 pm »
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.
« Last Edit: June 05, 2021, 11:17:54 pm by bplus »

Offline luke

  • Administrator
  • Seasoned Forum Regular
  • Posts: 324
    • View Profile
Re: Stringless Math
« Reply #4 on: June 05, 2021, 11:47:59 pm »
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:
Code: [Select]
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#

Code: [Select]
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.
« Last Edit: June 06, 2021, 12:46:25 am by luke »

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Stringless Math
« Reply #5 on: June 06, 2021, 01:23:54 am »
Impressive but not the print out here:
Code: QB64: [Select]
  1. Print result$, Len(result$)
  2.  

35660 digits confirmed on internet.
35660 and what do we see.PNG
* 35660 and what do we see.PNG (Filesize: 5.72 KB, Dimensions: 645x427, Views: 248)
« Last Edit: June 06, 2021, 01:25:08 am by bplus »

Offline luke

  • Administrator
  • Seasoned Forum Regular
  • Posts: 324
    • View Profile
Re: Stringless Math
« Reply #6 on: June 12, 2021, 09:26:20 am »
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 :)
« Last Edit: June 12, 2021, 10:36:53 am by luke »

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Stringless Math
« Reply #7 on: June 12, 2021, 11:04:33 am »
"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?

Offline luke

  • Administrator
  • Seasoned Forum Regular
  • Posts: 324
    • View Profile
Re: Stringless Math
« Reply #8 on: June 12, 2021, 11:36:44 am »
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.

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Stringless Math
« Reply #9 on: June 12, 2021, 01:02:58 pm »
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.

Offline luke

  • Administrator
  • Seasoned Forum Regular
  • Posts: 324
    • View Profile
Re: Stringless Math
« Reply #10 on: June 12, 2021, 08:54:40 pm »
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.

Offline Pete

  • Forum Resident
  • Posts: 2361
  • Cuz I sez so, varmint!
    • View Profile
Re: Stringless Math
« Reply #11 on: June 14, 2021, 10:37:06 pm »
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
Want to learn how to write code on cave walls? https://www.tapatalk.com/groups/qbasic/qbasic-f1/

Offline luke

  • Administrator
  • Seasoned Forum Regular
  • Posts: 324
    • View Profile
Re: Stringless Math
« Reply #12 on: June 14, 2021, 11:10:37 pm »
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.