64 bit arithmetic
Moderators: David Barker, Jerry Messina
64 bit arithmetic
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
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
- Posts: 139
- Joined: Tue Aug 10, 2010 5:19 pm
- Location: Ventura, CA
- Senacharim
- Posts: 139
- Joined: Tue Aug 10, 2010 5:19 pm
- Location: Ventura, CA
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?
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?
Surviving Member
Bermuda Triangle Battalion
from 2026 to 1992
Voted "Most likely to time travel"--Class of 2024.
Bermuda Triangle Battalion
from 2026 to 1992
Voted "Most likely to time travel"--Class of 2024.
-
- Swordfish Developer
- Posts: 1473
- Joined: Fri Jan 30, 2009 6:27 pm
- Location: US
Here's a routine adapted from some C code I had. Seems to work.
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.
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
I'll port them over and post them too when I get a chance.
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
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
- Swordfish Developer
- Posts: 1214
- Joined: Tue Oct 03, 2006 7:01 pm
- Location: Saltburn by the Sea, UK
- Contact:
-
- Swordfish Developer
- Posts: 1473
- Joined: Fri Jan 30, 2009 6:27 pm
- Location: US
Here are the other routines I mentioned.
I really haven't verified that the signed routines work, so you've been warned.
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
}