I've been writing some library functions to allow calculations (addition, subtraction, multiplication, division) with arbitrary length real (noninteger) numbers. The purpose (at the moment) is to calculate thousands of digit accuracy for a recursive equation like x(i) = x(i1)^2 + c. My functions are working well and am having no problem reaching the first 5000 digits of accuracy of the calculation, but of course as the number of digits increases, the calculation time is going up. My multiplication algorithm is an implementation of grade school multiplication, so the number of operations goes up with n^2, n= number of digits. I know there are C++ libraries out there for manipulating arbitrary length numbers with better algorithms like the Karatsuba. It sounds like the HIME library had this type of function. Does anyone know of Powerbasic source code out there that does this? I couldn't find an example with a search of the forums.
Announcement
Collapse
No announcement yet.
arbitrary length number math
Collapse
X

> an implementation of grade school multiplication
...which changes every 10 years or so. Seriously, I learned one way in the 60s, while kids in larger towns learned "New Math". That apparently fell out of favor while I was in college, and was replaced by a different technique that I never learned. My point being that one of the several methods might be less computationally intensive as the one you were taught.
If you are using strings, use an arbitrarylength BYTE array to store the digits, instead.
"Not my circus, not my monkeys."
Comment

I'm not using strings, the digits of the big numbers are stored as long integers (LONG) in an array, with the array index equal to the power of 10 the digit corresponds to.
It's not that my routines are "slow", but eventually i will hit the time limit i am willing to wait to get a result. The only way to improve that would be with better algorithms, like the Karatsuba for multiplication which seems be 85% faster than what my routine can do, just based on the number operations required. I was hoping someone may already have that done.
some performance numbers:
to do 1000 multiplications of.....
numbers with 2000 decimal places  takes 34 sec
numbers with 3000 decimal places  takes 61 sec
numbers with 4000 decimal places  takes 96 sec
numbers with 6000 decimal places  takes 181 sec
Eric  I think that "new math" rearranges how the kids handle the computation but from what i read there are still the same amount of multiplications to get the results as the old fashioned way.
A little more detail of how i store the big numbers is in the comments of this routine which takes a long number and multiplies it by a single digit:
Code:SUB real_x_single_digit(number1() AS LONG, multiplier AS LONG, result() AS LONG) 'Numbers are stored as an array of long integers, with the position in the array corresponding to the power of 10 associated with the digit 'The array goes from the minimum power of ten to the maximum power of 10 plus 1 more place for holding the sign of the number 'positive numbers are designated by +1 in the UBOUND position, negative numbers are designed by 1 in the UBOUND position 'example: minimum power of 10 is 10 'maximum power of 10 is 5 'LBOUND = 10 'UBOUND = 6 'the number to be stored is +452.47251 'array index: 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 'array value: 1 0 0 0 4 5 2 4 7 2 5 1 0 0 0 0 0 'the number to be stored is 39405.39405903 'array index: 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 'array value: 1 0 3 9 4 0 5 3 9 4 0 5 9 0 3 0 0 'number1() is number to be multiplied. 'multiplier = single digit number to multiply the big number by , 09 'result() = the answer to the multiplication LOCAL i AS LONG LOCAL tsum AS LONG LOCAL carryover AS LONG 'note global variables are GLOBAL min10, max10 as LONG 'step through each digit carryover = 0 'the first multiplication has no carryover from a prior operation FOR i = min10 TO max10 tsum = multiplier * number1(i) + carryover 'multiply the digits and add the carryover from the previous digit carryover = tsum\10 'carryover is the value of the 10's place of the sum result(i) = tsum MOD 10 'the digit is the 1's place value of the sum NEXT i 'if carryover >0 of the last digit then have overflowed available digits IF carryover >0 THEN PRINT "OVERFLOW ERROR  real_x_single_digit()" PRINT "Hit a key ..." WAITKEY$ EXIT SUB END IF 'set the sign of the result IF SGN(multiplier) = SGN(number1(max10+1)) THEN result(max10 + 1) = 1 'if both numbers being multiplied have the same sign, the result is a positive number ELSE result(max10 + 1) = 1 'else the result must be a negative number END IF END SUB
Code:SUB real_add_2_bignums(bignum1() AS LONG, bignum2() AS LONG, result() AS LONG, carryflag AS LONG) LOCAL i,j,k AS LONG LOCAL tsum, carryover AS LONG LOCAL overallmax AS LONG carryflag = 0 'deal with one or the other negative (subtraction) IF bignum1(max10+1) = 1 AND bignum2(max10 + 1) = 1 THEN 'print "case bignum1 = 1, bignum2 = 1" 'AB CALL real_subtract_2bignums(bignum1(), bignum2(), result()) EXIT SUB ELSEIF bignum1(max10+1) = 1 AND bignum2(max10 + 1) = 1 THEN 'PRINT "case bignum1 = 1, bignum2 = 1" 'BA CALL real_subtract_2bignums(bignum2(), bignum1(), result()) EXIT SUB END IF 'ok not subtraction, do addition FOR i = min10 TO max10 tsum = bignum1(i) + bignum2(i) + carryover carryover = tsum\10 result(i) = tsum MOD 10 NEXT i 'add in carryover if >0 IF carryover >0 THEN 'if true then exceeding allowable unless doing subtraction with 9's complements carryflag = 1 ' PRINT "OVERFLOW ERROR  adding 2 real BIGNUMS. Exceeded max digits." ' PRINT "Hit a key ..." ' WAITKEY$ ' EXIT SUB END IF 'deal with both positives or both negatives IF (bignum1(max10+1) + bignum2(max10+1)) = 2 THEN 'both postive result(max10 + 1) =1 ELSEIF (bignum1(max10+1) + bignum2(max10+1)) = 2 THEN 'both negative result(max10 + 1) = 1 END IF END SUB SUB real_subtract_2bignums(bignum1() AS LONG, bignum2() AS LONG, result() AS LONG) 'subtract bignum1  bignum2 = result LOCAL i,j,k AS LONG LOCAL carryflag AS LONG 'step 1 find 9'c complement of bignum2 DIM nines_bignum2(min10 TO max10 + 1) AS LONG DIM sum_nines2_bignum1(min10 TO max10 + 1) AS LONG FOR i = min10 TO max10 nines_bignum2(i) = 9  bignum2(i) NEXT i 'PRINT "nines of bignum2 = ";real_print_bignum(nines_bignum2()) 'WAITKEY$ 'now add bignum1 carryflag=0 CALL real_add_2_bignums(nines_bignum2(), bignum1(), sum_nines2_bignum1(), carryflag) 'PRINT "sum of bignum 1 and nines of bignum2 = ";real_print_bignum(sum_nines2_bignum1()) 'PRINT "sign holder =";sum_nines2_bignum1(max10+1) 'WAITKEY$ SELECT CASE carryflag CASE 0 'this is the case where the result is negative (bignum2>bignum1) FOR i = min10 TO max10 result(i) = 9  sum_nines2_bignum1(i) NEXT i result(max10 + 1) = 1 'set the negative EXIT SUB CASE 1 'this is the case where the result is positive 'leave the carryover in the ones position 'set to a positive number result(max10 + 1) = 1 carryflag = 0 CALL real_add_2_bignums(sum_nines2_bignum1(), addone(), result(), carryflag) EXIT SUB END SELECT END SUB
Comment

