extra/project-euler/255/255.factor

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