00001 /** This is the main header file of the SCS library, which defines the 00002 SCS data structure, and the functions that implement arithmetic on it. 00003 00004 @file scs.h 00005 00006 @author David Defour David.Defour@ens-lyon.fr 00007 @author Florent de Dinechin Florent.de.Dinechin@ens-lyon.fr 00008 00009 This file is part of the SCS library. 00010 00011 Copyright (C) 2002 David Defour and Florent de Dinechin 00012 00013 This library is free software; you can redistribute it and/or 00014 modify it under the terms of the GNU Lesser General Public 00015 License as published by the Free Software Foundation; either 00016 version 2.1 of the License, or (at your option) any later version. 00017 00018 This library is distributed in the hope that it will be useful, 00019 but WITHOUT ANY WARRANTY; without even the implied warranty of 00020 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00021 Lesser General Public License for more details. 00022 00023 You should have received a copy of the GNU Lesser General Public 00024 License along with this library; if not, write to the Free Software 00025 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00026 00027 */ 00028 00029 00030 00031 /* Avoid loading the header twice */ 00032 #ifndef INCLUDE_SCS 00033 #define INCLUDE_SCS 1 00034 00035 #ifndef DOXYGEN_SHOULD_SKIP_THIS /* because it is not very clean */ 00036 00037 #include <limits.h> 00038 00039 #include "scs_config.h" 00040 00041 #ifdef HAVE_GMP_H 00042 #include <gmp.h> 00043 #endif 00044 00045 #ifdef HAVE_MPFR_H 00046 #include <mpfr.h> 00047 #endif 00048 00049 00050 #endif /* DOXYGEN_SHOULD_SKIP_THIS */ 00051 00052 00053 00054 00055 00056 00057 /** @internal An union to cast floats into doubles or the other way round. For 00058 internal purpose only */ 00059 00060 typedef union { 00061 int i[2]; 00062 unsigned long long int l; 00063 double d; 00064 } scs_db_number; 00065 00066 00067 00068 00069 /* ****************************************************************** */ 00070 /**@name SCS data-types */ /**@{*/ 00071 00072 /** @struct scs 00073 The SCS data type. 00074 00075 An SCS number is a a floating-point number in base 2^32. 00076 00077 - Its mantissa is formed of SCS_NB_WORDS digits (currently 32 bits by default) 00078 00079 - Its exponent is a 32-bit integer 00080 00081 - It also has a sign field, and an exception field used to store and 00082 propagate IEEE-754 exceptions. 00083 00084 00085 The real number represented by a scs structure is equal to: 00086 @f$ 00087 \displaystyle 00088 \sum_{i=0}^{\mathtt{SCS\_NB\_WORDS}} 00089 2^{(\mathtt{index} -i)\mathtt{SCS\_NB\_BITS}} 00090 \times 00091 \mathtt{h\_word}[i] 00092 @f$ 00093 */ 00094 00095 /* 00096 (verbatim-mode formula for the above eqation:) the number represented by a 00097 SCS structure is : 00098 00099 i<SCS_NB_WORDS (index - i).SCS_NB_BITS 00100 sign . ( sum ( h_word[i] . 2^ ) 00101 i=0 00102 */ 00103 00104 struct scs { 00105 /** the digits, as 32 bits words */ 00106 unsigned int h_word[SCS_NB_WORDS]; 00107 /** Used to store Nan,+/-0, Inf, etc and then let the hardware handle them */ 00108 scs_db_number exception; 00109 /** This corresponds to the exponent in an FP format, but here we are 00110 in base 2^32 */ 00111 int index; 00112 /** The sign equals 1 or -1*/ 00113 int sign; 00114 }; 00115 00116 00117 typedef struct scs scs; 00118 00119 00120 00121 /** scs_ptr is a pointer on a SCS structure */ 00122 typedef struct scs * scs_ptr; 00123 00124 00125 00126 00127 /** scs_t is an array of one SCS struct to lighten syntax : you may 00128 declare a scs_t object, and pass it to the scs functions (which 00129 expect pointers) without using ampersands. 00130 */ 00131 typedef struct scs scs_t[1]; 00132 00133 /**@}*/ /* end doxygen group for SCS data-types */ 00134 00135 00136 00137 00138 00139 00140 00141 00142 /* ****************************************************************** */ 00143 /**@name Conversion and initialization functions */ /**@{*/ 00144 00145 /** Convert a SCS number to a double, rounding to the nearest */ 00146 void scs_get_d(double*, scs_ptr); 00147 00148 /** Convert a SCS number to a double, rounding towards minus infinity */ 00149 void scs_get_d_minf(double*, scs_ptr); 00150 00151 /** Convert a SCS number to a double, rounding towards plus infinity */ 00152 void scs_get_d_pinf(double*, scs_ptr); 00153 00154 /** Convert a SCS number to a double, rounding towards zero */ 00155 void scs_get_d_zero(double*, scs_ptr); 00156 00157 /** Convert a double into a SCS number (this is an exact operation) */ 00158 void scs_set_d(scs_ptr, double); 00159 00160 /** Convert an unsigned int into a SCS number (this is an exact operation) */ 00161 void scs_set_si(scs_ptr, signed int); 00162 00163 00164 /** Print out a SCS number. Sorry for the strange name, we are mimicking GMP */ 00165 void scs_get_std(scs_ptr); 00166 00167 00168 /** Copy a SCS number into another */ 00169 void scs_set(scs_ptr, scs_ptr); 00170 00171 00172 /** Set a SCS number to zero */ 00173 void scs_zero(scs_ptr); 00174 00175 00176 /** Generate a random SCS number. 00177 The index field of result will be between -expo_max and +expo_max. 00178 Example: to get a number in the double-precision floating-point range, 00179 expo_max should be smaller than 39. 00180 @warning No guarantee is made about the quality of the random algorithm 00181 used... */ 00182 void scs_rand(scs_ptr result, int expo_max); 00183 00184 /**@}*/ /* end doxygen group for conversion / initialisation functions*/ 00185 00186 00187 00188 00189 /* ****************************************************************** */ 00190 /**@name Addition and renormalisation functions */ /**@{*/ 00191 00192 00193 /** Addition of two SCS numbers. 00194 The arguments x, y and result may point to the same memory 00195 location. The result is a normalised SCS number. 00196 */ 00197 void scs_add(scs_ptr result, scs_ptr x, scs_ptr y); 00198 00199 00200 /** Subtraction of two SCS numbers. 00201 The arguments x, y and result may point to the same memory 00202 location. The result is a normalised SCS number. 00203 */ 00204 void scs_sub(scs_ptr result, scs_ptr x, scs_ptr y); 00205 00206 00207 /** Addition without renormalisation, to be used for adding many 00208 numbers. 00209 @warning In case of a cancellation, severe loss of precision could 00210 happen. Safe if the numbers are of the same sign. 00211 */ 00212 void scs_add_no_renorm(scs_ptr result, scs_ptr x, scs_ptr y); 00213 00214 00215 /** Renormalisation (to be used after several scs_add_no_renorm). 00216 This function removes the carry from each digit, and also shifts the 00217 digits in case of a cancellation (so that if result != 0 then its 00218 first digit is non-zero) 00219 00220 @warning THIS FUNCTION HAS NEVER BEEN PROPERLY TESTED and is 00221 currently unused in the library: instead, specific renormalisation 00222 steps are fused within the code of the operations which require it. 00223 */ 00224 00225 void scs_renorm(scs_ptr); 00226 00227 00228 /** Renormalisation assuming no cancellation. This function is useful 00229 for example when adding many numbers of the same sign */ 00230 void scs_renorm_no_cancel_check(scs_ptr); 00231 00232 /**@}*/ /* end doxygen group for addition and normalisation functions*/ 00233 00234 00235 00236 00237 /* ****************************************************************** */ 00238 /**@name Multiplication functions */ /**@{*/ 00239 00240 /** Multiplication of two SCS numbers. The arguments x, y and result 00241 may point to the same memory location. The result is a normalised SCS 00242 number. 00243 */ 00244 void scs_mul(scs_ptr result, const scs_ptr x, const scs_ptr y); 00245 00246 /** Multiplication of a SCS with an unsigned integer; result is 00247 returned in x. */ 00248 void scs_mul_ui(scs_ptr, const unsigned int); 00249 00250 /** Square. Result is normalised */ 00251 void scs_square(scs_ptr result, scs_ptr x); 00252 00253 /** Fused multiply-and-add (ab+c); Result is normalised 00254 \warning This function has not been tested thoroughly */ 00255 void scs_fma(scs_ptr result, scs_ptr a, scs_ptr b, scs_ptr c); 00256 00257 /**@}*/ /* end doxygen group for Multiplication functions*/ 00258 00259 00260 00261 00262 00263 /* ****************************************************************** */ 00264 /**@name Divisions */ /**@{*/ 00265 00266 /** SCS inverse. 00267 Stores 1/x in result. Result is normalised 00268 00269 @warning This function is known not to work for most precisions: it 00270 performs a fixed number of Newton-Raphson iterations (two), starting 00271 with a FP number (53 bits), so provides roughly 210 bits of 00272 precision. It should be modified to perform more iterations if more 00273 precision is needed. 00274 */ 00275 void scs_inv(scs_ptr result, scs_ptr x); 00276 00277 /** SCS division. Computes x/y. Result is normalised 00278 @warning This function is known not to work for most precisions: it 00279 performs a fixed number of Newton-Raphson iterations (two), starting 00280 with a FP number (53 bits), so provides roughly 210 bits of 00281 precision. It should be modified to perform more iterations if more 00282 precision is needed. 00283 */ 00284 void scs_div(scs_ptr result, scs_ptr x, scs_ptr y); 00285 00286 /**@}*/ /* end doxygen group for division functions*/ 00287 00288 00289 00290 00291 00292 /* ****************************************************************** */ 00293 /**@name Functions for testing purpose */ /**@{*/ 00294 00295 00296 #ifdef HAVE_LIBGMP 00297 /** Convert a SCS number into a GMP MPF (multiple precision, 00298 floating-point) number. Should be exact if the target number has 00299 more precision than the SCS number, otherwise the rounding is 00300 unspecified (the conversion uses MPF functions) */ 00301 void scs_get_mpf(scs_ptr, mpf_t); 00302 #endif 00303 00304 00305 00306 #ifdef HAVE_MPFR_H 00307 /** Convert a SCS number into a MPFR (multiple precision, 00308 floating-point) number. Should be exact if the target number has 00309 more precision than the SCS number, otherwise should be correctly 00310 rounded (the conversion uses MPFR functions). Not heavily tested 00311 though */ 00312 void scs_get_mpfr(scs_ptr, mpfr_t); 00313 #endif 00314 00315 00316 /**@}*/ /* end doxygen group for functions for testing purpose */ 00317 00318 #endif /* INCLUDE_SCS */ 00319 00320 00321 00322 00323