equation_solutions.sf

#!/usr/bin/ruby

# a(n) is the first prime p such that, with q the next prime, p + q^2 is 10^n times a prime.
# https://oeis.org/A352848

# Known terms:
#   2, 409, 25819, 101119, 3796711, 4160119, 264073519, 2310648079, 165231073519, 9671986711

# Trying to solve the following equation:
#
#   p + (p+c)^2 = 10^n * t
#
# where n is fixed and c = {10, 22, 40, 52, 60, 70}

# The solutions are given by solving the following quadratic modular equation:
#   a*x^2 + b*x + c == 0 (mod m)
#
# where:
#     a = 1
#     b = 1
#     c = -m - d    # where d is the even difference between p and q
#     m = 10^n

# The modular quadratic equation can be solved by finding all the solutions to `t` satisfying:
#   t^2 == D (mod 4*a*m)
#
# where D = b^2 - 4*a*c
#
# then the solutions to `x` are given by:
#   x_1 = (-b + t)/(2*a)
#   x_2 = (-b - t)/(2*a)

# a(7) <= 2310648079
# a(4) <= 1018890648079
# a(8) <= 2225520648079
# a(9) <= 2542620648079
# a(9) <= 2373701060119
# a(9) <= 5929131073519
# a(10) <= 18300671986711

# New bounds:
#   a(16) <= 3288779373987568759
#   a(17) <= 73421283931459964959
#   a(18) <= 3895482305490671986711
#   a(19) <= 78624665472066701060119
#   a(20) <= 1105336353410163225058471
#   a(21) <= 1538369277033631620648079
#   a(22) <= 89010754482305490671986711
#   a(23) <= 177815369277033631620648079
#   a(24) <= 17851980754482305490671986711
#   a(25) <= 69587976634261881432233925799
#   a(26) <= 125171915369277033631620648079
#   a(27) <= 31594928719850170531982385757759

func solve_modular_quadratic(a,b,c,m) {

    var D = (b**2 - 4*a*c).mod(m)

    var solutions = []

    sqrtmod_all(D, m).each {|t|
        for u in (-b + t, -b - t) {
            var x = (u/(2*a))%m

            if ((a*x*x + b*x + c) % m == 0) {
                solutions << x
            }
        }
    }

    return solutions.sort.uniq
}

var MIN_BOUNDS = Hash()

func a(n, d=10, max=nil) {

    var f = 10**n

    #var X = Poly(1)
    #var (a, b, c) = ((X+d)**2 + X -> coeffs.map { .tail }.flip...)

    #for x in (solve_modular_quadratic(a, b, -f + c, f)) {
    #for x in (modular_quadratic_formula(a, b, -f + c, f)) {
    #for x in (modular_quadratic_formula(1, 2*d + 1, d**2, f)) {
    for x in (solve_modular_quadratic(1, 2*d + 1, d**2, f) + modular_quadratic_formula(1, 2*d + 1, d**2, f) -> sort.uniq) {

        defined(max) && (x > max) && break

        if (x.is_even) {
            say ":: Found even x = #{x}. Skipping..."
            next
        }

        say ":: Searching with x = #{x}"

        for k in (0 .. Inf) {

            var p = (f*k + x)
            #var t = (f * n**2 + (2*x + 1)*n + ((x*x + x + 8) / f))

            defined(max) && (p > max) && break

            all_prime(p+d, p) || next

            var u = (p + p.next_prime**2)
            u.remdiv(10).is_prime || next

            var v = u.valuation(10)

            say "a(#{v}) <= #{p}"
            MIN_BOUNDS{v} := p
            MIN_BOUNDS{v} = min(MIN_BOUNDS{v}, p)

            if (v == n) {
                max := p
                max = min(p, max)
            }
        }
    }

    return max
}

var n = 13
var max = 2735490671986711

for d in (2..2000 `by` 2) {
    say "Checking: d = #{d}"
    max = a(n, d, max)
}

say "Minimum found bounds:"
say MIN_BOUNDS

__END__

Confirmed terms:

a(10) = 18300671986711
a(11) = 154590671986711
a(12) = 2237199609971479
a(13) = 2735490671986711
a(14) = 193086838131073519
a(15) = 1529978199609971479
a(16) = 3288779373987568759

Confirmed terms, assuming that Cramer's conjecture is true:

a(17) = 73421283931459964959
a(18) = 3895482305490671986711

For n >= 1, a(n) has the form k * 10^n + x, for some k >= 0, where x is a solution to the modular quadratic equation x^2 + (2*d + 1)*x + d^2 == 0 (mod 10^n), where d = q-p.


2, 409, 25819, 101119, 3796711, 4160119, 264073519, 2310648079, 165231073519, 9671986711, 18300671986711, 154590671986711, 2237199609971479, 2735490671986711, 193086838131073519, 1529978199609971479, 3288779373987568759




__END__