Formal verification of numerical subroutines

» Back to forum topics
  8 posts   Feedicon  
Replies: 7 - Last Post: March 13, 2010 05:16
by: Dmitry Nadezhin
« Older Topic » Newer Topic
showing 1 - 8 of 8
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 01, 2010 10:27 by Dmitry Nadezhin
» Permalink
На этом топике я буду рассказывать о моих попытках формальной верификации некоторых подпрограмм из различных вычислительных библиотек. Я надеюсь, что дело дойдет и до самой библиотеки JInterval, но начать я думаю не с ней.

Под формальной верификацией я понимаю следующее.
Пусть имеется спецификация того, что должна выдавать некоторая подпрограмма.
Рассматривается реализация этой подпрограммы, выполненная на неком языке - на псевдокоде, на языке программирования или на машинном языке. Язык должен иметь строгое определение семантики, которая может восприниматься как аксиомы некоторой машинной системой проверки доказательств. Формулируется теорема "после выполнения подпрограммы по семантическим правилам языка выдается результат, соответствующий спецификации", эта теорема доказывается, доказательство проверяется машинной системой проверки доказательств.

Вот пример спецификации подпрограммы вычисления натурального логарифма. Спецификация основывается на математическом понятии "натуральный логарифм" и кроме этого содержит требования к качеству аппроксимации и резульаты в специальных случаях.
/**
* Returns the natural logarithm (base <i>e</i>) of a {@code double}
* value. Special cases:
* <ul><li>If the argument is NaN or less than zero, then the result
* is NaN.
* <li>If the argument is positive infinity, then the result is
* positive infinity.
* <li>If the argument is positive zero or negative zero, then the
* result is negative infinity.</ul>
*
* <p>The computed result must be within 1 ulp of the exact result.
* Results must be semi-monotonic.
*
* @param a a value
* @return the value ln&nbsp;{@code a}, the natural logarithm of
* {@code a}.
*/
public static double log(double a);

Вот реализация подпрограммы для маленького диапазона на псевдокоде - входном языке системы gappa.

@rnd = float< ieee_64, ne >;
f = rnd(a - 1);
R = rnd(rnd(f*f)*rnd(rnd(0.5)-rnd(rnd(0.33333333333333333)*f)));
result = rnd(f - R);

Вот теорема о точности.
{ a in [0x0.fffffp0,0x1.00000ffffffffp0] -> | result -/ log(a) | <= 1b-53 }

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


------------------------
On this topic for I shall tell about my attempts to formal verification of some routines from various computational libraries. I hope that I shall prove subroutines form JInterval project simetimes, but I shall start from other libraries.

I understand the formal verification in such a way.
Suppose that we have a specification of what some subroutine must return.
We consider the implementation of this subroutine, written in a some language - in pseudocode, in programming language, or in machine language. The language should have a rigorous definition of semantics, which can be perceived as an axiom of some computer proof checker. A theorem is formulated "The result of execution of the subrouting by semantic rules of the language conforms to the specification", this theorem is proved, and the proof is checked by proof-checker.

Here is an example specification of subrouting to compute the natural logarithm The specification is based on the mathematical concept of "natural logarithm"; and in addition contains requirements for the quality of approximation and rezulaty in special cases.

/**
* Returns the natural logarithm (base <i>e</i>) of a {@code double}
* value. Special cases:
* <ul><li>If the argument is NaN or less than zero, then the result
* is NaN.
* <li>If the argument is positive infinity, then the result is
* positive infinity.
* <li>If the argument is positive zero or negative zero, then the
* result is negative infinity.</ul>
*
* <p>The computed result must be within 1 ulp of the exact result.
* Results must be semi-monotonic.
*
* @param a a value
* @return the value ln&nbsp;{@code a}, the natural logarithm of
* {@code a}.
*/
public static double log(double a);

