extra/project-euler/255/255.factor

   1 ! Copyright (c) 2009 Jon Harper.
   2 ! See https://factorcode.org/license.txt for BSD license.
   3 USING: kernel math math.functions project-euler.common ;
   4 IN: project-euler.255
   5
   6 ! https://projecteuler.net/problem=255
   7
   8 ! DESCRIPTION
   9 ! -----------
  10
  11 ! We define the rounded-square-root of a positive integer n as
  12 ! the square root of n rounded to the nearest integer.
  13
  14 ! The following procedure (essentially Heron's method adapted to
  15 ! integer arithmetic) finds the rounded-square-root of n:
  16
  17 !     Let d be the number of digits of the number n.
  18 !     If d is odd, set x_(0) = 2×10^((d-1)⁄2).
  19 !     If d is even, set x_(0) = 7×10^((d-2)⁄2).
  20
  21 !     Repeat: [see URL for figure ]
  22
  23 !     until x_(k+1) = x_(k).
  24
  25 ! As an example, let us find the rounded-square-root of n =
  26 ! 4321. n has 4 digits, so x_(0) = 7×10^((4-2)⁄2) = 70.
  27
  28 !     [ see URL for figure ]
  29
  30 ! Since x_(2) = x_(1), we stop here.
  31
  32 ! So, after just two iterations, we have found that the
  33 ! rounded-square-root of 4321 is 66 (the actual square root is
  34 ! 65.7343137…).
  35
  36 ! The number of iterations required when using this method is
  37 ! surprisingly low. For example, we can find the
  38 ! rounded-square-root of a 5-digit integer (10,000 ≤ n ≤ 99,999)
  39 ! with an average of 3.2102888889 iterations (the average value
  40 ! was rounded to 10 decimal places).
  41
  42 ! Using the procedure described above, what is the average
  43 ! number of iterations required to find the rounded-square-root
  44 ! of a 14-digit number (10^(13) ≤ n < 10^(14))? Give your answer
  45 ! rounded to 10 decimal places.
  46
  47 ! Note: The symbols ⌊x⌋ and ⌈x⌉ represent the floor function and
  48 ! ceiling function respectively.
  49
  50 ! SOLUTION
  51 ! --------
  52
  53 <PRIVATE
  54
  55 ! same as produce, but outputs the sum instead of the sequence of results
  56 : produce-sum ( id pred quot -- sum )
  57     [ 0 ] 2dip [ [ dip swap ] curry ] [ [ dip + ] curry ] bi* while ; inline
  58
  59 : x0 ( i -- x0 )
  60     number-length dup even?
  61     [ 2 - 2 / 10 swap ^ 7 * ]
  62     [ 1 - 2 / 10 swap ^ 2 * ] if ;
  63
  64 : ⌈a/b⌉  ( a b -- ⌈a/b⌉ )
  65     [ 1 - + ] keep /i ;
  66
  67 : xk+1 ( n xk -- xk+1 )
  68     [ ⌈a/b⌉ ] keep + 2 /i ;
  69
  70 : next-multiple ( a multiple -- next )
  71     [ [ 1 - ] dip /i 1 + ] keep * ;
  72
  73 DEFER: iteration#
  74 ! Gives the number of iterations when xk+1 has the same value for all a<=i<=n
  75 :: (iteration#) ( i xi a b -- # )
  76     a xi xk+1 dup xi =
  77     [ drop i b a - 1 + * ]
  78     [ i 1 + swap a b iteration# ] if ;
  79
  80 ! Gives the number of iterations in the general case by breaking into intervals
  81 ! in which xk+1 is the same.
  82 :: iteration# ( i xi a b -- # )
  83     a
  84     a xi next-multiple
  85     [ dup b < ]
  86     [
  87         ! set up the values for the next iteration
  88         [ nip [ 1 + ] [ xi + ] bi ] 2keep
  89         ! set up the arguments for (iteration#)
  90         [ i xi ] 2dip (iteration#)
  91     ] produce-sum
  92     ! deal with the last numbers
  93     [ drop b [ i xi ] 2dip (iteration#) ] dip
  94     + ;
  95
  96 : (euler255) ( a b -- answer )
  97     [ 10^ ] bi@ 1 -
  98     [ [ drop x0 1 swap ] 2keep iteration# ] 2keep
  99     swap - 1 + /f ;
 100
 101 PRIVATE>
 102
 103 : euler255 ( -- answer )
 104     13 14 (euler255) 10 nth-place ;
 105
 106 ! [ euler255 ] gc time
 107 ! Running time: 37.468911341 seconds
 108
 109 SOLUTION: euler255