Границы типов чисел в Common LISP и стеке перетекают в GHCI

Давайте определим промежуточную функцию Common Lisp:

(defun area (a b c)
  (let ((s (/ (+ a b c) 2)))
    (sqrt (* s (- s a) (- s b) (- s c)))))

Ваши тесты дают:

CL-USER> (area 333 333 333)
48016.344

CL-USER> (area 3333 3333 3333)
4810290.0

Во втором случае должно быть понятно, что и потолок, и пол равны. Это не относится к Haskell, где второй тест с 3333 возвращает:

4810290.040910754

Плавающая точка

В Common Lisp значение, из которого мы извлекаем квадратный корень:

370222244442963/16 

Это потому, что вычисления производятся с рациональными числами. До этого момента точность максимальна. Однако SQRT может возвращать либо рациональный, если возможно, либо приблизительный результат. Как особый случай, результат может быть целым числом в некоторых реализациях, как указал Райнер Джосвиг в комментарии. Это имеет смысл, потому что и integer, и ratio являются непересекающимися подтипами рационального типа. Но, как показывает ваша проблема, некоторые квадратные корни иррациональны (например, √2), и в этом случае CL ​​может возвращать число с плавающей запятой, приближающееся к значению (или сложное число с плавающей запятой).

Соответствующий раздел, касающийся чисел с плавающей запятой и математических функций, приведен в разделе 12.1.3.3 Правило заменяемости чисел с плавающей запятой. Короче говоря, результат преобразуется в single-float, когда вы вычисляете квадратный корень, что приводит к некоторой потере точности. Чтобы иметь двойника, вы должны быть более явными:

(defun area (a b c)
   (let ((s (/ (+ a b c) 2)))
     (sqrt (float (* s (- s a) (- s b) (- s c)) 0d0))))

Я также мог бы использовать (coerce ... 'double-float), но здесь я решил вызвать функцию преобразования FLOAT. Необязательный второй аргумент представляет собой прототип с плавающей запятой, т.е. значение целевого типа. Выше это 0d0, двойной поплавок. Вы также можете использовать 0l0 для длинных двойников или 0s0 для коротких. Этот параметр полезен, если вы хотите иметь ту же точность, что и входное число с плавающей запятой, но его можно использовать и с литералами, как в примере. Точное значение типов short, single, double или long float определяется реализацией, но они должны учитывать некоторые правила. Текущие реализации обычно дают большую точность, чем требуемый минимум.

CL-USER> (area 3333 3333 3333)
4810290.040910754d0

Теперь, если бы я хотел проверить, является ли результат целочисленным, я бы обрезал число с плавающей запятой и посмотрел, равно ли второе возвращаемое значение, остаток, нулю.

CL-USER> (zerop (nth-value 1 (truncate 4810290.040910754d0)))
NIL

Произвольной точности

Обратите внимание, что независимо от языка реализации (Haskell, CL или другой) подход будет давать неверные результаты для некоторых входных данных, учитывая то, как представлены числа с плавающей запятой. Действительно, та же проблема, что и с CL, может возникнуть для некоторых входных данных с более точными числами с плавающей запятой, где результат будет очень близок к целому числу. Вам может понадобиться другой математический подход или что-то вроде MPFR для вычислений с плавающей запятой произвольной точности. SBCL поставляется с sb-mpfr:

CL-USER> (require :sb-mpfr)
("SB-MPFR" "SB-GMP")

CL-USER> (in-package :sb-mpfr)
#<PACKAGE "SB-MPFR">

А потом:

SB-MPFR> (with-precision 256
           (sqrt (coerce 370222244442963/16 'mpfr-float)))
.4810290040910754427104204965311207243133723228299086361205561385039201180068712e+7
-1

Я отвечу на ваш второй вопрос, я не уверен насчет первого. В Haskell, поскольку это ленивый язык, когда вы используете хвостовую рекурсию с параметром-аккумулятором, может иметь место «накопление переходников». Преобразователь — это выражение, которое приостановлено и еще не вычислено. Возьмем более простой пример, суммируя все числа от 0 до n:

tri :: Int -> Int -> Int
tri 0 accum = accum
tri n accum = tri (n-1) (accum + n)

Если мы проследим оценку, мы увидим, что происходит:

tri 3 0
  = tri (3-1) (0+3)
  = tri 2 (0+3)
  = tri (2-1) ((0+3)+2)
  = tri 1 ((0+3)+2)
  = tri (1-1) (((0+3)+2)+1)
  = tri 0 (((0+3)+2)+1)
  = ((0+3)+2)+1     -- here is where ghc uses the C stack
  = (0+3)+2        (+1) on stack
  = 0+3            (+2) (+1) on stack
  = 0              (+3) (+2) (+1) on stack
  = 3              (+2) (+1) on stack
  = 5              (+1) on stack
  = 6             

Это, конечно, упрощение, но это интуиция, которая может помочь вам понять как переполнение стека, так и утечку пространства, вызванную накоплением транков. GHC оценивает преобразователь только тогда, когда это необходимо. Мы спрашиваем, равно ли значение n 0 каждый раз до tri, поэтому в этом параметре нет наращивания переходов, но никому не нужно знать значение accum до самого конца, что может быть действительно огромным переходом, как вы можете видеть из пример. При оценке этого огромного блока стек может переполниться.

Решение состоит в том, чтобы заставить tri оценить accum раньше. Обычно это делается с помощью BangPattern (но можно сделать и с seq, если вам не нравятся расширения).

{-# LANGUAGE BangPatterns #-}
tri :: Int -> Int -> Int
tri 0 !accum = accum
tri n !accum = tri (n-1) (accum + n)

! перед accum означает "оценить этот параметр в момент сопоставления с шаблоном" (хотя шаблону технически не нужно знать его значение). Затем мы получаем эту оценочную трассировку:

tri 3 0
  = tri (3-1) (0+3)
  = tri 2 3          -- we evaluate 0+3 because of the bang pattern
  = tri (2-1) (3+2)
  = tri 1 5
  = tri (1-1) (5+1)
  = tri 0 6
  = 6

Надеюсь, это поможет.

О стиле:

(if (predicate? ...) t nil)

просто

(predicate? ...)

Вы проверяете со своим IF, если T равно T, а затем возвращаете T. Но T это уже T, так что можно просто вернуть.