Search code examples
vb.netalgorithmperformancebigintegercode-structure

Is there an algorithm, to find values ​of a polynomial with big integers, quickly without loops?


For example, if I want to find

1085912312763120759250776993188102125849391224162 = a^9+b^9+c^9+d the code needs to brings

a=3456

b=78525

c=217423

d=215478

I do not need specific values, only that they comply with the fact that a, b and c have 6 digits at most and d is as small as possible.

Is there a quick way to find it?

I appreciate any help you can give me.

I have tried with nested loops but it is extremely slow and the code gets stuck.

Any help in VB or other code would be appreciated. I think the structure is more important than the language in this case

Imports System.Numerics
Public Class Form1
    Private Sub Button1_Click(sender As Object, e As EventArgs) Handles Button1.Click
        Dim Value As BigInteger = BigInteger.Parse("1085912312763120759250776993188102125849391224162")
        Dim powResult As BigInteger
        Dim dResult As BigInteger
        Dim a As Integer
        Dim b As Integer
        Dim c As Integer
        Dim d As Integer

        For i = 1 To 999999
            For j = 1 To 999999
                For k = 1 To 999999
                    powResult = BigInteger.Add(BigInteger.Add(BigInteger.Pow(i, 9), BigInteger.Pow(j, 9)), BigInteger.Pow(k, 9))
                    dResult = BigInteger.Subtract(Value, powResult)
                    If Len(dResult.ToString) <= 6 Then
                        a = i
                        b = j
                        c = k
                        d = dResult
                        RichTextBox1.Text = a & " , " & b & " , " & c & " , " & d
                        Exit For
                        Exit For
                        Exit For
                    End If
                Next
            Next
        Next
    End Sub
End Class

UPDATE

I wrote the code in vb. But with this code, a is correct, b is correct but c is incorrect, and the result is incorrect.

a^9 + b^9 + c^9 + d is a number bigger than the initial value.

The code should brings

a= 217423

b= 78525

c= 3456

d= 215478

Total Value is ok= 1085912312763120759250776993188102125849391224162

but code brings

a= 217423

b= 78525

c= 65957

d= 70333722607339201875244531009974

Total Value is bigger and not equal=1085935936469985777155428248430866412402362281319

Whats i need to change in the code to make c= 3456 and d= 215478?

the code is

Imports System.Numerics
Public Class Form1
    Private Function pow9(x As BigInteger) As BigInteger
        Dim y As BigInteger
        y = x * x  ' x^2
        y *= y   ' x^4
        y *= y   ' x^8
        y *= x   ' x^9
        Return y
    End Function

    Private Sub Button1_Click(sender As Object, e As EventArgs) Handles Button1.Click
        Dim a, b, c, d, D2, n As BigInteger
        Dim aa, bb, cc, dd, ae As BigInteger
        D2 = BigInteger.Parse("1085912312763120759250776993188102125849391224162")
        'first solution so a is maximal
        d = D2
        'a = BigIntegerSqrt(D2)
        'RichTextBox1.Text = a.ToString
        For a = 1 << ((Convert.ToInt32(Math.Ceiling(BigInteger.Log(d, 2))) + 8) / 9) To a > 0 Step -1
            If (pow9(a) <= d) Then
                d -= pow9(a)
                Exit For
            End If
        Next
        For b = 1 << ((Convert.ToInt32(Math.Ceiling(BigInteger.Log(d, 2))) + 8) / 9) To b > 0 Step -1
            If (pow9(b) <= d) Then
                d -= pow9(b)
                Exit For
            End If
        Next
        For c = 1 << ((Convert.ToInt32(Math.Ceiling(BigInteger.Log(d, 2))) + 8) / 9) To c > 0 Step -1
            If (pow9(c) <= d) Then
                d -= pow9(c)
                Exit For
            End If
        Next

        ' minimize d
        aa = a
        bb = b
        cc = c
        dd = d
        If (aa < 10) Then
            ae = 0
        Else
            ae = aa - 10
        End If

        For a = aa - 1 To a > ae Step -1 'a goes down few iterations
            d = D2 - pow9(a)
            For n = 1 << ((Convert.ToInt32(Math.Ceiling(BigInteger.Log(d, 2))) + 8) / 9) To b < n 'b goes up
                If (pow9(b) >= d) Then
                    b = b - 1
                    d -= pow9(b)
                    Exit For
                End If
            Next
            For c = 1 << ((Convert.ToInt32(Math.Ceiling(BigInteger.Log(d, 2))) + 8) / 9) To c > 0 Step -1 'c must be search fully
                If pow9(c) <= d Then
                    d -= pow9(c)
                    Exit For
                End If
            Next
            If d < dd Then 'remember better solution
                aa = a
                bb = b
                cc = c
                dd = d
            End If
            If a < ae Then
                Exit For
            End If
        Next
        a = aa
        b = bb
        c = cc
        d = dd
        ' a,b,c,d is the result
        RichTextBox1.Text = D2.ToString
        Dim Sum As BigInteger
        Dim a9 As BigInteger
        Dim b9 As BigInteger
        Dim c9 As BigInteger
        a9 = BigInteger.Pow(a, 9)
        b9 = BigInteger.Pow(b, 9)
        c9 = BigInteger.Pow(c, 9)
        Sum = BigInteger.Add(BigInteger.Add(BigInteger.Add(a9, b9), c9), d)
        RichTextBox2.Text = Sum.ToString
        Dim Subst As BigInteger
        Subst = BigInteger.Subtract(Sum, D2)
        RichTextBox3.Text = Subst.ToString
    End Sub
