64 bit arithmetic

[matherp]

Hi all

I need to multiply two unsigned 32 bit integers to get a 64 bit result and then divide the answer by a 32 bit unsigned integer to give a 32 bit result. The nature of the numbers are such that the answer will always be in the 32 bit range.

I'd appreciate any help on how to get started or pointers to any multiple precision libraries that may exist in basic.

Thanks

Peter

[Senacharim]

Were I doing it, I'd pass them to a function which temporarily puts the values into 64 bit unsigned variables, perform the math, and then pass back the result in 32 bit format.

[matherp]

Unless I'm missing something swordfish doesn't have a 64bit datatype - that is precisely the issue

[Senacharim]

Ah, in that case I'd suggest scaled integer math.

Or, if the inputs will always fall within a certain expected range a lookup-table would be faster depending upon your desired degree of accuracy.

Can you give some more specific examples of the code/operation you're attempting to perform/solve?

[matherp]

I need to calculate N=(B * 2^32)/125000000 where B is any integer between 1 and 40,000,000. In this case a lookup won't work.

[Jerry Messina]

Here's a routine adapted from some C code I had. Seems to work.

Code: Select all

(*
/*
   SF port of C code adapted from int64.c

   original code:
   Copyright (c) Stuff, February 2013

   Use of this program, for any purpose, is granted by the author,
   Harry W (1and0), as long as this copyright notice is included in
   the source code or any source code derived from this program.
   The user assumes all responsibility for using this code.
*/
*)
include "system.bas"

public type
    uint8_t = byte,
    uint16_t = word,
    uint32_t = longword,
    int16_t  = integer,
    int32_t  = longint

public structure uint64_t
    b(8) as byte

    L as b(0).asword    // uint16_t
    H as b(2).asword    // uint16_t
    U as b(4).asword    // uint16_t
    M as b(6).asword    // uint16_t

    LW as b(0).aslongword   // uint32_t
    HW as b(4).aslongword   // uint32_t
    MW as b(2).aslongword   // uint32_t
end structure

(*
//--------------------------------------------------------------
  Integer 64/32-bit Unsigned Division With Remainder

    quotient  = dividend / divisor
    remainder = dividend % divisor

    where
        dividend = quotient * divisor + remainder

  Arguments:
    dividend  = 64-bit unsigned integer to be divided
    divisor   = 32-bit unsigned integer being divided by
    remainder = pointer to 32-bit unsigned integer holding the remainder

  Return Values:
    The function returns a 64-bit unsigned integer of the quotient and
    a pointer to a 32-bit unsigned integer of the remainder.

  Method:
  * This division algorithm is known as long division, the same one we learned
    in elementary school.  It shifts gradually from the left to the right end
    of the dividend, subtracting the largest possible multiple of the divisor
    at each stage; the multiples become the digits of the quotient, and the
    final difference is the remainder.  The binary implementation is to
    repeatedly double the quotient.  When the divisor is less than or equal to
    the dividend, subtract the divisor from the dividend and set the least
    significant bit of the quotient to '1'.  Then halve the divisor and repeat
    the shift-and-subtract until 64 bits are computed.  The remaining dividend
    is the remainder.
//--------------------------------------------------------------
*)
public function div6432u(byval dividend as uint64_t, byval divisor as uint32_t, byref remainder as uint32_t) as uint64_t
    dim
        num_bits as uint8_t,
        divisorL as uint32_t,
        quotient as uint64_t

    quotient.HW = 0
    quotient.LW = 0
    remainder = 0

    if (divisor = 0) then
        result = quotient
        exit
    endif

    num_bits = sizeof(uint32_t)*8 + 1
    while (divisor.bits(31) = 0)      // divisor & (1UL<<31)
        divisor = divisor << 1
        num_bits = num_bits+1
    end while
    divisorL = 0

    repeat
        quotient.HW = quotient.HW << 1
        if (quotient.LW.bits(31) = 1) then        // quotient.LW & (1UL<<31)
            quotient.HW.bits(0) = 1
        endif
        quotient.LW = quotient.LW << 1

        if ((divisor < dividend.HW) or
            ((divisor = dividend.HW) and (divisorL <= dividend.LW))) then
            dividend.HW = dividend.HW - divisor
            if (dividend.LW < divisorL) then
                dividend.HW = dividend.HW-1
            endif
            dividend.LW = dividend.LW - divisorL
            quotient.LW.bits(0) = 1
        endif

        divisorL = divisorL >> 1
        if (divisor.bits(0) = 1) then
            divisorL.bits(31) = 1
        endif
        divisor = divisor >> 1

        num_bits = num_bits - 1
    until (num_bits = 0)

    remainder = dividend.LW
    result = quotient
end function

//
//

dim B as uint32_t
dim N as uint64_t

dim remainder as uint32_t
dim b64 as uint64_t

//
// compute N = (B * 2^32)/125000000
// where B = 1 to 40000000
// (same as N = B * 34.359738368)
//
main:

    for B = 0 to 40000000
        b64.HW = B          // b64 = (B * 2^32)
        b64.LW = 0
        N = div6432u(b64, 125000000, remainder)
        nop()
    next

There are other routines, such as 32x32 multiply, but you don't need them for that equation.

I'll port them over and post them too when I get a chance.

[matherp]

Jerry

Many thanks - it works perfectly. I can even round the result based on the remainder.

The calculation is that required to set the frequency for an Analog Devices AD9850 digital synthesizer.

I'd built a version on picaxe but limited to 32 bit resolution which gives about 1:10000 resolution but with your code I can redo and get full resolution.


I, and I'm sure others, would certainly be interested in other 64bit routines if you have the time.

Thanks again and best regards

Peter

[David Barker]

A very useful piece of code - thanks Jerry!

[Jerry Messina]

Here are the other routines I mentioned.

Code: Select all

(*
/*
   SF port of C code adapted from int64.c

   original code:
   Copyright (c) Stuff, February 2013

   Use of this program, for any purpose, is granted by the author,
   Harry W (1and0), as long as this copyright notice is included in
   the source code or any source code derived from this program.
   The user assumes all responsibility for using this code.
*/
*)
include "system.bas"

public type
    uint8_t = byte,
    uint16_t = word,
    int16_t  = integer,
    uint32_t = longword,
    int32_t  = longint

// 64-bit structure
// can hold both signed and unsigned types
public structure uint64_t
    b(8) as byte

    L as b(0).asword    // uint16_t
    H as b(2).asword    // uint16_t
    U as b(4).asword    // uint16_t
    M as b(6).asword    // uint16_t

    LW as b(0).aslongword   // uint32_t
    HW as b(4).aslongword   // uint32_t
    MW as b(2).aslongword   // uint32_t
end structure

(*
//--------------------------------------------------------------
  Integer 32x32-bit Unsigned Multiplication With 64-bit Product

    product = multiplicand * multiplier

  Arguments:
    multiplicand = 32-bit unsigned integer to be multiplied
    multiplier   = 32-bit unsigned integer of multiples

  Return Value:
    The function returns a 64-bit unsigned integer of the product.

  Method:
  * This binary multiplication algorithm is known as peasant multiplication,
    which is basically a shift-and-add algorithm.  For each '1' bit in the
    multiplier, shift the multiplicand an appropriate amount and then sum the
    shifted values.  The implementation is to repeatedly halve the multiplier,
    ignoring the remainder, and repeatedly double the multiplicand.  For each
    multiplier that is odd, add the corresponding multiplicand to obtain the
    product.  Repeat the shift-and-add until the multiplier is zero.
//--------------------------------------------------------------
*)
public function mul3232u(byval multiplicand as uint32_t, byval multiplier as uint32_t) as uint64_t
    dim multiplicandH as uint32_t,
        product as uint64_t

    multiplicandH = 0
    product.HW = 0
    product.LW = 0
    while (multiplier <> 0)
        if (multiplier.bits(0) = 1) then              // multiplier & 1
            product.HW = product.HW + multiplicandH
            product.LW = product.LW + multiplicand
            if (product.LW < multiplicand) then
                product.HW = product.HW + 1
            endif
        endif

        multiplicandH = multiplicandH << 1

        if (multiplicand.bits(31) = 1) then             // multiplicand & (1UL<<31)
            multiplicandH.bits(0) = 1                   // multiplicandH or $01
        endif
        multiplicand = multiplicand << 1

        multiplier = multiplier >> 1
    end while
    result = product
end function

(*
//--------------------------------------------------------------
  Integer 64/32-bit Unsigned Division With Remainder

    quotient  = dividend / divisor
    remainder = dividend % divisor

    where
        dividend = quotient * divisor + remainder

  Arguments:
    dividend  = 64-bit unsigned integer to be divided
    divisor   = 32-bit unsigned integer being divided by
    remainder = pointer to 32-bit unsigned integer holding the remainder

  Return Values:
    The function returns a 64-bit unsigned integer of the quotient and
    a pointer to a 32-bit unsigned integer of the remainder.

  Method:
  * This division algorithm is known as long division, the same one we learned
    in elementary school.  It shifts gradually from the left to the right end
    of the dividend, subtracting the largest possible multiple of the divisor
    at each stage; the multiples become the digits of the quotient, and the
    final difference is the remainder.  The binary implementation is to
    repeatedly double the quotient.  When the divisor is less than or equal to
    the dividend, subtract the divisor from the dividend and set the least
    significant bit of the quotient to '1'.  Then halve the divisor and repeat
    the shift-and-subtract until 64 bits are computed.  The remaining dividend
    is the remainder.
//--------------------------------------------------------------
*)
public function div6432u(byval dividend as uint64_t, byval divisor as uint32_t, byref remainder as uint32_t) as uint64_t
    dim
        num_bits as uint8_t,
        divisorL as uint32_t,
        quotient as uint64_t

    quotient.HW = 0
    quotient.LW = 0
    remainder = 0

    if (divisor = 0) then
        result = quotient
        exit
    endif

    num_bits = sizeof(uint32_t)*8 + 1
    while (divisor.bits(31) = 0)                // divisor & (1UL<<31)
        divisor = divisor << 1
        num_bits = num_bits+1
    end while
    divisorL = 0

    repeat
        quotient.HW = quotient.HW << 1
        if (quotient.LW.bits(31) = 1) then      // quotient.LW & (1UL<<31)
            quotient.HW.bits(0) = 1
        endif
        quotient.LW = quotient.LW << 1

        if ((divisor < dividend.HW) or
            ((divisor = dividend.HW) and (divisorL <= dividend.LW))) then
            dividend.HW = dividend.HW - divisor
            if (dividend.LW < divisorL) then
                dividend.HW = dividend.HW-1
            endif
            dividend.LW = dividend.LW - divisorL
            quotient.LW.bits(0) = 1
        endif

        divisorL = divisorL >> 1
        if (divisor.bits(0) = 1) then
            divisorL.bits(31) = 1
        endif
        divisor = divisor >> 1

        num_bits = num_bits - 1
    until (num_bits = 0)

    remainder = dividend.LW
    result = quotient