base 10^n mp arithmetic suffers due to slow shifting but conversion to/from string is easy whereas base 2^n arithmetic is faster but the conversion is slower and more complex
I am curious about your division implementation
[edit]
what are the specs of your PC?as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

Lenovo Win 10 laptop
OS Name Microsoft Windows 10 Pro
System Type x64based PC
System SKU LENOVO_MT_20W6_BU_Think_FM_ThinkPad P15s Gen 2i
Processor 11th Gen Intel(R) Core(TM) i71185G7 @ 3.00GHz, 2995 Mhz, 4 Core(s), 8 Logical Processor(s)
32G ram
1tera byte SSD
I actually haven't gotten around to division as i didn't need it quite yet. Here's the multiplication SUB
Code:SUB real_x_2_bignums(bignum1() AS LONG, bignum2() AS LONG, result() AS LONG) LOCAL i,j,k AS LONG LOCAL tsum, carryover AS LONG LOCAL carryflag, posn, single_line_mult_posn AS LONG DIM runningsum(min10 TO max10+1) AS LONG DIM single_line_mult(min10 TO max10+1) AS LONG 'notes: if max10 = 6, that is 7 digits. cannot exceed sqrt(1e7) for number size or get overflow (for squares) 'can multiply if value of a x b < 1e7 'multiplying 1000 x 1000 decimal numbers gives a result with 2000 decimals. Since i am restrcting numbers to a fixed 'number of decimals, this allows me to not multiply the entire numbers together, and cut the required number of multiplcations 'by 1/2 'signs IF bignum1(max10+1) = bignum2(max10+1) THEN result(max10 + 1) =1 'both negative or both positive yields positive ELSE result(max10 + 1) = 1 'one positive and one negative yields negative END IF posn = 0 'start stepping throught the digits of num2 FOR j = min10 TO max10/2 'if digit to be multiplied in bignum2(j) = 0 then iterate . This saves 10% of all required multiplications by ignoring multiply by zero and 'just skips to the next non zero digit IF bignum2(j) = 0 THEN DECR posn ITERATE FOR END IF single_line_mult_posn = min10 FOR k = posn TO max10/2 'this line is taking advangate of not needing multiply out the entire numbers IF k<min10 THEN single_line_mult(single_line_mult_posn) = 0 'PRINT k;" k<min10 ";" single_line_mult_posn=";single_line_mult_posn INCR single_line_mult_posn ITERATE FOR END IF tsum = bignum1(k) * bignum2(j) + carryover carryover = tsum\10 single_line_mult(single_line_mult_posn) = tsum MOD 10 'PRINT k;" k>=min10 ";" single_line_mult_posn=";single_line_mult_posn INCR single_line_mult_posn NEXT k 'PRINT real_print_bignum(single_line_mult()), bignum2(j),j IF carryover > 0 THEN PRINT "OVERFLOW in MULTIPLY 2 BIGNUMS" EXIT SUB END IF 'add to running sum 'real_add_2_bignums(nines_bignum2(), bignum1(), sum_nines2_bignum1(), carryflag) carryflag = 0 'PRINT "Current digit of bignum2 being being multiplied, (position, value) ",j, bignum2(j) 'PRINT "current value of running sum ";real_print_bignum(result()) 'PRINT "current value of line to be added to running sum ";real_print_bignum(single_line_mult()) real_add_2_bignums(result(), single_line_mult(), result(), carryflag) 'PRINT "new value of running sum after addition ";real_print_bignum(result()) DECR posn 'PRINT 'waitkey$ NEXT j 'waitkey$ 'deal with signs END SUB
Comment

thank you Larry, my PC Processor Intel(R) Core(TM) i99900K CPU @ 3.60GHz, 3600 Mhz, 8 Core(s), 16 Logical Processor(s)
128 GB ram
in your opening post you said noninteger bignums but in your code I see only biginteger functions, is floatingpoint to be implemented later?
I see, fixed point
I did my fiddling in another basic and started out as mp floatingpoint, so in my code multiplication of two mp numbers of the same precision does not double it's length, so 1000 mults of 6000 digits is done in less than 5 seconds in 32bit
looks like you are multiplying a digit at a time, you can speed thing up a lot by multiplying 8 digits at a timeas the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

I agree with Johan's suggestion of raising the base to 2^n, and have worked onandoff on base 2^16 in DWORD arrays. For unsigned integers that is the largest possible base in DWORD because products &h80000000 to &hFFFE0001 (which is &hFFFF * &hFFFF) require all 32 bits. Haven't posted anything because the operation set isn't complete. I need division and square root for prime number work, where negatives are not needed. The upper 16 bits are for "carry/borrow".
Maybe with base 2^8 in your LONG arrays negatives and fractions can be worked out. The number of CPU ops will be like 25 times less than with base 10. Conversion to/from decimal is at the beginning and/or end of all other calculations, if at all.
A page I started for posting code when a "roundtuit" occurs:
http://www.yarkerdsyc.info/D_Notebo...tic/index.html
The links to arithmetic ops are still dead. Two support functions have code posted.
Last note, HIME became unavailable years ago because it included cryptography routines, and the return was not worth the effort to renew the export license.
Cheers,Dale
Comment