Below is implementation of a subprogram for a small band of pseudo - input language system [url = http://gappa.gforge.inria.fr/] gappa [/ url].

@rnd = float< ieee_64, ne >;
f = rnd(a - 1);
R = rnd(rnd(f*f)*rnd(rnd(0.5)-rnd(rnd(0.33333333333333333)*f)));
result = rnd(f - R);

This is the theorem about the subroutine
{ a in [0x0.fffffp0,0x1.00000ffffffffp0] -> | result -/ log(a) | <= 1b-53 }

This topic will be a kind of blog. I will learn to prove and put the notes about the learing here. I will write in Russian, and I will optionally translate some posts into English. <br>
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 01, 2010 12:42 by Dmitry Nadezhin
» Permalink
Выбор проверяльщика доказательств

Я выбрал Coq, потому что он был разработан в INRIA, а там я нашел еще несколько тулов, которые подходят для локазательств.
Gappa
Why
Krakatoa

Coq - поверяльщик доказательств с интерактивной средой для их разработки. Язык у него довольно сложный.
Gappa - простая системя для оценки ошибок округления в вычислениях. Внутри она реализована на интервалах,
но может использоваться для проверки любых вычислений как интервальных, так и неинтервальных.

Я долго и медленно изучал Coq, а с Gappa начал работать сразу. Gappa выводит переводит свое доказательства на язык Coq для проверки.

Для сборки Gappa потребовались:
boost библиотека (без установки, только .hpp файлы)ю
gmp
mpfr

Со сборкой Coq под OpenSolaris возни было больше, а в Linux-ах он кажется уже есть в репозиториях в виде пакета.
Чтобы coq смог проверить доказательство, выведенное из Gappa потребуется еще coq-библиотека gappalib-coq .
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 03, 2010 04:15 by Dmitry Nadezhin
» Permalink
Пробую Gappa для оценки констант


В файле java.lang.Math определены double константы E и PI,
являющиеся аппроксимациями иррациональных математических констант e и π.
Комментарии содержат спецификации этих констант.
----
/**
* The {@code double} value that is closer than any other to
* <i>e</i>, the base of the natural logarithms.
*/
public static final double E = 2.7182818284590452354;

/**
* The {@code double} value that is closer than any other to
* <i>pi</i>, the ratio of the circumference of a circle to its
* diameter.
*/
public static final double PI = 3.14159265358979323846;
----

Посмотрим, как Gappa может помочь в проверке того , что константы удовлетворяют спецификации.

Используем разложение в бесконечный ряд
e = 1 + 1/1! + 1/2! + 1/3! + 1/4! + 1/5! + ... .
Если возьмем несколько n=2 первых членов (1 + 1/1!),
то остаточный ряд ограничивается
0 < 1/2! + 1/3! + 1/4! + 1/5! + ... = 1/2! * (1 +1/3+1/(3*4) + 1/(3*4*5) + ... ) <
1/2! * (1 + (1/3) + (1/3)^2 + (1/3)^3 + ...) = 1/2! * 1/(1 - 1/3) = 1/2!*3/2 = 1/2!*(1 + 1/2)

То есть остаток ряда есть 1/n!*(1+r/n), где r - некое число между 0 и 1.
Запишем это на Gappa

---- e0.g
e = 1 + 1/1*(1 + 1/2*(1 + r/2));

{ r in [0,1] -> e in ? }
-----
Сначала идут определения. Действительное e определяется через действительное r .
Потом идет предикат, который хотим доказать. В нашем случае он импликация вида
"Если r в интервале [0.1], то e в интервале [?,?]".

Запускаем Gappa, она отвечает:
----
Results for r in [0, 1]:
e in [5b-1 {2.5, 2^(1.32193)}, 11b-2 {2.75, 2^(1.45943)}]
-----
То есть e в интервале
[ 5*2^(-1) , 11*2^(-2) ] = [5/2 , 11/4] = [2.5 , 2.75] ~ [ 2^(1.32193) , 2^(1.45943)]

Увеличим количество членов ряда
---- e1.g
e = 1 +
1/1*(1 +
1/2*(1 +
1/3*(1 +
r/3)));

{ r in [0,1] -> e in ? }
----
Results for r in [0, 1]:
e in [384307168202282325b-57 {2.66667, 2^(1.41504)}, 196156783769914937b-56 {2.72222, 2^(1.44478)}]
----


Запишем теперь на Gappa double константу Math.E .
Введем сокращение "rnd" для операции округления к ближайшему double
Math.E - это результат округления числа, заданного десятичной константой, к double .
Попросим Gappa выдать оценки Math_E, e и их разности
---- e2.g
@rnd = float< ieee_64, ne >;

Math_E = rnd(2.7182818284590452354e0);

e = 1 +
1/1*(1 +
1/2*(1 +
1/3*(1 +
r/3)));

{ r in [0,1] -> Math_E in ? /\ e in ? /\ Math_E - e in ? }
----
Results for r in [0, 1]:
Math_E in [6121026514868073b-51 {2.71828, 2^(1.44269)}, 6121026514868073b-51 {2.71828, 2^(1.44269)}]
Math_E - e in [-283935294136601b-56 {-0.00394039, -2^(-7.98744)}, 7438528749274347b-57 {0.0516152, 2^(-4.27606)}]
e in [384307168202282325b-57 {2.66667, 2^(1.41504)}, 196156783769914937b-56 {2.72222, 2^(1.44478)}]
----
Math_E - это точное двоичное число, для e и для разности получена оценка.


Возьмем тепер столько членов ряда, чтобы оценить
с точностью до Unit in the Last Place (ULP)
Для числа 2 <= e < 4 ULP равен 4*2^(-53)=2^(-51) .

---- e3.g
@rnd = float< ieee_64, ne >;

Math_E = rnd(2.7182818284590452354e0);

e = 1 +
1/1*(1 +
1/2*(1 +
1/3*(1 +
1/4*(1 +
1/5*(1 +
1/6*(1 +
1/7*(1 +
1/8*(1 +
1/9*(1 +
1/10*(1 +
1/11*(1 +
1/12*(1 +
1/13*(1 +
1/14*(1 +
1/15*(1 +
1/16*(1 +
1/17*(1 +
1/18*(1 +
r/18))))))))))))))))));

{ r in [0,1] -> Math_E in ? /\ e in ? /\ Math_E - e in ? /\ Math_E - rnd(e) in ? /\ (Math_E - e) * 1b51 in ? }
----
Results for r in [0, 1]:
(Math_E - e) * 1b51 in [-21b-6 {-0.328125, -2^(-1.60768)}, -39b-7 {-0.304688, -2^(-1.7146)}]
Math_E in [6121026514868073b-51 {2.71828, 2^(1.44269)}, 6121026514868073b-51 {2.71828, 2^(1.44269)}]
Math_E - e in [-21b-57 {-1.45717e-16, -2^(-52.6077)}, -39b-58 {-1.35308e-16, -2^(-52.7146)}]
Math_E - float<53,-1074,ne>(e) in [0, 0]
e in [783491393903113383b-58 {2.71828, 2^(1.44269)}, 391745696951556693b-57 {2.71828, 2^(1.44269)}]
----
Math_E - e находится в диапазоне [-033,-0.30] ULP . Это значит, что Math_E есть ближайшее double к e .
По другому это записывается Math_e - rnd(e) in [0,0] .
Кроме того, это значит, что точное e закличено между Math.E и Math.nextUp(Math.E) .

Аналогично для числа π .
Используем простое разложение
π/4 = 4*arctg(1/5) - arctg(1/239)

---- pi.g
@rnd = float< ieee_64, ne >;

Math_PI = rnd(3.14159265358979323846e0);

sqr_1_5 = 1/(5*5);
atan_1_5 = 1/5*(
1/1 - sqr_1_5*(
1/3 - sqr_1_5*(
1/5 - sqr_1_5*(
1/7 - sqr_1_5*(
1/9 - sqr_1_5*(
1/11 - sqr_1_5*(
1/13 - sqr_1_5*(
1/15 - sqr_1_5*(
1/17 - sqr_1_5*(
1/19 - sqr_1_5*(
1/21 - sqr_1_5*(
1/23 - sqr_1_5*(
1/25 * r1)))))))))))));

sqr_1_239 = 1/(239*239);
atan_1_239 = 1/239*(
1/1 - sqr_1_239*(
1/3 - sqr_1_239*(
1/5 - sqr_1_239*(
1/7 - sqr_1_239*(
1/9 * r2)))));

pi = 4*(4*atan_1_5 - atan_1_239);

{ r1 in [0,1] /\ r2 in [0,1] -> Math_PI in ? /\ pi in ? /\ Math_PI - pi in ? /\ Math_PI - rnd(pi) in ? /\ (Math_PI - pi) * 1b51 in ? }
-----
Results for r1 in [0, 1] and r2 in [0, 1]:
(Math_PI - pi) * 1b51 in [-37b-7 {-0.289062, -2^(-1.79055)}, -1b-2 {-0.25, -2^(-2)}]
Math_PI in [884279719003555b-48 {3.14159, 2^(1.6515)}, 884279719003555b-48 {3.14159, 2^(1.6515)}]
Math_PI - float<53,-1074,ne>(pi) in [0, 0]
Math_PI - pi in [-37b-58 {-1.2837e-16, -2^(-52.7905)}, -1b-53 {-1.11022e-16, -2^(-53)}]
pi in [28296951008113761b-53 {3.14159, 2^(1.6515)}, 905502432259640357b-58 {3.14159, 2^(1.6515)}]
b
---

Math.PI - π находится в диапазоне [-0.29,-0.25] ULP , то есть π заключено между Math.PI и Math.nextUp(Math.PI)

Gappa файлы можно найти на SVN
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 03, 2010 09:22 by Dmitry Nadezhin
» Permalink
com.kenai.jinterval.stdfun.RealStdFun


Интервальные оценки e и π, полученные в предудущем посте,
позволют узнать, что аппроксимации в классе com.kenai.jinterval.stdfun.java_math.StdFun
более точны, чем было ранее известно. Напомню, что это за класс.

Этот класс реализает спецификацию, прописанную в базовом абстрактном классе
com.kenai.jinterval.stdfun.RealStdFun .

Класс RealStdFun описывает, что должна делать вычислялка элеметарных функций
(пока только двух констант). Она будет иметь конкретные подклассы, вычисляющие
элементарные функции разными способами.

Рассмотрим методы, вычисляющие e. Они выдают приближенный результат ea.
Абсолютной ошибкой результата называется (ea - e).

Методы предъявляют разные требования к знаку и величине ошибки.

Binary eLower(int maxErrExp) -2^maxErrExp <= ea - e <= 0
Binary eUpper(int maxErrExp) 0 <= ea - e <= +2^maxErrExp
Binary eNear (int maxErrExp) -2^maxErrExp <= ea - e <= +2^maxErrExp

Класс com.kenai.jinterval.stdfun.RealStdFun - пример подкласса,
частично реализующего это спецификацию.
Он выдает результат double-аппроксимацию числа e, если она допустима по заданной
точности и вываливается в противном случае.

Используя нтервльные оценки из предудущего поста можно немного увеличить диапазон точности, в которых эти приближения еще верны.
Про нижнюю границу eLower теперь знаем что ее точность не хуже 2^(-52) ( было 2^(-50) ).
Про верхнюю границу eUpper теперь знаем что ее точность не хуже 2^(-51) ( было 2^(-50) ).

Аналогично по границам piLower и piUpper.
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 06, 2010 09:04 by Dmitry Nadezhin
» Permalink
Разглядываю Coq-файл, сгенерированный Gappa'ой


Gappa находит доказательство некоторых фактов о вычислении. Однако она не настаивает, чтобы мы верили ее доказательству. Gappa может вывести сформулированную теорему вместе с ее доказательством в одну из систем проверки доказательств: Coq или HOL Light.

Посмотрю, какой файл для Coq выдаст Gappa для простого входного файла.
---- e0.g
e = 1 + 1/1*(1 + 1/2*(1 + r/2));

{ r in [0,1] -> e in ? }
----

Командная строка:
$ gappa e0.g -Bcoq >e0.v

выдает довольно длинный файл:
---- e0.v
Require Import Gappa_library.
Section Generated_by_Gappa.
Definition f1 := Float2 (1) (0).
Definition f2 := Float2 (0) (0).
Definition i1 := makepairF f2 f1.
Variable _r : R.
Notation p1 := (BND _r i1). (* BND(r, [0, 1]) *)
Definition f3 := Float2 (11) (-2).
Definition f4 := Float2 (5) (-1).
Definition i2 := makepairF f4 f3.
Notation r8 := (Float1 (2)).
Notation r7 := ((_r / r8)%R).
Notation r9 := (Float1 (1)).
Notation r6 := ((r9 + r7)%R).
Notation r10 := ((r9 / r8)%R).
Notation r5 := ((r10 * r6)%R).
Notation r4 := ((r9 + r5)%R).
Notation r11 := ((r9 / r9)%R).
Notation r3 := ((r11 * r4)%R).
Notation _e := ((r9 + r3)%R).
Notation p2 := (BND _e i2). (* BND(e, [2.5, 2.75]) *)
Definition i3 := makepairF f1 f1.
Notation p3 := (BND r9 i3). (* BND(1, [1, 1]) *)
Lemma t1 : p3.
apply constant1 ; finalize.
Qed.
Lemma l3 : p3. (* BND(1, [1, 1]) *)
apply t1.
Qed.
Definition f5 := Float2 (7) (-2).
Definition f6 := Float2 (3) (-1).
Definition i4 := makepairF f6 f5.
Notation p4 := (BND r3 i4). (* BND(1 / 1 * (1 + 1 / 2 * (1 + r / 2)), [1.5, 1.75]) *)
Notation p5 := (BND r11 i3). (* BND(1 / 1, [1, 1]) *)
Notation p6 := (NZR r9). (* NZR(1) *)
Notation p7 := (ABS r9 i3). (* ABS(1, [1, 1]) *)
Lemma t2 : p3 -> p7.
intros h0.
apply abs_of_bnd_p with (1 := h0) ; finalize.
Qed.
Lemma l7 : p7. (* ABS(1, [1, 1]) *)
assert (h0 : p3). apply l3.
apply t2. exact h0.
Qed.
Lemma t3 : p7 -> p6.
intros h0.
apply nzr_of_abs with (1 := h0) ; finalize.
Qed.
Lemma l6 : p6. (* NZR(1) *)
assert (h0 : p7). apply l7.
apply t3. exact h0.
Qed.
Lemma t4 : p6 -> p5.
intros h0.
apply div_refl with (1 := h0) ; finalize.
Qed.
Lemma l5 : p5. (* BND(1 / 1, [1, 1]) *)
assert (h0 : p6). apply l6.
apply t4. exact h0.
Qed.
Notation p8 := (BND r4 i4). (* BND(1 + 1 / 2 * (1 + r / 2), [1.5, 1.75]) *)
Definition f7 := Float2 (3) (-2).
Definition f8 := Float2 (1) (-1).
Definition i5 := makepairF f8 f7.
Notation p9 := (BND r5 i5). (* BND(1 / 2 * (1 + r / 2), [0.5, 0.75]) *)
Definition i6 := makepairF f8 f8.
Notation p10 := (BND r10 i6). (* BND(1 / 2, [0.5, 0.5]) *)
Definition f9 := Float2 (1) (1).
Definition i7 := makepairF f9 f9.
Notation p11 := (BND r8 i7). (* BND(2, [2, 2]) *)
Lemma t5 : p11.
apply constant1 ; finalize.
Qed.
Lemma l11 : p11. (* BND(2, [2, 2]) *)
apply t5.
Qed.
Lemma t6 : p3 -> p11 -> p10.
intros h0 h1.
apply div_pp with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l10 : p10. (* BND(1 / 2, [0.5, 0.5]) *)
assert (h0 : p3). apply l3.
assert (h1 : p11). apply l11.
apply t6. exact h0. exact h1.
Qed.
Definition i8 := makepairF f1 f6.
Notation p12 := (BND r6 i8). (* BND(1 + r / 2, [1, 1.5]) *)
Definition i9 := makepairF f2 f8.
Notation p13 := (BND r7 i9). (* BND(r / 2, [0, 0.5]) *)
Lemma t7 : p1 -> p11 -> p13.
intros h0 h1.
apply div_pp with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l13 : p1 -> p13. (* BND(r / 2, [0, 0.5]) *)
intros h0.
assert (h1 : p11). apply l11.
apply t7. exact h0. exact h1.
Qed.
Lemma t8 : p3 -> p13 -> p12.
intros h0 h1.
apply add with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l12 : p1 -> p12. (* BND(1 + r / 2, [1, 1.5]) *)
intros h0.
assert (h1 : p3). apply l3.
assert (h2 : p13). apply l13. exact h0.
apply t8. exact h1. exact h2.
Qed.
Lemma t9 : p10 -> p12 -> p9.
intros h0 h1.
apply mul_pp with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l9 : p1 -> p9. (* BND(1 / 2 * (1 + r / 2), [0.5, 0.75]) *)
intros h0.
assert (h1 : p10). apply l10.
assert (h2 : p12). apply l12. exact h0.
apply t9. exact h1. exact h2.
Qed.
Lemma t10 : p3 -> p9 -> p8.
intros h0 h1.
apply add with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l8 : p1 -> p8. (* BND(1 + 1 / 2 * (1 + r / 2), [1.5, 1.75]) *)
intros h0.
assert (h1 : p3). apply l3.
assert (h2 : p9). apply l9. exact h0.
apply t10. exact h1. exact h2.
Qed.
Lemma t11 : p5 -> p8 -> p4.
intros h0 h1.
apply mul_pp with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l4 : p1 -> p4. (* BND(1 / 1 * (1 + 1 / 2 * (1 + r / 2)), [1.5, 1.75]) *)
intros h0.
assert (h1 : p5). apply l5.
assert (h2 : p8). apply l8. exact h0.
apply t11. exact h1. exact h2.
Qed.
Lemma t12 : p3 -> p4 -> p2.
intros h0 h1.
apply add with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l2 : p1 -> p2. (* BND(e, [2.5, 2.75]) *)
intros h0.
assert (h1 : p3). apply l3.
assert (h2 : p4). apply l4. exact h0.
apply t12. exact h1. exact h2.
Qed.
Lemma l1 : p1 -> p2. (* BND(e, [2.5, 2.75]) *)
intros h0.
apply l2. exact h0.
Qed.
End Generated_by_Gappa.
----

Первая строка этого файла импортирует библиотеку поддержки Gappa на языке Coq .
Require Import Gappa_library.
Грубо можно сказать, что это библиотека, определяющая интервальную арифметику для Coq.

---- Gappa_library.v
Require Export Rdefinitions.
Require Export Rbasic_fun.
Require Export R_sqrt.
Require Export Bool.
Require Export Gappa_definitions.
Require Export Gappa_pred_bnd.
Require Export Gappa_pred_abs.
Require Export Gappa_pred_fixflt.
Require Export Gappa_pred_nzr.
Require Export Gappa_pred_rel.
Require Export Gappa_rewriting.
Require Export Gappa_round_def.
Require Export Gappa_fixed.
Require Export Gappa_float.
Require Export Gappa_user.
Require Export Gappa_proxy.

Ltac finalize := vm_cast_no_check (refl_equal true).
Ltac next_interval t h :=
apply t with (1 := h) ; [ finalize |
clear h ; intro h ; simpl in h ; generalize h ; clear h ].
----

Ипортируя Gappa_library, мы импортируем тем самым еще ряд модулей.
Часть из них из стандартной библиотеки Coq: Rdefinitions, Rbasic_fun, R_sqrt, Bool. Они содержат основные теоремы о Real числах.

Следующий файл - Gappa_definitions.
---- Gappa_definitions.v
Require Import Reals.
Require Import ZArith.
Require Import Gappa_integer.

Section Gappa_definitions.

Record float2 : Set := Float2 {Fnum : Z; Fexp : Z}.
Coercion float2R (x : float2) := F2R 2 (Fnum x) (Fexp x).
Record float10 : Set := Float10 { Fnum10 : Z ; Fexp10 : Z }.
Coercion float10R (x : float10) := F2R 10 (Fnum10 x) (Fexp10 x).

Record FF: Set := makepairF { lower : float2 ; upper : float2 }.

Definition BND (x : R) (xi : FF) :=
(lower xi <= x <= upper xi)%R.
Definition ABS (x : R) (xi : FF) :=
(0 <= lower xi)%R /\ (lower xi <= Rabs x <= upper xi)%R.
Definition REL (x1 x2 : R) (xi : FF) :=
exists x : R, (lower xi <= x <= upper xi)%R /\ (x1 = x2 * (1 + x))%R.
Definition FIX (x : R) (n : Z) :=
exists f : float2, float2R f = x /\ (n <= Fexp f)%Z.
Definition FLT (x : R) (n : positive) :=
exists f : float2, float2R f = x /\ (Zabs (Fnum f) < Zpower_pos 2 n)%Z.
Definition NZR (x : R) := (x <> 0)%R.

Definition contradiction := forall P : Prop, P.

End Gappa_definitions.
----
Определяются типы float2 и float10 - плавающие числа с снованиями 2 и 10.
Каждый из типов состоеит из двух целых полей с мантиссой и с экспонентой.
Функции-конструкторы этих типов - Float2 и Float10 .
(Float2 5 -2) создает плавающее число, значение которого 5*2^(-2)=1.25 .
Функции-селекторы полей: Fnum, Fexp и Fnum10, Fexp10 .

Определяется тип FF, который собственно и представляет интервал.
Функция-конструктор - makepair2.
Функции-селекторы - lower и upper .
(makepair2 (Float2 1 0) (Float2 1 1)) создает интервал [1.0 , 2.0]

Определены также 5 предикатов
(BND x xi) - действительное "x" содержится в интервале с плавающими границами "xi"
(ABS x xi) - абсолютная величина действительного "x" содержится в неотрицательном интервале "xi"
(REL approx exact xi) - относительная погрешность, с которой действительное "approx" приближает
действительное "exact", содержится в "xi".
(FIX x n) - действительное "x" есть кратное числа 2^n
(FLT x n) - действительное "x" представимо в виде действительного числа с длиной мантиссы "n" и неограниченной экспонентой
(NZR x) - действительное "x" отлично от нуля

В e0.v выражение для e разбито на подвыражения. Все они имеют действительный тип.
----
Variable _r : R.
Notation r8 := (Float1 (2)). (* 2 *)
Notation r7 := ((_r / r8)%R). (* r/2 *)
Notation r9 := (Float1 (1)). (* 1 *)
Notation r6 := ((r9 + r7)%R). (* 1 + r/2 *)
Notation r10 := ((r9 / r8)%R). (* 1/2 *)
Notation r5 := ((r10 * r6)%R). (* 1/2*(1 + r/2) *)
Notation r4 := ((r9 + r5)%R). (* 1 + 1/2*(1 + r/2) *)
Notation r11 := ((r9 / r9)%R). (* 1/1 *)
Notation r3 := ((r11 * r4)%R). (* 1/1*(1 + 1/2*(1 + r/2)) *)
Notation _e := ((r9 + r3)%R). (* 1 + 1/1*(1 + 1/2*(1 + r/2)) *)
----
В e0.v определены также плаващие и интервальные константы
----
Definition f1 := Float2 (1) (0). (* 1.0 *)
Definition f2 := Float2 (0) (0). (* 0.0 *)
Definition i1 := makepairF f2 f1. (* [0.0, 1.0] *)
Definition f3 := Float2 (11) (-2). (* 2.75 *)
Definition f4 := Float2 (5) (-1). (* 2.5 *)
Definition i2 := makepairF f4 f3. (* [2.5,2.75] *)
Definition i3 := makepairF f1 f1. (* [1.0, 1.0] *)
Definition f5 := Float2 (7) (-2). (* 1.75 *)
Definition f6 := Float2 (3) (-1). (* 1.5 *)
Definition i4 := makepairF f6 f5. (* [1.5, 1.75] *)
Definition f7 := Float2 (3) (-2). (* 0.75 *)
Definition f8 := Float2 (1) (-1). (* 0.5 *)
Definition i5 := makepairF f8 f7. (* [0.5, 0.75] *)
Definition i6 := makepairF f8 f8. (* [0.5, 0.5] *)
Definition f9 := Float2 (1) (1). (* 2.0 *)
Definition i7 := makepairF f9 f9. (* [2.0, 2.0] *)
Definition i8 := makepairF f1 f6. (* [1.0, 1.5] *)
Definition i9 := makepairF f2 f8. (* [0.0, 0.5] *)
----

Из предыдущих определений выражены элементарные предикаты
----
Notation p1 := (BND _r i1). (* BND(r, [0, 1]) *)
Notation p11 := (BND r8 i7). (* BND(2, [2, 2]) *)
Notation p13 := (BND r7 i9). (* BND(r / 2, [0, 0.5]) *)
Notation p3 := (BND r9 i3). (* BND(1, [1, 1]) *)
Notation p7 := (ABS r9 i3). (* ABS(1, [1, 1]) *)
Notation p6 := (NZR r9). (* NZR(1) *)
Notation p12 := (BND r6 i8). (* BND(1 + r / 2, [1, 1.5]) *)
Notation p10 := (BND r10 i6). (* BND(1 / 2, [0.5, 0.5]) *)
Notation p9 := (BND r5 i5). (* BND(1 / 2 * (1 + r / 2), [0.5, 0.75]) *)
Notation p8 := (BND r4 i4). (* BND(1 + 1 / 2 * (1 + r / 2), [1.5, 1.75]) *)
Notation p5 := (BND r11 i3). (* BND(1 / 1, [1, 1]) *)
Notation p4 := (BND r3 i4). (* BND(1 / 1 * (1 + 1 / 2 * (1 + r / 2)), [1.5, 1.75]) *)
Notation p2 := (BND _e i2). (* BND(e, [2.5, 2.75]) *)
----

Используя все предыдущее, в e0.v сформулированы и доказаны ряд лемм,
последняя из которых имеет вид:

Lemma l1 : p1 -> p2. (* BND(e, [2.5, 2.75]) *)

Раскрыв все обозначения, что фактически это означает импликацию:

r in [0.0, 1.0] -> 1 + 1/1*(1 + 1/2*(1 + r/2)) in [2.5, 2.75]

Это и есть основной результат e0.v .

Разглядывание других лемм из e0.v и их доказательств отложу до следующего поста.
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 09, 2010 11:15 by Dmitry Nadezhin
» Permalink
Разглядываю доказательство e0.v в Coq


Посмотрю, что Gappa подготовила для Coq. Простая вещь - интервальное оценивание рационального числа 1/2 :.
------
Require Import Gappa_library.

Notation r9 := (Float1 (1)).
Notation r8 := (Float1 (2)).
Notation r10 := ((r9 / r8)%R).

Definition f1 := Float2 (1) (0).
Definition i3 := makepairF f1 f1.
Definition f9 := Float2 (1) (1).
Definition i7 := makepairF f9 f9.
Definition f8 := Float2 (1) (-1).
Definition i6 := makepairF f8 f8.

Notation p3 := (BND r9 i3). (* BND(1, [1, 1]) *)
Notation p11 := (BND r8 i7). (* BND(2, [2, 2]) *)
Notation p10 := (BND r10 i6). (* BND(1 / 2, [0.5, 0.5]) *)

Lemma t1 : p3.
apply constant1 ; finalize.
Qed.
Lemma l3 : p3. (* BND(1, [1, 1]) *)
apply t1.
Qed.
Lemma t5 : p11.
apply constant1 ; finalize.
Qed.
Lemma l11 : p11. (* BND(2, [2, 2]) *)
apply t5.
Qed.
Lemma t6 : p3 -> p11 -> p10.
intros h0 h1.
apply div_pp with (1 := h0) (2 := h1) ; finalize.
Qed.
Lemma l10 : p10. (* BND(1 / 2, [0.5, 0.5]) *)
assert (h0 : p3). apply l3.
assert (h1 : p11). apply l11.
apply t6. exact h0.exact h1.
Qed.
----------
Запускаю coqtop и подаю ему объявления перед Lemma t1. Он проглатывает.
Теперь подаю формулировку леммы:
--------------------------- Ввод
Lemma t1 : p3.
apply constant1.
finalize.
Qed.
Print t1.
------------------------ Вывод.
Coq < Lemma t1 : p3.
1 subgoal

============================
p3

t1 < apply constant1.
1 subgoal

============================
constant2_helper (Float2 1 0) i3 = true

t1 < finalize.
Proof completed.

t1 < Qed.

t1 is defined

Coq < Print t1.
t1 =
constant1 1 i3 (refl_equal true<:constant2_helper (Float2 1 0) i3 = true)
: p3

Coq <
--------------------------------------
Доказательство в Coq похоже на постооение дерева вывода в грамматике.
Нетерминалами являются логические формулы, которые надо доказать.
Есть правила вывода, которые говорят - чтобы доказать такую-то логическую
формулу, достаточно доказать некоторые другие логические формулы - потомки.
Построение дерева вывода проходит интерактивно от корня к листьям.

В нашем случае корень - формула p3 , которая утверждает, что целая константа "1" содержится в интервале [0b0,1b0] = [0.0,1.0].

На первом шаге с помощи тактики apply ставим в корень правило constant1.
Тактика - что-то вроде макрокоманды, которая помогает подобрать нужное правило.
Правилом может быть ранее доказанная лемма.
Правило constant1 имеет вид
------
Coq < Check constant1.
constant1
: forall (x : Z) (zi : FF),
constant2_helper (Float2 x 0) zi = true -> BND (Float1 x) zi
-------
Это правило доказывает формулу вида BND (Float1 x) zi . У правила два параметра-терма: целое число "x" и интервал "zi". У правила одна формула-потомок
constant2_helper (Float2 x 0) zi = true.
После того, как мы поместили в корень правило constant1, текущей недоказанной формулой становится формула
constant2_helper (Float2 1 0) i3 = true
Видим, что формальные параметры правила "x" и "zi" заменились на фактические термы "1" и "i3".
Следующий шаг вызывает тактику "finalize", которая просто проверяет, что "1" находится между границами интервала "i3", и ставит в формулу правило-заглушку без формул-потомков
refl_equal true<:constant2_helper (Float2 1 0) i3 = true.

В результате мы поставили в корень "p3" законченное дерево вывода
constant1 1 i3 (refl_equal true<:constant2_helper (Float2 1 0) i3 = true)
Дерево вывода записано в префиксной форме.
Сначала имя правило верхнего уровня "constant1", затем его параметры "1" и "i3", затем в скобках поддерево вывода.

Команда Qed перепроверяет дерево вывода на корректность и закрывает доказательство леммы.

Если мы уже знаем дерево вывода, то мы можем поставить его в корень в неизменном виде тактикой exact (дерево вывода). Эта возможность полезна, когда мы строим дерево вывода не интерактивно, а генерим Coq-скрипт в другой программе (типа Gappa) только для того, чтобы перепроверить его в Coq.
-----------------------
Coq < Lemma t1_exact : p3 .
1 subgoal

============================
p3

t1_exact < exact (constant1 1 i3 (refl_equal true<:constant2_helper (Float2 1 0) i3 = true)).
Proof completed.

t1_exact < Qed.
exact
(constant1 1 i3 (refl_equal true<:constant2_helper (Float2 1 0) i3 = true)).

t1_exact is defined

Coq <
-----------------------

Есть чуть более умная команда refine, которая так же как exact подставляет дерево вывода. В ней можно опустить ряд параметров если их можно восстановить из других параметров.
-------------------
Coq < Lemma t1_refine : p3.
1 subgoal

============================
p3

t1_refine < refine (constant1 _ _ (refl_equal true<:constant2_helper (Float2 1 0) i3 = true)).
Proof completed.

t1_refine < Qed.
refine
(constant1 _ _ (refl_equal true<:constant2_helper (Float2 1 0) i3 = true)).

t1_refine is defined

Coq <
-------------------
Coq < Lemma t1_refine_incomplete : p3.
1 subgoal

============================
p3

t1_refine_incomplete < refine (constant1 _ _ _).
1 subgoal

============================
constant2_helper (Float2 1 0) i3 = true

t1_refine_incomplete < finalize.
Proof completed.

t1_refine_incomplete < Qed.
refine (constant1 _ _ _).
finalize.

t1_refine_incomplete is defined

Coq <
-------------------------

И, наконец, в командах exact и refine можно указать неполное дерево вывода.
------------------
Coq < Lemma t1_refine_incomplete : p3.
1 subgoal

============================
p3

t1_refine_incomplete < refine (constant1 _ _ _).
1 subgoal

============================
constant2_helper (Float2 1 0) i3 = true

t1_refine_incomplete < finalize.
Proof completed.

t1_refine_incomplete < Qed.
refine (constant1 _ _ _).
finalize.

t1_refine_incomplete is defined

Coq <
-----------------------

Доказанная лемма становится новым правилом вывода
----------
Coq < Check t1.
t1
: p3
-----------
То, что у нее нет параметров и потомков, означает, что она может поставить в корень с форумлой p3 законченной дерево вывода. Лемма
Lemma l3 : p3. (* BND(1, [1, 1]) *)
apply t1.
Qed.
зачем-то так и делает
---------------------------
Coq < Lemma l3 : p3. (* BND(1, [1, 1]) *)
1 subgoal

============================
p3

l3 < apply t1.
Proof completed.

l3 < Qed.
apply t1.

l3 is defined

Coq < Print l3.
l3 = t1
: p3

Coq <
----------------------------
Дерево, построенное леммой l3 - это в точности то же дерево, что и в лемме t1.



Рассмотрем теперь построение дерева вывода для леммы t6, которая доказывает формулу
(BND 1 [1.0,1.0]) -> (BND 2 [2.0,2.0]) -> (BND (1/2) [0.5,0.5]) .Знак "->" - это знак логической имликации.
Смысл нашей формулу "если доказано, что (BND 1 ..) и (BND 2 ..), то верно что (BND (1/2 ...).
-------------------
Coq < Lemma t6 : p3 -> p11 -> p10.
1 subgoal

============================
p3 -> p11 -> p10

t6 < intros h0 h1.
1 subgoal

h0 : p3
h1 : p11
============================
p10

t6 < apply div_pp.
Toplevel input, characters 0-12:
> apply div_pp.
> ^^^^^^^^^^^^
Error: Unable to find an instance for the variables xi, yi.

t6 < apply div_pp with (1 := h0) (2 := h1).
1 subgoal

h0 : p3
h1 : p11
============================
div_pp_helper i3 i7 i6 = true

t6 < finalize.
Proof completed.

t6 < Qed.
intros h0 h1.
apply div_pp with (1 := h0) (2 := h1).
finalize.

t6 is defined

Coq < Print t6.
t6 =
fun (h0 : p3) (h1 : p11) =>
div_pp r9 r8 i3 i7 i6 h0 h1 (refl_equal true<:div_pp_helper i3 i7 i6 = true)
: p3 -> p11 -> p10

Coq < Check t6.
t6
: p3 -> p11 -> p10

Coq <
-------------------
Эта лемма строит неполный вывод формулы p10 по правилу вывода div_pp.
------------
Coq < Check div_pp.
div_pp
: forall (x y : R) (xi yi zi : FF),
BND x xi ->
BND y yi -> div_pp_helper xi yi zi = true -> BND (x / y) zi
-----------
У правила div_pp три формулы потомка. Третьий мы подсоединили в лемме (с правилом div_pp_helper),
а первые два еще не достроены.
Таким образом, эта лемма строит новое правило вывода t6 с двумя формулами-потомками.
Посмотрим, как это правило будет использоваться в доказательстве леммы
Lemma l10 : p10. (* BND(1 / 2, [0.5, 0.5]) *)
-----------------------------------
Coq < Lemma l10 : p10. (* BND(1 / 2, [0.5, 0.5]) *)
1 subgoal

============================
p10

l10 < assert (h0 : p3).
2 subgoals

============================
p3

subgoal 2 is:
p10

l10 < apply l3.
1 subgoal

h0 : p3
============================
p10

l10 < assert (h1 : p11).
2 subgoals

h0 : p3
============================
p11

subgoal 2 is:
p10

l10 < apply l11.
1 subgoal

h0 : p3
h1 : p11
============================
p10

l10 < apply t6.
2 subgoals

h0 : p3
h1 : p11
============================
p3

subgoal 2 is:
p11

l10 < exact h0.
1 subgoal

h0 : p3
h1 : p11
============================
p11

l10 < exact h1.
Proof completed.

l10 < Qed.
assert p3 as h0.
apply l3.

assert p11 as h1.
apply l11.

apply t6.
exact h0.

exact h1.

l10 is defined

Coq < Print l10.
l10 = let h0 := l3 in let h1 := l11 in t6 h0 h1
: p10

Coq < Check l10.
l10
: p10

Coq <
-----------------------------------
Лемма l10 построила полное дерево вывода формулы p10, которое получается пристыковкой к неполному дереву из
леммы t6 нужных поддеревьев из лемм l3 и l11 .

Все это кажется избыточным и действительно, при желании построить доказательство леммы l10 , которое сразу строит полное правило вывода
--------------------------------
Coq < Lemma l10_full : p10 . (* BND(1 / 2, [0.5, 0.5]) *)
1 subgoal

============================
p10

l10_full < refine (div_pp _ _ _ _ _ l3 l11 _).
1 subgoal

============================
div_pp_helper i3 i7 i6 = true

l10_full < finalize.
Proof completed.

l10_full < Qed.
refine (div_pp _ _ _ _ _ l3 l11 _).
finalize.

l10_full is defined

Coq < Print l10_full.
l10_full =
div_pp r9 r8 i3 i7 i6 l3 l11 (refl_equal true<:div_pp_helper i3 i7 i6 = true)
: p10

Coq < Check l10_full.
l10_full
: p10

Coq <
--------------------------------

И даже совсем не использующее других лемм

-------------------------------
Coq < Lemma l10_very_full : p10 . (* BND(1 / 2, [0.5, 0.5]) *)
1 subgoal

============================
p10

l10_very_full < refine (div_pp _ _ _ _ _
l10_very_full < (constant1 _ (makepairF (Float2 1 0) (Float2 1 0)) _)
l10_very_full < (constant1 _ (makepairF (Float2 1 1) (Float2 1 1)) _)
l10_very_full < _); finalize.
Proof completed.

l10_very_full < Qed.
refine
(div_pp _ _ _ _ _ (constant1 _ (makepairF (Float2 1 0) (Float2 1 0)) _)
(constant1 _ (makepairF (Float2 1 1) (Float2 1 1)) _) _); finalize.

l10_very_full is defined

Coq < Print l10_very_full.
l10_very_full =
div_pp r9 r8 (makepairF (Float2 1 0) (Float2 1 0))
(makepairF (Float2 1 1) (Float2 1 1)) i6
(constant1 1 (makepairF (Float2 1 0) (Float2 1 0))
(refl_equal true
<:constant2_helper (Float2 1 0) (makepairF (Float2 1 0) (Float2 1 0)) =
true))
(constant1 2 (makepairF (Float2 1 1) (Float2 1 1))
(refl_equal true
<:constant2_helper (Float2 2 0) (makepairF (Float2 1 1) (Float2 1 1)) =
true))
(refl_equal true
<:div_pp_helper (makepairF (Float2 1 0) (Float2 1 0))
(makepairF (Float2 1 1) (Float2 1 1)) i6 = true)
: p10

Coq < Check l10_very_full.
l10_very_full
: p10

Coq <
------------------------------


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

Ее результат e.v можно найти на SVN.

В этих упражнения я немного понял, что такое Gappa и Coq и в следующий раз попробую сформулировать задачу верификации подпрограммы java.lang.StrictMath.log() .
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 10, 2010 11:13 by Dmitry Nadezhin
» Permalink
Verification of double StrictMath.log(double x)


I choose the static method java.lang.StrictMath.log(double x) as an object of verification.

JDK contains two classes for computation of elementary functions - Math and StrictMath .
Math is specified in terms of approximation quality (accuracy of the returned result and monotonicity of the
method).
The specification of StrictMath requires the same result as cirtain published C-algorithm.

Я выбираю в качестве объекта верификации подпрограмму натурального логарифма из библиотеки java.lang.StrictMath .

В JDK есть два класса для вычисления элементарных функций - Math и StrictMath.
Math специфицирован в терминах качества аппроксимации (ошибка результата и монотонность).
Спецификация StrictMath требует, чтобы результат совпадал с точностью до бита с результатом, выдаваемой
известной C-библиотекой.

--------------- Math.java
The class {@code Math} contains methods for performing basic
numeric operations such as the elementary exponential, logarithm,
square root, and trigonometric functions.

<p>Unlike some of the numeric methods of class
{@code StrictMath}, all implementations of the equivalent
functions of class {@code Math} are not defined to return the
bit-for-bit same results. This relaxation permits
better-performing implementations where strict reproducibility is
not required.

<p>By default many of the {@code Math} methods simply call
the equivalent method in {@code StrictMath} for their
implementation. Code generators are encouraged to use
platform-specific native libraries or microprocessor instructions,
where available, to provide higher-performance implementations of
{@code Math} methods. Such higher-performance
implementations still must conform to the specification for
{@code Math}.

<p>The quality of implementation specifications concern two
properties, accuracy of the returned result and monotonicity of the
method. Accuracy of the floating-point {@code Math} methods
is measured in terms of <i>ulps</i>, units in the last place. For
a given floating-point format, an ulp of a specific real number
value is the distance between the two floating-point values
bracketing that numerical value. When discussing the accuracy of a
method as a whole rather than at a specific argument, the number of
ulps cited is for the worst-case error at any argument. If a
method always has an error less than 0.5 ulps, the method always
returns the floating-point number nearest the exact result; such a
method is <i>correctly rounded</i>. A correctly rounded method is
generally the best a floating-point approximation can be; however,
it is impractical for many floating-point methods to be correctly
rounded. Instead, for the {@code Math} class, a larger error
bound of 1 or 2 ulps is allowed for certain methods. Informally,
with a 1 ulp error bound, when the exact result is a representable
number, the exact result should be returned as the computed result;
otherwise, either of the two floating-point values which bracket
the exact result may be returned. For exact results large in
magnitude, one of the endpoints of the bracket may be infinite.
Besides accuracy at individual arguments, maintaining proper
relations between the method at different arguments is also
important. Therefore, most methods with more than 0.5 ulp errors
are required to be <i>semi-monotonic</i>: whenever the mathematical
function is non-decreasing, so is the floating-point approximation,
likewise, whenever the mathematical function is non-increasing, so
is the floating-point approximation. Not all approximations that
have 1 ulp accuracy will automatically meet the monotonicity
requirements.
------------
------------ StrictMath.java
The class {@code StrictMath} contains methods for performing basic
numeric operations such as the elementary exponential, logarithm,
square root, and trigonometric functions.

<p>To help ensure portability of Java programs, the definitions of
some of the numeric functions in this package require that they
produce the same results as certain published algorithms. These
algorithms are available from the well-known network library
{@code netlib} as the package "Freely Distributable Math
Library," <a
href="ftp://ftp.netlib.org/fdlibm.tar">{@code fdlibm}</a>. These
algorithms, which are written in the C programming language, are
then to be understood as executed with all floating-point
operations following the rules of Java floating-point arithmetic.
--------------

Here is the specification of log() method
Вот как выглядит спецификация функции log()
---------- Math.java
/**
* Returns the natural logarithm (base <i>e</i>) of a {@code double}
* value. Special cases:
* <ul><li>If the argument is NaN or less than zero, then the result
* is NaN.
* <li>If the argument is positive infinity, then the result is
* positive infinity.
* <li>If the argument is positive zero or negative zero, then the
* result is negative infinity.</ul>
*
* <p>The computed result must be within 1 ulp of the exact result.
* Results must be semi-monotonic.
*
* @param a a value
* @return the value ln&nbsp;{@code a}, the natural logarithm of
* {@code a}.
*/
public static double log(double a) {
return StrictMath.log(a); // default impl. delegates to StrictMath
}
--------------

Here is the implementation of log() in the reference C-library
Вот как реализована функция log() в этой библиотеке.
--------------

/* @(#)e_log.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/

/* __ieee754_log(x)
* Return the logrithm of x
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/

#include "fdlibm.h"

#ifdef __STDC__
static const double
#else
static double
#endif
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */

static double zero = 0.0;

#ifdef __STDC__
double __ieee754_log(double x)
#else
double __ieee754_log(x)
double x;
#endif
{
double hfsq,f,s,z,R,w,t1,t2,dk;
int k,hx,i,j;
unsigned lx;

hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */

k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
__HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
k += (i>>20);
f = x-1.0;
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
if(f==zero) if(k==0) return zero; else {dkFrowndouble)k;
return dk*ln2_hi+dk*ln2_lo;}
R = f*f*(0.5-0.33333333333333333*f);
if(k==0) return f-R; else {dkFrowndouble)k;
return dk*ln2_hi-((R-dk*ln2_lo)-f);}
}
s = f/(2.0+f);
dk = (double)k;
z = s*s;
i = hx-0x6147a;
w = z*z;
j = 0x6b851-hx;
t1= w*(Lg2+w*(Lg4+w*Lg6));
t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
i |= j;
R = t2+t1;
if(i>0) {
hfsq=0.5*f*f;
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
} else {
if(k==0) return f-s*(f-R); else
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}
--------------

