extra/project-euler/255/255.factor

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