in my first implementation I used base 10^8 I then tried base 2^32 which is not 100% ironed out, it still has wrinkles in it but for the most part it's functional
btw, 1000 mults of 6000 digit numbers is less than a second in base 2^32Last edited by Johan Klassen; 13 May 2022, 04:36 AM.as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

For doing low level arithmetic, you really need to learn some ASM. It'll make your job a LOT easier.
The CPU has more appropriate instructions built in to do what you need and it'll be a hundred times faster.
Here's something I wrote years ago to do square roots to many digits.
Here, it'll do the square root of a number to 30,000 digits in about a second.
Code:'PBCC program #REGISTER NONE 'calculate square roots of arbitrary precision 'The number is stored in the array n() and the answer appears in the array a() 'The initial number must be stored as 32 bit values which looks a bit odd when written as bytes,eg, for the sqrt(12345678) use 'n(0)=&h01020304 and n(1)=&h05060708 'the answer appears 1 digit in each array elememt of a(). 'the user is left to sort out the position of the decimal point. FUNCTION PBMAIN t##=TIMER p&=30000 'number of digits to calculate DIM a(p&) AS BYTE, c(p&) AS BYTE, n(p&/4) AS DWORD DIM carry AS LONG,digit AS LONG 'Set n() to the number you want the square root of 'n(0)=&h01020304 :n(1)=&h05060708 'This is the number 12345678 n(0)=&h00020000 'this is the number 2 'n(0)=&h01000000 'this is the number 10 c(2)=1 'the initial subtractor np&=VARPTR(n(0)) 'pointers to the start iof the number.. cp&=VARPTR(c(0)) '..the subtractor ap&=VARPTR(a(0)) '.. and the answer digit =1 'the digit I'm calculating !pushad 'save all registers !mov edi,cp& 'get pointers to the 2 main arrays !mov esi,np& mainlp: 'short cut, do a check of the most significant digit to see if we really need to do the whole subtraction !mov eax,[esi] !sub eax,[edi] !jc skp3 'subtract will fail, so don't waste time doing it! 'subtract !mov ecx,digit 'start the full subtraction, digit/4 words worth !sar ecx,2 !clc 'must be clear at start lp1: 'do the subtraction 4 bytes at a time for speed !mov eax,[esi+ecx*4] 'get part of number !sbb eax,[edi+ecx*4] 'subtract the subtractor !mov ebx,eax 'keep a copy !and eax, &h80808080 'check for individual byte overflows !shr eax,7 !imul eax,&h0f6 'correct for overflows (F6=9) !sub ebx,eax !mov [esi+ecx*4],ebx 'store result !shl eax,1 'keep the carry for the next word !dec ecx 'done them all? !jns lp1 'no, go back for next word !jnc elsejmp 'if carry is clear at the end then the subtract succeeded 'must have gone negative so this digit is finished 'restore n() to what it was before the subtract by adding the number back in.. !mov ecx,digit 'do digit/4 words !sar ecx,2 !clc 'must be clear at start lp2: !mov eax,[esi+ecx*4] 'get a word !adc eax,[edi+ecx*4] 'add in the right bit !add eax, &hf6f6f6f6 'check for overflow between digits !mov ebx,eax !and eax, &h80808080 'correct for overflows !shr eax,7 !imul eax,&h0f6 !sub ebx,eax !mov [esi+ecx*4],ebx 'store result !not eax 'set carry as required (set if eax bit 31 is clear) !shl eax,1 !dec ecx 'all done? !jns lp2 'no, go back for next word skp3: 'update the subtractor !mov eax,digit 'dec current digit !add eax,edi !inc digit 'update digit pointer !xor eax,3 !dec byte ptr [eax] !mov eax,digit 'and inc next digit !add eax,edi !xor eax,3 !inc byte ptr [eax] 'and shift the remaining number !mov ecx,digit 'digit/4 words to do !sar ecx,2 !mov edx,0 !inc ecx shiftlp: !mov eax,[esi+edx*4+4] 'get the next word !shld [esi+edx*4],eax,8 'shift in 1 byte from this word !inc edx 'inc pointer to word !dec ecx 'loop counter !jns shiftlp 'Done them all? !jmp endifjmp 'yes, skip to next digit elsejmp: 'we get here if a subtract succeeded, so we update the answer and the subtractor accordingly !mov ecx,digit 'get which digit to update !add ecx,edi 'offset into subtractor array !xor ecx,3 'invert last 2 bits so BSWAP isn't needed everywhere !mov al,[ecx] 'get current value of this digit !add al,2 'add 2 !cmp al,9 'make sure it didn't overflow !jle skp4 'no, we're OK !mov byte ptr [ecx], 1 'overflow, reset this digit to 1 !mov ecx,digit 'calculate location of next digit !add ecx,edi 'offset into subtractor array !dec ecx 'one digit before !xor ecx,3 'compensate for bytes being swapped !inc byte ptr [ecx] 'update digit !jmp skp5 'finished skp4: !mov [ecx],al skp5: 'update answer !mov eax,ap& 'get pointer to answer array !add eax,digit 'get offset into array of current digit !inc byte ptr [eax] 'increment current digit endifjmp: 'LOOP UNTIL digit =p& !mov eax,p& 'the number of digits we want to calculate !cmp eax,digit 'have we done them all? !jne mainlp 'no, go back for another !popad 'yes, restore all registers t## = TIMERt## 'extract the digits to print them FOR r& = 1 TO p& '256 a$=a$+TRIM$((STR$(a(r&)))) NEXT PRINT "result=" + a$ PRINT "" PRINT "calculation took ";t##;" secs." WAITKEY$ END FUNCTION
Comment

