equation_solutions.sf

#!/usr/bin/ruby

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

# Attempt to find integer solutions to the equation:
#   p^2 + p + c = t * 10^k

# See also:
#   https://www.wolframalpha.com/input?i=p%5E2+%2B+p+%2B+8+%3D+t+*+10%5E10

# p^2 + p + c = t * 10^n
# p^2 + p - 10000000000 t + 8 = 0

# b^2 - 4ac
# a*x^2 + b*x + c = 0
# p^2 + p + (-10000000000*t + 8)

# 1 - 4*(-10000000000*1 + 8)
# -10000000000 floor((x^2 + x + 8)/10000000000) + x^2 + x + 8 = 0

#var x = 2632873031
#var f = 10000000000
#var n = 0

# How to find integer solutions to p and t in the following equation (?):
#   p^2 + p + 8 = t * 10^10
#   p^2 + p + (-10^10*t + 8) = 0

# Solution examples:
#   p = 10^10 n + 1810486343,
#   t = 10^10 n^2 + 3620972687 n + 327786080, n element Z
#
#   p = 10^10 n + 2632873031,
#   t = 10^10 n^2 + 5265746063 n + 693202040, n element Z

# The equation:
#   p^2 + p + 8 = t * 10^k
#
# has solutions of the form:
#
#   p = (10^k * n + x)
#   t = (10^k * n^2 + (2*x + 1)*n + ((x*x + x + 8) / 10^k))
#
# for positive values of x and n.

# a(7) <= 5440486343
# a(8) <= 460510486343
# a(9) <= 2244105527413
# a(9) <= 7144810486343
# a(9) <= 15157702304929
# a(10) <= 24331810486343
# a(10) <= 88324702304929

# New bounds:
#   a(10) <= 17012632873031
#   a(14) <= 99466287423954043
#   a(15) <= 1917321601810486343
#   a(16) <= 6091565756519625353
#   a(17) <= 379430283012423635659
#   a(18) <= 1857717470295105527413
#   a(19) <= 20115100846994093616419
#   a(20) <= 116117030954234923769461
#   a(21) <= 4649486272754125989703301
#   a(22) <= 6060465438321601810486343
#   a(23) <= 7385439332929272937836433
#   a(24) <= 8316715055212188813502662611
#   a(25) <= 29901250515623818002025296521
#   a(26) <= 82180030043271641984960850223

# How to efficiently find the smallest integer value x, such that n is also an integer?
#   (x^2 + x + 8) / 10^8 = n
#   (x^2 + x + 8) = 10^8 * n
#   (x^2 + x + 8) - 10^8 * n = 0
#   x^2 + x + (-10^8*n + 8) = 0

# Answer: by solving the 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)

# n = 10000000000000 m^2 -  3203620972687 m +  256579683416, x = 1601810486343 - 10000000000000 m, m element Z
# n = 10000000000000 m^2 - 18405265746063 m + 8468845179580, x = 9202632873031 - 10000000000000 m, m element Z

# :: Searching with x = 10486343 and k = 1099634

# See also:
#   https://en.wikipedia.org/wiki/Binary_quadratic_form

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=8, max=nil) {

    var f = 10**n

    #for x in (solve_modular_quadratic(1, 1, d, f)) {
    for x in (modular_quadratic_formula(1, 1, d, f) + solve_modular_quadratic(1, 1, d, f) -> sort.uniq) {
    #for x in (modular_quadratic_formula(1, 1, d, f)) {

        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 + p.next_prime)

            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 = 16
var max = 6091565756519625353

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) = 17012632873031
a(11) = 43502632873031
a(12) = 2322601810486343
a(13) = 5470654702304929
a(14) = 99466287423954043
a(15) = 1917321601810486343
a(16) = 6091565756519625353

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

a(17) = 379430283012423635659
a(18) = 1857717470295105527413

2, 523, 2243, 39419, 763031, 37427413, 594527413, 5440486343, 1619625353, 35960850223, 17012632873031, 43502632873031, 2322601810486343, 5470654702304929, 99466287423954043, 1917321601810486343, 6091565756519625353, 379430283012423635659, 1857717470295105527413

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 + x + d == 0 (mod 10^n), where d = q-p.

__END__

var k = 10
var n = 0
var f = 10**7

var good_x = %n[

31, 3031, 13656, 126968, 486343,
126968
486343
513656
873031
1126968
1486343
1513656
1873031
486343,
2873031,

].sort.uniq

var seen = Set()

#~ a(1) <= 701
#~ a(2) <= 650761
#~ a(3) <= 1762031
#~ a(4) <= 763031
#~ a(5) <= 565186343
#~ a(6) <= 4603486343
#~ a(7) <= 5440486343

for x in (good_x) {

    for k in (1..100) {

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

        if (t.is_int) {

            #say ("Found: ", [n, x, p, t])
            #next

            for n in (1..1e5) {

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

                if (p.is_prime && t.is_prime) {

                    var t = (p*p + p.next_prime)
                    var v = valuation(t, 10)

                    if (!seen.has(v) && is_prime(t / 10**v)) {
                        seen << v
                        say "a(#{v}) <= #{p}"
                    }
                }
            }
        }
    }

    #say [p, t]
    #var v = ((x*x + x + 8)/f)
    #say [x, quadratic_formula(f, f, 8*f - f*((p**2 + p + 8)))]
    #say [x, quadratic_formula(f, f, 8*f - f*(t * f))]
    #next
}

__END__

#~ say bsearch_le(2632873031/5, 2632873031*5, {|x|

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

    #~ say "p = #{p}"

    #~ t *= t.de

    #~ x = iquadratic_formula(1, 1, -t * f)
    #~ p = (f * n + x)

    #~ var u = (p*p + p + (-f * t + 8))

    #~ #say [p, t]
    #~ #say [p**2 + p + 8, t * 10**k,  (p**2 + p + 8) / (t * 10**k) -> float]
    #~ #(p**2 + p + 8) - (t * 10**k) <=> 0

    #~ u <=> 0
#~ })

#~ __END__

say bsearch_le(2632873031 / 5, 5*2632873031, {|x|

    var p = (f*n + x)
    var t = (f*n*n + (2*x + 1)*n + floor((x*x + x + 8)/f))
    var u = (p*p + p + (-f*t + 8))

    assert_eq(((f*n + x)**2 + (f*n + x) + (-f*(f*n*n + (2*x + 1)*n + ((x*x + x + 8)/f)) + 8)), 0)
    assert_eq(-f * ((x*x + x + 8) / f) + (x*x + x + 8), 0)

    say [p, t, x, u]

    #-f * floor((x*x + x + 8) / f) + (x*x + x + 8) <=> 0

    u <=> 0
    #-10000000000 * ceil((x**2 + 20000000001*x + 8)/10000000000) + x**2 + 20000000001*x - u <=> -8

    #(x*x + x + 8)/f <=> round(x*x / f)

    #((p**2 + p + 8)/t <=> f)
    #((p**2 + p + 8)/f <=> t)
    #((p*p + p + (-f*t + 8)) <=> 0)
})

#p = 10000000000 n + 2632873031, t = 10000000000 n^2 + 5265746063 n + 693202040, n element Z

#p^2 - (t * 10^10) = -(p + 8)