I'm going to solve the following tasks:
Я собираюсь решить следующие задачи

1) to write functional equivalent of this algorithm in the Coq input language.
1) Написать функциональный эквивалент этого алгоритма на языке Coq .

2) To prove that Coq-algorithm satisfies quality requirements:
- accuracy of 1 ULP in whole argument range;
- monotonicity of the method;
2) Доказать, что этот функциональный эквивалент удовлетворяет требованиям:
- ошибка результата не более 1 ULP на всем диапазоне аргумента;
- результат монотонно возрастает от аргумента

3) To prove that Coq algorithm returs the same results as the reference C algorithm. This is either by
- to use algorithm in C notation and to use somebody's definition of C language in Coq language;
- to compile the algorithm from C to object file of some processor and to create formal description
of processor architecture subset that is used in the object file.
Currently I prefer the second approach with MIPS architecture.
MIPS is one of the simplest architecture and I like the NestedVM project
3) Доказать, что алгоритм на языке Coq всегда дает те же результаты, что и исходный алгоритм.
Для этого придется
- либо использовать алгоритм на C кем-то написанное определение языка C на языке Coq .
- либо странслировать алгоритм из C в объектный файл какого-то процесора и построить формальное описание используемое в объектном файл подмножество системы команд этого процессора.
На данный момент я склоняюсь ко второму варианту с трансляцией в систему команд MIPS.
MIPS я выбираю, потому что это одна из самых простейших систем команд и потому-что мне понравился проект
NestedVM.

4) If the algorithm will be implemented in Java by this time, the to prove that JVM running this algorithm will return the same results as the Coq algorithm.
Again two variants:
- to use the file StrictMath.java and somebody's description of Java language in the Coq language;
- to use the file StrictMath.class and to define Java Virtual Machine in the Coq language.
4) Если к тому времени этот алгоритм будет реализован на языке Java, то доказать что алгоритм, исполняемый на JVM, будет давать те же результаты, что и алгоритм, записанный на Coq.
Тут опять два варианта:
- использовать файл StrictMath.java и кем-то написанное определения языка Java на Coq;
- использовать файл StrictMath.class и определение Java Virtual Machine на языке Coq .
 
Dmitry Nadezhin
Nadezhin__small
Posted: March 13, 2010 05:16 by Dmitry Nadezhin
» Permalink
Write algorithm from e_log2.c in Gappa and Coq


I rewrote a part of algorithm from FDLIBM e_log2.c in Gappa
------------------ log1.g
@rnd = float< ieee_64, ne >;

x = rnd(x_);
dk = int<ne>(k_);

ln2_hi = rnd(6.93147180369123816490e-01); # /* 3fe62e42 fee00000 */
ln2_lo = rnd(1.90821492927058770002e-10); # /* 3dea39ef 35793c76 */
Lg1 = rnd(6.666666666666735130e-01); # /* 3FE55555 55555593 */
Lg2 = rnd(3.999999999940941908e-01); # /* 3FD99999 9997FA04 */
Lg3 = rnd(2.857142874366239149e-01); # /* 3FD24924 94229359 */
Lg4 = rnd(2.222219843214978396e-01); # /* 3FCC71C5 1D8E78AF */
Lg5 = rnd(1.818357216161805012e-01); # /* 3FC74664 96CB03DE */
Lg6 = rnd(1.531383769920937332e-01); # /* 3FC39A09 D078C69F */
Lg7 = rnd(1.479819860511658591e-01); # /* 3FC2F112 DF3E5244 */