hello Paul
I have found that most of the time gcc with optimization level 2 or 3 does just as well or better than hand crafted asm, some basic compilers translate code to C and then take advantage of the advanced optimization that gcc/g++ offer
granted there are asm functions like addwithcarry that don't exist in C (there may be intrinsic functions available for some) but gcc optimization boils the inefficient C into efficient asm
I just compared the time of your square root to the time of my mpfloat square root and mine was faster, but I must say that in this case the algorithm used probably was the determining factor, I am guessing that your program extracted the square root longhand whereas I used the NewtonRaphson algorithmas the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

Thank you all for the very interesting posts. I have really just started to tinker around in this area of very large number math, can you please clarify some of the terminology you are using? I don't understand when you talk about 2^8 or 2^32 base. does this refer to the number of bits in each array value? and what does each array value represent? one digit of the base 10 number? 4 digits of the base 10 number?
It also sounds like you are doing math in base 2 (binary?). are you saving the decimal values as binary numbers to do operations faster, or converting the entire long number to a binary number?
I think my approach of saving the base 10 digits in LONG arrays has the advantage of simplicity and familiarity, but clearly it is not very fast. I chose LONG vs INTEGER because the PB help said that using LONG numbers is faster and more efficient than any other number. i actually did confirm that, when i use all INTEGER calculations it is 2x more time than LONG.
Comment

