Announcement

Collapse
No announcement yet.

arbitrary length number math

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • arbitrary length number math

    I've been writing some library functions to allow calculations (addition, subtraction, multiplication, division) with arbitrary length real (non-integer) numbers. The purpose (at the moment) is to calculate thousands of digit accuracy for a recursive equation like x(i) = x(i-1)^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.

  • #2
    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.

    Comment


    • #3
      > 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.

      "Not my circus, not my monkeys."

      Comment


      • #4
        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

        Comment


        • #5
          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


          • #6
            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

            Comment


            • #7
              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
              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


              • #8
                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,
                Dale

                Comment


                • #9
                  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
                  Last 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


                  • #10
                    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

                    Comment


                    • #11
                      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
                      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


                      • #12
                        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


                        • #13
                          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
                          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


                          • #14
                            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


                            • #15
                              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
                              Last 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


                              • #16
                                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
                                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


                                • #17
                                  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
                                  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


                                  • #18
                                    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,
                                    Dale

                                    Comment


                                    • #19
                                      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


                                      • #20
                                        I may try to translate the decimal mp-float to PB
                                        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

                                        Working...
                                        X