c2 = rnd(2.0);
c0_5 = rnd(0.5);
c0_333 = rnd(0.33333333333333333);

f rnd= x - 1;
R_1 rnd= f*f*(c0_5 - f*c0_333);
ln_1 rnd= dk*ln2_hi - ((R_1 - dk*ln2_lo) - f);
s rnd= f/(c2+f);
z rnd= s*s;
w rnd= z*z;
t1 rnd= w*(Lg2 + w*(Lg4 + w*Lg6));
t2 rnd= z*(Lg1 + w*(Lg3 + w*(Lg5 + w*Lg7)));
R rnd= t2 + t1;
ln_2 rnd= dk*ln2_hi - ((s*(f - R) - dk*ln2_lo) - f);

hfsq rnd= c0_5*f*f;
ln_3 rnd= dk*ln2_hi - ((hfsq - (s*(hfsq + R) + dk*ln2_lo)) - f);

# xk = x*2^k
# x in [ 0x1.6a09c00000000p-1 , 0x1.6a09bffffffffp0 ]
#
# x in [ 0x1.6a09c00000000p-1 , 0x1.6b851ffffffffp-1] -> ln_3
# x in [ 0x1.6b85200000000p-1 , 0x1.ffffdffffffffp-1] -> ln_2
# x in [ 0x1.ffffe00000000p-1 , 0x1.00000ffffffffp0 ] -> ln_1
# x in [ 0x1.0000100000000p0 , 0x1.61479ffffffffp0 ] -> ln_2
# x in [ 0x1.6147a00000000p0 , 0x1.6a09bffffffffp0 ] -> ln_3

