x минус x не 0 в R

Question

это первый раз, когда я задаю вопрос здесь, поэтому надеюсь, что вы поймете мою проблему.

Вещь - это то, что я хочу сделать свой собственный fft(), не используя данный в R. Пока это хорошо работает для серии, такой как seq (1,5).

Но для c (1,1) произошло что-то странное. Насколько я мог указать на это, кажется, что x - x не является 0 в этом случае. Вот строки кода:

series <- c(1,1)     # defining the Serie 
     nr_samples <- length(series)  # getting the length 

################# 
# Calculating the harmonic frequncy 
################# 

     harmonic <- seq(0,(nr_samples-1)) 
     harmonic <- 2*pi*harmonic 
     harmonic <- harmonic/nr_samples 

################# 
# Exponential funktion needed for summing up 
################# 

     exponential <- function(index, omega){ 

      result <- exp(-((0+1i)*omega*index)) 

      return(result) 

     } 

################# 
# The sum for calculating the fft 
################# 

     my_fft <- function(time_series, omega){ 

      nr_samples <- length(time_series) 
      summand <- 0 

    # In the next loop the mistakes Happens  
    # While running this loop for harmonic[2] 
    # the rseult should be 0 because 
    # summand = - exp_factor 
    # The result for summand + exp_factor 
    # is 0-1.22464679914735e-16i 

    for (i in 1:nr_samples){ 

      exp_factor <- exponential((i-1), omega)    
      summand <- summand + time_series[i]*exp_factor  
      print(paste("Summand", summand, "Exp", exp_factor)) 
      }              

      return(summand)          
     } 


    transform <- sapply(harmonic, function(x){my_fft(series,x)}) 
    fft_transform <- fft(series) 
    df <- data.frame(transform, fft_transform) 
    print(df) 

Может кто-нибудь сказать мне, почему слагаемым + exp_factor для гармоники [2] не равно нулю ??

Answer 1

Это обычно называют FAQ 7.31, который говорит:

Единственные цифры, которые могут быть представлены именно в числовом типа R являются целые числа и дроби, знаменатель которой является степенью 2. Другие номера должны быть закруглены (обычно) 53 двоичных разряда. В результате два числа с плавающей запятой не будут надежно равны, если они не были вычислены по тому же алгоритму, и не всегда даже тогда. Например

R> a <- sqrt(2) 
R> a * a == 2 
[1] FALSE 
R> a * a - 2 
[1] 4.440892e-16 

Функция all.equal() сравнивает два объекта с использованием числового допуска .Machine $ double.eps^0,5. Если вы хотите гораздо большей точности, чем это, вам нужно будет внимательно рассмотреть распространение ошибок.

Для получения дополнительной информации см., Например, Дэвид Голдберг (1991), «Что каждый компьютерный ученый должен знать о арифметике с плавающей точкой», ACM Computing Surveys, 23/1, 5-48, также доступны через http://www.validlab.com/goldberg/paper.pdf.

Цитирую из «Элементов Programming Style» по Керниган и Plauger:

10,0 раз 0.1 вряд ли когда-либо 1,0.

(Конец цитаты)

Бумага Голдберг легендарна, и вы можете прочитать его. Это свойство все вычисления с плавающей запятой и не относятся к R.