1

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

Mencey
  • 13
  • 4
  • So, your loop will run up to (10^6)^3 = 10^18 iterations? And each iteration will create a string from a number to check that the number of *digits* is fewer than some quantity, instead of just checking if the *numeric value* is less than some number? Oh, and it's all done on your UI > – djv Mar 25 '22 at 19:21
  • I would make `a` as big as possible without going over. Then do the same with `b` and `c`. Whatever's left is `d`. You could then search around that solution for better solutions. For example, reduce `a` by 1, and recompute `b` and `c`. Or reduce `b` by 1, and recompute `c`. – user3386109 Mar 25 '22 at 19:24
  • Brute force is no good. This is further away from a programming problem than it is to a PhD thesis. Look into [How to find all roots of complex polynomials by Newton’s method, John Hubbard, Dierk Schleicher, Scott Sutherland](http://pi.math.cornell.edu/~hubbard/NewtonInventiones.pdf) to got an idea of what you're dealing with. Perhaps if you had much, much smaller constraints, you could use brute force, but not for this. – djv Mar 25 '22 at 19:29
  • You could also heavily decimate your loops, and find close approximations, saving the curves of the independent variables, then fit those curves each to a 9th order approximation, and iterate the process, bounding subsequent steps' variables to progressively smaller ranges. Use least-squares approximation for "goodness" checks. Does this seem like too much? – djv Mar 25 '22 at 19:33
  • Note `1.08591231276312+048` is about 160 bit integer. – chux - Reinstate Monica Mar 25 '22 at 19:55
  • You are talking about the roots of multivariate polynomials and it's a very sophisticated topic let alone being 9th degree and 4 variables are involved in general. I have read even 3rd degree of two varibale polynomial is a very hard candy. What problem exactly you want to solve? – Redu Mar 25 '22 at 21:23
  • You *can* do this with brute force. You only need two nested loops to select `a` and `b`. Given `a` and `b`, there is only one value for `c`. That value is `root9(N - a^9 - b^9)`, where `root9(x)` returns the largest integer less than or equal to the ninth root of `x`. Given `abc`, it's easy to compute `d`. The next improvement is to reduce the search range for `a` and `b`. If we assume, without loss of generality, that `a >= b >= c`, then the minimum value of `a` is `root9(N/3)`, and the maximum value of `a` is `root9(N)`. So the range of `a` is [192440, 217425]. Python solves in 20 minutes. – user3386109 Mar 26 '22 at 04:09
  • can be any of `a,b,c,d` zero ? or they must be all non zero? – Spektre Mar 26 '22 at 09:32
  • They can be zero, but a, b and c should have 6 digits at most, and d needs to be the minimum posible. All of them cant be zero at same time – Mencey Mar 26 '22 at 14:45
  • @Mencey then the solution idea in my answer is correct (ignoring bigint math leads to `O(N^2)` against your `O(N^3)`). The only uncertain thing is the limit of final iterations your case took 2 iterations (~250ms) but some that I tried did need 19 iterations ... I set 50 (there surely is some math formula for it but too lazy to dig in deep into it...) just to be more or less sure (~1.6 sec on my 10+ years old PC) ... I even used my own bigint lib so its not very optimized. With state of the art lib and x64 bit code it should be much faster than that... – Spektre Mar 26 '22 at 15:20
  • @Mencey I added `[Edit1]` into my answer ... I added binary search where I could and did the full a tuning search so the result is numerically stable (do not miss possible solution) ... Its much faster than before and even with the full search is less than a second for your input ... – Spektre Mar 27 '22 at 09:36
  • Its amazing. Thanks you very much. I very appreciatte your time to help me. And to everyone helps with ideas and possible solutions. Thanks a lot to everyone. – Mencey Mar 27 '22 at 13:12

2 Answers2

2

[Update]

The below code is an attempt to solve a problem like OP's, yet I erred in reading it.

The below is for 1085912312763120759250776993188102125849391224162 = a^9+b^9+c^9+d^9+e and to minimize e.

Just became too excite about OP's interesting conundrum and read too quick.

I review this more later.


OP's approach is O(N*N*N*N) - slow

Below is a O(N*N*log(N)) one.

Algorithm

Let N = 1,000,000. (Looks like 250,000 is good enough for OP's sum of 1.0859e48.)

Define 160+ wide integer math routines.

Define type: pow9

int x,y,
int160least_t z

Form array pow9 a[N*N] populated with x, y, x^9 + y^9, for every x,y in the [1...N] range.

Sort array on z.

Cost so far O(N*N*log(N).

For array elements indexed [0... N*N/2] do a binary search for another array element such that the sum is 1085912312763120759250776993188102125849391224162

Sum closest is the answer.

Time: O(N*N*log(N))

Space: O(N*N)


Easy to start with FP math and then later get a better answer with crafter extended integer math.

Try with smaller N and total sum targets to iron out implementation issues.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • Space `O(N*N)` works out to be `1e6 * 1e6 * 160/8 = 20 terabytes`. That's a lot of space. – user3386109 Mar 25 '22 at 20:26
  • 1
    @user3386109 Yes, certainly pushing what could be implemented today. Using `double` instead of `int160_t` brings its down to 8 TB. Possible today with a large drive. I suspect with more thought applied, we could eliminate redundant combinations (easy to gain at least 50% reduction.). Further (1e6)^9 is far larger than OP's 1.09e48, so `N` could be smaller. Many ideas to apply. Simply something to get OP started. – chux - Reinstate Monica Mar 25 '22 at 20:33
  • @user3386109 As long as we stick with positive numbers , looks like `N < 270,000` is required. I suspecting that leads to just a few 100G of space. – chux - Reinstate Monica Mar 25 '22 at 20:37
  • Fair enough, a few space saving optimizations probably could get it down to a size that's manageable with a few terabytes of HDD. – user3386109 Mar 25 '22 at 20:37
  • 2
    In cryptography (but applicable for discrete math in general) this algorithm is sometimes called a "meet-in-the-middle" algorithm. – President James K. Polk Mar 25 '22 at 22:30
  • Mencey, For a sum like 1.0e30, smashed together code took 1s for an good answer. A sum like 1.0e34, code took 15s for an good answer. A sum like 1.0e38, code took 133s for an good answer. I suspect 1.1e48 would take ~100,000s. IAC, its a start. With improvements, do-able in reasonably time. – chux - Reinstate Monica Mar 25 '22 at 22:31
  • @chux-ReinstateMonica can you please check my answer if there is any more room for improvement? Discrete polynomial math on bigints is not my thing ... – Spektre Mar 27 '22 at 09:44
  • 1
    @Spektre Your approach certainly and quickly results in a smallish `d`. I do not think it results in absolutely the smallest possible `d` as the optimal solution may have a smaller `a`. But I already facepalmed on this question and need to look at it anew. – chux - Reinstate Monica Mar 27 '22 at 11:16
  • @chux-ReinstateMonica I was afraid of that with my previous version ... however with the full search of `a` I am doing now should that not happen unless some weird solution has `a,b,c` relatively close to each other (so the fine tuning might skip it) ... maybe that can be searched differently just to test its the case ... need some more thinking about it ... and yes this is interesting and fun question that is why I answered too :) – Spektre Mar 27 '22 at 11:21
  • IMO, the `pow(big_int, 9)` is a distraction. I see the task now as more of a [Knapsack problem](https://en.wikipedia.org/wiki/Knapsack_problem) of 3 items with many (1 million) potential "weights". – chux - Reinstate Monica Mar 27 '22 at 11:27
  • Some insight: If `a` is the largest of the 3, certainly `a` is restricted to the range [pow(1.1e48/3, 1/9) ... pow(1.1e48, 1/9)], `b` (next largest) to the range [0...pow(1.1e48 - a^9, 1/9)] and `c` to range [0...pow(1.1e48 - a^9 - b^9, 1/9)]. I still suspect, at worst, a O(n * n * log(n)) effort where `n` is pow(1.1e48, 1/9). – chux - Reinstate Monica Mar 27 '22 at 11:45
  • @chux-ReinstateMonica the `1<<((d.bits()+8)/9)` with decresing remainder `d` is the same in case for binary search, for linear one the 9th root is slightly better but nut by much ... its the same equation ported to log2, exp2 on integer for simplicity and speed – Spektre Mar 27 '22 at 13:22
1

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...

Spektre
  • 49,595
  • 11
  • 110
  • 380
  • I am trying to write the code that you have provided me, but I have some doubts. Although I am not fluent in c++, I am trying to do it in that language. When I put typedef uint<8> bigint; it gives me an error. On the other hand, how would you implement the bits() function and the isnonzero() function? since they also give me an error. Thank you very much – Mencey Mar 28 '22 at 10:13
  • @Mencey ignore that line of code and just exchange all the `bigint` with yours `BigInteger` ... btw `typedef` just creates a datatype with name `bigint` that in real is `uint<8>` just for my convenience so I can change the `<8>` in one place in case I want to change the bit width of mine unsigned integer number ... – Spektre Mar 28 '22 at 10:39
  • Perfect. Now I don't have that error anymore. But the bits() function and the isnonzero() function still give me an error. How could I solve it? Thank you very much – Mencey Mar 28 '22 at 14:55
  • @Mencey if you do not have `x.bits()` you can use `log2(x)` if you do not have `log2(x)` you can use `log(x)/log(2)` however you have to be careful with integer rounding so maybe better would be `(log(x)*3401)>>10` ... as I mentioned `x.isnonzero()` is the same as `x!=0` ... however You should look at your biginteger its most likely you got `bits()` there maybe with slightly different name... You can even implement it with simple bitshift right until x is zero and the result is the number of shifts ... so simple for loop – Spektre Mar 28 '22 at 16:01
  • spektre, i added the code translated to vb to de initial question, because n brings error, in the line for (n=1<<((d.bits()+8)/9),b++;b – Mencey Mar 29 '22 at 15:19
  • @Mencey try to use `n = (1 as BigInteger) << ((d.bits()+8)/9)` or use some variable of type BigInteger loaded with `1` to ensure the result is BigInteger – Spektre Mar 29 '22 at 15:29
  • Spektre, now the code runs ok, but c is bigger and d is bigger. And the total value adding all the values is bigger than D. a^9+b^9+c^9+d is bigger than D. Something goes wrong. needs to be a= 217423 b= 78525 c= 3456 d= 215478 but code brings a= 217423 b= 78525 c= 65957 d= 70333722607339201875244531009974 In the initial question, i added an update with the entire code. Thanks very much – Mencey Mar 29 '22 at 19:29
  • @Mencey the code looks OK to me it makes no sense that `c` is bigger ... The only thing I can think of is maybe you got somewhere off by 1 bug and `d` is negative at some point can you check it ? – Spektre Mar 30 '22 at 08:34