x1 = x - 1;
ln = x1*(1/1 - x1*(1/2 - x1*(1/3 - x1*r/4)));

ln_1a rnd= f - R_1;
ln_2a rnd= f + s*(R-f);
ln_3a rnd= f + (s*(hfsq+R) - hfsq);

R_1_ = x1*x1*(1/2 - x1*(1/3 - x1*r/4));
ln_ = x1 - R_1_;

threshold = fixed<-20, dn> (sqrt(2));

{ (dk in [0,0] /\ x in [ 0x1.0000000000001p0 , 0x1.0000000000001p0 ] /\ r in [0,1] -> ln_1a in ? /\ ln_1a - ln_ in ? /\ (ln_1a - ln_)*1b105 in ? ) /\
(dk in [0,0] /\ x in [ 0x1.00000ffffffffp0 , 0x1.00000ffffffffp0 ] /\ r in [0,1] -> ln_1a in ? /\ ln_1a - ln_ in ? /\ (ln_1a - ln_)*1b73 in ? ) /\
(dk in [0,0] /\ x in [ 0x1.0000000000001p0 , 0x1.00000ffffffffp0 ] /\ r in [0,1] -> (R_1 - R_1_)/x1 in ?)
}

(R_1 - R_1_)/x1 -> (R_1 - R_1_)/R_1_*x1*(1/2 - x1*(1/3 - x1*r/4));
-----------------
Gappa code doesn't include algorithm reduction.
Also it doesn't include argument range splittion.
There are expressions ln_1, ln_2, ln_3 which encodes C algorithm for different ranges.
And there is a comment that says which expression should be used in every range.

