Announcement

**Rod Macia** · 12 May 2022, 10:51 AM

Do you care to share one of your functions, we may spot something, instead of guessing.

I'm assuming you are storing the digits in a string? are you concatenating, concatenation is notoriously slow, using PB Stringbuilder is much much faster.

**Eric Pearson** · 12 May 2022, 12:42 PM

> 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 arbitrary-length BYTE array to store the digits, instead.

**Larry Chickola** · 12 May 2022, 03:03 PM

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 , 0-9
'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

Here's the code for adding or subtracting two big numbers that are stored in my format. Subtraction is done by 9's complement addition, which was very interesting to learn.

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"
'A-B
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"
'B-A
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

**Johan Klassen** · 12 May 2022, 03:45 PM

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?

**Larry Chickola** · 12 May 2022, 04:08 PM

Lenovo Win 10 laptop
OS Name Microsoft Windows 10 Pro
System Type x64-based PC
System SKU LENOVO_MT_20W6_BU_Think_FM_ThinkPad P15s Gen 2i
Processor 11th Gen Intel(R) Core(TM) i7-1185G7 @ 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

**Johan Klassen** · 12 May 2022, 04:31 PM

thank you Larry, my PC Processor Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz, 3600 Mhz, 8 Core(s), 16 Logical Processor(s)
128 GB ram
~~in your opening post you said non-integer bignums but in your code I see only big-integer functions, is floating-point to be implemented later?~~
I see, fixed point
I did my fiddling in another basic and started out as mp floating-point, 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 32-bit
looks like you are multiplying a digit at a time, you can speed thing up a lot by multiplying 8 digits at a time

**Dale Yarker** · 12 May 2022, 09:01 PM

I agree with Johan's suggestion of raising the base to 2^n, and have worked on-and-off 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 "round-tu-it" occurs:
http://www.yarker-dsyc.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,

**Johan Klassen** · 12 May 2022, 09:53 PM

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^32

**Paul Dixon** · 13 May 2022, 06:53 AM

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## = TIMER-t##

'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

**Johan Klassen** · 13 May 2022, 09:45 AM

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 add-with-carry 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 mp-float 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 Newton-Raphson algorithm

**Larry Chickola** · 13 May 2022, 12:12 PM

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.

**Johan Klassen** · 13 May 2022, 12:55 PM

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 stack

**Dale Yarker** · 13 May 2022, 03:10 PM

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,

**Johan Klassen** · 13 May 2022, 03:55 PM

hello Dale

I use unsigned long (32-bit), but for multiplication I use unsigned 64-bit integer to hold the intermediate product, I post here the multiplication loop but it's not in PB
fac1, fac2 are structures of the mp-float
fac3 is an an array of unsigned long 2 times the size of the input mantissa
digit, carry, prod are 64-bit 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

note: I scan for trailing 0's in fac1 and fac2 and set den and num to point to the least-non-zero element
also swap fac1 with fac2 if num<den, also swap den with num, this way the smaller number is at the bottom

**Johan Klassen** · 13 May 2022, 04:03 PM

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 addition

**Johan Klassen** · 13 May 2022, 07:45 PM

some references/documents about multiprecision arithmetic

A Multiple-Precision Division Algorithm http://dmsmith.lmu.build/MComp1996.pdf

Simple Algorithm for Arbitrary-Precision 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...rithmetic-demo

BigInteger by Richard/Stephan https://github.com/stephanbrunker/big_integer

**Dale Yarker** · 14 May 2022, 01:11 AM

Unsigned long would be PB DWORD.

Stucture would be TYPE/END TYPE. The name "mp-float" 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

Cheers,

**Larry Chickola** · 14 May 2022, 08:34 AM

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.

**Johan Klassen** · 14 May 2022, 03:11 PM

I may try to translate the decimal mp-float to PB

Announcement

arbitrary length number math

arbitrary length number math

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment