extra/project-euler/066/066.factor

   1 ! Copyright (C) 2023 Giftpflanze.
   2 ! See https://factorcode.org/license.txt for BSD license.
   3 USING: kernel math math.functions math.order
   4 project-euler.common ranges sequences ;
   5 IN: project-euler.066
   6
   7 ! https://projecteuler.net/problem=66
   8
   9 ! DESCRIPTION
  10 ! -----------
  11
  12 ! Consider quadratic Diophantine equations of the form:
  13 ! x² - Dy² = 1
  14
  15 ! For example, when D = 13, the minimal solution in x is
  16 ! 649² - 13 × 180² = 1.
  17
  18 ! It can be assumed that there are no solutions in positive
  19 ! integers when D is square.
  20
  21 ! By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we
  22 ! obtain the following:
  23 ! 3² - 2 × 2² = 1
  24 ! 2² - 3 × 1² = 1
  25 ! 9² - 5 × 4² = 1
  26 ! 5² - 6 × 2² = 1
  27 ! 8² - 7 × 3² = 1
  28
  29 ! Hence, by considering minimal solutions in x for D <= 7, the
  30 ! largest x is obtained when D = 5.
  31
  32 ! Find the value of D <= 1000 in minimal solutions of x for
  33 ! which the largest value of x is obtained.
  34
  35
  36 ! SOLUTION
  37 ! --------
  38
  39 ! https://www.isibang.ac.in/~sury/chakravala.pdf
  40 ! N = D
  41 ! x0 = p0 = [sqrt(N)]
  42 ! q0 = 1
  43 ! m0 = p0² - N
  44 ! x' ≡ -x (mod |m|), x' < sqrt(N) < x'+|m|
  45 ! p' = (px'+Nq)/|m|
  46 ! q' = (p+x'q)/|m|
  47 ! m' = (x'²-N)/m
  48
  49 :: chakravala ( n x p q m -- n x' p' q' m' )
  50     n sqrt ceiling >integer :> upper-bound
  51     upper-bound m abs - :> lower-bound
  52     x neg m abs rem :> reminder
  53     lower-bound reminder - :> distance
  54     reminder distance m abs / ceiling 0 max m abs * + :> x'
  55     n
  56     x'
  57     p x' * n q * + m abs /
  58     p x' q * + m abs /
  59     x' sq n - m / ;
  60
  61 :: minimal-x ( D -- x )
  62     D sqrt floor >integer :> p0
  63     p0 sq D - :> m0
  64     D p0 p0 1 m0
  65     [ dup 1 = ] [ chakravala ] do until 2drop 2nip ;
  66
  67 : euler066 ( -- D )
  68     1000 [1..b] [ perfect-square? ] reject
  69     [ minimal-x ] supremum-by ;
  70
  71 SOLUTION: euler066