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