1 /* $OpenBSD: header.h,v 1.1 2009/04/09 01:24:43 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2009 Gaston H. Gonnet <gonnet@inf.ethz.ch> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* C program to test the accuracy of a given function. 20 This will be done in terms of ulps. A ulp is a "unit in 21 the last place". No program can guarantee results better 22 than 0.5 ulp. As a matter of fact, 0.5 ulp could require 23 infinite computation. So very good implementations 24 usually achieve 0.51 ulp. 0.6 ulp is still very good, 1 ulp 25 starts to get a bit sloppy. 26 27 The model of computation is the standard one, the result of f(x) 28 is the correctly rounded value of f(x) when x is assumed to be 29 an exact rational. 30 31 This program selects values which are difficult to compute to 32 measure the quality of the computation. 33 34 Hardware which has floating point registers with extra bits of 35 precision (e.g. Pentium) usually do much better than others. 36 This function computes the accuracy for both cases. */ 37 38 /* We represent double precision floating point numbers x as a 39 pair of an integer and a base 2 exponent. x = mant*2^expo, 40 or x = scalb(mant,expo). 2^52 <= mant < 2^53. This is to 41 avoid sloppy compilers which do not transcribe correctly 42 some input doubles */ 43 44 /* scalb( F(arg[i]), val_e[i] ) == val[i]+eps[i] 45 to about twice precision. 46 The scaling is done so that val[i] is in ulps 47 (val is consequently an integer) 48 and to avoid underflows in certain cases. */ 49 50 /* Shell sort, needs the definition of KEY(x) 51 From: Gonnet & Baeza, Hanbook of Algorithms and Data structures, 4.1.4. 52 Can sort any array r from r[lo] to r[hi] in place. 53 Uses a temporary variable t, of the same type as r[]. 54 The comparison is done between KEY(r[i]) vs KEY(r[j]), 55 so the macro KEY(x) has to be defined appropriately */ 56 #define SORT(r,lo,up,t) \ 57 { int SORTd, SORTi, SORTj; \ 58 for ( SORTd=(up)-(lo)+1; SORTd>1; ) { \ 59 SORTd = SORTd<5 ? 1 : (5*SORTd-1)/11; \ 60 for ( SORTi=(up)-SORTd; SORTi>=(lo); SORTi-- ) { \ 61 (t) = (r)[SORTi]; \ 62 for ( SORTj=SORTi+SORTd; SORTj<=(up) && KEY((t)) > KEY((r)[SORTj]); SORTj+=SORTd ) \ 63 (r)[SORTj-SORTd] = (r)[SORTj]; \ 64 (r)[SORTj-SORTd] = (t); \ 65 } \ 66 }} 67 68 double scalb( double x, double n ); 69