[seqfan] Re: More terms for needed for A086748, Numbers m such that when C(2k, k) == 1 (mod m) then k is necessarily even.

[Richard J. Mathar]

Here a Maple support program for A086748.
It shows within a few minutes that all m < 2000
that are in A086748 have a prime divisors 3 or 5 (or both).




# Determine binomial(m,n) modulo the prime p via the Lucas theorem
# @param m the numerator in the binomial
# @param n the denominator in the binomial
# @param p a prime number (validity not tested for speed here!)
# @return binomial(m,n) mod p.
LucasT := proc(m,n,p::integer)
        local a,mp,np,i ;
        if m < n then
                1 ;
        else
                a :=1 ;
                mp := convert(m,base,p) ;
                np := convert(n,base,p) ;
                for i from 1 to nops(np) do
                        a := modp(a* binomial(op(i,mp),op(i,np)),p) ;
                end do:
                a ;
        end if;
end proc:

# Check whether the number m is in A086748
# Returns true if m is in A086748, else false.
# For odd m where the smallest prime divisors is >=7 ("7-rough")
# the program prints an exclamation mark and an odd "certifikate" k
# which solves binomial(2k,k) =1 mod m, i.e. in these cases
# there is a proof that m is not in the sequence.
isA086748 := proc(m)
        local p,k ;
        if type(m,'even') then
                # by Wang comment all entries are odd
                false ;
        elif modp(m,3) = 0 or modp(m,5) = 0 then
                # by Wang comment all entries which are odd multiples of 3 or 5
                # are in the sequence
                true; 
        else
                # p is the smallest prime divisor of m.
                # At this branch in the program (p=2,3,5 alrdy excluded)
                # this means p >= 7.
                p := min(op(numtheory[factorset](m))) ;
                # Determine whether an odd k exists that solves binomial(2k,k) = 1 mod m.
                # (This may be an endless loop if the search fails.)
                # We want to solve binomial(2k,k) = 1+x*m for any x. If p is a divisor of m,
                # a solution k also yields binomial(2k,k) = 1+x'*p, i.e., solving modulo m
                # is sufficient to solve modulo p. In turn, solving modulo p is
                # necessary to solve modulo m. Because the Lucas theorem of
                # computing binomials is modulo primes and fast, we use LucasT
                # as a quick filter before running again the test with the full m.
                for k from 1 by 2 do
                        if LucasT(2*k,k,p) = 1 then
                                if modp( binomial(2*k,k),m) = 1 then
                                        # Echo that we found an odd k which solves binomia(2k,k) = 1 mod m
                                        printf("! %d %d\n",m,k) ;
                                        return false ;
                                end if;
                        end if;
                end do:
        end if;
end proc:

# loop through odd m
for m from 1 by 2 do
        if modp(m,3) = 0 or modp(m,5) = 0 then
		# if m divisible by 3 or 5: not interesting....
                ;
        else
		# TBD: by the Chinese remainder theorem we
		# may only be interested in prime m here (?!)
                isA086748(m) ;
        end if;
end do: