Por el día del matemático, conociendo algo de las entrañas de R

                                                   
                                      88          
                                ,d    88          
                                88    88          
88,dPYba,,adPYba,  ,adPPYYba, MM88MMM 88,dPPYba,  
88P'   "88"    "8a ""     `Y8   88    88P'    "8a 
88      88      88 ,adPPPPP88   88    88       88 
88      88      88 88,    ,88   88,   88       88 
88      88      88 `"8bbdP"Y8   "Y888 88       88   
                                                   

 

Por el día del matemático, conociendo algo de las entrañas de R

    Hoy, por el día del matemático, quise por fin resolver una duda que me surgió hace algunas semanas: ¿se puede llegar a la unidad, dividiendo a la mitad iterativamente? Es decir: 0.5 + 0.25 + 0.126 + 0.0625 + …
    El razonamiento nos dice que no es posible llegar, y la contraparte de la pregunta, llegar a 0 dividiendo por mitad desde 1, parece más intuitiva, no hay un número que dividido entre dos sea 0. Al intentar hacerlo en R, el código es sencillo para programar la división por mitades, sin embargo, existe una característica, de las entrañas de R que nos impide llegar más profundo: los doubles (números decimales) están codificados con 53 bits por lo cual no llegamos más allá de la iteración 54:

temp = 0
for(i in 2^seq(1, 60, by = 1))  { 
  temp = temp + 1/i
  print(sprintf("%.40f", temp))
  flush.console()
  }
# nos da (el print es el mismo loop pero usando foreach:  
  result.50 "1.999999999999999111821580299874767661094665527343750000000000"
  result.51 "1.999999999999999555910790149937383830547332763671875000000000"
  result.52 "1.999999999999999777955395074968691915273666381835937500000000"
  result.53 "2.000000000000000000000000000000000000000000000000000000000000"
  result.54 "2.000000000000000000000000000000000000000000000000000000000000"
  result.55 "2.000000000000000000000000000000000000000000000000000000000000"

Entonces, tenemos que recurrir a la GNU Multiple precision floating-point reliably, que permite ir más al abismo infinitesimal, hagámoslo hasta 1/(2^1000):

library(Rmpfr)
temp = mpfr(0, 1000)
for(i in 2^seq(1, 1000, by = 1))  { 
  temp = temp + mpfr(1, 120)/i
  print(temp)
  flush.console()
}

# la iteración 1000:
1 'mpfr' number of precision  1000   bits 
[1] 0.99999999999999999978315956550289911319850943982601165771484374999999

99999999999999999999999999999999999999999999999999999999999999999999999999

99999999999999999999999999999999999999999999999999999999999999999999999999

99999999999999999999999999999999999999999999999999999999999999999999999999

999999999991

datamarindo.blogspot.mx, CCbySA, 2019-05-31