I am just an amateur but what I mean by base 10^8 and 2^32 is that for 10^8 each element holds 8 decimal digits and 2^32 holds 32 binary digits per array element, in the first case the arithmetic is done in base 10 but with 8 digits at a time and in the second case it's done in base 2 with 32 bits at a time
using a long array is good in my opinion but depending on how the array is allocated you may be restricted to a certain memory amount, for example if the array uses the stackas the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

Decimal is base 10. First digit holds 0 to 9. One more gives a 1 in second digit. Each digit has a "weight" of 10 times the previous digit. Or, first digit's weight is 1 * 10^0. (0 is exponent), second digit's weight is 1 * 10^1.
Base 2^8 is actually base 256 (2 to eighth power (takes 8 bits) First digit holds 0 to 255. One more gives 1 in the next digit which is 256, or 1 * 256^1. Saying base 2^8 makes clearer that it takes 8 bits to store, and for larger bases is shorter to write. For example base 2^32 can also be written as base 4294967296.
Doing base 10 arithmetic in LONG arrays is okay because 9 * 9 = 81, 81 fits in a LONG then the 8 can be carried to the next LONG. It is inefficient because 9 only uses 4 bits of 31 in a LONG (sign bit is 32nd), and 6 of combinations of the 4 bits aren't used.

Johan, for base 2^32 what integer type are you using? Because the largest base 2^32 number times the largest base 2^32 number takes 64 bits for the product. PB's QUAD is 63 bits and 1 sign bit; many possible products overflow a QUAD.
(4294967295 * 4294967295 = 18446744065119617025 in binary is:
1111111111111111111111111111111000000000000000000000000000000001)
That is why I use base 2^16 for unsigned in DWORDs, the largest possible product will not overflow.

Assembly does make sense ultimately, but making sure "of all the details" in BASIC is a good way to start.
Cheers,Dale
Comment

hello Dale
I use unsigned long (32bit), but for multiplication I use unsigned 64bit integer to hold the intermediate product, I post here the multiplication loop but it's not in PB
fac1, fac2 are structures of the mpfloat
fac3 is an an array of unsigned long 2 times the size of the input mantissa
digit, carry, prod are 64bit unsigned integers
Code:For j=den To 0 Step 1 carry=0 digit=fac2.mantissa[j] For i=num To 0 Step 1 prod=fac3(i+j+1)+digit*fac1.mantissa[i]+carry carry=prod Shr 32'\&H100000000 fac3(i+j+1)=prod '(prod mod &H100000000) Next fac3(j)=carry Next For i=0 To dwords fac1.mantissa[i]=fac3(i) Next
also swap fac1 with fac2 if num<den, also swap den with num, this way the smaller number is at the bottomLast edited by Johan Klassen; 13 May 2022, 07:48 PM.as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

btw, in my testing/debugging I used gcc __builtin_umulll_overflow and __builtin_uaddll_overflow to test for overflow but in all the tests that I did overflow never occurred so I switched to normal multiplication and additionas the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

some references/documents about multiprecision arithmetic
A MultiplePrecision Division Algorithm http://dmsmith.lmu.build/MComp1996.pdf
Simple Algorithm for ArbitraryPrecision Integer Division https://youtu.be/6bpLYxk9TUQ
Arithmetic Operations Beyond Floating Point Number Precision https://ia800403.us.archive.org/10/i.../1009.5911.pdf
Arbitrary Precision Arithmetic https://www.youtube.com/playlist?lis...UNYwKick1LxiWR
Git repo https://github.com/nickelcarbide/arb...rithmeticdemo
BigInteger by Richard/Stephan https://github.com/stephanbrunker/big_integer
as the salmon fish is compelled to go back to it's birthplace to spawn, so comes a day when man is compelled to go back to it's source.. GOD
Comment

Unsigned long would be PB DWORD.
Stucture would be TYPE/END TYPE. The name "mpfloat" sounds like it contains a DOUBLE (52 bit fraction), but use seems DWORD. And, what is/are other member(s)? Only mantissa used in shown code.
PB has signed QUAD, no unsigned 64 bit type. (could be called QWORD if it existed ) So with PB we are 1 bit short of complete base 2^32.
((another reason to look at assembly  ! MUL puts product in EDX (high 32) and EAX (low 32)))
". . . fac3 is an an array of unsigned long 2 times the size of the input mantissa . . ." Do you mean 2 times the number of digits? I have found that upper bound of product array needs to be upper bound of first multiplier plus upper bound of second multiplier plus one if arrays are zero based indexed.
I need to get back to higher priority projects. But first some UNTESTED PB multiply code. Thanks for sharing yours!
Code:' function MultiplyVLU alias "Multiply_VLU" _ (byref Mult1() as dword, _ 'multiplicand byref Mult2() as dword, _ 'the multiplier byref Prod() as dword, _ 'Product export as long '                                       local Mult1HDE, Mult2HDE, ProdHDE as long 'High Digit Exponents local CurMult1DE, CurMult2DE, as long 'Current digit being multiplied local CurProdDE as long 'Current Product destination digit local ProdNext as long 'Carry destination digit Mult1HDE = ubound(Mult1()) Mult2HDE = ubound(Mult2()) '                   check for undimensioned multiplicands  if Mult1HDE = 1 or Mult2HDE = 1 or then function = 1 exit function end if '                 leading 0 fix, and multiply by 0 shortcut  if Mult1(Mult1HDE) = 0) or Mult2(Mult2HDE) = 0 then if Mult1HDE = 0 or Mult2HDE = 0 then 'is multiply by 0 redim Prod(0) exit function else 'is a leading 0 if Mult1(Mult1HDE) = 0 then 'trim pTBTrim0 = varptr(Mult1(0)) pTBTrimH = varptr(Mult1(Mult1HDE)) pTBTrimHDE = varptr(Mult1HDE) gosub FindTrueHDE redim preserve Mult1(Mult1HDE) end if if Mult2(Mult2HDE) = 0 then 'trim pTBTrim0 = varptr(Mult2(0)) pTBTrimH = varptr(Mult2(Mult2HDE) pTBTrimHDE = varptr(Mult2HDE) gosub FindTrueHDE redim preserve Mult2(Mult2HDE) end if end if end if '                                 Multiply  ProdHDE = Multp1HDE + Mult2HDE + 1 ' . . . . . . prepare product variables . . redim Prod(ProductHDE) for CurMult2DE = 0 to Mult2HDE ' . . . . . . . . . . . . . multiplier loop . . for CurMult1DE = 0 to Mult1HDE ' . . . . . . . . . . . multipicand loop . . CurProdDE = CurMult1DE + CurMult2DE 'calc target product digit TmpDigProd = Mult1(CurMult1DE) * Mult2(CurMult2DE) 'multiply digits CarryDig = TmpDigProdDig ProdDig and= &h0000FFFF??? shift right CarryDig, 16 Prod(CurProdDE) += TmpProdDig Prod(CurProdDE + 1) += CarryDig ' carry for CurCarryDE = CurProdDE to CurProdDE + Mult1HDE next next next '                                 Exit Trim  'The product can only have one, or no, leading zero digit. if Prod(ProdHDE) = 0 then redim preserve Prod(ProdHDE  1) '                               Subroutine(s)  exit function 'Don't "fall" into subroutine when done. FindTrueHDE: for pCurTBTrim = pTBTrimH  4 to pTBTrim0 step 4 poke long, pTBTrimHDE, (peek(long, pTBTrimHDE)  1) if peek(long, pCurTBTrim) then return end if next return end function
Dale
Comment

Dale and Johan, thank you for the ideas and reading materials. I will post my routines and some usage examples in the source code forum soon. Then I will start going for more speed, for me it probably makes sense to try working with 8 digit base 10 number chunks before trying other base number systems.
Comment
Comment