 |
 |
| View previous topic :: View next topic |
| Author |
Message |
stma
Joined: 16 Feb 2004
Posts: 26
|
| 32 bit multiplication |
 Posted: Wed Nov 08, 2006 5:00 am |
 |
|
Hi,
Can anyone point me in the right direction for 32 bit multiplication help.
Obviously any 32bit multiplication will result in a 64 bit result.
Does some form of scaling need to be performed to produce a sensible result or is it split into two 32bit variables.
Thanks
Steve |
|
 |
Ttelmah
Guest
|
|
| Posted: Wed Nov 08, 2006 5:51 am |
|
|
CCS, does not include any support for 64bit.
The existing 32bit multiplication routines, will simply overflow if the result is larger than 32bits.
You need to look at some of the standard libraries from MicroChip for example, and translate them yourself into a form that will work on CCS. What you would probably do is declare a structure containing 8 bytes, and 'typedef' this to be a int64. I did this in the past for an int24 type.
Best Wishes |
|
|
stma
Joined: 16 Feb 2004
Posts: 26
|
|
| Posted: Wed Nov 08, 2006 7:59 am |
|
|
Thanks very much for the reply. Will have a look at doing it that way.
Just done some scribblings and need only to multiply two 18 bit numbers initially.
Has anyone any ideas on scaling.
Im not a mathematician so keep it simple!
Thanks
Steve |
|
|
Ken Johnson
Joined: 23 Mar 2006
Posts: 197
Location: Lewisburg, WV
|
|
| Posted: Wed Nov 08, 2006 10:01 am |
|
|
To multiply 2 18-bit numbers, I'd just do c = a * b, where a, b, c are 32-bit. Unsigned would be better than signed, if they're always positive numbers.
Ken |
|
|
Ttelmah
Guest
|
|
| Posted: Wed Nov 08, 2006 4:44 pm |
|
|
Just for anyone who wants to start thinking about multiplying 32bit by 32bit, to give a 64bit result, I have generated the following:
| Code: |
union bits64 {
int8 b[8];
int32 w[2];
};
typedef union bits64 int64;
int64 add64x64(int64 a,int64 b) {
//simple 64 bit addition
#asm
movf b.b[0],0
addwf a.b[0],1
movf b.b[1],0
addwfc a.b[1],1
movf b.b[2],0
addwfc a.b[2],1
movf b.b[3],0
addwfc a.b[3],1
movf b.b[4],0
addwfc a.b[4],1
movf b.b[5],0
addwfc a.b[5],1
movf b.b[5],0
addwfc a.b[6],1
movf b.b[6],0
addwfc a.b[6],1
movf b.b[7],0
addwfc a.b[7],1
#endasm
return a;
}
int64 mult32x32(int32 a,int32 b) {
int64 temp1,temp2;
int8 i;
//zero output
for (i=0;i<8;i++)
temp2.b[i]=0;
//Routine to add two 64bit values stored LSB first
//Not trying to break any records for efficiency, so all in C, and not using
//hardware multiply
//start by moving the int32, into the low half of the int64
temp1.w[0]=b;
//and clear the top half
for (i=4;i<8;i++)
temp1.b[i]=0;
//Now need to work through all 32 bits in the 'a' variable
for (i=0;i<32;i++) {
//if source bit in 'a' is one, perform addition
if (bit_test(a,i))
temp2=add64x64(temp1,temp2);
//rotate the second value here
shift_left(&temp1,8,0);
}
}
|
Called using:
test=mult32x32(123456789,456789);
Where test is another int64 variable, correctly generates the byte pattern:
F9 6F A5 2E 4A 33
Which is the right value for 56393703190521.
Best Wishes |
|
|
Ttelmah
Guest
|
|
| Posted: Thu Nov 09, 2006 3:08 am |
|
|
| Ken Johnson wrote: | To multiply 2 18-bit numbers, I'd just do c = a * b, where a, b, c are 32-bit. Unsigned would be better than signed, if they're always positive numbers.
Ken |
Unfortunately, that doesn't solve the problem.
Multiplying 2 18bit numbers, in the 'worst case', potentially gives a 36bit result. Won't fit into a int32.
Best Wishes |
|
|
Ken Johnson
Joined: 23 Mar 2006
Posts: 197
Location: Lewisburg, WV
|
|
| Posted: Thu Nov 09, 2006 8:07 am |
|
|
18 plus 18 is more than 32 ! right, of course.
I had a funny feeling when I posted that, but it didn't sink in until sometime last night, just before I went to sleep . . .
Sorry 'bout that!
Ken |
|
|
Ttelmah
Guest
|
|
| Posted: Thu Nov 09, 2006 8:24 am |
|
|
It is one of those things that is not obvious, till you think about it, and is made 'worse' by having as standard a '32bit' arithmetic library, which just throws away the bits about 32...
Generally, for most applications, 32bits is 'big enough', but it is a problem sometimes waiting to sneak out and catch you!.
I wrote 32bit libraries years ago, for the Z80, and for those, had to deal with generating a potentially 64bit result for the multiplication. In fact the 'flowchart', is what I implemented in the code I posted. This could be made fairly efficient, if the 'scratch' was made global, so the value didn't have to keep being passed backwards/forwards to the addition routine, but even as posted, it'll work quite well. Not using the hardware multiply function, makes it portable to 16x, and 18x chips.
The 'extra' needed, is the division routine, since this will be required to print the result if required. Not terribly hard, but requires the generation of a 64bit 'compare' as well.
Best Wishes |
|
|
SherpaDoug
Joined: 07 Sep 2003
Posts: 1640
Location: Cape Cod Mass USA
|
|
| Posted: Thu Nov 09, 2006 1:39 pm |
|
|
I am curious, what is the application that rerquires more than 32 bits of resolution? I can maybe imagine some things in astronomy, nuclear physics, and maybe international finances, but none of them would be handled by a PIC! Are you sure you couldn't make do with an 18X14 multiply?
_________________
The search for better is endless. Instead simply find very good and get the job done. |
|
|
Ttelmah
Guest
|
|
| Posted: Thu Nov 09, 2006 5:00 pm |
|
|
Ongoing to my earlier post, though not 'elegant', I have assembled a set of basic routines for 64bit unsigned integer arithmetic.
| Code: |
union bits64 {
int8 b[8];
int32 w[2];
};
typedef union bits64 int64;
int64 scratch64;
int1 BORROW=FALSE;
#bit CARRY=0xFD8.0
//Macro to add two 64bit values. Result into the first
#define M_add64x64(M_a,M_b) \
#asm\
movf M_b.b[0],0\
addwf M_a.b[0],1\
movf M_b.b[1],0\
addwfc M_a.b[1],1\
movf M_b.b[2],0\
addwfc M_a.b[2],1\
movf M_b.b[3],0\
addwfc M_a.b[3],1\
movf M_b.b[4],0\
addwfc M_a.b[4],1\
movf M_b.b[5],0\
addwfc M_a.b[5],1\
movf M_b.b[6],0\
addwfc M_a.b[6],1\
movf M_b.b[7],0\
addwfc M_a.b[7],1\
#endasm
//Macro as above, to subtract two 64bit values.
#define M_sub64x64(M_a,M_b) \
#asm\
movf M_b.b[0],0\
subwf M_a.b[0],1\
movf M_b.b[1],0\
subwfb M_a.b[1],1\
movf M_b.b[2],0\
subwfb M_a.b[2],1\
movf M_b.b[3],0\
subwfb M_a.b[3],1\
movf M_b.b[4],0\
subwfb M_a.b[4],1\
movf M_b.b[5],0\
subwfb M_a.b[5],1\
movf M_b.b[6],0\
subwfb M_a.b[6],1\
movf M_b.b[7],0\
subwfb M_a.b[7],1\
#endasm
//Macros to shift a 64bit value
#define M_shiftleft64(x) shift_left(&x,8,0)
#define M_shiftright64(x) shift_right(&x,8,0)
//Macro to zero a 64bit value
#define M_zero64(x) x.w[0]=0L;x.w[1]=0L
#define M_iszero64(x) ((x.w[0]==0L)&&(x.w[1]==0L))
#define bit64_set(x,i) x.b[i>>3]|=(1<<(i&7))
int64 cast32x64(int32 a) {
//routine to convert a 32bit int to 64bit
int64 temp;
int8 i;
temp.w[0]=a;
//and clear the top half
temp.w[1]=0L;
return temp;
}
int64 I_A_LT_B(int64 a,int64 b) {
//Internal routine to _compare_ two 64 bit values
//Actually performs subtraction, and returns this, with the 'borrow'
//flag in the global variable 'BORROW'. Performs a-b, hence flag is
//set if A less than B
M_sub64x64(a,b);
if (CARRY==0) BORROW=TRUE;
else BORROW=FALSE;
return a;
}
int64 add64x64(int64 a,int64 b) {
//simple 64 bit addition call - returns (a+b)
M_add64x64(a,b);
return a;
}
int64 sub64x64(int64 a,int64 b) {
//64 bit subtraction call as above - returns (a-b)
M_sub64x64(a,b);
return a;
}
int64 mult32x32(int32 a,int32 b) {
//Routine to multiply two 32bit values with a 64bit result
int64 temp2;
int8 i;
//zero output
M_zero64(temp2);
//Not trying to break any records for efficiency, so all in C, and not using
//hardware multiply
//start by moving the int32, into the low half of the int64
scratch64=cast32x64(b);
//Now need to work through all 32 bits in the 'a' variable
for (i=0;i<32;i++) {
//if source bit in 'a' is one, perform addition
if (bit_test(a,i))
M_add64x64(temp2,scratch64);
//rotate the second value here
M_shiftleft64(scratch64);
}
return temp2;
}
int64 div64x64(int64 a,int64 b) {
//Routine to divide 'a' by 'b' in 64bit arithmetic
//returns with result, leaving remainder in the scratch64
int8 bitno=0,ctr;
int64 temp;
if (M_iszero64(b)) {
if (!M_iszero64(a)) {
//need maximum result
temp.w[0]=temp.w[1]=0xFFFFFFFF;
M_zero64(scratch64);
return temp;
}
//Else 0/0=1
M_zero64(temp);
M_zero64(scratch64);
temp.b[0]=1;
return temp;
//return 1
}
M_zero64(temp);
bitno=0;
//Now position divisor to suit find top bit in A, and rotate b till it's
//top bit is in the same position
if (a.w[1]!=0) {
for(ctr=31;!bit_test(a.w[1], ctr);ctr--) ;
//Now rotate b till it's top bit matches.
while (!bit_test(b.w[1], ctr)) {
M_shiftleft64(b);
++bitno;
}
}
else {
for(ctr=31;!bit_test(a.w[0], ctr);ctr--) ;
bitno=ctr;
//Now rotate b till it's top bit matches.
while (!bit_test(b.w[1], ctr)) {
M_shiftleft64(b);
++bitno;
}
}
//bitno now stores how far b had to rotate
while (bitno) {
//Loop for bitno bits, performing subtract, shift & test
scratch64=I_A_LT_B(a,b);
if (!BORROW) {
a=scratch64;
//Set output bit if no borrow
bit64_set(temp,bitno);
}
if (M_iszero64(a)) break;
M_shiftright64(b);
bitno--;
}
scratch64=a;
return temp;
}
|
These are only for the PIC18 (for the PIC16, the location of the borrow bit would need to be recoded, and a set of macros would need to be added to simulate the add with carry, and subtract with borrow instruction which doesn't exist on these chips). There is also a 'caveat' (which is one that exists fairly frequently in the CCS code), that the cast32x64 routine, should always put it's result into a temporary variable, and not be directly used inside a function call. If you code:
test=div64x64(test,cast32x64(45678L));
The compiler will overwrite the top bit of the value handed to the subroutine (Ihave no idea why!), while if you code:
temp=cast32x64(45678L);
test=div64x64(test,temp);
It works fine.
I have slightly tested these on half a dozen 'long' arithmetic operations, so am fairly sure they are close. 
Because the remainder is left in scratch64, from any division, doing decimal output is fairly easy. You can just repeatedly divide by ten, output the resulting digit, and then carry on using the remainder of the last sum.
The multiplication, is now a lot more efficient.
Best Wishes |
|
|
stma
Joined: 16 Feb 2004
Posts: 26
|
|
| Posted: Fri Nov 10, 2006 6:13 am |
|
|
Thanks to all for the replies.
Might get away with 18X14 (or less).
However the 32x32 to 64 routines will come in extremely useful.
On a similar subject, where is a good resource to look into floating point maths in ccs c. Starting from scratch!
Thanks again
Steve |
|
|
kkelis
Joined: 13 Dec 2006
Posts: 5
|
|
| Posted: Mon Mar 10, 2014 4:45 pm |
|
|
I know this is a very old thread, but hopefully i will get some help. Using the above routines provided by Ttelmah i am trying to make a multiplication and then divide the product. So far i was able to do addition,subtractiom, multiplication or a division but one at a time. When i try to divide the multiplication product i get 0 result or the LCD just goes blank. I dont know if this is too much for the PIC to handle...Any help appreciated
Here is my code:
| Code: |
#include <18f452.h>
#use delay(clock=40000000)
#fuses HS,NOWDT,NOLVP,PROTECT,PUT,BROWNOUT
#include "LCD.c"
union bits64 {
int8 b[8];
int32 w[2];
};
typedef union bits64 int64;
int64 scratch64;
int1 BORROW=FALSE;
#bit CARRY=0xFD8.0
//Macro to add two 64bit values. Result into the first
#define M_add64x64(M_a,M_b) \
#asm\
movf M_b.b[0],0\
addwf M_a.b[0],1\
movf M_b.b[1],0\
addwfc M_a.b[1],1\
movf M_b.b[2],0\
addwfc M_a.b[2],1\
movf M_b.b[3],0\
addwfc M_a.b[3],1\
movf M_b.b[4],0\
addwfc M_a.b[4],1\
movf M_b.b[5],0\
addwfc M_a.b[5],1\
movf M_b.b[6],0\
addwfc M_a.b[6],1\
movf M_b.b[7],0\
addwfc M_a.b[7],1\
#endasm
//Macro as above, to subtract two 64bit values.
#define M_sub64x64(M_a,M_b) \
#asm\
movf M_b.b[0],0\
subwf M_a.b[0],1\
movf M_b.b[1],0\
subwfb M_a.b[1],1\
movf M_b.b[2],0\
subwfb M_a.b[2],1\
movf M_b.b[3],0\
subwfb M_a.b[3],1\
movf M_b.b[4],0\
subwfb M_a.b[4],1\
movf M_b.b[5],0\
subwfb M_a.b[5],1\
movf M_b.b[6],0\
subwfb M_a.b[6],1\
movf M_b.b[7],0\
subwfb M_a.b[7],1\
#endasm
//Macros to shift a 64bit value
#define M_shiftleft64(x) shift_left(&x,8,0)
#define M_shiftright64(x) shift_right(&x,8,0)
//Macro to zero a 64bit value
#define M_zero64(x) x.w[0]=0L;x.w[1]=0L
#define M_iszero64(x) ((x.w[0]==0L)&&(x.w[1]==0L))
#define bit64_set(x,i) x.b[i>>3]|=(1<<(i&7))
int64 cast32x64(int32 a) {
//routine to convert a 32bit int to 64bit
int64 temp;
int8 i;
temp.w[0]=a;
//and clear the top half
temp.w[1]=0L;
return temp;
}
int64 I_A_LT_B(int64 a,int64 b) {
//Internal routine to _compare_ two 64 bit values
//Actually performs subtraction, and returns this, with the 'borrow'
//flag in the global variable 'BORROW'. Performs a-b, hence flag is
//set if A less than B
M_sub64x64(a,b);
if (CARRY==0) BORROW=TRUE;
else BORROW=FALSE;
return a;
}
int64 add64x64(int64 a,int64 b) {
//simple 64 bit addition call - returns (a+b)
M_add64x64(a,b);
return a;
}
int64 sub64x64(int64 a,int64 b) {
//64 bit subtraction call as above - returns (a-b)
M_sub64x64(a,b);
return a;
}
int64 mult32x32(int32 a,int32 b) {
//Routine to multiply two 32bit values with a 64bit result
int64 temp2;
int8 i;
//zero output
M_zero64(temp2);
//Not trying to break any records for efficiency, so all in C, and not using
//hardware multiply
//start by moving the int32, into the low half of the int64
scratch64=cast32x64(b);
//Now need to work through all 32 bits in the 'a' variable
for (i=0;i<32;i++) {
//if source bit in 'a' is one, perform addition
if (bit_test(a,i))
M_add64x64(temp2,scratch64);
//rotate the second value here
M_shiftleft64(scratch64);
}
return temp2;
}
int64 div64x64(int64 a,int64 b) {
//Routine to divide 'a' by 'b' in 64bit arithmetic
//returns with result, leaving remainder in the scratch64
int8 bitno=0,ctr;
int64 temp;
if (M_iszero64(b)) {
if (!M_iszero64(a)) {
//need maximum result
temp.w[0]=temp.w[1]=0xFFFFFFFF;
M_zero64(scratch64);
return temp;
}
//Else 0/0=1
M_zero64(temp);
M_zero64(scratch64);
temp.b[0]=1;
return temp;
//return 1
}
M_zero64(temp);
bitno=0;
//Now position divisor to suit find top bit in A, and rotate b till it's
//top bit is in the same position
if (a.w[1]!=0) {
for(ctr=31;!bit_test(a.w[1], ctr);ctr--) ;
//Now rotate b till it's top bit matches.
while (!bit_test(b.w[1], ctr)) {
M_shiftleft64(b);
++bitno;
}
}
else {
for(ctr=31;!bit_test(a.w[0], ctr);ctr--) ;
bitno=ctr;
//Now rotate b till it's top bit matches.
while (!bit_test(b.w[1], ctr)) {
M_shiftleft64(b);
++bitno;
}
}
//bitno now stores how far b had to rotate
while (bitno) {
//Loop for bitno bits, performing subtract, shift & test
scratch64=I_A_LT_B(a,b);
if (!BORROW) {
a=scratch64;
//Set output bit if no borrow
bit64_set(temp,bitno);
}
if (M_iszero64(a)) break;
M_shiftright64(b);
bitno--;
}
scratch64=a;
return temp;
}
//////////////////////////////////////////////////////////////////////
void main()
{
lcd_init();
for(;;)
{
int64 test1,test2,b;
int32 fout,a;
unsigned int32 d;
a=34359738368;
fout=20000000;
b=1000000000;
test1=mult32x32(fout,a);
//test1=150000000;
test2=div64x64(test1,b);
d=test2;
printf(lcd_putc,"\f%lu",d);
delay_ms(100);
}
} |
|
|
|
Ttelmah
Joined: 11 Mar 2010
Posts: 19515
|
|
| Posted: Tue Mar 11, 2014 3:44 am |
|
|
The reason it is hanging, is an 'a' typed where a 'b' should be.
| Code: |
else {
for(ctr=31;!bit_test(b.w[0], ctr);ctr--) ;
bitno=ctr;
//Now rotate b till it's top bit matches.
while (!bit_test(b.w[1], ctr)) {
M_shiftleft64(b);
++bitno;
}
}
|
for the 'else' section in the div64 routine.
As I said, 'close'. Thought I had checked this before posting...
However there are problems (it'll always give zero with your values).
Remember that the biggest number CCS can code directly, is an int32 at 4294967295. Your value into 'a', is above this. The compiler does not check for overflows like this....
Remember also that the compiler does not understand how to write 'into' an int64, so will not clear the high bytes for you. So to load a value, you have to write to the int32 'parts'. So (if you coded 'a' as an int64:
| Code: |
a.w[0]=0x0L;
a.w[1]=0x8;
|
Will write your 64bit value (0x800000000) into the int64.
While you can put in the smaller value into b, with:
| Code: |
b=cast32x64(1000000000);
|
The point is that the compiler knows nothing about 64 bit values, so does not know how to cast up/down to these. You need to use the supplied cast if you want to do this. Similarly for the returned values, you should access the 32bit values from the union to read these.
Best Wishes |
|
|
Douglas Kennedy
Joined: 07 Sep 2003
Posts: 755
Location: Florida
|
|
| Posted: Tue Mar 11, 2014 7:15 am |
|
|
Well ,
I'll echo a point Sherpadoug made. Why do you need it?
The most accurate physical calculation from Quantum Mechanics can be notated with just 12 significant decimal places. Quantum Mechanics underpins every physical mechanism and measurement.
Now Mathematics in abstraction can use notation of unlimited precision.
If you count all the grains of sand on all the world's beaches you may need integer notation with more than 12 decimal places. With 58 decimal digits you may be able to count all the atoms in the visible universe. 64 bit binary notation has equivalency to 19 decimal significant digit notation.
Often a quest for precision is triggered by the misconception that a result in binary notation that is translated into decimal notation is inaccurate.
The translational difference doesn't exist for integers but is a natural consequence for fractions and irrational numbers.
1/10 is accurate in decimal but in binary notation it can only be approximated just as 1/3 can only be approximated in the decimal system. |
|
|
asmboy
Joined: 20 Nov 2007
Posts: 2128
Location: albany ny
|
|
| Posted: Tue Mar 11, 2014 8:00 am |
|
|
| Quote: |
Why do you need it?
|
cuz it would be handy and faster than the way i do it now.
l example: equation for the ohmic value of
a synchronous(phase sensitive) detector"X,Y,Z"
instrument i'm working on right now:
ohms=(((fullscale_span_value*ADC_count)/2^16*pregain*postgain)*16bit_correction_ref)/frequency_based_correction_value
in an actual instrument the first term can easily be
fullscale*adc_count
100000*65535
and correction ref( in my particular case is 75500 and a typical corr_value is in the range of 50000 to 82000
depending on pre and post gain - fractional ohm results are possible,and all this handily overflows 32bits
thats why we have floats
 |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ 
Search
Register
Profile 
Log in to check your private messages 
Log in
