maitag / littleBigInt

little pure haxe cross-platform arbitrary-precision integer library
MIT License
15 stars 1 forks source link

on modolu operation, (in?)finite loop in mul #1

Closed neimanpinchas closed 4 months ago

neimanpinchas commented 5 months ago

When I am trying to run the following code, the code enters a long loop that didn't finish yet.

    var a=BigInt.fromHexString('D5D4921855F54CC4D0EC0E0CBEF8B582A11E959D54243145088DEF36C75886EEBAE84160B64E011F6B9AD301F1D3466E247FC2ACF7B5DF1798D72918C459CF0');

    var b=BigInt.fromHexString('7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFED');
    var c= a%b;
    trace(a,b,c);

What I've observed is that in the following method

    static function divModLong(a:BigInt, b:BigInt):{quotient:BigInt, remainder:BigInt} {        
        var e = b.length - 1;
        var r:BigInt;
        var x:LittleInt = b.get(e);
        var q:BigInt = divFast(a.splitHigh(e), x);

The q variable is being initialized to a huge BigInt with 81904 chunks, so this may take very long to process, I am not so good in math to debug it, but it does not seem to make sense, python and nodejs are doing it in a second.

maitag commented 5 months ago

mhm, that seems really a beasty bug somewhere into loop here: https://github.com/maitag/littleBigInt/blob/master/src/BigInt.hx#L651 ("r" did not shrink for that numbers!) or is it maybe some glitch in "divModLittle()"? I tested much around here https://try.haxe.org/#57394790 (line 775) but not found a solution yet -_-, (the algorithm i am using there is explained here: https://www.youtube.com/watch?v=6bpLYxk9TUQ)

neimanpinchas commented 5 months ago

Hi

Attaching a HTML explanation form the video author https://justinparrtech.com/JustinParr-Tech/an-algorithm-for-arbitrary-precision-integer-division/ (I do not look on videos for religious reasons).

I tried my best to find the bug, like mutating thr R instead of taking the abs value only on check (so the R/A is getting inverted as well, casuing q and qn to rise instead of going down), It wasn'st easy to compare to Justins Magic because of the left and right shifting.

After trying to demysify the magic, I found that his point is to figure out in which leaf of the binary tree (Sorry if I use wrong terminalogy) the correct qoutient is, this recalled me of an algorithm that I figured out in my teens based on the math that my parents tought me as a child, when I wasn't able to attend math in school, Ihave no idea today how these methods are called, what their sources are, they might be common they might be old Jewish traditions...

I think comments are not needed, the code talks for itself, if anything is not clear I will update.

 static function divModLong(a:BigInt, b:BigInt):{quotient:BigInt, remainder:BigInt} {
//PSN div 2006
  //static function divModLong(a:BigInt, b:BigInt):{quotient:BigInt, remainder:BigInt} {
    var abin=a.toBaseString(2);
    var bbin=b.toBaseString(2);
    if (bbin<abin){
      //impossible to express as int
            return {quotient:0, remainder:b};
    }

    trace(abin,bbin);
    var rema=abin.substr(abin.indexOf("1"));//take off leading 0 to save iterations
    trace(rema);
    var quotientbin="";
//first bit of qoutient is result of first bbin.length chrs;
    var divider_bin=abin.substr(0,bbin.length);
    var devider=BigInt.fromBinaryString(divider_bin);

    rema=rema.substr(bbin.length);

    while (true){
      trace(quotientbin,rema,divider_bin);
    //once we got another instance of b or >b we add another 1 to the out;
    if (devider>=b){
            quotientbin+="1";
      divider_bin=(BigInt.fromBinaryString(divider_bin)-b).toBinaryString();
    } else {
            quotientbin+="0";
    }

      if (rema.length==0){
                break;
      }

    divider_bin+=rema.charAt(0);
    rema=rema.substr(1);
    devider=BigInt.fromBinaryString(divider_bin);

    }
    //remainder is ourcollector that didn't reached b;
    var rem_out=if (divider_bin.length>1){
            BigInt.fromBinaryString(divider_bin);
    } else {
            BigInt.fromInt(0);
    }

    trace(a,b,{quotient:BigInt.fromBinaryString(quotientbin), remainder:rem_out},quotientbin);
    return {quotient:BigInt.fromBinaryString(quotientbin), remainder:rem_out};
  }

The maximum number of loops is abin.length-bbin.length, which is the same as Justins max when realizing that the Q length is naturally that same number +-1.

this method is very CPU efficient as integer substraction is cheaper than integer multiplication.

What could be enhanced is that instead of converting the numbers to binstrings, we could make methods that could give us bigints of bit subsets of a bigint, like substr is doing for a string.

I added a small scale example to your link, to verify that it is working.

What is your opinion regarding changing the method?

What is your opinion about changing the algorithm used?

maitag commented 5 months ago

Interesting algorithm, you can try to implement it by my little-chunks as alternative way. Instead of using Strings i am split all up there into 15bit Integer-chunks because its optimized for karatzuba multiplication. (so more easy to fetch 2 integers to multiplicate... can also optimize that array-handling later targetspecific) Was testing also "divModLittle()" again (so bigint divided by only a 15 bit chunky one) but found no bug yet so it seems its maybe into that algorithm from Justin or into my implementation. The bitshifting there is only for faster multiply/divide by "2" and "Qn" is reduced like:

Qn = Q + R / A    ,    Q = (Q + Qn) / 2    
->   Q = (Q + Q + R/A) / 2
->   Q = ( 2*Q + R/A) / 2

Wonder only why i am was using "-" instead of "+" there (to long ago as i am implemented it but i am think i did it for optimization .. doesn't work also by reverting this) really have no idea at now how to fix ... or why "r" is not shrinking for some rare cases ~^

neimanpinchas commented 5 months ago

All I need would be a a BigNum.substr(startbit,endbit) that would give me a bigint out of the this bigint, I am not sure how to use littlechunks, maybe I could try to apply a mask of (2^start)-1 to clear out the left bits, and shift right to clear the end bits,this would be more efficient.

I see that divFast(r.splitHigh(e), x); will return a positive int for a negetive Dividend, this might be intended but will break the algorithm by causing r to increase and not to shrink. I have uncovered it, after I changed the mutable r.setPositive() to a unmutable r.abs() within the while, this may explain why you had to use minus instead of plus in some cases. Also I think we should update r before checking for next iteration (not sure)

r = a - (q * b);
        do { 
      //trace(q);

            if (r != null) {
                q.shiftOneBitLeft();
        q = q + divFast(r.splitHigh(e), x); // problem into divModLittle?
                q.shiftOneBitRight();
        //r.setPositive();
                trace( r, q, divFast(r.splitHigh(e), x) ); // Why "r" is not shrinking ?
                stop++;
            } 
      r = a - (q * b);
        } while (r.abs() >= b && stop < 150);

Thanks for the insight about reduction, actually the mutilpication is useless, we could simplify it forther to.

Q = Q + ((R/A)/2);

or

var diff=divFast(r.splitHigh(e), x);
diff..shiftOneBitRight();
q=q+diff;

Could divFast or divModLittle be adjusted easily to support negative Q?

maitag commented 5 months ago

ah, maybe cos of the sign of new "r" ... mhm, maybe:

do {
    r = a - (q * b);
    if (r != null) {
        if (r.isPositive)
            q = q + divFast(r.splitHigh(e), x)/2;
        else
            q = q - divFast(r.splitHigh(e), x)/2;

        // same as r = r.abs() so no need to check into while-condition
        r.setPositive();
    }
} while (r >= b);

for "q" start value caluculation a signcheck is not need cos "a" and "b" are made positive already before calling the "divModLong()"

maitag commented 5 months ago

mhm, this also not helps.... more and more i am think it is some rounding problem here divFast(r.splitHigh(e), x) what only appears into rare cases of number ( why the original was work at 99% *lol ;:) So will try next time to check the "remainder" there also to round the "quotient" up or down ( best into that "divFast" helper function! )

-> no no no .. bad idea .... (^_^)

I am think it is a side-effect of the "splitHigh()"-optimization, because this one not really copy the whole chunk-arraydata and instead only create new object with new start/end-values on it. So it needs some temporary variable or a ".clone()" ... will check that later !

maitag commented 5 months ago

neimanpinchas, please feel free to add another of yours or the common-ones there (we can switch between "justins" [who ever that is ;:] and our own ones easy by haxecompiler-conditionals in code) Not sure i mean what i am was implemented there 3 years ago ... i m only r o u g h t ly remember [more &more] -> A N Y <- have to be found WAY into the unit-testcode! [4your numberalgorithms neimanpinchas ... https://www.youtube.com/watch?v=W8RE2NyAiJg&t=107s (^_^) \o/

maitag commented 5 months ago

What do you think is a fast and easy to ìmplement one ? http://maitag.de/semmi/stable-diffusion/possie_02/00023-618132019.png

maitag commented 5 months ago

math-proove is need // this time ~_s into A L G O R Y T H M s https://www.youtube.com/watch?v=hK8SgMXOi24

maitag commented 5 months ago

What can be the fastest way to fix ... mhm ... the 7-one? (until we not found the glitch into the justin-ones) -> i mean to keep that lib MATH->100% working again<- (shame on me /o\ into div there)

maitag commented 5 months ago
    o-o    o-o  o-o-o  o-o
   o   o  o        o      o
  o-o-o  o-o   o    o    o-o
 o      o     (_\    o      o
o      o-o     |\     o    o-o

https://www.youtube.com/watch?v=NpQBfvEf1WA

maitag commented 5 months ago

anyway, best would be to implement this one: https://en.wikipedia.org/wiki/Long_division#Algorithm_for_arbitrary_base (where base is our variable bitchunk -> https://github.com/maitag/littleBigInt/blob/master/src/LittleIntChunks.hx#L26 ) sry, can't spend more time at now (but will keep in mind or pls feel free to implement and do a pullR;:)~

maitag commented 5 months ago

neimanpinchas, it was a great stupidity of mine that I so easily trusted this J. Parr with his algorithm (without mathematical proof!) -> lets rewrite it <-

neimanpinchas commented 5 months ago

My local worktree is still a mess, I think I need to cleanup the unneceray in order to make a PR you will appriciate.

I've implemented substring without strings, making having use of bitwise And operator, causing a preformance bottle neck of slow array unshift in javascript, so I changed it to push and reverse, javascript GC is still my bottleneck, so I think to try to change the LittleChunk atom into something more allocation friendly (linked list of Uint8arrays?).

I think that I also have an alternative multiply algorithm that could defeat karatzuba in perf, and does not rely on divmod, it is based on the principal, that every multiplication result is the diffrance of the 2 squares S1=n1-n2/2 and S2 (n1+n2)/2, plus the fact that squaring in binary can be eaisly done by bit shifting in a maximum of nbits squared time over 2, befords law in binary would mean that most binary digits are zeros (I havn't seen this anywhere but I will write an article about it, when I will write the article that demystifys benfords law), pheraps I need to open another ticket about it.

I hav'nt checked the vids, but In a small community as haxe, all libs are appriciated, even if some edge cases are failing, if it is ever fixed thats also good, better later than never.

I will see this week what I could assamble.

maitag commented 5 months ago

The bug was in the multiplication and is still fixed now.
https://github.com/maitag/littleBigInt/blob/master/src/BigInt.hx#L534

Yes, that's how it is sometimes. I had a strange feeling right from the start because the values of the quotient and the remainder made such an unusual jump.

So, now I have to cool my brain down again and I'll do the optimization later. I already have a plan, but first I want to refactor the code and move the targetspecific chunk-list-handling into a new haxe type/file.

(^_^)