# Findm(n) Get min m s.t. n | m(m+1)
# Findm(n) returns [m,X,Y] with m = smallest number such that n divides m*(m+1),
# and n=X*Y, gcd(X,Y)=1, and either X|m, Y|m+1, or X|m+1, Y|m.

with(combinat); with(NumberTheory);
# Find max e such that p^e divides N:
powA:=proc(N,p) local c,w; w:=N; c:=0;
   while (w mod p) = 0 do c:=c+1; w:=w/p; od; c; end;
# Extended gcd: find d=GCD(a,b) and (s,t) such that s*a-t*b = +-1,
#  then return [a,b,d,s,t]
mygcd2:=proc(a,b) local d,s,t,m;
   d := igcdex(a,b,`s`,`t`);
   if s*t=0 then m:=max(a,b)-1; else m := min(abs(s)*a, abs(t)*b); fi;
   [a,b,d,s,t,m]; end;

Findm := proc(n)
local XX,YY,X,Y,m,mp,S,pflis,explis,nopf,i,j,t1;
global powA,mygcd2;
XX:=1; YY:=n;
if n=1 then return([1,XX,YY]); fi;
m:=n-1;
pflis:=convert(PrimeFactors(n),list);
nopf:=nops(pflis);
# lprint(pflis,nopf);
explis:=Array(1..nopf,0);
for i from 1 to nopf do
   explis[i]:=powA(n,pflis[i]); od:
# lprint(n,pflis,nopf,[seq(explis[j],j=1..nopf)]);
S:=subsets([seq(j,j=1..nopf)]);
while not S[finished] do
t1:=S[nextvalue]();
X:=1; Y:=1;
# lprint(t1);
for i from 1 to nopf do
  if member(i,t1) then X:=X*pflis[i]^explis[i]; else Y:=Y*pflis[i]^explis[i];  fi;
                     od;
mp:=mygcd2(X,Y)[6];
if mp < m then m := mp; XX:=X; YY:=Y; fi;
# lprint("X,Y", X,Y,mp,m);
end do:
[m,XX,YY];
end;
Findm(60); # should return [15, 4, 15]