End Class

Solution

  • In case a,b,c,d might be zero I got an Idea for fast and simple solution:

    First something better than brute force search of a^9 + d = x so that a is maximal (that ensures minimal d)...

    1. let d = 1085912312763120759250776993188102125849391224162

    2. find max value a such that a^9 <= d

      this is simple as we know 9th power will multiply the bitwidth of operand 9 times so the max value can be at most a <= 2^(log2(d)/9) Now just search all numbers from this number down to zero (decrementing) until its 9th power is less or equal to x. This value will be our a.

      Its still brute force search however from much better starting point so much less iterations are required.

    3. We also need to update d so let

      d = d - a^9
      

    Now just find b,c in the same way (using smaller and smaller remainder d)... these searches are not nested so they are fast ...

    b^9 <= d; d-=b^9;
    c^9 <= d; c-=b^9;
    

    To improve speed even more you can hardcode the 9th power using power by squaring ...

    This will be our initial solution (on mine setup it took ~200ms with 32*8 bits uints) with these results:

    x = 1085912312763120759250776993188102125849391224162
        1085912312763120759250776993188102125849391224162 (reference)
    a = 217425
    b = 65957
    c = 22886
    d = 39113777348346762582909125401671564
    

    Now we want to minimize d so simply decrement a and search b upwards until still a^9 + b^9 <= d is lower. Then search c as before and remember better solution. The a should be search downwards to meet b in the middle but as both a and bhave the same powers only few iterations might suffice (I used 50) from the first solution (but I have no proof of this its just my feeling). But still even if full range is used this has less complexity than yours as I have just 2 nested fors instead of yours 3 and they all are with lower ranges...

    Here small working C++ example (sorry do not code in BASIC for decades):

    //---------------------------------------------------------------------------
    typedef uint<8> bigint;
    //---------------------------------------------------------------------------
    bigint pow9(bigint &x)
        {
        bigint y;
        y=x*x;  // x^2
        y*=y;   // x^4
        y*=y;   // x^8
        y*=x;   // x^9
        return y;
        }
    //---------------------------------------------------------------------------
    void compute()
        {
        bigint a,b,c,d,D,n;
        bigint aa,bb,cc,dd,ae;
        D="1085912312763120759250776993188102125849391224162";
        // first solution so a is maximal
        d=D;
        for (a=1<<((d.bits()+8)/9);a>0;a--) if (pow9(a)<=d) break; d-=pow9(a);
        for (b=1<<((d.bits()+8)/9);b>0;b--) if (pow9(b)<=d) break; d-=pow9(b);
        for (c=1<<((d.bits()+8)/9);c>0;c--) if (pow9(c)<=d) break; d-=pow9(c);
    
        // minimize d
        aa=a; bb=b; cc=c; dd=d;
        if (aa<50) ae=0; else ae=aa-50;
        for (a=aa-1;a>ae;a--)       // a goes down few iterations
            {
            d=D-pow9(a);
            for (n=1<<((d.bits()+8)/9),b++;b<n;b++) if (pow9(b)>=d) break; b--; d-=pow9(b); // b goes up
            for (c=1<<((d.bits()+8)/9);c>0;c--) if (pow9(c)<=d) break; d-=pow9(c);          // c must be search fully
            if (d<dd)               // remember better solution
                {
                aa=a; bb=b; cc=c; dd=d;
                }
            }
        a=aa; b=bb; c=cc; d=dd; // a,b,c,d is the result
        }
    //-------------------------------------------------------------------------
    

    The function bits() just returns number of occupied bits (similar to log2 but much faster). Here final results:

    x = 1085912312763120759250776993188102125849391224162
        1085912312763120759250776993188102125849391224162 (reference)
    a = 217423
    b = 78525
    c = 3456
    d = 215478
    

    It took 1689.651 ms ... As you can see this is much faster than yours however I am not sure with the number of search iterations while fine tuning ais OK or it should be scaled by a/b or even full range down to (a+b)/2 which will be much slower than this...

    One last thing I did not bound a,b,c to 999999 so if you want it you just add if (a>999999) a=999999; statement after any a=1<<((d.bits()+8)/9)...

    [Edit1] adding binary search

    Ok now all the full searches for 9th root (except of the fine tunnig of a) can be done with binary search which will improve speed a lot more while ignoring bigint multiplication complexity leads to O(n.log(n)) against your O(n^3)... Here updated code (will full iteration of a while fitting so its safe):

    //---------------------------------------------------------------------------
    typedef uint<8> bigint;
    //---------------------------------------------------------------------------
    bigint pow9(bigint &x)
        {
        bigint y;
        y=x*x;  // x^2
        y*=y;   // x^4
        y*=y;   // x^8
        y*=x;   // x^9
        return y;
        }
    //---------------------------------------------------------------------------
    bigint binsearch_max_pow9(bigint &d)    // return biggest x, where x^9 <= d, and lower d by x^9
        {                                   // x = floor(d^(1/9)) , d = remainder
        bigint m,x;
        for (m=bigint(1)<<((d.bits()+8)/9),x=0;m.isnonzero();m>>=1)
         { x|=m; if (pow9(x)>d) x^=m; }
        d-=pow9(x);
        return x;
        }
    //---------------------------------------------------------------------------
    void compute()
        {
        bigint a,b,c,d,D,n;
        bigint aa,bb,cc,dd;
        D="1085912312763120759250776993188102125849391224162";
        // first solution so a is maximal
        d=D;
        a=binsearch_max_pow9(d);
        b=binsearch_max_pow9(d);
        c=binsearch_max_pow9(d);
        // minimize d
        aa=a; bb=b; cc=c; dd=d;
        for (a=aa-1;a>=b;a--)       // a goes down few iterations
            {
            d=D-pow9(a);
            for (n=1<<((d.bits()+8)/9),b++;b<n;b++) if (pow9(b)>=d) break; b--; d-=pow9(b); // b goes up
            c=binsearch_max_pow9(d);
            if (d<dd)               // remember better solution
                {
                aa=a; bb=b; cc=c; dd=d;
                }
            }
        a=aa; b=bb; c=cc; d=dd;     // a,b,c,d is the result
        }
    //-------------------------------------------------------------------------
    

    function m.isnonzero() is the same as m!=0 just faster... The results are the same as above code but the time duration is only 821 ms for full iteration of a which would be several thousands seconds with previous code.

    I think except using some polynomial discrete math trick I do not know of there is only one more thing to improve and that is to compute consequent pow9 without multiplication which will boost the speed a lot (as bigint multiplication is slowest operation by far) like I did in here:

    but I am too lazy to derive it...