6.3 abc-conjecture (12.9.96)

6.3.1 John B. Cosgrave

In the coming academic year I want to try to explain to my students the so-called ’abc-conjecture’, a conjecture (initially by Oesterle’, and refined by Masser) of the past decade which is rightly considered to be a truly profound one, with many deep consequences (see - for example - Dorian Goldfeld’s article "Beyond the Last Theorem" in the March/April ’96 issue of "The Sciences" (New York), or the article by Fields medallist Alan Baker in the recent centennial issue of the ’Mathematical Gazette’).

Because of the help that I got in connection with my square-free question, and the related ’product’ question, I am now in a position to experiment with the following programme, whose ’meaning’ is simply this: one is trying to find relatively prime values of ’a’ and ’b’ (i.e. igcd(a, b)=1) such that their sum ’c’ has the PROPERTY that the square-free part of a*b*c DIVIDED by ’c’ has a ’small’ value (in the procedure below, ’small’ means ’less than 1’):

 with(numtheory): # this is needed for 'factorset' in the following 
 
 small_Masser:=proc(n)  # the 'Masser' is explained below 
 local a, b;  # one doesn't need an extra local 'c' 
 for a to n do 
 for b to (a-1) do 
 if igcd(a,b)=1 and 
 convert(factorset(a*b*(a+b)),`*`)/(a+b) < 1 
 then print(a, b, convert(factorset(a*b*(a+b)),`*`)/(a+b)) 
 fi od od end:
 

Comment One: My computing facilities don’t allow me to let ’n’ get too large. n=500 took 1452 seconds on my 486. The last of the outputs - for n=500 - was:

                        343, 169, 91/256
 

the SIGNIFICANCE of which you should sense from:

                        343 + 169 = 512, 
namely:                 7^3 + 13^2 = 2^9.
 

They are not, of course, all as structured as that one. Here is another:

                        343, 32, 14/25 
 
which is related to:    343 + 32 = 375, 
 
namely:                 7^3 + 2^5 = 3*5^3
 

Comment Two: the ’c’ of the 'abc-conjecture' is the 'a+b' in the above.

Comment Three: It is a wonderful theorem of David Masser’s that IF the ’1’in line six of the above programme is replaced by ANY small positive quantity, THEN there will always be output PROVIDED one chooses ’n’ to be SUFFICIENTLY LARGE.

(If I had had ’.1’ in place of ’1’ when I executed 'small_Masser(500)', then there would have been no output.)

Comment Four: The ’abc-conjecture’ is ... (the margin is too small to contain it). By the way, one of the very many consequences of the conjecture is Fermat’s famous ’Last Theorem’ (and also K.F.Roth’s legendary theorem on rational approximations to algebraic numbers - for which Roth won the Fields medal).

If any of you are able to suggest ways of improving the speed of the above procedure, I would appreciate hearing from you.

Also, I would appreciate hearing from those of you who have fast computing facilities about the following: HOW BIG did you have to make ’n’ before you got output with the ’1’ in line six replaced by .1, or, better still, replaced by .001 (you probably won’t get an ’n’, though - by Masser’s Theorem - there has to be such an ’n’), or whatever limit you are able to push it to.

6.3.2 Robert Israel

This can be speeded up quite a bit, I think. Note that the test will fail if a+b is squarefree. Since most integers are squarefree, you could do the outer loop over a+b instead of a, checking first for squarefreeness and at the same time getting the factorset of a+b. Instead of factoring a*b*(a+b), take the union of the factors of a, b and a+b. This way, no actual factoring ends up being done in the inner loop, because factorset remembers its previous results. Then instead of dividing by a+b and checking if the resulting fraction is less than 1, avoid division and stay with integers by checking if the product of the factors is less than a+b. Here, then, is my version.

 with(numtheory): 
 new_Masser:=proc(n) 
 local a, b, s, f; 
 for s from 3 to 2*n-1 do 
    f:= factorset(s); 
    if convert(f, `*`) < s then 
      for a from floor(s/2)+1 to min(s-1,n) do 
         b:= s-a; 
         if igcd(a,b) = 1 and 
            convert(`union`(f, factorset(a), factorset(b)),`*`) < s 
         then print(a, b, convert(`union`(f, factorset(a), factorset(b)),`*`)/s) 
         fi 
      od 
    fi 
od 
end:
 

new_Masser(500) in Release 4 ran for me in 42.208 CPU seconds (admittedly, on a faster machine than a 486).

6.3.3 Andrei Broder

On my machine, the code below for \(n=500\), runs about 12 times faster than the original code. I guess that further optimizations are possible. In particular, on the math side, can b really get as large as a-1?

If you plan to show the students your code for various problems, may I suggest to consult a colleague in Computer Science that specializes in algorithms or algorithms engineering? They might help you write faster code and in the process teach the students some "tricks of the trade", such as pre-computation of values, cheap weak tests to avoid expensive exact tests, etc.

small_Masser_1:=proc(n) 
local a, b, i, fset, sqf, t; 
fset := array(1..2*n); 
for i to 2*n do fset[i] := factorset(i) od;  #  Precompute factorsets 
sqf := map(x -> convert(x,`*`), fset); # Precompute square free values 
    for a to n do 
        for b to (a-1) do 
            if igcd(a,b) = 1 and sqf[a]*sqf[b]/(a+b) < 1 then # Second tests avoids expensive test 
              t := convert(fset[a] union fset[b] union fset[a+b], `*`)/(a+b); 
              if t < 1 then print(a, b, t) fi; 
           fi 
        od 
    od 
end;
 

6.3.4 Andrei Broder

Dr. Israel’s code (new_Masser) is about 10% faster than mine (small_Masser_1). Fortunately our ideas can be combined as in comb_Masser (all codes below) which runs more than twice as fast.

BTW, changing the test sqf[a]*sqf[b]/(a+b) < 1 to sqf[a]*sqf[b] < (a+b) in small_Masser_1 (= small_Masser_2) reduced its running time by about 10%.

I never realized that integer comparison is so much cheaper than rationals comparison.

Approximate times (n=500) on my machine under Maple V.3: 
comb_Masser ~ 8.5 
small_Masser_2 ~ 18.1 
new_Masser ~ 18.6 
small_Masser_1 ~ 20.5 
small_Masser (original code) ~ 256.7
 
comb_Masser := 
proc(n) 
local a,b,s,i,fset,sqf; 
    fset := array(1 .. 2*n); 
    for i to 2*n do  fset[i] := factorset(i) od; 
    sqf := map(x -> convert(x,`*`),fset); 
    for s from 3 to 2*n-1 do 
        if sqf[s] < s then 
            for a from floor(1/2*s)+1 to min(s-1,n) do 
                b := s-a; 
                if igcd(a,b) = 1 and sqf[a]*sqf[b] < s and 
                    convert( 
                    (fset[a] union fset[b]) union fset[s],`*`) 
                     < s then 
                    print(a,b,convert( 
                       (fset[a] union fset[b]) union fset[s],`*` 
                       )/s) 
                fi 
            od 
        fi 
    od 
end;
 
small_Masser_2 := 
proc(n) 
local a,b,i,fset,sqf,t; 
    fset := array(1 .. 2*n); 
    for i to 2*n do  fset[i] := factorset(i) od; 
    sqf := map(x -> convert(x,`*`),fset); 
    for a to n do 
        for b to a-1 do 
            if igcd(a,b) = 1 and sqf[a]*sqf[b] < a+b then 
                t := convert( 
                    (fset[a] union fset[b]) union fset[a+b],`*`) 
                    /(a+b); 
                if t < 1 then print(a,b,t) fi 
            fi 
        od 
    od 
end;
 
