extra/project-euler/255/255.factor

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