Some gappa test code computes the algorithm for two particular input.
Also it computesa bounding for R_1 error on one of ranges.

Also I wrote a coq code
------------------ fdlibm_log2.v
Require Import Reals.
Require Import Gappa_library.
Open Local Scope R_scope.

Definition o(x: R) := (rounding_float roundNE 53 1074) x.

Definition threshold := Float2 370727 (-18). (* 0x1.6a09cp0 *)
Definition lower1 := Float2 1048575 (-20). (* 0x1.ffffep-1 *)
Definition upper1 := Float2 1048577 (-20). (* 0x1.00001p0 *)
Definition lower2 := Float2 744489 (-20). (* 0x1.6b852p-1 *)
Definition upper2 := Float2 723517 (-19). (* 0x1.6147ap0 *)

Definition ln2_hi := float2R (o(Float10 69314718036912381649 (-20))).
Definition ln2_lo := float2R (o(Float10 190821492927058770002 (-30))).
Definition Lg1 := float2R (o(Float10 666666666666673513 (-18))).
Definition Lg2 := float2R (o(Float10 3999999999940941908 (-19))).
Definition Lg3 := float2R (o(Float10 2857142874366239149 (-19))).
Definition Lg4 := float2R (o(Float10 2222219843214978396 (-19))).
Definition Lg5 := float2R (o(Float10 1818357216161805012 (-19))).
Definition Lg6 := float2R (o(Float10 1531383769920937332 (-19))).
Definition Lg7 := float2R (o(Float10 (1479819860511658591) (-19))).
Definition c2 := float2R (o(Float10 20 (-1))).
Definition c0_5 := float2R (o(Float10 5 (-1))).
Definition c0_333 := float2R (o(Float10 33333333333333333 (-17))).