new_Masser := 
proc(n) 
local a,b,s,f; 
    for s from 3 to 2*n-1 do 
        f := factorset(s); 
        if convert(f,`*`) < s then 
            for a from floor(1/2*s)+1 to min(s-1,n) do 
                b := s-a; 
                if igcd(a,b) = 1 and convert( 
                    `union`(f,factorset(a),factorset(b)),`*`) < 
                    s then 
                    print(a,b,convert( 
                        `union`(f,factorset(a),factorset(b)),`*` 
                        )/s) 
                fi 
            od 
        fi 
    od 
end;
 
small_Masser_1:=proc(n) 
local a, b, i, fset, sqf, t; 
fset := array(1..2*n); 
for i to 2*n do fset[i] := factorset(i) od;  #  Precompute factorsets 
sqf := map(x -> convert(x,`*`), fset); # Precompute square free values 
    for a to n do 
        for b to (a-1) do 
            if igcd(a,b) = 1 and sqf[a]*sqf[b]/(a+b) < 1 then # Second tests avoids expensive test 
              t := convert(fset[a] union fset[b] union fset[a+b], `*`)/(a+b); 
              if t < 1 then print(a, b, t) fi; 
           fi 
        od 
    od 
end;
 
small_Masser := proc(n) 
local a,b; 
    for a to n do 
        for b to a-1 do 
            if igcd(a,b) = 1 and 
                convert(factorset(a*b*(a+b)),`*`)/(a+b) < 1 then 
                print( 
                    a,b,convert(factorset(a*b*(a+b)),`*`)/(a+b)) 
            fi 
        od 
    od 
end;
 

6.3.5 Olaf Becken

> Also, I would appreciate hearing from those of you who have fast computing ...

I found this thread very interesting, especially, how the "collective intellect" works. Now I wish to add my own 5-cent-ideas (sp?).

From sqf[a]*sqf[b]*sqf[s] < s, a+b=s, a>b and igcd(a,b) = 1 one can prove that (in order of importance)

- sqf[a] < a  ( from sqf[a]*1*2 <= sqf[a]*sqf[b]*sqf[s] < s <= 2*a )

- 2*sqf[s] < s - igcd(sqf[a],sqf[s]) = 1 and igcd(sqf[b],sqf[s]) = 1

Using these conditions one can speed up the procedure. I tried

Masser_6th_version :=proc(n) 
local sqf,a,b,s; 
   sqf := array(1 .. 2*n-1,[1,2]); 
   for s from 3 to 2*n-1 do 
       sqf[s] := convert(factorset(s),`*`); 
       if 2*% < s then 
          for a from floor(s/2)+1 to min(s-1,n) do 
             if sqf[a] < a and sqf[a]*sqf[s] < s then 
                 b := s-a; 
                 if igcd(a,b) = 1 and sqf[a]*sqf[b]*sqf[s] < s then 
                        print(a,b,sqf[a]*sqf[b]*sqf[s]/s) 
                 fi 
             fi 
          od 
       fi 
   od 
end;
 

which needs (with \(n=500\)) 9.2 seconds on my computer whereas comb_Masser needs 15.7 seconds.

BTW, the lowest number I got for sqf[a]*sqf[b]*sqf[s]/s is 21/1024 ~ 0.0205 with 3^10+7^3=29*2^11.

6.3.6 Dave Reynolds

Okay, I’ve had a chance to analyze and I see what’s happening. Let \(sqf(x)\) represent the product of the factorset of \(x\).

Hence we are looking for small values of

   sqf(a)*sqf(b)*sqf(a+b)/(a+b) 
 
                    n 
let b = 1, let a = 3  - 1 
 
                 n                      n            n 
let f(n)  = sqf(3  - 1) * sqf(1) * sqf(3  - 1 + 1)/(3  - 1 + 1) 
 
                 n                 n 
          = sqf(3  - 1) * 1 * 3 / 3 
 
                 n         (n - 1) 
          = sqf(3  - 1) / 3

Now we analyze f(2n) from above.

                 2n        (2n - 1) 
    f(2n) = sqf(3  - 1) / 3 
 
        2n         n       n 
since  3  - 1  = (3  - 1)(3  + 1) 
 
                  n       n          (2n - 1) 
    f(2n) = sqf((3  - 1)(3  + 1)) / 3 
       n         n 
Now,  3 - 1 and 3 + 1 may share common factors.

Since both expressions are even, we know they at least share 2 as a common factor.

Hence

          n       n                 n              n 
    sqf((3  - 1)(3  + 1))  <=  sqf(3  - 1) * sqf((3  + 1) / 2)

Therefore

                  n         (n - 1)               n              n 
    f(2n) <= sqf(3  - 1) / 3           *   sqf(((3  + 1) / 2) / 3 
 
                             n              n 
    f(2n) <= f(n)  *  sqf(((3  + 1) / 2) / 3

Since

          n                              n 
    sqf((3 + 1) / 2) is at most  1/2 *  3 + 1, 
 
          n 
         3 + 1 
    and  -----  is very close to 1 for large n, 
           n 
          3

we are basically saying that \(f(2n)\) is less than or equal to roughly 1/2 * f(n). We can always get lucky and obtain even more common factors and make a leap, but halving each time we double n is good enough.

We have a way of calculating arbitrarily small values for \(f(n)\). This demonstrates:

 
                              n, f(n), evalf(f(n)) 
 
 
                                  22 
                              5, ----, .2716049383 
                                  81 
 
                                 1342 
                            10, -----, .06818066352 
                                19683 
 
                               7924510 
                         20, ----------, .006818181816 
                             1162261467 
 
                          13815528930746510 
                     40, -------------------, .003409090909 
                         4052555153018976267 
 
                  83982289439969274611410914889990510 
            80, --------------------------------------, .001704545455 
                49269609804781974438694403402127765867 

6.3.7 Dave Reynolds

I started following this thread very late, so I don’t know if these comments apply.

Factoring is very time consuming. A simple "sieve" technique can generate the sqf array much more efficiently. This is a times savings for moderate ’n’.

However, for fairly large ’n’, the O(n^2) behavior of the secondary loop overshadows all else.

Compare with Masser_6th_version:

daver1 := proc(n) 
local sqf,a,b,s,m; 
   m := 2*n-1; 
   sqf := array(1 .. m); 
   for s to m do 
       sqf[s] := 1; 
   od; 
   for s from 2 to m do 
       if sqf[s] = 1 then 
           a := s; 
           while a <= m do 
               sqf[a] := sqf[a] * s; 
               a := a + s; 
           od 
       fi 
   od; 
   for s from 3 to 2*n-1 do 
       if 2*sqf[s] < s then 
          for a from floor(s/2)+1 to min(s-1,n) do 
             if sqf[a] < a and sqf[a]*sqf[s] < s then 
                 b := s-a; 
                 if igcd(a,b) = 1 and sqf[a]*sqf[b]*sqf[s] < s then 
                        print(a,b,sqf[a]*sqf[b]*sqf[s]/s) 
                 fi 
             fi 
          od 
       fi 
   od 
end;
 

It is also easier to convert to C code.

#include <stdio.h> 
 
unsigned long gcd(m, n) 
unsigned long m, n; 
{ 
    unsigned long r; 
 
    do { 
    r = m % n; 
    m = n; 
    n = r; 
    } while (r != 0); 
 
    return(m); 
} 
 
#define n 5160 
#define m (2 * (n) - 1) 
 
unsigned long sqf[m + 1];  /* C arrays are zero based; ignore 1st element */ 
 
nt main() 
{ 
    unsigned long a, b, s, t, u, v; 
    double best, d; 
 
    for (s = 1; s <= m; s++) 
    sqf[s] = 1; 
 
    for (s = 2; s <= m; s++) 
    if (sqf[s] == 1) 
        for (a = s; a <= m; a += s) 
        sqf[a] *= s; 
 
    best = 1.0; 
    for (s = 3; s <= m; s++) { 
    if ((2 * sqf[s]) < s) { 
        t = s - 1; 
        if (t > n) 
        t = n; 
        for (a = s / 2 + 1;  a <= t; a++) { 
            if ((sqf[a] < a) && ((u = sqf[a] * sqf[s]) < s)) { 
                b = s - a; 
                if ((gcd(a, b) == 1) && ((v = (u * sqf[b])) < s)) { 
                d = (double) v / (double) s; 
                if (d < best) { 
                    printf("%ld, %ld, %ld/%ld = %f\n", a, b, v, s, d); 
                    best = d; 
                   } 
                } 
           } 
        } 
      } 
    } 
}
 

Output of above program:

8, 1, 6/9 = 0.666667 
63, 1, 42/64 = 0.656250 
49, 32, 42/81 = 0.518519 
80, 1, 30/81 = 0.370370 
125, 3, 30/128 = 0.234375 
512, 1, 114/513 = 0.222222 
1024, 5, 210/1029 = 0.204082 
2187, 10, 390/2197 = 0.177515 
2400, 1, 210/2401 = 0.087464 
4374, 1, 210/4375 = 0.048000 
32768, 37, 1110/32805 = 0.033836 
59049, 343, 1218/59392 = 0.020508 
255879, 121, 4290/256000 = 0.016758 
(The next line is invalid. The result of 32 bit overflow!) 234375, 165941, 2546/400316 = 0.006360
 

Anyone with 64 bit integers wish to take it further?

Instead of exhaustive search, I decided to quickly check special cases.

             x        y 
        if  a  + b = c   and gcd(a, b) = 1 then 
 
                                     (y-1) 
           the result is   a*sqf(b)/c
 

If y is sufficiently large, this could be very small.

I decided to simplify this further to:

                     x                   x 
         t + 1 = base     Hence t = (base  - 1) 
 
                               (x-1) 
   This results in  sqf(t)/base
 

Here is some simple Maple code to play with.

with(numtheory): 
sqf:=x->convert(factorset(x),`*`); 
try:=proc(base,max_power) 
  local a,s,t,best; 
  best:=1: 
  s:=base: 
  for i from 2 to max_power do 
    olds:=s: 
    s:=s*base; 
    a:=s-1; 
    t:=sqf(a)/olds: 
    if (t < best) then 
      print(i, t, evalf(t)): 
      best:=t: 
    fi 
  od 
end;
 
> try(3,50); 
 
                              2, 2/3, .6666666667 
 
                                  10 
                              4, ----, .3703703704 
                                  27 
 
                                  22 
                              5, ----, .2716049383 
                                  81 
 
                                  410 
                              8, ----, .1874714220 
                                 2187 
 
                                 1342 
                            10, -----, .06818066352 
                                19683 
 
                               7924510 
                         20, ----------, .006818181816 
                             1162261467 
 
                          13815528930746510 
                     40, -------------------, .003409090909 
                         4052555153018976267 
 
bytes used=1000096, alloc=655240, time=1.17
 
> ifactor(3^10-1); 
 
                                   3     2 
                                (2)  (11)  (61) 
 
> ifactor(3^20-1); 
 
                             4    2     2 
                          (2)  (5)  (11)  (61) (1181) 
 
> ifactor(3^40-1); 
 
                     5    2     2 
                  (2)  (5)  (11)  (41) (61) (1181) (42521761) 
 
#  Hmmmm... this is peculiar! 
 
> ifactor(3^80-1); 
 
    6    2     2 
 (2)  (5)  (11)  (17) (41) (61) (193) (1181) (128653413121) (42521761) (14401) 
 
 
#  Is this a coincidence? 
 
> sqf(3^80-1)/3^79; 
 
                       83982289439969274611410914889990510 
                     -------------------------------------- 
                     49269609804781974438694403402127765867 
 
> evalf(%); 
 
                                 .001704545455