end function

(*
//--------------------------------------------------------------
  Integer 32x32-bit Signed Multiplication With 64-bit Product

    product = multiplicand * multiplier

  Arguments:
    multiplicand = 32-bit signed integer to be multiplied
    multiplier   = 32-bit signed integer of multiples

  Return Value:
    The function returns a 64-bit signed integer of the product.
//--------------------------------------------------------------
*)
public function mul3232s(byval multiplicand as int32_t, byval multiplier as int32_t) as uint64_t
    dim product as uint64_t,
        sign as uint8_t

    sign = 0

    if (multiplier < 0) then
        multiplier = -multiplier
        sign = sign or $01
    endif
    if (multiplicand < 0) then
        multiplicand = -multiplicand
        sign = sign + 1
    endif

    product = mul3232u(multiplicand, multiplier)

    if ((sign and $01) <> 0) then
        product.HW = not product.HW
        product.LW = not product.LW
        product.LW = product.LW + 1
        if (product.LW = 0) then
            product.HW = product.HW + 1
        endif
    endif
    result = product
end function

(*
//--------------------------------------------------------------
  Integer 64/32-bit Signed Division With Remainder

    quotient  = dividend / divisor
    remainder = dividend % divisor

    where
        dividend = quotient * divisor + remainder

  Arguments:
    dividend  = 64-bit signed integer to be divided
    divisor   = 32-bit signed integer being divided by
    remainder = pointer to 32-bit signed integer holding the remainder

  Return Values:
    The function returns a 64-bit signed integer of the quotient and
    a pointer to a 32-bit signed integer of the remainder.
//--------------------------------------------------------------
*)
public function div6432s(byval dividend as uint64_t, byval divisor as int32_t, byref remainder as int32_t) as uint64_t
    dim quotient as uint64_t,
        sign as uint8_t

    dim rem as uint32_t     // temp unsigned remainder. see below

    sign = 0

    if (divisor < 0) then
        divisor = -divisor
        sign = sign or $02
    endif
    if (dividend.HW < 0) then
        dividend.HW = not dividend.HW
        dividend.LW = not dividend.LW
        dividend.LW = dividend.LW + 1
        if (dividend.LW = 0) then
            dividend.HW = dividend.HW + 1
        endif
        sign = sign + 3
    endif

    // we need to pass div6432u a ref to a uint32_t, but remainder is an int32_t
    // since we can't cast a reference, pass it a local uint32_t and we'll
    // copy them afterwards
    quotient = div6432u(dividend, divisor, rem)
    remainder = int32_t(rem)        // copy unsigned -> signed var

    if ((sign and $01) <> 0) then
        remainder = -remainder
    endif

    if ((sign and $02) <> 0) then
        quotient.HW = not quotient.HW
        quotient.LW = not quotient.LW
        quotient.LW = quotient.LW + 1
        if (quotient.LW = 0) then
            quotient.HW = quotient.HW + 1
        endif
    endif
    result = quotient
end function

//
//--------------------------------------------------------------
// examples
//--------------------------------------------------------------
//

dim x as uint32_t,
    y as uint32_t,
    remain as uint32_t

dim xi as int32_t,
    yi as int32_t,
    ri as int32_t

dim p as uint64_t,
    q as uint64_t

    x = 1
    y = 1
    p = mul3232u(x, y)
    q = div6432u(p, x, remain)

    xi = 1
    yi = 1
    p = mul3232s(xi, yi)
    q = div6432s(p, xi, ri)

{
// more examples
dim B as uint32_t
dim N as uint64_t

dim remainder as uint32_t
dim b64 as uint64_t

//
// compute N = (B * 2^32)/125000000
// where B = 1 to 40000000
// (same as N = B * 34.359738368)
//
    for B = 0 to 40000000
        b64.HW = B          // b64 = (B * 2^32)
        b64.LW = 0
        N = div6432u(b64, 125000000, remainder)
        nop()
    next
}

I really haven't verified that the signed routines work, so you've been warned.