Record posfloat2 : Type := mkposfloat2 {pos :> float2; cond_pos : Gappa_dyadic.Fpos pos = true }.

Definition log2_reduce (x: posfloat2) : float2*Z :=
let k0 := Zplus (Zlogarithm.N_digits (Fnum x)) (Fexp x) in
let k := if Gappa_dyadic.Flt2 x (Float2 (Fnum threshold) (Fexp threshold + k0)) then k0 else Zpred k0 in
(Float2 (Fnum x) (Fexp x - k), k).

Definition fdlibm_log2 (x: posfloat2 ) : float2 :=
let (x',k) := log2_reduce x in
let f := o(x' - 1) in
let dk := IZR k in
if (Gappa_dyadic.Fle2 lower1 x') && (Gappa_dyadic.Flt2 x' upper1) then
let R := o(o(f*f) * o(c0_5 - o(f*c0_333)))
in o(o(dk*ln2_hi) - o(o(R - o(dk*ln2_lo)) - f))
else
let s := o(f/o(c2+f)) in
let z := o(s*s) in
let w := o(z*z) in
let t1 := o(w*o(Lg2 + o(w*o(Lg4 + o(w*Lg6))))) in
let t2 := o(z*o(Lg1 + o(w*o(Lg3 + o(w*o(Lg5 + o(w*Lg7))))))) in
let R := o(t2 + t1) in
if (Gappa_dyadic.Fle2 lower2 x') && (Gappa_dyadic.Flt2 x' upper2) then
o(o(dk*ln2_hi) - o(o(o(s*o(f - R)) - o(dk*ln2_lo)) - f))
else
let hfsq := o(o(c0_5*f)*f) in
o(o(dk*ln2_hi) - o(o(hfsq - o(o(s*o(hfsq + R)) + o(dk*ln2_lo))) - f))
.
------------
This algorithm does the argument reduction and case splitting.
However it handles only finite positive arguments.

The next step will be to prove some facts about Coq algorithm combining Coq and Gappa tool.
showing 1 - 8 of 8
Replies: 7 - Last Post: March 13, 2010 05:16
by: Dmitry Nadezhin
« Older Topic » Newer Topic
  • Mysql
  • Glassfish
  • Jruby
  • Rails
  • Nblogo
Terms of Use; Privacy Policy;
© 2010, Oracle Corporation and/or its affiliates
(revision 20120518.3c65429)
 
 
Close
loading
Please Confirm
Close