Stephen Paulger / ttmath

Library for big numbers from http://www.ttmath.org/

Dependents:   PIDHeater82 Conceptcontroller_v_1_0 AlarmClockApp COG4050_adxl355_tilt ... more

TTMath is a small library which allows one to perform arithmetic operations with big unsigned integer, big signed integer and big floating point numbers. It provides standard mathematical operations like adding, subtracting, multiplying, dividing.

TTMath is BSD Licensed (new/modified BSD)

For more information about ttmath see http://www.ttmath.org/

ttmath.h@0:04a9f72bbca7, 2013-07-30 (annotated)

Committer:
stevep
Date:
Tue Jul 30 18:43:48 2013 +0000
Revision:
0:04a9f72bbca7
v0.9.3 of ttmath

Who changed what in which revision?

UserRevisionLine numberNew contents of line
stevep 0:04a9f72bbca7 1 /*
stevep 0:04a9f72bbca7 2 * This file is a part of TTMath Bignum Library
stevep 0:04a9f72bbca7 3 * and is distributed under the (new) BSD licence.
stevep 0:04a9f72bbca7 4 * Author: Tomasz Sowa <t.sowa@ttmath.org>
stevep 0:04a9f72bbca7 5 */
stevep 0:04a9f72bbca7 6
stevep 0:04a9f72bbca7 7 /*
stevep 0:04a9f72bbca7 8 * Copyright (c) 2006-2012, Tomasz Sowa
stevep 0:04a9f72bbca7 9 * All rights reserved.
stevep 0:04a9f72bbca7 10 *
stevep 0:04a9f72bbca7 11 * Redistribution and use in source and binary forms, with or without
stevep 0:04a9f72bbca7 12 * modification, are permitted provided that the following conditions are met:
stevep 0:04a9f72bbca7 13 *
stevep 0:04a9f72bbca7 14 * * Redistributions of source code must retain the above copyright notice,
stevep 0:04a9f72bbca7 15 * this list of conditions and the following disclaimer.
stevep 0:04a9f72bbca7 16 *
stevep 0:04a9f72bbca7 17 * * Redistributions in binary form must reproduce the above copyright
stevep 0:04a9f72bbca7 18 * notice, this list of conditions and the following disclaimer in the
stevep 0:04a9f72bbca7 19 * documentation and/or other materials provided with the distribution.
stevep 0:04a9f72bbca7 20 *
stevep 0:04a9f72bbca7 21 * * Neither the name Tomasz Sowa nor the names of contributors to this
stevep 0:04a9f72bbca7 22 * project may be used to endorse or promote products derived
stevep 0:04a9f72bbca7 23 * from this software without specific prior written permission.
stevep 0:04a9f72bbca7 24 *
stevep 0:04a9f72bbca7 25 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
stevep 0:04a9f72bbca7 26 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
stevep 0:04a9f72bbca7 27 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
stevep 0:04a9f72bbca7 28 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
stevep 0:04a9f72bbca7 29 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
stevep 0:04a9f72bbca7 30 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
stevep 0:04a9f72bbca7 31 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
stevep 0:04a9f72bbca7 32 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
stevep 0:04a9f72bbca7 33 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
stevep 0:04a9f72bbca7 34 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
stevep 0:04a9f72bbca7 35 * THE POSSIBILITY OF SUCH DAMAGE.
stevep 0:04a9f72bbca7 36 */
stevep 0:04a9f72bbca7 37
stevep 0:04a9f72bbca7 38
stevep 0:04a9f72bbca7 39
stevep 0:04a9f72bbca7 40 #ifndef headerfilettmathmathtt
stevep 0:04a9f72bbca7 41 #define headerfilettmathmathtt
stevep 0:04a9f72bbca7 42
stevep 0:04a9f72bbca7 43 /*!
stevep 0:04a9f72bbca7 44 \file ttmath.h
stevep 0:04a9f72bbca7 45 \brief Mathematics functions.
stevep 0:04a9f72bbca7 46 */
stevep 0:04a9f72bbca7 47
stevep 0:04a9f72bbca7 48 #ifdef _MSC_VER
stevep 0:04a9f72bbca7 49 //warning C4127: conditional expression is constant
stevep 0:04a9f72bbca7 50 #pragma warning( disable: 4127 )
stevep 0:04a9f72bbca7 51 //warning C4702: unreachable code
stevep 0:04a9f72bbca7 52 #pragma warning( disable: 4702 )
stevep 0:04a9f72bbca7 53 //warning C4800: forcing value to bool 'true' or 'false' (performance warning)
stevep 0:04a9f72bbca7 54 #pragma warning( disable: 4800 )
stevep 0:04a9f72bbca7 55 #endif
stevep 0:04a9f72bbca7 56
stevep 0:04a9f72bbca7 57
stevep 0:04a9f72bbca7 58 #include "ttmathbig.h"
stevep 0:04a9f72bbca7 59 #include "ttmathobjects.h"
stevep 0:04a9f72bbca7 60
stevep 0:04a9f72bbca7 61
stevep 0:04a9f72bbca7 62 namespace ttmath
stevep 0:04a9f72bbca7 63 {
stevep 0:04a9f72bbca7 64 /*
stevep 0:04a9f72bbca7 65 *
stevep 0:04a9f72bbca7 66 * functions defined here are used only with Big<> types
stevep 0:04a9f72bbca7 67 *
stevep 0:04a9f72bbca7 68 *
stevep 0:04a9f72bbca7 69 */
stevep 0:04a9f72bbca7 70
stevep 0:04a9f72bbca7 71
stevep 0:04a9f72bbca7 72 /*
stevep 0:04a9f72bbca7 73 *
stevep 0:04a9f72bbca7 74 * functions for rounding
stevep 0:04a9f72bbca7 75 *
stevep 0:04a9f72bbca7 76 *
stevep 0:04a9f72bbca7 77 */
stevep 0:04a9f72bbca7 78
stevep 0:04a9f72bbca7 79
stevep 0:04a9f72bbca7 80 /*!
stevep 0:04a9f72bbca7 81 this function skips the fraction from x
stevep 0:04a9f72bbca7 82 e.g 2.2 = 2
stevep 0:04a9f72bbca7 83 2.7 = 2
stevep 0:04a9f72bbca7 84 -2.2 = 2
stevep 0:04a9f72bbca7 85 -2.7 = 2
stevep 0:04a9f72bbca7 86 */
stevep 0:04a9f72bbca7 87 template<class ValueType>
stevep 0:04a9f72bbca7 88 ValueType SkipFraction(const ValueType & x)
stevep 0:04a9f72bbca7 89 {
stevep 0:04a9f72bbca7 90 ValueType result( x );
stevep 0:04a9f72bbca7 91 result.SkipFraction();
stevep 0:04a9f72bbca7 92
stevep 0:04a9f72bbca7 93 return result;
stevep 0:04a9f72bbca7 94 }
stevep 0:04a9f72bbca7 95
stevep 0:04a9f72bbca7 96
stevep 0:04a9f72bbca7 97 /*!
stevep 0:04a9f72bbca7 98 this function rounds to the nearest integer value
stevep 0:04a9f72bbca7 99 e.g 2.2 = 2
stevep 0:04a9f72bbca7 100 2.7 = 3
stevep 0:04a9f72bbca7 101 -2.2 = -2
stevep 0:04a9f72bbca7 102 -2.7 = -3
stevep 0:04a9f72bbca7 103 */
stevep 0:04a9f72bbca7 104 template<class ValueType>
stevep 0:04a9f72bbca7 105 ValueType Round(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 106 {
stevep 0:04a9f72bbca7 107 if( x.IsNan() )
stevep 0:04a9f72bbca7 108 {
stevep 0:04a9f72bbca7 109 if( err )
stevep 0:04a9f72bbca7 110 *err = err_improper_argument;
stevep 0:04a9f72bbca7 111
stevep 0:04a9f72bbca7 112 return x; // NaN
stevep 0:04a9f72bbca7 113 }
stevep 0:04a9f72bbca7 114
stevep 0:04a9f72bbca7 115 ValueType result( x );
stevep 0:04a9f72bbca7 116 uint c = result.Round();
stevep 0:04a9f72bbca7 117
stevep 0:04a9f72bbca7 118 if( err )
stevep 0:04a9f72bbca7 119 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 120
stevep 0:04a9f72bbca7 121 return result;
stevep 0:04a9f72bbca7 122 }
stevep 0:04a9f72bbca7 123
stevep 0:04a9f72bbca7 124
stevep 0:04a9f72bbca7 125
stevep 0:04a9f72bbca7 126 /*!
stevep 0:04a9f72bbca7 127 this function returns a value representing the smallest integer
stevep 0:04a9f72bbca7 128 that is greater than or equal to x
stevep 0:04a9f72bbca7 129
stevep 0:04a9f72bbca7 130 Ceil(-3.7) = -3
stevep 0:04a9f72bbca7 131 Ceil(-3.1) = -3
stevep 0:04a9f72bbca7 132 Ceil(-3.0) = -3
stevep 0:04a9f72bbca7 133 Ceil(4.0) = 4
stevep 0:04a9f72bbca7 134 Ceil(4.2) = 5
stevep 0:04a9f72bbca7 135 Ceil(4.8) = 5
stevep 0:04a9f72bbca7 136 */
stevep 0:04a9f72bbca7 137 template<class ValueType>
stevep 0:04a9f72bbca7 138 ValueType Ceil(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 139 {
stevep 0:04a9f72bbca7 140 if( x.IsNan() )
stevep 0:04a9f72bbca7 141 {
stevep 0:04a9f72bbca7 142 if( err )
stevep 0:04a9f72bbca7 143 *err = err_improper_argument;
stevep 0:04a9f72bbca7 144
stevep 0:04a9f72bbca7 145 return x; // NaN
stevep 0:04a9f72bbca7 146 }
stevep 0:04a9f72bbca7 147
stevep 0:04a9f72bbca7 148 ValueType result(x);
stevep 0:04a9f72bbca7 149 uint c = 0;
stevep 0:04a9f72bbca7 150
stevep 0:04a9f72bbca7 151 result.SkipFraction();
stevep 0:04a9f72bbca7 152
stevep 0:04a9f72bbca7 153 if( result != x )
stevep 0:04a9f72bbca7 154 {
stevep 0:04a9f72bbca7 155 // x is with fraction
stevep 0:04a9f72bbca7 156 // if x is negative we don't have to do anything
stevep 0:04a9f72bbca7 157 if( !x.IsSign() )
stevep 0:04a9f72bbca7 158 {
stevep 0:04a9f72bbca7 159 ValueType one;
stevep 0:04a9f72bbca7 160 one.SetOne();
stevep 0:04a9f72bbca7 161
stevep 0:04a9f72bbca7 162 c += result.Add(one);
stevep 0:04a9f72bbca7 163 }
stevep 0:04a9f72bbca7 164 }
stevep 0:04a9f72bbca7 165
stevep 0:04a9f72bbca7 166 if( err )
stevep 0:04a9f72bbca7 167 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 168
stevep 0:04a9f72bbca7 169 return result;
stevep 0:04a9f72bbca7 170 }
stevep 0:04a9f72bbca7 171
stevep 0:04a9f72bbca7 172
stevep 0:04a9f72bbca7 173 /*!
stevep 0:04a9f72bbca7 174 this function returns a value representing the largest integer
stevep 0:04a9f72bbca7 175 that is less than or equal to x
stevep 0:04a9f72bbca7 176
stevep 0:04a9f72bbca7 177 Floor(-3.6) = -4
stevep 0:04a9f72bbca7 178 Floor(-3.1) = -4
stevep 0:04a9f72bbca7 179 Floor(-3) = -3
stevep 0:04a9f72bbca7 180 Floor(2) = 2
stevep 0:04a9f72bbca7 181 Floor(2.3) = 2
stevep 0:04a9f72bbca7 182 Floor(2.8) = 2
stevep 0:04a9f72bbca7 183 */
stevep 0:04a9f72bbca7 184 template<class ValueType>
stevep 0:04a9f72bbca7 185 ValueType Floor(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 186 {
stevep 0:04a9f72bbca7 187 if( x.IsNan() )
stevep 0:04a9f72bbca7 188 {
stevep 0:04a9f72bbca7 189 if( err )
stevep 0:04a9f72bbca7 190 *err = err_improper_argument;
stevep 0:04a9f72bbca7 191
stevep 0:04a9f72bbca7 192 return x; // NaN
stevep 0:04a9f72bbca7 193 }
stevep 0:04a9f72bbca7 194
stevep 0:04a9f72bbca7 195 ValueType result(x);
stevep 0:04a9f72bbca7 196 uint c = 0;
stevep 0:04a9f72bbca7 197
stevep 0:04a9f72bbca7 198 result.SkipFraction();
stevep 0:04a9f72bbca7 199
stevep 0:04a9f72bbca7 200 if( result != x )
stevep 0:04a9f72bbca7 201 {
stevep 0:04a9f72bbca7 202 // x is with fraction
stevep 0:04a9f72bbca7 203 // if x is positive we don't have to do anything
stevep 0:04a9f72bbca7 204 if( x.IsSign() )
stevep 0:04a9f72bbca7 205 {
stevep 0:04a9f72bbca7 206 ValueType one;
stevep 0:04a9f72bbca7 207 one.SetOne();
stevep 0:04a9f72bbca7 208
stevep 0:04a9f72bbca7 209 c += result.Sub(one);
stevep 0:04a9f72bbca7 210 }
stevep 0:04a9f72bbca7 211 }
stevep 0:04a9f72bbca7 212
stevep 0:04a9f72bbca7 213 if( err )
stevep 0:04a9f72bbca7 214 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 215
stevep 0:04a9f72bbca7 216 return result;
stevep 0:04a9f72bbca7 217 }
stevep 0:04a9f72bbca7 218
stevep 0:04a9f72bbca7 219
stevep 0:04a9f72bbca7 220
stevep 0:04a9f72bbca7 221 /*
stevep 0:04a9f72bbca7 222 *
stevep 0:04a9f72bbca7 223 * logarithms and the exponent
stevep 0:04a9f72bbca7 224 *
stevep 0:04a9f72bbca7 225 *
stevep 0:04a9f72bbca7 226 */
stevep 0:04a9f72bbca7 227
stevep 0:04a9f72bbca7 228
stevep 0:04a9f72bbca7 229 /*!
stevep 0:04a9f72bbca7 230 this function calculates the natural logarithm (logarithm with the base 'e')
stevep 0:04a9f72bbca7 231 */
stevep 0:04a9f72bbca7 232 template<class ValueType>
stevep 0:04a9f72bbca7 233 ValueType Ln(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 234 {
stevep 0:04a9f72bbca7 235 if( x.IsNan() )
stevep 0:04a9f72bbca7 236 {
stevep 0:04a9f72bbca7 237 if( err )
stevep 0:04a9f72bbca7 238 *err = err_improper_argument;
stevep 0:04a9f72bbca7 239
stevep 0:04a9f72bbca7 240 return x; // NaN
stevep 0:04a9f72bbca7 241 }
stevep 0:04a9f72bbca7 242
stevep 0:04a9f72bbca7 243 ValueType result;
stevep 0:04a9f72bbca7 244 uint state = result.Ln(x);
stevep 0:04a9f72bbca7 245
stevep 0:04a9f72bbca7 246 if( err )
stevep 0:04a9f72bbca7 247 {
stevep 0:04a9f72bbca7 248 switch( state )
stevep 0:04a9f72bbca7 249 {
stevep 0:04a9f72bbca7 250 case 0:
stevep 0:04a9f72bbca7 251 *err = err_ok;
stevep 0:04a9f72bbca7 252 break;
stevep 0:04a9f72bbca7 253 case 1:
stevep 0:04a9f72bbca7 254 *err = err_overflow;
stevep 0:04a9f72bbca7 255 break;
stevep 0:04a9f72bbca7 256 case 2:
stevep 0:04a9f72bbca7 257 *err = err_improper_argument;
stevep 0:04a9f72bbca7 258 break;
stevep 0:04a9f72bbca7 259 default:
stevep 0:04a9f72bbca7 260 *err = err_internal_error;
stevep 0:04a9f72bbca7 261 break;
stevep 0:04a9f72bbca7 262 }
stevep 0:04a9f72bbca7 263 }
stevep 0:04a9f72bbca7 264
stevep 0:04a9f72bbca7 265
stevep 0:04a9f72bbca7 266 return result;
stevep 0:04a9f72bbca7 267 }
stevep 0:04a9f72bbca7 268
stevep 0:04a9f72bbca7 269
stevep 0:04a9f72bbca7 270 /*!
stevep 0:04a9f72bbca7 271 this function calculates the logarithm
stevep 0:04a9f72bbca7 272 */
stevep 0:04a9f72bbca7 273 template<class ValueType>
stevep 0:04a9f72bbca7 274 ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 275 {
stevep 0:04a9f72bbca7 276 if( x.IsNan() )
stevep 0:04a9f72bbca7 277 {
stevep 0:04a9f72bbca7 278 if( err ) *err = err_improper_argument;
stevep 0:04a9f72bbca7 279 return x;
stevep 0:04a9f72bbca7 280 }
stevep 0:04a9f72bbca7 281
stevep 0:04a9f72bbca7 282 if( base.IsNan() )
stevep 0:04a9f72bbca7 283 {
stevep 0:04a9f72bbca7 284 if( err ) *err = err_improper_argument;
stevep 0:04a9f72bbca7 285 return base;
stevep 0:04a9f72bbca7 286 }
stevep 0:04a9f72bbca7 287
stevep 0:04a9f72bbca7 288 ValueType result;
stevep 0:04a9f72bbca7 289 uint state = result.Log(x, base);
stevep 0:04a9f72bbca7 290
stevep 0:04a9f72bbca7 291 if( err )
stevep 0:04a9f72bbca7 292 {
stevep 0:04a9f72bbca7 293 switch( state )
stevep 0:04a9f72bbca7 294 {
stevep 0:04a9f72bbca7 295 case 0:
stevep 0:04a9f72bbca7 296 *err = err_ok;
stevep 0:04a9f72bbca7 297 break;
stevep 0:04a9f72bbca7 298 case 1:
stevep 0:04a9f72bbca7 299 *err = err_overflow;
stevep 0:04a9f72bbca7 300 break;
stevep 0:04a9f72bbca7 301 case 2:
stevep 0:04a9f72bbca7 302 case 3:
stevep 0:04a9f72bbca7 303 *err = err_improper_argument;
stevep 0:04a9f72bbca7 304 break;
stevep 0:04a9f72bbca7 305 default:
stevep 0:04a9f72bbca7 306 *err = err_internal_error;
stevep 0:04a9f72bbca7 307 break;
stevep 0:04a9f72bbca7 308 }
stevep 0:04a9f72bbca7 309 }
stevep 0:04a9f72bbca7 310
stevep 0:04a9f72bbca7 311 return result;
stevep 0:04a9f72bbca7 312 }
stevep 0:04a9f72bbca7 313
stevep 0:04a9f72bbca7 314
stevep 0:04a9f72bbca7 315 /*!
stevep 0:04a9f72bbca7 316 this function calculates the expression e^x
stevep 0:04a9f72bbca7 317 */
stevep 0:04a9f72bbca7 318 template<class ValueType>
stevep 0:04a9f72bbca7 319 ValueType Exp(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 320 {
stevep 0:04a9f72bbca7 321 if( x.IsNan() )
stevep 0:04a9f72bbca7 322 {
stevep 0:04a9f72bbca7 323 if( err )
stevep 0:04a9f72bbca7 324 *err = err_improper_argument;
stevep 0:04a9f72bbca7 325
stevep 0:04a9f72bbca7 326 return x; // NaN
stevep 0:04a9f72bbca7 327 }
stevep 0:04a9f72bbca7 328
stevep 0:04a9f72bbca7 329 ValueType result;
stevep 0:04a9f72bbca7 330 uint c = result.Exp(x);
stevep 0:04a9f72bbca7 331
stevep 0:04a9f72bbca7 332 if( err )
stevep 0:04a9f72bbca7 333 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 334
stevep 0:04a9f72bbca7 335 return result;
stevep 0:04a9f72bbca7 336 }
stevep 0:04a9f72bbca7 337
stevep 0:04a9f72bbca7 338
stevep 0:04a9f72bbca7 339 /*!
stevep 0:04a9f72bbca7 340 *
stevep 0:04a9f72bbca7 341 * trigonometric functions
stevep 0:04a9f72bbca7 342 *
stevep 0:04a9f72bbca7 343 */
stevep 0:04a9f72bbca7 344
stevep 0:04a9f72bbca7 345
stevep 0:04a9f72bbca7 346 /*
stevep 0:04a9f72bbca7 347 this namespace consists of auxiliary functions
stevep 0:04a9f72bbca7 348 (something like 'private' in a class)
stevep 0:04a9f72bbca7 349 */
stevep 0:04a9f72bbca7 350 namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 351 {
stevep 0:04a9f72bbca7 352
stevep 0:04a9f72bbca7 353 /*!
stevep 0:04a9f72bbca7 354 an auxiliary function for calculating the Sine
stevep 0:04a9f72bbca7 355 (you don't have to call this function)
stevep 0:04a9f72bbca7 356 */
stevep 0:04a9f72bbca7 357 template<class ValueType>
stevep 0:04a9f72bbca7 358 uint PrepareSin(ValueType & x, bool & change_sign)
stevep 0:04a9f72bbca7 359 {
stevep 0:04a9f72bbca7 360 ValueType temp;
stevep 0:04a9f72bbca7 361
stevep 0:04a9f72bbca7 362 change_sign = false;
stevep 0:04a9f72bbca7 363
stevep 0:04a9f72bbca7 364 if( x.IsSign() )
stevep 0:04a9f72bbca7 365 {
stevep 0:04a9f72bbca7 366 // we're using the formula 'sin(-x) = -sin(x)'
stevep 0:04a9f72bbca7 367 change_sign = !change_sign;
stevep 0:04a9f72bbca7 368 x.ChangeSign();
stevep 0:04a9f72bbca7 369 }
stevep 0:04a9f72bbca7 370
stevep 0:04a9f72bbca7 371 // we're reducing the period 2*PI
stevep 0:04a9f72bbca7 372 // (for big values there'll always be zero)
stevep 0:04a9f72bbca7 373 temp.Set2Pi();
stevep 0:04a9f72bbca7 374
stevep 0:04a9f72bbca7 375 if( x.Mod(temp) )
stevep 0:04a9f72bbca7 376 return 1;
stevep 0:04a9f72bbca7 377
stevep 0:04a9f72bbca7 378
stevep 0:04a9f72bbca7 379 // we're setting 'x' as being in the range of <0, 0.5PI>
stevep 0:04a9f72bbca7 380
stevep 0:04a9f72bbca7 381 temp.SetPi();
stevep 0:04a9f72bbca7 382
stevep 0:04a9f72bbca7 383 if( x > temp )
stevep 0:04a9f72bbca7 384 {
stevep 0:04a9f72bbca7 385 // x is in (pi, 2*pi>
stevep 0:04a9f72bbca7 386 x.Sub( temp );
stevep 0:04a9f72bbca7 387 change_sign = !change_sign;
stevep 0:04a9f72bbca7 388 }
stevep 0:04a9f72bbca7 389
stevep 0:04a9f72bbca7 390 temp.Set05Pi();
stevep 0:04a9f72bbca7 391
stevep 0:04a9f72bbca7 392 if( x > temp )
stevep 0:04a9f72bbca7 393 {
stevep 0:04a9f72bbca7 394 // x is in (0.5pi, pi>
stevep 0:04a9f72bbca7 395 x.Sub( temp );
stevep 0:04a9f72bbca7 396 x = temp - x;
stevep 0:04a9f72bbca7 397 }
stevep 0:04a9f72bbca7 398
stevep 0:04a9f72bbca7 399 return 0;
stevep 0:04a9f72bbca7 400 }
stevep 0:04a9f72bbca7 401
stevep 0:04a9f72bbca7 402
stevep 0:04a9f72bbca7 403 /*!
stevep 0:04a9f72bbca7 404 an auxiliary function for calculating the Sine
stevep 0:04a9f72bbca7 405 (you don't have to call this function)
stevep 0:04a9f72bbca7 406
stevep 0:04a9f72bbca7 407 it returns Sin(x) where 'x' is from <0, PI/2>
stevep 0:04a9f72bbca7 408 we're calculating the Sin with using Taylor series in zero or PI/2
stevep 0:04a9f72bbca7 409 (depending on which point of these two points is nearer to the 'x')
stevep 0:04a9f72bbca7 410
stevep 0:04a9f72bbca7 411 Taylor series:
stevep 0:04a9f72bbca7 412 sin(x) = sin(a) + cos(a)*(x-a)/(1!)
stevep 0:04a9f72bbca7 413 - sin(a)*((x-a)^2)/(2!) - cos(a)*((x-a)^3)/(3!)
stevep 0:04a9f72bbca7 414 + sin(a)*((x-a)^4)/(4!) + ...
stevep 0:04a9f72bbca7 415
stevep 0:04a9f72bbca7 416 when a=0 it'll be:
stevep 0:04a9f72bbca7 417 sin(x) = (x)/(1!) - (x^3)/(3!) + (x^5)/(5!) - (x^7)/(7!) + (x^9)/(9!) ...
stevep 0:04a9f72bbca7 418
stevep 0:04a9f72bbca7 419 and when a=PI/2:
stevep 0:04a9f72bbca7 420 sin(x) = 1 - ((x-PI/2)^2)/(2!) + ((x-PI/2)^4)/(4!) - ((x-PI/2)^6)/(6!) ...
stevep 0:04a9f72bbca7 421 */
stevep 0:04a9f72bbca7 422 template<class ValueType>
stevep 0:04a9f72bbca7 423 ValueType Sin0pi05(const ValueType & x)
stevep 0:04a9f72bbca7 424 {
stevep 0:04a9f72bbca7 425 ValueType result;
stevep 0:04a9f72bbca7 426 ValueType numerator, denominator;
stevep 0:04a9f72bbca7 427 ValueType d_numerator, d_denominator;
stevep 0:04a9f72bbca7 428 ValueType one, temp, old_result;
stevep 0:04a9f72bbca7 429
stevep 0:04a9f72bbca7 430 // temp = pi/4
stevep 0:04a9f72bbca7 431 temp.Set05Pi();
stevep 0:04a9f72bbca7 432 temp.exponent.SubOne();
stevep 0:04a9f72bbca7 433
stevep 0:04a9f72bbca7 434 one.SetOne();
stevep 0:04a9f72bbca7 435
stevep 0:04a9f72bbca7 436 if( x < temp )
stevep 0:04a9f72bbca7 437 {
stevep 0:04a9f72bbca7 438 // we're using the Taylor series with a=0
stevep 0:04a9f72bbca7 439 result = x;
stevep 0:04a9f72bbca7 440 numerator = x;
stevep 0:04a9f72bbca7 441 denominator = one;
stevep 0:04a9f72bbca7 442
stevep 0:04a9f72bbca7 443 // d_numerator = x^2
stevep 0:04a9f72bbca7 444 d_numerator = x;
stevep 0:04a9f72bbca7 445 d_numerator.Mul(x);
stevep 0:04a9f72bbca7 446
stevep 0:04a9f72bbca7 447 d_denominator = 2;
stevep 0:04a9f72bbca7 448 }
stevep 0:04a9f72bbca7 449 else
stevep 0:04a9f72bbca7 450 {
stevep 0:04a9f72bbca7 451 // we're using the Taylor series with a=PI/2
stevep 0:04a9f72bbca7 452 result = one;
stevep 0:04a9f72bbca7 453 numerator = one;
stevep 0:04a9f72bbca7 454 denominator = one;
stevep 0:04a9f72bbca7 455
stevep 0:04a9f72bbca7 456 // d_numerator = (x-pi/2)^2
stevep 0:04a9f72bbca7 457 ValueType pi05;
stevep 0:04a9f72bbca7 458 pi05.Set05Pi();
stevep 0:04a9f72bbca7 459
stevep 0:04a9f72bbca7 460 temp = x;
stevep 0:04a9f72bbca7 461 temp.Sub( pi05 );
stevep 0:04a9f72bbca7 462 d_numerator = temp;
stevep 0:04a9f72bbca7 463 d_numerator.Mul( temp );
stevep 0:04a9f72bbca7 464
stevep 0:04a9f72bbca7 465 d_denominator = one;
stevep 0:04a9f72bbca7 466 }
stevep 0:04a9f72bbca7 467
stevep 0:04a9f72bbca7 468 uint c = 0;
stevep 0:04a9f72bbca7 469 bool addition = false;
stevep 0:04a9f72bbca7 470
stevep 0:04a9f72bbca7 471 old_result = result;
stevep 0:04a9f72bbca7 472 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
stevep 0:04a9f72bbca7 473 {
stevep 0:04a9f72bbca7 474 // we're starting from a second part of the formula
stevep 0:04a9f72bbca7 475 c += numerator. Mul( d_numerator );
stevep 0:04a9f72bbca7 476 c += denominator. Mul( d_denominator );
stevep 0:04a9f72bbca7 477 c += d_denominator.Add( one );
stevep 0:04a9f72bbca7 478 c += denominator. Mul( d_denominator );
stevep 0:04a9f72bbca7 479 c += d_denominator.Add( one );
stevep 0:04a9f72bbca7 480 temp = numerator;
stevep 0:04a9f72bbca7 481 c += temp.Div(denominator);
stevep 0:04a9f72bbca7 482
stevep 0:04a9f72bbca7 483 if( c )
stevep 0:04a9f72bbca7 484 // Sin is from <-1,1> and cannot make an overflow
stevep 0:04a9f72bbca7 485 // but the carry can be from the Taylor series
stevep 0:04a9f72bbca7 486 // (then we only break our calculations)
stevep 0:04a9f72bbca7 487 break;
stevep 0:04a9f72bbca7 488
stevep 0:04a9f72bbca7 489 if( addition )
stevep 0:04a9f72bbca7 490 result.Add( temp );
stevep 0:04a9f72bbca7 491 else
stevep 0:04a9f72bbca7 492 result.Sub( temp );
stevep 0:04a9f72bbca7 493
stevep 0:04a9f72bbca7 494
stevep 0:04a9f72bbca7 495 addition = !addition;
stevep 0:04a9f72bbca7 496
stevep 0:04a9f72bbca7 497 // we're testing whether the result has changed after adding
stevep 0:04a9f72bbca7 498 // the next part of the Taylor formula, if not we end the loop
stevep 0:04a9f72bbca7 499 // (it means 'x' is zero or 'x' is PI/2 or this part of the formula
stevep 0:04a9f72bbca7 500 // is too small)
stevep 0:04a9f72bbca7 501 if( result == old_result )
stevep 0:04a9f72bbca7 502 break;
stevep 0:04a9f72bbca7 503
stevep 0:04a9f72bbca7 504 old_result = result;
stevep 0:04a9f72bbca7 505 }
stevep 0:04a9f72bbca7 506
stevep 0:04a9f72bbca7 507 return result;
stevep 0:04a9f72bbca7 508 }
stevep 0:04a9f72bbca7 509
stevep 0:04a9f72bbca7 510 } // namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 511
stevep 0:04a9f72bbca7 512
stevep 0:04a9f72bbca7 513
stevep 0:04a9f72bbca7 514 /*!
stevep 0:04a9f72bbca7 515 this function calculates the Sine
stevep 0:04a9f72bbca7 516 */
stevep 0:04a9f72bbca7 517 template<class ValueType>
stevep 0:04a9f72bbca7 518 ValueType Sin(ValueType x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 519 {
stevep 0:04a9f72bbca7 520 using namespace auxiliaryfunctions;
stevep 0:04a9f72bbca7 521
stevep 0:04a9f72bbca7 522 ValueType one, result;
stevep 0:04a9f72bbca7 523 bool change_sign;
stevep 0:04a9f72bbca7 524
stevep 0:04a9f72bbca7 525 if( x.IsNan() )
stevep 0:04a9f72bbca7 526 {
stevep 0:04a9f72bbca7 527 if( err )
stevep 0:04a9f72bbca7 528 *err = err_improper_argument;
stevep 0:04a9f72bbca7 529
stevep 0:04a9f72bbca7 530 return x;
stevep 0:04a9f72bbca7 531 }
stevep 0:04a9f72bbca7 532
stevep 0:04a9f72bbca7 533 if( err )
stevep 0:04a9f72bbca7 534 *err = err_ok;
stevep 0:04a9f72bbca7 535
stevep 0:04a9f72bbca7 536 if( PrepareSin( x, change_sign ) )
stevep 0:04a9f72bbca7 537 {
stevep 0:04a9f72bbca7 538 // x is too big, we cannnot reduce the 2*PI period
stevep 0:04a9f72bbca7 539 // prior to version 0.8.5 the result was zero
stevep 0:04a9f72bbca7 540
stevep 0:04a9f72bbca7 541 // result has NaN flag set by default
stevep 0:04a9f72bbca7 542
stevep 0:04a9f72bbca7 543 if( err )
stevep 0:04a9f72bbca7 544 *err = err_overflow; // maybe another error code? err_improper_argument?
stevep 0:04a9f72bbca7 545
stevep 0:04a9f72bbca7 546 return result; // NaN is set by default
stevep 0:04a9f72bbca7 547 }
stevep 0:04a9f72bbca7 548
stevep 0:04a9f72bbca7 549 result = Sin0pi05( x );
stevep 0:04a9f72bbca7 550
stevep 0:04a9f72bbca7 551 one.SetOne();
stevep 0:04a9f72bbca7 552
stevep 0:04a9f72bbca7 553 // after calculations there can be small distortions in the result
stevep 0:04a9f72bbca7 554 if( result > one )
stevep 0:04a9f72bbca7 555 result = one;
stevep 0:04a9f72bbca7 556 else
stevep 0:04a9f72bbca7 557 if( result.IsSign() )
stevep 0:04a9f72bbca7 558 // we've calculated the sin from <0, pi/2> and the result
stevep 0:04a9f72bbca7 559 // should be positive
stevep 0:04a9f72bbca7 560 result.SetZero();
stevep 0:04a9f72bbca7 561
stevep 0:04a9f72bbca7 562 if( change_sign )
stevep 0:04a9f72bbca7 563 result.ChangeSign();
stevep 0:04a9f72bbca7 564
stevep 0:04a9f72bbca7 565 return result;
stevep 0:04a9f72bbca7 566 }
stevep 0:04a9f72bbca7 567
stevep 0:04a9f72bbca7 568
stevep 0:04a9f72bbca7 569 /*!
stevep 0:04a9f72bbca7 570 this function calulates the Cosine
stevep 0:04a9f72bbca7 571 we're using the formula cos(x) = sin(x + PI/2)
stevep 0:04a9f72bbca7 572 */
stevep 0:04a9f72bbca7 573 template<class ValueType>
stevep 0:04a9f72bbca7 574 ValueType Cos(ValueType x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 575 {
stevep 0:04a9f72bbca7 576 if( x.IsNan() )
stevep 0:04a9f72bbca7 577 {
stevep 0:04a9f72bbca7 578 if( err )
stevep 0:04a9f72bbca7 579 *err = err_improper_argument;
stevep 0:04a9f72bbca7 580
stevep 0:04a9f72bbca7 581 return x; // NaN
stevep 0:04a9f72bbca7 582 }
stevep 0:04a9f72bbca7 583
stevep 0:04a9f72bbca7 584 ValueType pi05;
stevep 0:04a9f72bbca7 585 pi05.Set05Pi();
stevep 0:04a9f72bbca7 586
stevep 0:04a9f72bbca7 587 uint c = x.Add( pi05 );
stevep 0:04a9f72bbca7 588
stevep 0:04a9f72bbca7 589 if( c )
stevep 0:04a9f72bbca7 590 {
stevep 0:04a9f72bbca7 591 if( err )
stevep 0:04a9f72bbca7 592 *err = err_overflow;
stevep 0:04a9f72bbca7 593
stevep 0:04a9f72bbca7 594 return ValueType(); // result is undefined (NaN is set by default)
stevep 0:04a9f72bbca7 595 }
stevep 0:04a9f72bbca7 596
stevep 0:04a9f72bbca7 597 return Sin(x, err);
stevep 0:04a9f72bbca7 598 }
stevep 0:04a9f72bbca7 599
stevep 0:04a9f72bbca7 600
stevep 0:04a9f72bbca7 601 /*!
stevep 0:04a9f72bbca7 602 this function calulates the Tangent
stevep 0:04a9f72bbca7 603 we're using the formula tan(x) = sin(x) / cos(x)
stevep 0:04a9f72bbca7 604
stevep 0:04a9f72bbca7 605 it takes more time than calculating the Tan directly
stevep 0:04a9f72bbca7 606 from for example Taylor series but should be a bit preciser
stevep 0:04a9f72bbca7 607 because Tan receives its values from -infinity to +infinity
stevep 0:04a9f72bbca7 608 and when we calculate it from any series then we can make
stevep 0:04a9f72bbca7 609 a greater mistake than calculating 'sin/cos'
stevep 0:04a9f72bbca7 610 */
stevep 0:04a9f72bbca7 611 template<class ValueType>
stevep 0:04a9f72bbca7 612 ValueType Tan(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 613 {
stevep 0:04a9f72bbca7 614 ValueType result = Cos(x, err);
stevep 0:04a9f72bbca7 615
stevep 0:04a9f72bbca7 616 if( err && *err != err_ok )
stevep 0:04a9f72bbca7 617 return result;
stevep 0:04a9f72bbca7 618
stevep 0:04a9f72bbca7 619 if( result.IsZero() )
stevep 0:04a9f72bbca7 620 {
stevep 0:04a9f72bbca7 621 if( err )
stevep 0:04a9f72bbca7 622 *err = err_improper_argument;
stevep 0:04a9f72bbca7 623
stevep 0:04a9f72bbca7 624 result.SetNan();
stevep 0:04a9f72bbca7 625
stevep 0:04a9f72bbca7 626 return result;
stevep 0:04a9f72bbca7 627 }
stevep 0:04a9f72bbca7 628
stevep 0:04a9f72bbca7 629 return Sin(x, err) / result;
stevep 0:04a9f72bbca7 630 }
stevep 0:04a9f72bbca7 631
stevep 0:04a9f72bbca7 632
stevep 0:04a9f72bbca7 633 /*!
stevep 0:04a9f72bbca7 634 this function calulates the Tangent
stevep 0:04a9f72bbca7 635 look at the description of Tan(...)
stevep 0:04a9f72bbca7 636
stevep 0:04a9f72bbca7 637 (the abbreviation of Tangent can be 'tg' as well)
stevep 0:04a9f72bbca7 638 */
stevep 0:04a9f72bbca7 639 template<class ValueType>
stevep 0:04a9f72bbca7 640 ValueType Tg(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 641 {
stevep 0:04a9f72bbca7 642 return Tan(x, err);
stevep 0:04a9f72bbca7 643 }
stevep 0:04a9f72bbca7 644
stevep 0:04a9f72bbca7 645
stevep 0:04a9f72bbca7 646 /*!
stevep 0:04a9f72bbca7 647 this function calulates the Cotangent
stevep 0:04a9f72bbca7 648 we're using the formula tan(x) = cos(x) / sin(x)
stevep 0:04a9f72bbca7 649
stevep 0:04a9f72bbca7 650 (why do we make it in this way?
stevep 0:04a9f72bbca7 651 look at information in Tan() function)
stevep 0:04a9f72bbca7 652 */
stevep 0:04a9f72bbca7 653 template<class ValueType>
stevep 0:04a9f72bbca7 654 ValueType Cot(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 655 {
stevep 0:04a9f72bbca7 656 ValueType result = Sin(x, err);
stevep 0:04a9f72bbca7 657
stevep 0:04a9f72bbca7 658 if( err && *err != err_ok )
stevep 0:04a9f72bbca7 659 return result;
stevep 0:04a9f72bbca7 660
stevep 0:04a9f72bbca7 661 if( result.IsZero() )
stevep 0:04a9f72bbca7 662 {
stevep 0:04a9f72bbca7 663 if( err )
stevep 0:04a9f72bbca7 664 *err = err_improper_argument;
stevep 0:04a9f72bbca7 665
stevep 0:04a9f72bbca7 666 result.SetNan();
stevep 0:04a9f72bbca7 667
stevep 0:04a9f72bbca7 668 return result;
stevep 0:04a9f72bbca7 669 }
stevep 0:04a9f72bbca7 670
stevep 0:04a9f72bbca7 671 return Cos(x, err) / result;
stevep 0:04a9f72bbca7 672 }
stevep 0:04a9f72bbca7 673
stevep 0:04a9f72bbca7 674
stevep 0:04a9f72bbca7 675 /*!
stevep 0:04a9f72bbca7 676 this function calulates the Cotangent
stevep 0:04a9f72bbca7 677 look at the description of Cot(...)
stevep 0:04a9f72bbca7 678
stevep 0:04a9f72bbca7 679 (the abbreviation of Cotangent can be 'ctg' as well)
stevep 0:04a9f72bbca7 680 */
stevep 0:04a9f72bbca7 681 template<class ValueType>
stevep 0:04a9f72bbca7 682 ValueType Ctg(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 683 {
stevep 0:04a9f72bbca7 684 return Cot(x, err);
stevep 0:04a9f72bbca7 685 }
stevep 0:04a9f72bbca7 686
stevep 0:04a9f72bbca7 687
stevep 0:04a9f72bbca7 688 /*
stevep 0:04a9f72bbca7 689 *
stevep 0:04a9f72bbca7 690 * inverse trigonometric functions
stevep 0:04a9f72bbca7 691 *
stevep 0:04a9f72bbca7 692 *
stevep 0:04a9f72bbca7 693 */
stevep 0:04a9f72bbca7 694
stevep 0:04a9f72bbca7 695 namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 696 {
stevep 0:04a9f72bbca7 697
stevep 0:04a9f72bbca7 698 /*!
stevep 0:04a9f72bbca7 699 an auxiliary function for calculating the Arc Sine
stevep 0:04a9f72bbca7 700
stevep 0:04a9f72bbca7 701 we're calculating asin from the following formula:
stevep 0:04a9f72bbca7 702 asin(x) = x + (1*x^3)/(2*3) + (1*3*x^5)/(2*4*5) + (1*3*5*x^7)/(2*4*6*7) + ...
stevep 0:04a9f72bbca7 703 where abs(x) <= 1
stevep 0:04a9f72bbca7 704
stevep 0:04a9f72bbca7 705 we're using this formula when x is from <0, 1/2>
stevep 0:04a9f72bbca7 706 */
stevep 0:04a9f72bbca7 707 template<class ValueType>
stevep 0:04a9f72bbca7 708 ValueType ASin_0(const ValueType & x)
stevep 0:04a9f72bbca7 709 {
stevep 0:04a9f72bbca7 710 ValueType nominator, denominator, nominator_add, nominator_x, denominator_add, denominator_x;
stevep 0:04a9f72bbca7 711 ValueType two, result(x), x2(x);
stevep 0:04a9f72bbca7 712 ValueType nominator_temp, denominator_temp, old_result = result;
stevep 0:04a9f72bbca7 713 uint c = 0;
stevep 0:04a9f72bbca7 714
stevep 0:04a9f72bbca7 715 x2.Mul(x);
stevep 0:04a9f72bbca7 716 two = 2;
stevep 0:04a9f72bbca7 717
stevep 0:04a9f72bbca7 718 nominator.SetOne();
stevep 0:04a9f72bbca7 719 denominator = two;
stevep 0:04a9f72bbca7 720 nominator_add = nominator;
stevep 0:04a9f72bbca7 721 denominator_add = denominator;
stevep 0:04a9f72bbca7 722 nominator_x = x;
stevep 0:04a9f72bbca7 723 denominator_x = 3;
stevep 0:04a9f72bbca7 724
stevep 0:04a9f72bbca7 725 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
stevep 0:04a9f72bbca7 726 {
stevep 0:04a9f72bbca7 727 c += nominator_x.Mul(x2);
stevep 0:04a9f72bbca7 728 nominator_temp = nominator_x;
stevep 0:04a9f72bbca7 729 c += nominator_temp.Mul(nominator);
stevep 0:04a9f72bbca7 730 denominator_temp = denominator;
stevep 0:04a9f72bbca7 731 c += denominator_temp.Mul(denominator_x);
stevep 0:04a9f72bbca7 732 c += nominator_temp.Div(denominator_temp);
stevep 0:04a9f72bbca7 733
stevep 0:04a9f72bbca7 734 // if there is a carry somewhere we only break the calculating
stevep 0:04a9f72bbca7 735 // the result should be ok -- it's from <-pi/2, pi/2>
stevep 0:04a9f72bbca7 736 if( c )
stevep 0:04a9f72bbca7 737 break;
stevep 0:04a9f72bbca7 738
stevep 0:04a9f72bbca7 739 result.Add(nominator_temp);
stevep 0:04a9f72bbca7 740
stevep 0:04a9f72bbca7 741 if( result == old_result )
stevep 0:04a9f72bbca7 742 // there's no sense to calculate more
stevep 0:04a9f72bbca7 743 break;
stevep 0:04a9f72bbca7 744
stevep 0:04a9f72bbca7 745 old_result = result;
stevep 0:04a9f72bbca7 746
stevep 0:04a9f72bbca7 747
stevep 0:04a9f72bbca7 748 c += nominator_add.Add(two);
stevep 0:04a9f72bbca7 749 c += denominator_add.Add(two);
stevep 0:04a9f72bbca7 750 c += nominator.Mul(nominator_add);
stevep 0:04a9f72bbca7 751 c += denominator.Mul(denominator_add);
stevep 0:04a9f72bbca7 752 c += denominator_x.Add(two);
stevep 0:04a9f72bbca7 753 }
stevep 0:04a9f72bbca7 754
stevep 0:04a9f72bbca7 755 return result;
stevep 0:04a9f72bbca7 756 }
stevep 0:04a9f72bbca7 757
stevep 0:04a9f72bbca7 758
stevep 0:04a9f72bbca7 759
stevep 0:04a9f72bbca7 760 /*!
stevep 0:04a9f72bbca7 761 an auxiliary function for calculating the Arc Sine
stevep 0:04a9f72bbca7 762
stevep 0:04a9f72bbca7 763 we're calculating asin from the following formula:
stevep 0:04a9f72bbca7 764 asin(x) = pi/2 - sqrt(2)*sqrt(1-x) * asin_temp
stevep 0:04a9f72bbca7 765 asin_temp = 1 + (1*(1-x))/((2*3)*(2)) + (1*3*(1-x)^2)/((2*4*5)*(4)) + (1*3*5*(1-x)^3)/((2*4*6*7)*(8)) + ...
stevep 0:04a9f72bbca7 766
stevep 0:04a9f72bbca7 767 where abs(x) <= 1
stevep 0:04a9f72bbca7 768
stevep 0:04a9f72bbca7 769 we're using this formula when x is from (1/2, 1>
stevep 0:04a9f72bbca7 770 */
stevep 0:04a9f72bbca7 771 template<class ValueType>
stevep 0:04a9f72bbca7 772 ValueType ASin_1(const ValueType & x)
stevep 0:04a9f72bbca7 773 {
stevep 0:04a9f72bbca7 774 ValueType nominator, denominator, nominator_add, nominator_x, nominator_x_add, denominator_add, denominator_x;
stevep 0:04a9f72bbca7 775 ValueType denominator2;
stevep 0:04a9f72bbca7 776 ValueType one, two, result;
stevep 0:04a9f72bbca7 777 ValueType nominator_temp, denominator_temp, old_result;
stevep 0:04a9f72bbca7 778 uint c = 0;
stevep 0:04a9f72bbca7 779
stevep 0:04a9f72bbca7 780 two = 2;
stevep 0:04a9f72bbca7 781
stevep 0:04a9f72bbca7 782 one.SetOne();
stevep 0:04a9f72bbca7 783 nominator = one;
stevep 0:04a9f72bbca7 784 result = one;
stevep 0:04a9f72bbca7 785 old_result = result;
stevep 0:04a9f72bbca7 786 denominator = two;
stevep 0:04a9f72bbca7 787 nominator_add = nominator;
stevep 0:04a9f72bbca7 788 denominator_add = denominator;
stevep 0:04a9f72bbca7 789 nominator_x = one;
stevep 0:04a9f72bbca7 790 nominator_x.Sub(x);
stevep 0:04a9f72bbca7 791 nominator_x_add = nominator_x;
stevep 0:04a9f72bbca7 792 denominator_x = 3;
stevep 0:04a9f72bbca7 793 denominator2 = two;
stevep 0:04a9f72bbca7 794
stevep 0:04a9f72bbca7 795
stevep 0:04a9f72bbca7 796 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
stevep 0:04a9f72bbca7 797 {
stevep 0:04a9f72bbca7 798 nominator_temp = nominator_x;
stevep 0:04a9f72bbca7 799 c += nominator_temp.Mul(nominator);
stevep 0:04a9f72bbca7 800 denominator_temp = denominator;
stevep 0:04a9f72bbca7 801 c += denominator_temp.Mul(denominator_x);
stevep 0:04a9f72bbca7 802 c += denominator_temp.Mul(denominator2);
stevep 0:04a9f72bbca7 803 c += nominator_temp.Div(denominator_temp);
stevep 0:04a9f72bbca7 804
stevep 0:04a9f72bbca7 805 // if there is a carry somewhere we only break the calculating
stevep 0:04a9f72bbca7 806 // the result should be ok -- it's from <-pi/2, pi/2>
stevep 0:04a9f72bbca7 807 if( c )
stevep 0:04a9f72bbca7 808 break;
stevep 0:04a9f72bbca7 809
stevep 0:04a9f72bbca7 810 result.Add(nominator_temp);
stevep 0:04a9f72bbca7 811
stevep 0:04a9f72bbca7 812 if( result == old_result )
stevep 0:04a9f72bbca7 813 // there's no sense to calculate more
stevep 0:04a9f72bbca7 814 break;
stevep 0:04a9f72bbca7 815
stevep 0:04a9f72bbca7 816 old_result = result;
stevep 0:04a9f72bbca7 817
stevep 0:04a9f72bbca7 818 c += nominator_x.Mul(nominator_x_add);
stevep 0:04a9f72bbca7 819 c += nominator_add.Add(two);
stevep 0:04a9f72bbca7 820 c += denominator_add.Add(two);
stevep 0:04a9f72bbca7 821 c += nominator.Mul(nominator_add);
stevep 0:04a9f72bbca7 822 c += denominator.Mul(denominator_add);
stevep 0:04a9f72bbca7 823 c += denominator_x.Add(two);
stevep 0:04a9f72bbca7 824 c += denominator2.Mul(two);
stevep 0:04a9f72bbca7 825 }
stevep 0:04a9f72bbca7 826
stevep 0:04a9f72bbca7 827
stevep 0:04a9f72bbca7 828 nominator_x_add.exponent.AddOne(); // *2
stevep 0:04a9f72bbca7 829 one.exponent.SubOne(); // =0.5
stevep 0:04a9f72bbca7 830 nominator_x_add.Pow(one); // =sqrt(nominator_x_add)
stevep 0:04a9f72bbca7 831 result.Mul(nominator_x_add);
stevep 0:04a9f72bbca7 832
stevep 0:04a9f72bbca7 833 one.Set05Pi();
stevep 0:04a9f72bbca7 834 one.Sub(result);
stevep 0:04a9f72bbca7 835
stevep 0:04a9f72bbca7 836 return one;
stevep 0:04a9f72bbca7 837 }
stevep 0:04a9f72bbca7 838
stevep 0:04a9f72bbca7 839
stevep 0:04a9f72bbca7 840 } // namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 841
stevep 0:04a9f72bbca7 842
stevep 0:04a9f72bbca7 843 /*!
stevep 0:04a9f72bbca7 844 this function calculates the Arc Sine
stevep 0:04a9f72bbca7 845 x is from <-1,1>
stevep 0:04a9f72bbca7 846 */
stevep 0:04a9f72bbca7 847 template<class ValueType>
stevep 0:04a9f72bbca7 848 ValueType ASin(ValueType x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 849 {
stevep 0:04a9f72bbca7 850 using namespace auxiliaryfunctions;
stevep 0:04a9f72bbca7 851
stevep 0:04a9f72bbca7 852 ValueType result, one;
stevep 0:04a9f72bbca7 853 one.SetOne();
stevep 0:04a9f72bbca7 854 bool change_sign = false;
stevep 0:04a9f72bbca7 855
stevep 0:04a9f72bbca7 856 if( x.IsNan() )
stevep 0:04a9f72bbca7 857 {
stevep 0:04a9f72bbca7 858 if( err )
stevep 0:04a9f72bbca7 859 *err = err_improper_argument;
stevep 0:04a9f72bbca7 860
stevep 0:04a9f72bbca7 861 return x;
stevep 0:04a9f72bbca7 862 }
stevep 0:04a9f72bbca7 863
stevep 0:04a9f72bbca7 864 if( x.GreaterWithoutSignThan(one) )
stevep 0:04a9f72bbca7 865 {
stevep 0:04a9f72bbca7 866 if( err )
stevep 0:04a9f72bbca7 867 *err = err_improper_argument;
stevep 0:04a9f72bbca7 868
stevep 0:04a9f72bbca7 869 return result; // NaN is set by default
stevep 0:04a9f72bbca7 870 }
stevep 0:04a9f72bbca7 871
stevep 0:04a9f72bbca7 872 if( x.IsSign() )
stevep 0:04a9f72bbca7 873 {
stevep 0:04a9f72bbca7 874 change_sign = true;
stevep 0:04a9f72bbca7 875 x.Abs();
stevep 0:04a9f72bbca7 876 }
stevep 0:04a9f72bbca7 877
stevep 0:04a9f72bbca7 878 one.exponent.SubOne(); // =0.5
stevep 0:04a9f72bbca7 879
stevep 0:04a9f72bbca7 880 // asin(-x) = -asin(x)
stevep 0:04a9f72bbca7 881 if( x.GreaterWithoutSignThan(one) )
stevep 0:04a9f72bbca7 882 result = ASin_1(x);
stevep 0:04a9f72bbca7 883 else
stevep 0:04a9f72bbca7 884 result = ASin_0(x);
stevep 0:04a9f72bbca7 885
stevep 0:04a9f72bbca7 886 if( change_sign )
stevep 0:04a9f72bbca7 887 result.ChangeSign();
stevep 0:04a9f72bbca7 888
stevep 0:04a9f72bbca7 889 if( err )
stevep 0:04a9f72bbca7 890 *err = err_ok;
stevep 0:04a9f72bbca7 891
stevep 0:04a9f72bbca7 892 return result;
stevep 0:04a9f72bbca7 893 }
stevep 0:04a9f72bbca7 894
stevep 0:04a9f72bbca7 895
stevep 0:04a9f72bbca7 896 /*!
stevep 0:04a9f72bbca7 897 this function calculates the Arc Cosine
stevep 0:04a9f72bbca7 898
stevep 0:04a9f72bbca7 899 we're using the formula:
stevep 0:04a9f72bbca7 900 acos(x) = pi/2 - asin(x)
stevep 0:04a9f72bbca7 901 */
stevep 0:04a9f72bbca7 902 template<class ValueType>
stevep 0:04a9f72bbca7 903 ValueType ACos(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 904 {
stevep 0:04a9f72bbca7 905 ValueType temp;
stevep 0:04a9f72bbca7 906
stevep 0:04a9f72bbca7 907 temp.Set05Pi();
stevep 0:04a9f72bbca7 908 temp.Sub(ASin(x, err));
stevep 0:04a9f72bbca7 909
stevep 0:04a9f72bbca7 910 return temp;
stevep 0:04a9f72bbca7 911 }
stevep 0:04a9f72bbca7 912
stevep 0:04a9f72bbca7 913
stevep 0:04a9f72bbca7 914
stevep 0:04a9f72bbca7 915 namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 916 {
stevep 0:04a9f72bbca7 917
stevep 0:04a9f72bbca7 918 /*!
stevep 0:04a9f72bbca7 919 an auxiliary function for calculating the Arc Tangent
stevep 0:04a9f72bbca7 920
stevep 0:04a9f72bbca7 921 arc tan (x) where x is in <0; 0.5)
stevep 0:04a9f72bbca7 922 (x can be in (-0.5 ; 0.5) too)
stevep 0:04a9f72bbca7 923
stevep 0:04a9f72bbca7 924 we're using the Taylor series expanded in zero:
stevep 0:04a9f72bbca7 925 atan(x) = x - (x^3)/3 + (x^5)/5 - (x^7)/7 + ...
stevep 0:04a9f72bbca7 926 */
stevep 0:04a9f72bbca7 927 template<class ValueType>
stevep 0:04a9f72bbca7 928 ValueType ATan0(const ValueType & x)
stevep 0:04a9f72bbca7 929 {
stevep 0:04a9f72bbca7 930 ValueType nominator, denominator, nominator_add, denominator_add, temp;
stevep 0:04a9f72bbca7 931 ValueType result, old_result;
stevep 0:04a9f72bbca7 932 bool adding = false;
stevep 0:04a9f72bbca7 933 uint c = 0;
stevep 0:04a9f72bbca7 934
stevep 0:04a9f72bbca7 935 result = x;
stevep 0:04a9f72bbca7 936 old_result = result;
stevep 0:04a9f72bbca7 937 nominator = x;
stevep 0:04a9f72bbca7 938 nominator_add = x;
stevep 0:04a9f72bbca7 939 nominator_add.Mul(x);
stevep 0:04a9f72bbca7 940
stevep 0:04a9f72bbca7 941 denominator.SetOne();
stevep 0:04a9f72bbca7 942 denominator_add = 2;
stevep 0:04a9f72bbca7 943
stevep 0:04a9f72bbca7 944 for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
stevep 0:04a9f72bbca7 945 {
stevep 0:04a9f72bbca7 946 c += nominator.Mul(nominator_add);
stevep 0:04a9f72bbca7 947 c += denominator.Add(denominator_add);
stevep 0:04a9f72bbca7 948
stevep 0:04a9f72bbca7 949 temp = nominator;
stevep 0:04a9f72bbca7 950 c += temp.Div(denominator);
stevep 0:04a9f72bbca7 951
stevep 0:04a9f72bbca7 952 if( c )
stevep 0:04a9f72bbca7 953 // the result should be ok
stevep 0:04a9f72bbca7 954 break;
stevep 0:04a9f72bbca7 955
stevep 0:04a9f72bbca7 956 if( adding )
stevep 0:04a9f72bbca7 957 result.Add(temp);
stevep 0:04a9f72bbca7 958 else
stevep 0:04a9f72bbca7 959 result.Sub(temp);
stevep 0:04a9f72bbca7 960
stevep 0:04a9f72bbca7 961 if( result == old_result )
stevep 0:04a9f72bbca7 962 // there's no sense to calculate more
stevep 0:04a9f72bbca7 963 break;
stevep 0:04a9f72bbca7 964
stevep 0:04a9f72bbca7 965 old_result = result;
stevep 0:04a9f72bbca7 966 adding = !adding;
stevep 0:04a9f72bbca7 967 }
stevep 0:04a9f72bbca7 968
stevep 0:04a9f72bbca7 969 return result;
stevep 0:04a9f72bbca7 970 }
stevep 0:04a9f72bbca7 971
stevep 0:04a9f72bbca7 972
stevep 0:04a9f72bbca7 973 /*!
stevep 0:04a9f72bbca7 974 an auxiliary function for calculating the Arc Tangent
stevep 0:04a9f72bbca7 975
stevep 0:04a9f72bbca7 976 where x is in <0 ; 1>
stevep 0:04a9f72bbca7 977 */
stevep 0:04a9f72bbca7 978 template<class ValueType>
stevep 0:04a9f72bbca7 979 ValueType ATan01(const ValueType & x)
stevep 0:04a9f72bbca7 980 {
stevep 0:04a9f72bbca7 981 ValueType half;
stevep 0:04a9f72bbca7 982 half.Set05();
stevep 0:04a9f72bbca7 983
stevep 0:04a9f72bbca7 984 /*
stevep 0:04a9f72bbca7 985 it would be better if we chose about sqrt(2)-1=0.41... instead of 0.5 here
stevep 0:04a9f72bbca7 986
stevep 0:04a9f72bbca7 987 because as you can see below:
stevep 0:04a9f72bbca7 988 when x = sqrt(2)-1
stevep 0:04a9f72bbca7 989 abs(x) = abs( (x-1)/(1+x) )
stevep 0:04a9f72bbca7 990 so when we're calculating values around x
stevep 0:04a9f72bbca7 991 then they will be better converged to each other
stevep 0:04a9f72bbca7 992
stevep 0:04a9f72bbca7 993 for example if we have x=0.4999 then during calculating ATan0(0.4999)
stevep 0:04a9f72bbca7 994 we have to make about 141 iterations but when we have x=0.5
stevep 0:04a9f72bbca7 995 then during calculating ATan0( (x-1)/(1+x) ) we have to make
stevep 0:04a9f72bbca7 996 only about 89 iterations (both for Big<3,9>)
stevep 0:04a9f72bbca7 997
stevep 0:04a9f72bbca7 998 in the future this 0.5 can be changed
stevep 0:04a9f72bbca7 999 */
stevep 0:04a9f72bbca7 1000 if( x.SmallerWithoutSignThan(half) )
stevep 0:04a9f72bbca7 1001 return ATan0(x);
stevep 0:04a9f72bbca7 1002
stevep 0:04a9f72bbca7 1003
stevep 0:04a9f72bbca7 1004 /*
stevep 0:04a9f72bbca7 1005 x>=0.5 and x<=1
stevep 0:04a9f72bbca7 1006 (x can be even smaller than 0.5)
stevep 0:04a9f72bbca7 1007
stevep 0:04a9f72bbca7 1008 y = atac(x)
stevep 0:04a9f72bbca7 1009 x = tan(y)
stevep 0:04a9f72bbca7 1010
stevep 0:04a9f72bbca7 1011 tan(y-b) = (tan(y)-tab(b)) / (1+tan(y)*tan(b))
stevep 0:04a9f72bbca7 1012 y-b = atan( (tan(y)-tab(b)) / (1+tan(y)*tan(b)) )
stevep 0:04a9f72bbca7 1013 y = b + atan( (x-tab(b)) / (1+x*tan(b)) )
stevep 0:04a9f72bbca7 1014
stevep 0:04a9f72bbca7 1015 let b = pi/4
stevep 0:04a9f72bbca7 1016 tan(b) = tan(pi/4) = 1
stevep 0:04a9f72bbca7 1017 y = pi/4 + atan( (x-1)/(1+x) )
stevep 0:04a9f72bbca7 1018
stevep 0:04a9f72bbca7 1019 so
stevep 0:04a9f72bbca7 1020 atac(x) = pi/4 + atan( (x-1)/(1+x) )
stevep 0:04a9f72bbca7 1021 when x->1 (x converges to 1) the (x-1)/(1+x) -> 0
stevep 0:04a9f72bbca7 1022 and we can use ATan0() function here
stevep 0:04a9f72bbca7 1023 */
stevep 0:04a9f72bbca7 1024
stevep 0:04a9f72bbca7 1025 ValueType n(x),d(x),one,result;
stevep 0:04a9f72bbca7 1026
stevep 0:04a9f72bbca7 1027 one.SetOne();
stevep 0:04a9f72bbca7 1028 n.Sub(one);
stevep 0:04a9f72bbca7 1029 d.Add(one);
stevep 0:04a9f72bbca7 1030 n.Div(d);
stevep 0:04a9f72bbca7 1031
stevep 0:04a9f72bbca7 1032 result = ATan0(n);
stevep 0:04a9f72bbca7 1033
stevep 0:04a9f72bbca7 1034 n.Set05Pi();
stevep 0:04a9f72bbca7 1035 n.exponent.SubOne(); // =pi/4
stevep 0:04a9f72bbca7 1036 result.Add(n);
stevep 0:04a9f72bbca7 1037
stevep 0:04a9f72bbca7 1038 return result;
stevep 0:04a9f72bbca7 1039 }
stevep 0:04a9f72bbca7 1040
stevep 0:04a9f72bbca7 1041
stevep 0:04a9f72bbca7 1042 /*!
stevep 0:04a9f72bbca7 1043 an auxiliary function for calculating the Arc Tangent
stevep 0:04a9f72bbca7 1044 where x > 1
stevep 0:04a9f72bbca7 1045
stevep 0:04a9f72bbca7 1046 we're using the formula:
stevep 0:04a9f72bbca7 1047 atan(x) = pi/2 - atan(1/x) for x>0
stevep 0:04a9f72bbca7 1048 */
stevep 0:04a9f72bbca7 1049 template<class ValueType>
stevep 0:04a9f72bbca7 1050 ValueType ATanGreaterThanPlusOne(const ValueType & x)
stevep 0:04a9f72bbca7 1051 {
stevep 0:04a9f72bbca7 1052 ValueType temp, atan;
stevep 0:04a9f72bbca7 1053
stevep 0:04a9f72bbca7 1054 temp.SetOne();
stevep 0:04a9f72bbca7 1055
stevep 0:04a9f72bbca7 1056 if( temp.Div(x) )
stevep 0:04a9f72bbca7 1057 {
stevep 0:04a9f72bbca7 1058 // if there was a carry here that means x is very big
stevep 0:04a9f72bbca7 1059 // and atan(1/x) fast converged to 0
stevep 0:04a9f72bbca7 1060 atan.SetZero();
stevep 0:04a9f72bbca7 1061 }
stevep 0:04a9f72bbca7 1062 else
stevep 0:04a9f72bbca7 1063 atan = ATan01(temp);
stevep 0:04a9f72bbca7 1064
stevep 0:04a9f72bbca7 1065 temp.Set05Pi();
stevep 0:04a9f72bbca7 1066 temp.Sub(atan);
stevep 0:04a9f72bbca7 1067
stevep 0:04a9f72bbca7 1068 return temp;
stevep 0:04a9f72bbca7 1069 }
stevep 0:04a9f72bbca7 1070
stevep 0:04a9f72bbca7 1071 } // namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 1072
stevep 0:04a9f72bbca7 1073
stevep 0:04a9f72bbca7 1074 /*!
stevep 0:04a9f72bbca7 1075 this function calculates the Arc Tangent
stevep 0:04a9f72bbca7 1076 */
stevep 0:04a9f72bbca7 1077 template<class ValueType>
stevep 0:04a9f72bbca7 1078 ValueType ATan(ValueType x)
stevep 0:04a9f72bbca7 1079 {
stevep 0:04a9f72bbca7 1080 using namespace auxiliaryfunctions;
stevep 0:04a9f72bbca7 1081
stevep 0:04a9f72bbca7 1082 ValueType one, result;
stevep 0:04a9f72bbca7 1083 one.SetOne();
stevep 0:04a9f72bbca7 1084 bool change_sign = false;
stevep 0:04a9f72bbca7 1085
stevep 0:04a9f72bbca7 1086 if( x.IsNan() )
stevep 0:04a9f72bbca7 1087 return x;
stevep 0:04a9f72bbca7 1088
stevep 0:04a9f72bbca7 1089 // if x is negative we're using the formula:
stevep 0:04a9f72bbca7 1090 // atan(-x) = -atan(x)
stevep 0:04a9f72bbca7 1091 if( x.IsSign() )
stevep 0:04a9f72bbca7 1092 {
stevep 0:04a9f72bbca7 1093 change_sign = true;
stevep 0:04a9f72bbca7 1094 x.Abs();
stevep 0:04a9f72bbca7 1095 }
stevep 0:04a9f72bbca7 1096
stevep 0:04a9f72bbca7 1097 if( x.GreaterWithoutSignThan(one) )
stevep 0:04a9f72bbca7 1098 result = ATanGreaterThanPlusOne(x);
stevep 0:04a9f72bbca7 1099 else
stevep 0:04a9f72bbca7 1100 result = ATan01(x);
stevep 0:04a9f72bbca7 1101
stevep 0:04a9f72bbca7 1102 if( change_sign )
stevep 0:04a9f72bbca7 1103 result.ChangeSign();
stevep 0:04a9f72bbca7 1104
stevep 0:04a9f72bbca7 1105 return result;
stevep 0:04a9f72bbca7 1106 }
stevep 0:04a9f72bbca7 1107
stevep 0:04a9f72bbca7 1108
stevep 0:04a9f72bbca7 1109 /*!
stevep 0:04a9f72bbca7 1110 this function calculates the Arc Tangent
stevep 0:04a9f72bbca7 1111 look at the description of ATan(...)
stevep 0:04a9f72bbca7 1112
stevep 0:04a9f72bbca7 1113 (the abbreviation of Arc Tangent can be 'atg' as well)
stevep 0:04a9f72bbca7 1114 */
stevep 0:04a9f72bbca7 1115 template<class ValueType>
stevep 0:04a9f72bbca7 1116 ValueType ATg(const ValueType & x)
stevep 0:04a9f72bbca7 1117 {
stevep 0:04a9f72bbca7 1118 return ATan(x);
stevep 0:04a9f72bbca7 1119 }
stevep 0:04a9f72bbca7 1120
stevep 0:04a9f72bbca7 1121
stevep 0:04a9f72bbca7 1122 /*!
stevep 0:04a9f72bbca7 1123 this function calculates the Arc Cotangent
stevep 0:04a9f72bbca7 1124
stevep 0:04a9f72bbca7 1125 we're using the formula:
stevep 0:04a9f72bbca7 1126 actan(x) = pi/2 - atan(x)
stevep 0:04a9f72bbca7 1127 */
stevep 0:04a9f72bbca7 1128 template<class ValueType>
stevep 0:04a9f72bbca7 1129 ValueType ACot(const ValueType & x)
stevep 0:04a9f72bbca7 1130 {
stevep 0:04a9f72bbca7 1131 ValueType result;
stevep 0:04a9f72bbca7 1132
stevep 0:04a9f72bbca7 1133 result.Set05Pi();
stevep 0:04a9f72bbca7 1134 result.Sub(ATan(x));
stevep 0:04a9f72bbca7 1135
stevep 0:04a9f72bbca7 1136 return result;
stevep 0:04a9f72bbca7 1137 }
stevep 0:04a9f72bbca7 1138
stevep 0:04a9f72bbca7 1139
stevep 0:04a9f72bbca7 1140 /*!
stevep 0:04a9f72bbca7 1141 this function calculates the Arc Cotangent
stevep 0:04a9f72bbca7 1142 look at the description of ACot(...)
stevep 0:04a9f72bbca7 1143
stevep 0:04a9f72bbca7 1144 (the abbreviation of Arc Cotangent can be 'actg' as well)
stevep 0:04a9f72bbca7 1145 */
stevep 0:04a9f72bbca7 1146 template<class ValueType>
stevep 0:04a9f72bbca7 1147 ValueType ACtg(const ValueType & x)
stevep 0:04a9f72bbca7 1148 {
stevep 0:04a9f72bbca7 1149 return ACot(x);
stevep 0:04a9f72bbca7 1150 }
stevep 0:04a9f72bbca7 1151
stevep 0:04a9f72bbca7 1152
stevep 0:04a9f72bbca7 1153 /*
stevep 0:04a9f72bbca7 1154 *
stevep 0:04a9f72bbca7 1155 * hyperbolic functions
stevep 0:04a9f72bbca7 1156 *
stevep 0:04a9f72bbca7 1157 *
stevep 0:04a9f72bbca7 1158 */
stevep 0:04a9f72bbca7 1159
stevep 0:04a9f72bbca7 1160
stevep 0:04a9f72bbca7 1161 /*!
stevep 0:04a9f72bbca7 1162 this function calculates the Hyperbolic Sine
stevep 0:04a9f72bbca7 1163
stevep 0:04a9f72bbca7 1164 we're using the formula sinh(x)= ( e^x - e^(-x) ) / 2
stevep 0:04a9f72bbca7 1165 */
stevep 0:04a9f72bbca7 1166 template<class ValueType>
stevep 0:04a9f72bbca7 1167 ValueType Sinh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1168 {
stevep 0:04a9f72bbca7 1169 if( x.IsNan() )
stevep 0:04a9f72bbca7 1170 {
stevep 0:04a9f72bbca7 1171 if( err )
stevep 0:04a9f72bbca7 1172 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1173
stevep 0:04a9f72bbca7 1174 return x; // NaN
stevep 0:04a9f72bbca7 1175 }
stevep 0:04a9f72bbca7 1176
stevep 0:04a9f72bbca7 1177 ValueType ex, emx;
stevep 0:04a9f72bbca7 1178 uint c = 0;
stevep 0:04a9f72bbca7 1179
stevep 0:04a9f72bbca7 1180 c += ex.Exp(x);
stevep 0:04a9f72bbca7 1181 c += emx.Exp(-x);
stevep 0:04a9f72bbca7 1182
stevep 0:04a9f72bbca7 1183 c += ex.Sub(emx);
stevep 0:04a9f72bbca7 1184 c += ex.exponent.SubOne();
stevep 0:04a9f72bbca7 1185
stevep 0:04a9f72bbca7 1186 if( err )
stevep 0:04a9f72bbca7 1187 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1188
stevep 0:04a9f72bbca7 1189 return ex;
stevep 0:04a9f72bbca7 1190 }
stevep 0:04a9f72bbca7 1191
stevep 0:04a9f72bbca7 1192
stevep 0:04a9f72bbca7 1193 /*!
stevep 0:04a9f72bbca7 1194 this function calculates the Hyperbolic Cosine
stevep 0:04a9f72bbca7 1195
stevep 0:04a9f72bbca7 1196 we're using the formula cosh(x)= ( e^x + e^(-x) ) / 2
stevep 0:04a9f72bbca7 1197 */
stevep 0:04a9f72bbca7 1198 template<class ValueType>
stevep 0:04a9f72bbca7 1199 ValueType Cosh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1200 {
stevep 0:04a9f72bbca7 1201 if( x.IsNan() )
stevep 0:04a9f72bbca7 1202 {
stevep 0:04a9f72bbca7 1203 if( err )
stevep 0:04a9f72bbca7 1204 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1205
stevep 0:04a9f72bbca7 1206 return x; // NaN
stevep 0:04a9f72bbca7 1207 }
stevep 0:04a9f72bbca7 1208
stevep 0:04a9f72bbca7 1209 ValueType ex, emx;
stevep 0:04a9f72bbca7 1210 uint c = 0;
stevep 0:04a9f72bbca7 1211
stevep 0:04a9f72bbca7 1212 c += ex.Exp(x);
stevep 0:04a9f72bbca7 1213 c += emx.Exp(-x);
stevep 0:04a9f72bbca7 1214
stevep 0:04a9f72bbca7 1215 c += ex.Add(emx);
stevep 0:04a9f72bbca7 1216 c += ex.exponent.SubOne();
stevep 0:04a9f72bbca7 1217
stevep 0:04a9f72bbca7 1218 if( err )
stevep 0:04a9f72bbca7 1219 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1220
stevep 0:04a9f72bbca7 1221 return ex;
stevep 0:04a9f72bbca7 1222 }
stevep 0:04a9f72bbca7 1223
stevep 0:04a9f72bbca7 1224
stevep 0:04a9f72bbca7 1225 /*!
stevep 0:04a9f72bbca7 1226 this function calculates the Hyperbolic Tangent
stevep 0:04a9f72bbca7 1227
stevep 0:04a9f72bbca7 1228 we're using the formula tanh(x)= ( e^x - e^(-x) ) / ( e^x + e^(-x) )
stevep 0:04a9f72bbca7 1229 */
stevep 0:04a9f72bbca7 1230 template<class ValueType>
stevep 0:04a9f72bbca7 1231 ValueType Tanh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1232 {
stevep 0:04a9f72bbca7 1233 if( x.IsNan() )
stevep 0:04a9f72bbca7 1234 {
stevep 0:04a9f72bbca7 1235 if( err )
stevep 0:04a9f72bbca7 1236 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1237
stevep 0:04a9f72bbca7 1238 return x; // NaN
stevep 0:04a9f72bbca7 1239 }
stevep 0:04a9f72bbca7 1240
stevep 0:04a9f72bbca7 1241 ValueType ex, emx, nominator, denominator;
stevep 0:04a9f72bbca7 1242 uint c = 0;
stevep 0:04a9f72bbca7 1243
stevep 0:04a9f72bbca7 1244 c += ex.Exp(x);
stevep 0:04a9f72bbca7 1245 c += emx.Exp(-x);
stevep 0:04a9f72bbca7 1246
stevep 0:04a9f72bbca7 1247 nominator = ex;
stevep 0:04a9f72bbca7 1248 c += nominator.Sub(emx);
stevep 0:04a9f72bbca7 1249 denominator = ex;
stevep 0:04a9f72bbca7 1250 c += denominator.Add(emx);
stevep 0:04a9f72bbca7 1251
stevep 0:04a9f72bbca7 1252 c += nominator.Div(denominator);
stevep 0:04a9f72bbca7 1253
stevep 0:04a9f72bbca7 1254 if( err )
stevep 0:04a9f72bbca7 1255 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1256
stevep 0:04a9f72bbca7 1257 return nominator;
stevep 0:04a9f72bbca7 1258 }
stevep 0:04a9f72bbca7 1259
stevep 0:04a9f72bbca7 1260
stevep 0:04a9f72bbca7 1261 /*!
stevep 0:04a9f72bbca7 1262 this function calculates the Hyperbolic Tangent
stevep 0:04a9f72bbca7 1263 look at the description of Tanh(...)
stevep 0:04a9f72bbca7 1264
stevep 0:04a9f72bbca7 1265 (the abbreviation of Hyperbolic Tangent can be 'tgh' as well)
stevep 0:04a9f72bbca7 1266 */
stevep 0:04a9f72bbca7 1267 template<class ValueType>
stevep 0:04a9f72bbca7 1268 ValueType Tgh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1269 {
stevep 0:04a9f72bbca7 1270 return Tanh(x, err);
stevep 0:04a9f72bbca7 1271 }
stevep 0:04a9f72bbca7 1272
stevep 0:04a9f72bbca7 1273 /*!
stevep 0:04a9f72bbca7 1274 this function calculates the Hyperbolic Cotangent
stevep 0:04a9f72bbca7 1275
stevep 0:04a9f72bbca7 1276 we're using the formula coth(x)= ( e^x + e^(-x) ) / ( e^x - e^(-x) )
stevep 0:04a9f72bbca7 1277 */
stevep 0:04a9f72bbca7 1278 template<class ValueType>
stevep 0:04a9f72bbca7 1279 ValueType Coth(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1280 {
stevep 0:04a9f72bbca7 1281 if( x.IsNan() )
stevep 0:04a9f72bbca7 1282 {
stevep 0:04a9f72bbca7 1283 if( err )
stevep 0:04a9f72bbca7 1284 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1285
stevep 0:04a9f72bbca7 1286 return x; // NaN
stevep 0:04a9f72bbca7 1287 }
stevep 0:04a9f72bbca7 1288
stevep 0:04a9f72bbca7 1289 if( x.IsZero() )
stevep 0:04a9f72bbca7 1290 {
stevep 0:04a9f72bbca7 1291 if( err )
stevep 0:04a9f72bbca7 1292 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1293
stevep 0:04a9f72bbca7 1294 return ValueType(); // NaN is set by default
stevep 0:04a9f72bbca7 1295 }
stevep 0:04a9f72bbca7 1296
stevep 0:04a9f72bbca7 1297 ValueType ex, emx, nominator, denominator;
stevep 0:04a9f72bbca7 1298 uint c = 0;
stevep 0:04a9f72bbca7 1299
stevep 0:04a9f72bbca7 1300 c += ex.Exp(x);
stevep 0:04a9f72bbca7 1301 c += emx.Exp(-x);
stevep 0:04a9f72bbca7 1302
stevep 0:04a9f72bbca7 1303 nominator = ex;
stevep 0:04a9f72bbca7 1304 c += nominator.Add(emx);
stevep 0:04a9f72bbca7 1305 denominator = ex;
stevep 0:04a9f72bbca7 1306 c += denominator.Sub(emx);
stevep 0:04a9f72bbca7 1307
stevep 0:04a9f72bbca7 1308 c += nominator.Div(denominator);
stevep 0:04a9f72bbca7 1309
stevep 0:04a9f72bbca7 1310 if( err )
stevep 0:04a9f72bbca7 1311 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1312
stevep 0:04a9f72bbca7 1313 return nominator;
stevep 0:04a9f72bbca7 1314 }
stevep 0:04a9f72bbca7 1315
stevep 0:04a9f72bbca7 1316
stevep 0:04a9f72bbca7 1317 /*!
stevep 0:04a9f72bbca7 1318 this function calculates the Hyperbolic Cotangent
stevep 0:04a9f72bbca7 1319 look at the description of Coth(...)
stevep 0:04a9f72bbca7 1320
stevep 0:04a9f72bbca7 1321 (the abbreviation of Hyperbolic Cotangent can be 'ctgh' as well)
stevep 0:04a9f72bbca7 1322 */
stevep 0:04a9f72bbca7 1323 template<class ValueType>
stevep 0:04a9f72bbca7 1324 ValueType Ctgh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1325 {
stevep 0:04a9f72bbca7 1326 return Coth(x, err);
stevep 0:04a9f72bbca7 1327 }
stevep 0:04a9f72bbca7 1328
stevep 0:04a9f72bbca7 1329
stevep 0:04a9f72bbca7 1330 /*
stevep 0:04a9f72bbca7 1331 *
stevep 0:04a9f72bbca7 1332 * inverse hyperbolic functions
stevep 0:04a9f72bbca7 1333 *
stevep 0:04a9f72bbca7 1334 *
stevep 0:04a9f72bbca7 1335 */
stevep 0:04a9f72bbca7 1336
stevep 0:04a9f72bbca7 1337
stevep 0:04a9f72bbca7 1338 /*!
stevep 0:04a9f72bbca7 1339 inverse hyperbolic sine
stevep 0:04a9f72bbca7 1340
stevep 0:04a9f72bbca7 1341 asinh(x) = ln( x + sqrt(x^2 + 1) )
stevep 0:04a9f72bbca7 1342 */
stevep 0:04a9f72bbca7 1343 template<class ValueType>
stevep 0:04a9f72bbca7 1344 ValueType ASinh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1345 {
stevep 0:04a9f72bbca7 1346 if( x.IsNan() )
stevep 0:04a9f72bbca7 1347 {
stevep 0:04a9f72bbca7 1348 if( err )
stevep 0:04a9f72bbca7 1349 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1350
stevep 0:04a9f72bbca7 1351 return x; // NaN
stevep 0:04a9f72bbca7 1352 }
stevep 0:04a9f72bbca7 1353
stevep 0:04a9f72bbca7 1354 ValueType xx(x), one, result;
stevep 0:04a9f72bbca7 1355 uint c = 0;
stevep 0:04a9f72bbca7 1356 one.SetOne();
stevep 0:04a9f72bbca7 1357
stevep 0:04a9f72bbca7 1358 c += xx.Mul(x);
stevep 0:04a9f72bbca7 1359 c += xx.Add(one);
stevep 0:04a9f72bbca7 1360 one.exponent.SubOne(); // one=0.5
stevep 0:04a9f72bbca7 1361 // xx is >= 1
stevep 0:04a9f72bbca7 1362 c += xx.PowFrac(one); // xx=sqrt(xx)
stevep 0:04a9f72bbca7 1363 c += xx.Add(x);
stevep 0:04a9f72bbca7 1364 c += result.Ln(xx); // xx > 0
stevep 0:04a9f72bbca7 1365
stevep 0:04a9f72bbca7 1366 // here can only be a carry
stevep 0:04a9f72bbca7 1367 if( err )
stevep 0:04a9f72bbca7 1368 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1369
stevep 0:04a9f72bbca7 1370 return result;
stevep 0:04a9f72bbca7 1371 }
stevep 0:04a9f72bbca7 1372
stevep 0:04a9f72bbca7 1373
stevep 0:04a9f72bbca7 1374 /*!
stevep 0:04a9f72bbca7 1375 inverse hyperbolic cosine
stevep 0:04a9f72bbca7 1376
stevep 0:04a9f72bbca7 1377 acosh(x) = ln( x + sqrt(x^2 - 1) ) x in <1, infinity)
stevep 0:04a9f72bbca7 1378 */
stevep 0:04a9f72bbca7 1379 template<class ValueType>
stevep 0:04a9f72bbca7 1380 ValueType ACosh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1381 {
stevep 0:04a9f72bbca7 1382 if( x.IsNan() )
stevep 0:04a9f72bbca7 1383 {
stevep 0:04a9f72bbca7 1384 if( err )
stevep 0:04a9f72bbca7 1385 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1386
stevep 0:04a9f72bbca7 1387 return x; // NaN
stevep 0:04a9f72bbca7 1388 }
stevep 0:04a9f72bbca7 1389
stevep 0:04a9f72bbca7 1390 ValueType xx(x), one, result;
stevep 0:04a9f72bbca7 1391 uint c = 0;
stevep 0:04a9f72bbca7 1392 one.SetOne();
stevep 0:04a9f72bbca7 1393
stevep 0:04a9f72bbca7 1394 if( x < one )
stevep 0:04a9f72bbca7 1395 {
stevep 0:04a9f72bbca7 1396 if( err )
stevep 0:04a9f72bbca7 1397 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1398
stevep 0:04a9f72bbca7 1399 return result; // NaN is set by default
stevep 0:04a9f72bbca7 1400 }
stevep 0:04a9f72bbca7 1401
stevep 0:04a9f72bbca7 1402 c += xx.Mul(x);
stevep 0:04a9f72bbca7 1403 c += xx.Sub(one);
stevep 0:04a9f72bbca7 1404 // xx is >= 0
stevep 0:04a9f72bbca7 1405 // we can't call a PowFrac when the 'x' is zero
stevep 0:04a9f72bbca7 1406 // if x is 0 the sqrt(0) is 0
stevep 0:04a9f72bbca7 1407 if( !xx.IsZero() )
stevep 0:04a9f72bbca7 1408 {
stevep 0:04a9f72bbca7 1409 one.exponent.SubOne(); // one=0.5
stevep 0:04a9f72bbca7 1410 c += xx.PowFrac(one); // xx=sqrt(xx)
stevep 0:04a9f72bbca7 1411 }
stevep 0:04a9f72bbca7 1412 c += xx.Add(x);
stevep 0:04a9f72bbca7 1413 c += result.Ln(xx); // xx >= 1
stevep 0:04a9f72bbca7 1414
stevep 0:04a9f72bbca7 1415 // here can only be a carry
stevep 0:04a9f72bbca7 1416 if( err )
stevep 0:04a9f72bbca7 1417 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1418
stevep 0:04a9f72bbca7 1419 return result;
stevep 0:04a9f72bbca7 1420 }
stevep 0:04a9f72bbca7 1421
stevep 0:04a9f72bbca7 1422
stevep 0:04a9f72bbca7 1423 /*!
stevep 0:04a9f72bbca7 1424 inverse hyperbolic tangent
stevep 0:04a9f72bbca7 1425
stevep 0:04a9f72bbca7 1426 atanh(x) = 0.5 * ln( (1+x) / (1-x) ) x in (-1, 1)
stevep 0:04a9f72bbca7 1427 */
stevep 0:04a9f72bbca7 1428 template<class ValueType>
stevep 0:04a9f72bbca7 1429 ValueType ATanh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1430 {
stevep 0:04a9f72bbca7 1431 if( x.IsNan() )
stevep 0:04a9f72bbca7 1432 {
stevep 0:04a9f72bbca7 1433 if( err )
stevep 0:04a9f72bbca7 1434 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1435
stevep 0:04a9f72bbca7 1436 return x; // NaN
stevep 0:04a9f72bbca7 1437 }
stevep 0:04a9f72bbca7 1438
stevep 0:04a9f72bbca7 1439 ValueType nominator(x), denominator, one, result;
stevep 0:04a9f72bbca7 1440 uint c = 0;
stevep 0:04a9f72bbca7 1441 one.SetOne();
stevep 0:04a9f72bbca7 1442
stevep 0:04a9f72bbca7 1443 if( !x.SmallerWithoutSignThan(one) )
stevep 0:04a9f72bbca7 1444 {
stevep 0:04a9f72bbca7 1445 if( err )
stevep 0:04a9f72bbca7 1446 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1447
stevep 0:04a9f72bbca7 1448 return result; // NaN is set by default
stevep 0:04a9f72bbca7 1449 }
stevep 0:04a9f72bbca7 1450
stevep 0:04a9f72bbca7 1451 c += nominator.Add(one);
stevep 0:04a9f72bbca7 1452 denominator = one;
stevep 0:04a9f72bbca7 1453 c += denominator.Sub(x);
stevep 0:04a9f72bbca7 1454 c += nominator.Div(denominator);
stevep 0:04a9f72bbca7 1455 c += result.Ln(nominator);
stevep 0:04a9f72bbca7 1456 c += result.exponent.SubOne();
stevep 0:04a9f72bbca7 1457
stevep 0:04a9f72bbca7 1458 // here can only be a carry
stevep 0:04a9f72bbca7 1459 if( err )
stevep 0:04a9f72bbca7 1460 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1461
stevep 0:04a9f72bbca7 1462 return result;
stevep 0:04a9f72bbca7 1463 }
stevep 0:04a9f72bbca7 1464
stevep 0:04a9f72bbca7 1465
stevep 0:04a9f72bbca7 1466 /*!
stevep 0:04a9f72bbca7 1467 inverse hyperbolic tantent
stevep 0:04a9f72bbca7 1468 */
stevep 0:04a9f72bbca7 1469 template<class ValueType>
stevep 0:04a9f72bbca7 1470 ValueType ATgh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1471 {
stevep 0:04a9f72bbca7 1472 return ATanh(x, err);
stevep 0:04a9f72bbca7 1473 }
stevep 0:04a9f72bbca7 1474
stevep 0:04a9f72bbca7 1475
stevep 0:04a9f72bbca7 1476 /*!
stevep 0:04a9f72bbca7 1477 inverse hyperbolic cotangent
stevep 0:04a9f72bbca7 1478
stevep 0:04a9f72bbca7 1479 acoth(x) = 0.5 * ln( (x+1) / (x-1) ) x in (-infinity, -1) or (1, infinity)
stevep 0:04a9f72bbca7 1480 */
stevep 0:04a9f72bbca7 1481 template<class ValueType>
stevep 0:04a9f72bbca7 1482 ValueType ACoth(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1483 {
stevep 0:04a9f72bbca7 1484 if( x.IsNan() )
stevep 0:04a9f72bbca7 1485 {
stevep 0:04a9f72bbca7 1486 if( err )
stevep 0:04a9f72bbca7 1487 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1488
stevep 0:04a9f72bbca7 1489 return x; // NaN
stevep 0:04a9f72bbca7 1490 }
stevep 0:04a9f72bbca7 1491
stevep 0:04a9f72bbca7 1492 ValueType nominator(x), denominator(x), one, result;
stevep 0:04a9f72bbca7 1493 uint c = 0;
stevep 0:04a9f72bbca7 1494 one.SetOne();
stevep 0:04a9f72bbca7 1495
stevep 0:04a9f72bbca7 1496 if( !x.GreaterWithoutSignThan(one) )
stevep 0:04a9f72bbca7 1497 {
stevep 0:04a9f72bbca7 1498 if( err )
stevep 0:04a9f72bbca7 1499 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1500
stevep 0:04a9f72bbca7 1501 return result; // NaN is set by default
stevep 0:04a9f72bbca7 1502 }
stevep 0:04a9f72bbca7 1503
stevep 0:04a9f72bbca7 1504 c += nominator.Add(one);
stevep 0:04a9f72bbca7 1505 c += denominator.Sub(one);
stevep 0:04a9f72bbca7 1506 c += nominator.Div(denominator);
stevep 0:04a9f72bbca7 1507 c += result.Ln(nominator);
stevep 0:04a9f72bbca7 1508 c += result.exponent.SubOne();
stevep 0:04a9f72bbca7 1509
stevep 0:04a9f72bbca7 1510 // here can only be a carry
stevep 0:04a9f72bbca7 1511 if( err )
stevep 0:04a9f72bbca7 1512 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1513
stevep 0:04a9f72bbca7 1514 return result;
stevep 0:04a9f72bbca7 1515 }
stevep 0:04a9f72bbca7 1516
stevep 0:04a9f72bbca7 1517
stevep 0:04a9f72bbca7 1518 /*!
stevep 0:04a9f72bbca7 1519 inverse hyperbolic cotantent
stevep 0:04a9f72bbca7 1520 */
stevep 0:04a9f72bbca7 1521 template<class ValueType>
stevep 0:04a9f72bbca7 1522 ValueType ACtgh(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1523 {
stevep 0:04a9f72bbca7 1524 return ACoth(x, err);
stevep 0:04a9f72bbca7 1525 }
stevep 0:04a9f72bbca7 1526
stevep 0:04a9f72bbca7 1527
stevep 0:04a9f72bbca7 1528
stevep 0:04a9f72bbca7 1529
stevep 0:04a9f72bbca7 1530
stevep 0:04a9f72bbca7 1531 /*
stevep 0:04a9f72bbca7 1532 *
stevep 0:04a9f72bbca7 1533 * functions for converting between degrees, radians and gradians
stevep 0:04a9f72bbca7 1534 *
stevep 0:04a9f72bbca7 1535 *
stevep 0:04a9f72bbca7 1536 */
stevep 0:04a9f72bbca7 1537
stevep 0:04a9f72bbca7 1538
stevep 0:04a9f72bbca7 1539 /*!
stevep 0:04a9f72bbca7 1540 this function converts degrees to radians
stevep 0:04a9f72bbca7 1541
stevep 0:04a9f72bbca7 1542 it returns: x * pi / 180
stevep 0:04a9f72bbca7 1543 */
stevep 0:04a9f72bbca7 1544 template<class ValueType>
stevep 0:04a9f72bbca7 1545 ValueType DegToRad(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1546 {
stevep 0:04a9f72bbca7 1547 ValueType result, temp;
stevep 0:04a9f72bbca7 1548 uint c = 0;
stevep 0:04a9f72bbca7 1549
stevep 0:04a9f72bbca7 1550 if( x.IsNan() )
stevep 0:04a9f72bbca7 1551 {
stevep 0:04a9f72bbca7 1552 if( err )
stevep 0:04a9f72bbca7 1553 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1554
stevep 0:04a9f72bbca7 1555 return x;
stevep 0:04a9f72bbca7 1556 }
stevep 0:04a9f72bbca7 1557
stevep 0:04a9f72bbca7 1558 result = x;
stevep 0:04a9f72bbca7 1559
stevep 0:04a9f72bbca7 1560 // it is better to make division first and then multiplication
stevep 0:04a9f72bbca7 1561 // the result is more accurate especially when x is: 90,180,270 or 360
stevep 0:04a9f72bbca7 1562 temp = 180;
stevep 0:04a9f72bbca7 1563 c += result.Div(temp);
stevep 0:04a9f72bbca7 1564
stevep 0:04a9f72bbca7 1565 temp.SetPi();
stevep 0:04a9f72bbca7 1566 c += result.Mul(temp);
stevep 0:04a9f72bbca7 1567
stevep 0:04a9f72bbca7 1568 if( err )
stevep 0:04a9f72bbca7 1569 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1570
stevep 0:04a9f72bbca7 1571 return result;
stevep 0:04a9f72bbca7 1572 }
stevep 0:04a9f72bbca7 1573
stevep 0:04a9f72bbca7 1574
stevep 0:04a9f72bbca7 1575 /*!
stevep 0:04a9f72bbca7 1576 this function converts radians to degrees
stevep 0:04a9f72bbca7 1577
stevep 0:04a9f72bbca7 1578 it returns: x * 180 / pi
stevep 0:04a9f72bbca7 1579 */
stevep 0:04a9f72bbca7 1580 template<class ValueType>
stevep 0:04a9f72bbca7 1581 ValueType RadToDeg(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1582 {
stevep 0:04a9f72bbca7 1583 ValueType result, delimiter;
stevep 0:04a9f72bbca7 1584 uint c = 0;
stevep 0:04a9f72bbca7 1585
stevep 0:04a9f72bbca7 1586 if( x.IsNan() )
stevep 0:04a9f72bbca7 1587 {
stevep 0:04a9f72bbca7 1588 if( err )
stevep 0:04a9f72bbca7 1589 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1590
stevep 0:04a9f72bbca7 1591 return x;
stevep 0:04a9f72bbca7 1592 }
stevep 0:04a9f72bbca7 1593
stevep 0:04a9f72bbca7 1594 result = 180;
stevep 0:04a9f72bbca7 1595 c += result.Mul(x);
stevep 0:04a9f72bbca7 1596
stevep 0:04a9f72bbca7 1597 delimiter.SetPi();
stevep 0:04a9f72bbca7 1598 c += result.Div(delimiter);
stevep 0:04a9f72bbca7 1599
stevep 0:04a9f72bbca7 1600 if( err )
stevep 0:04a9f72bbca7 1601 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1602
stevep 0:04a9f72bbca7 1603 return result;
stevep 0:04a9f72bbca7 1604 }
stevep 0:04a9f72bbca7 1605
stevep 0:04a9f72bbca7 1606
stevep 0:04a9f72bbca7 1607 /*!
stevep 0:04a9f72bbca7 1608 this function converts degrees in the long format into one value
stevep 0:04a9f72bbca7 1609
stevep 0:04a9f72bbca7 1610 long format: (degrees, minutes, seconds)
stevep 0:04a9f72bbca7 1611 minutes and seconds must be greater than or equal zero
stevep 0:04a9f72bbca7 1612
stevep 0:04a9f72bbca7 1613 result:
stevep 0:04a9f72bbca7 1614 if d>=0 : result= d + ((s/60)+m)/60
stevep 0:04a9f72bbca7 1615 if d<0 : result= d - ((s/60)+m)/60
stevep 0:04a9f72bbca7 1616
stevep 0:04a9f72bbca7 1617 ((s/60)+m)/60 = (s+60*m)/3600 (second version is faster because
stevep 0:04a9f72bbca7 1618 there's only one division)
stevep 0:04a9f72bbca7 1619
stevep 0:04a9f72bbca7 1620 for example:
stevep 0:04a9f72bbca7 1621 DegToDeg(10, 30, 0) = 10.5
stevep 0:04a9f72bbca7 1622 DegToDeg(10, 24, 35.6)=10.4098(8)
stevep 0:04a9f72bbca7 1623 */
stevep 0:04a9f72bbca7 1624 template<class ValueType>
stevep 0:04a9f72bbca7 1625 ValueType DegToDeg( const ValueType & d, const ValueType & m, const ValueType & s,
stevep 0:04a9f72bbca7 1626 ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1627 {
stevep 0:04a9f72bbca7 1628 ValueType delimiter, multipler;
stevep 0:04a9f72bbca7 1629 uint c = 0;
stevep 0:04a9f72bbca7 1630
stevep 0:04a9f72bbca7 1631 if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() )
stevep 0:04a9f72bbca7 1632 {
stevep 0:04a9f72bbca7 1633 if( err )
stevep 0:04a9f72bbca7 1634 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1635
stevep 0:04a9f72bbca7 1636 delimiter.SetZeroNan(); // not needed, only to get rid of GCC warning about an uninitialized variable
stevep 0:04a9f72bbca7 1637
stevep 0:04a9f72bbca7 1638 return delimiter;
stevep 0:04a9f72bbca7 1639 }
stevep 0:04a9f72bbca7 1640
stevep 0:04a9f72bbca7 1641 multipler = 60;
stevep 0:04a9f72bbca7 1642 delimiter = 3600;
stevep 0:04a9f72bbca7 1643
stevep 0:04a9f72bbca7 1644 c += multipler.Mul(m);
stevep 0:04a9f72bbca7 1645 c += multipler.Add(s);
stevep 0:04a9f72bbca7 1646 c += multipler.Div(delimiter);
stevep 0:04a9f72bbca7 1647
stevep 0:04a9f72bbca7 1648 if( d.IsSign() )
stevep 0:04a9f72bbca7 1649 multipler.ChangeSign();
stevep 0:04a9f72bbca7 1650
stevep 0:04a9f72bbca7 1651 c += multipler.Add(d);
stevep 0:04a9f72bbca7 1652
stevep 0:04a9f72bbca7 1653 if( err )
stevep 0:04a9f72bbca7 1654 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1655
stevep 0:04a9f72bbca7 1656 return multipler;
stevep 0:04a9f72bbca7 1657 }
stevep 0:04a9f72bbca7 1658
stevep 0:04a9f72bbca7 1659
stevep 0:04a9f72bbca7 1660 /*!
stevep 0:04a9f72bbca7 1661 this function converts degrees in the long format to radians
stevep 0:04a9f72bbca7 1662 */
stevep 0:04a9f72bbca7 1663 template<class ValueType>
stevep 0:04a9f72bbca7 1664 ValueType DegToRad( const ValueType & d, const ValueType & m, const ValueType & s,
stevep 0:04a9f72bbca7 1665 ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1666 {
stevep 0:04a9f72bbca7 1667 ValueType temp_deg = DegToDeg(d,m,s,err);
stevep 0:04a9f72bbca7 1668
stevep 0:04a9f72bbca7 1669 if( err && *err!=err_ok )
stevep 0:04a9f72bbca7 1670 return temp_deg;
stevep 0:04a9f72bbca7 1671
stevep 0:04a9f72bbca7 1672 return DegToRad(temp_deg, err);
stevep 0:04a9f72bbca7 1673 }
stevep 0:04a9f72bbca7 1674
stevep 0:04a9f72bbca7 1675
stevep 0:04a9f72bbca7 1676 /*!
stevep 0:04a9f72bbca7 1677 this function converts gradians to radians
stevep 0:04a9f72bbca7 1678
stevep 0:04a9f72bbca7 1679 it returns: x * pi / 200
stevep 0:04a9f72bbca7 1680 */
stevep 0:04a9f72bbca7 1681 template<class ValueType>
stevep 0:04a9f72bbca7 1682 ValueType GradToRad(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1683 {
stevep 0:04a9f72bbca7 1684 ValueType result, temp;
stevep 0:04a9f72bbca7 1685 uint c = 0;
stevep 0:04a9f72bbca7 1686
stevep 0:04a9f72bbca7 1687 if( x.IsNan() )
stevep 0:04a9f72bbca7 1688 {
stevep 0:04a9f72bbca7 1689 if( err )
stevep 0:04a9f72bbca7 1690 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1691
stevep 0:04a9f72bbca7 1692 return x;
stevep 0:04a9f72bbca7 1693 }
stevep 0:04a9f72bbca7 1694
stevep 0:04a9f72bbca7 1695 result = x;
stevep 0:04a9f72bbca7 1696
stevep 0:04a9f72bbca7 1697 // it is better to make division first and then multiplication
stevep 0:04a9f72bbca7 1698 // the result is more accurate especially when x is: 100,200,300 or 400
stevep 0:04a9f72bbca7 1699 temp = 200;
stevep 0:04a9f72bbca7 1700 c += result.Div(temp);
stevep 0:04a9f72bbca7 1701
stevep 0:04a9f72bbca7 1702 temp.SetPi();
stevep 0:04a9f72bbca7 1703 c += result.Mul(temp);
stevep 0:04a9f72bbca7 1704
stevep 0:04a9f72bbca7 1705 if( err )
stevep 0:04a9f72bbca7 1706 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1707
stevep 0:04a9f72bbca7 1708 return result;
stevep 0:04a9f72bbca7 1709 }
stevep 0:04a9f72bbca7 1710
stevep 0:04a9f72bbca7 1711
stevep 0:04a9f72bbca7 1712 /*!
stevep 0:04a9f72bbca7 1713 this function converts radians to gradians
stevep 0:04a9f72bbca7 1714
stevep 0:04a9f72bbca7 1715 it returns: x * 200 / pi
stevep 0:04a9f72bbca7 1716 */
stevep 0:04a9f72bbca7 1717 template<class ValueType>
stevep 0:04a9f72bbca7 1718 ValueType RadToGrad(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1719 {
stevep 0:04a9f72bbca7 1720 ValueType result, delimiter;
stevep 0:04a9f72bbca7 1721 uint c = 0;
stevep 0:04a9f72bbca7 1722
stevep 0:04a9f72bbca7 1723 if( x.IsNan() )
stevep 0:04a9f72bbca7 1724 {
stevep 0:04a9f72bbca7 1725 if( err )
stevep 0:04a9f72bbca7 1726 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1727
stevep 0:04a9f72bbca7 1728 return x;
stevep 0:04a9f72bbca7 1729 }
stevep 0:04a9f72bbca7 1730
stevep 0:04a9f72bbca7 1731 result = 200;
stevep 0:04a9f72bbca7 1732 c += result.Mul(x);
stevep 0:04a9f72bbca7 1733
stevep 0:04a9f72bbca7 1734 delimiter.SetPi();
stevep 0:04a9f72bbca7 1735 c += result.Div(delimiter);
stevep 0:04a9f72bbca7 1736
stevep 0:04a9f72bbca7 1737 if( err )
stevep 0:04a9f72bbca7 1738 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1739
stevep 0:04a9f72bbca7 1740 return result;
stevep 0:04a9f72bbca7 1741 }
stevep 0:04a9f72bbca7 1742
stevep 0:04a9f72bbca7 1743
stevep 0:04a9f72bbca7 1744 /*!
stevep 0:04a9f72bbca7 1745 this function converts degrees to gradians
stevep 0:04a9f72bbca7 1746
stevep 0:04a9f72bbca7 1747 it returns: x * 200 / 180
stevep 0:04a9f72bbca7 1748 */
stevep 0:04a9f72bbca7 1749 template<class ValueType>
stevep 0:04a9f72bbca7 1750 ValueType DegToGrad(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1751 {
stevep 0:04a9f72bbca7 1752 ValueType result, temp;
stevep 0:04a9f72bbca7 1753 uint c = 0;
stevep 0:04a9f72bbca7 1754
stevep 0:04a9f72bbca7 1755 if( x.IsNan() )
stevep 0:04a9f72bbca7 1756 {
stevep 0:04a9f72bbca7 1757 if( err )
stevep 0:04a9f72bbca7 1758 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1759
stevep 0:04a9f72bbca7 1760 return x;
stevep 0:04a9f72bbca7 1761 }
stevep 0:04a9f72bbca7 1762
stevep 0:04a9f72bbca7 1763 result = x;
stevep 0:04a9f72bbca7 1764
stevep 0:04a9f72bbca7 1765 temp = 200;
stevep 0:04a9f72bbca7 1766 c += result.Mul(temp);
stevep 0:04a9f72bbca7 1767
stevep 0:04a9f72bbca7 1768 temp = 180;
stevep 0:04a9f72bbca7 1769 c += result.Div(temp);
stevep 0:04a9f72bbca7 1770
stevep 0:04a9f72bbca7 1771 if( err )
stevep 0:04a9f72bbca7 1772 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1773
stevep 0:04a9f72bbca7 1774 return result;
stevep 0:04a9f72bbca7 1775 }
stevep 0:04a9f72bbca7 1776
stevep 0:04a9f72bbca7 1777
stevep 0:04a9f72bbca7 1778 /*!
stevep 0:04a9f72bbca7 1779 this function converts degrees in the long format to gradians
stevep 0:04a9f72bbca7 1780 */
stevep 0:04a9f72bbca7 1781 template<class ValueType>
stevep 0:04a9f72bbca7 1782 ValueType DegToGrad( const ValueType & d, const ValueType & m, const ValueType & s,
stevep 0:04a9f72bbca7 1783 ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1784 {
stevep 0:04a9f72bbca7 1785 ValueType temp_deg = DegToDeg(d,m,s,err);
stevep 0:04a9f72bbca7 1786
stevep 0:04a9f72bbca7 1787 if( err && *err!=err_ok )
stevep 0:04a9f72bbca7 1788 return temp_deg;
stevep 0:04a9f72bbca7 1789
stevep 0:04a9f72bbca7 1790 return DegToGrad(temp_deg, err);
stevep 0:04a9f72bbca7 1791 }
stevep 0:04a9f72bbca7 1792
stevep 0:04a9f72bbca7 1793
stevep 0:04a9f72bbca7 1794 /*!
stevep 0:04a9f72bbca7 1795 this function converts degrees to gradians
stevep 0:04a9f72bbca7 1796
stevep 0:04a9f72bbca7 1797 it returns: x * 180 / 200
stevep 0:04a9f72bbca7 1798 */
stevep 0:04a9f72bbca7 1799 template<class ValueType>
stevep 0:04a9f72bbca7 1800 ValueType GradToDeg(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1801 {
stevep 0:04a9f72bbca7 1802 ValueType result, temp;
stevep 0:04a9f72bbca7 1803 uint c = 0;
stevep 0:04a9f72bbca7 1804
stevep 0:04a9f72bbca7 1805 if( x.IsNan() )
stevep 0:04a9f72bbca7 1806 {
stevep 0:04a9f72bbca7 1807 if( err )
stevep 0:04a9f72bbca7 1808 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1809
stevep 0:04a9f72bbca7 1810 return x;
stevep 0:04a9f72bbca7 1811 }
stevep 0:04a9f72bbca7 1812
stevep 0:04a9f72bbca7 1813 result = x;
stevep 0:04a9f72bbca7 1814
stevep 0:04a9f72bbca7 1815 temp = 180;
stevep 0:04a9f72bbca7 1816 c += result.Mul(temp);
stevep 0:04a9f72bbca7 1817
stevep 0:04a9f72bbca7 1818 temp = 200;
stevep 0:04a9f72bbca7 1819 c += result.Div(temp);
stevep 0:04a9f72bbca7 1820
stevep 0:04a9f72bbca7 1821 if( err )
stevep 0:04a9f72bbca7 1822 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1823
stevep 0:04a9f72bbca7 1824 return result;
stevep 0:04a9f72bbca7 1825 }
stevep 0:04a9f72bbca7 1826
stevep 0:04a9f72bbca7 1827
stevep 0:04a9f72bbca7 1828
stevep 0:04a9f72bbca7 1829
stevep 0:04a9f72bbca7 1830 /*
stevep 0:04a9f72bbca7 1831 *
stevep 0:04a9f72bbca7 1832 * another functions
stevep 0:04a9f72bbca7 1833 *
stevep 0:04a9f72bbca7 1834 *
stevep 0:04a9f72bbca7 1835 */
stevep 0:04a9f72bbca7 1836
stevep 0:04a9f72bbca7 1837
stevep 0:04a9f72bbca7 1838 /*!
stevep 0:04a9f72bbca7 1839 this function calculates the square root
stevep 0:04a9f72bbca7 1840
stevep 0:04a9f72bbca7 1841 Sqrt(9) = 3
stevep 0:04a9f72bbca7 1842 */
stevep 0:04a9f72bbca7 1843 template<class ValueType>
stevep 0:04a9f72bbca7 1844 ValueType Sqrt(ValueType x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 1845 {
stevep 0:04a9f72bbca7 1846 if( x.IsNan() || x.IsSign() )
stevep 0:04a9f72bbca7 1847 {
stevep 0:04a9f72bbca7 1848 if( err )
stevep 0:04a9f72bbca7 1849 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1850
stevep 0:04a9f72bbca7 1851 x.SetNan();
stevep 0:04a9f72bbca7 1852
stevep 0:04a9f72bbca7 1853 return x;
stevep 0:04a9f72bbca7 1854 }
stevep 0:04a9f72bbca7 1855
stevep 0:04a9f72bbca7 1856 uint c = x.Sqrt();
stevep 0:04a9f72bbca7 1857
stevep 0:04a9f72bbca7 1858 if( err )
stevep 0:04a9f72bbca7 1859 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 1860
stevep 0:04a9f72bbca7 1861 return x;
stevep 0:04a9f72bbca7 1862 }
stevep 0:04a9f72bbca7 1863
stevep 0:04a9f72bbca7 1864
stevep 0:04a9f72bbca7 1865
stevep 0:04a9f72bbca7 1866 namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 1867 {
stevep 0:04a9f72bbca7 1868
stevep 0:04a9f72bbca7 1869 template<class ValueType>
stevep 0:04a9f72bbca7 1870 bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err)
stevep 0:04a9f72bbca7 1871 {
stevep 0:04a9f72bbca7 1872 if( index.IsSign() )
stevep 0:04a9f72bbca7 1873 {
stevep 0:04a9f72bbca7 1874 // index cannot be negative
stevep 0:04a9f72bbca7 1875 if( err )
stevep 0:04a9f72bbca7 1876 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1877
stevep 0:04a9f72bbca7 1878 x.SetNan();
stevep 0:04a9f72bbca7 1879
stevep 0:04a9f72bbca7 1880 return true;
stevep 0:04a9f72bbca7 1881 }
stevep 0:04a9f72bbca7 1882
stevep 0:04a9f72bbca7 1883 return false;
stevep 0:04a9f72bbca7 1884 }
stevep 0:04a9f72bbca7 1885
stevep 0:04a9f72bbca7 1886
stevep 0:04a9f72bbca7 1887 template<class ValueType>
stevep 0:04a9f72bbca7 1888 bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err)
stevep 0:04a9f72bbca7 1889 {
stevep 0:04a9f72bbca7 1890 if( index.IsZero() )
stevep 0:04a9f72bbca7 1891 {
stevep 0:04a9f72bbca7 1892 if( x.IsZero() )
stevep 0:04a9f72bbca7 1893 {
stevep 0:04a9f72bbca7 1894 // there isn't root(0;0) - we assume it's not defined
stevep 0:04a9f72bbca7 1895 if( err )
stevep 0:04a9f72bbca7 1896 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1897
stevep 0:04a9f72bbca7 1898 x.SetNan();
stevep 0:04a9f72bbca7 1899
stevep 0:04a9f72bbca7 1900 return true;
stevep 0:04a9f72bbca7 1901 }
stevep 0:04a9f72bbca7 1902
stevep 0:04a9f72bbca7 1903 // root(x;0) is 1 (if x!=0)
stevep 0:04a9f72bbca7 1904 x.SetOne();
stevep 0:04a9f72bbca7 1905
stevep 0:04a9f72bbca7 1906 if( err )
stevep 0:04a9f72bbca7 1907 *err = err_ok;
stevep 0:04a9f72bbca7 1908
stevep 0:04a9f72bbca7 1909 return true;
stevep 0:04a9f72bbca7 1910 }
stevep 0:04a9f72bbca7 1911
stevep 0:04a9f72bbca7 1912 return false;
stevep 0:04a9f72bbca7 1913 }
stevep 0:04a9f72bbca7 1914
stevep 0:04a9f72bbca7 1915
stevep 0:04a9f72bbca7 1916 template<class ValueType>
stevep 0:04a9f72bbca7 1917 bool RootCheckIndexOne(const ValueType & index, ErrorCode * err)
stevep 0:04a9f72bbca7 1918 {
stevep 0:04a9f72bbca7 1919 ValueType one;
stevep 0:04a9f72bbca7 1920 one.SetOne();
stevep 0:04a9f72bbca7 1921
stevep 0:04a9f72bbca7 1922 if( index == one )
stevep 0:04a9f72bbca7 1923 {
stevep 0:04a9f72bbca7 1924 //root(x;1) is x
stevep 0:04a9f72bbca7 1925 // we do it because if we used the PowFrac function
stevep 0:04a9f72bbca7 1926 // we would lose the precision
stevep 0:04a9f72bbca7 1927 if( err )
stevep 0:04a9f72bbca7 1928 *err = err_ok;
stevep 0:04a9f72bbca7 1929
stevep 0:04a9f72bbca7 1930 return true;
stevep 0:04a9f72bbca7 1931 }
stevep 0:04a9f72bbca7 1932
stevep 0:04a9f72bbca7 1933 return false;
stevep 0:04a9f72bbca7 1934 }
stevep 0:04a9f72bbca7 1935
stevep 0:04a9f72bbca7 1936
stevep 0:04a9f72bbca7 1937 template<class ValueType>
stevep 0:04a9f72bbca7 1938 bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err)
stevep 0:04a9f72bbca7 1939 {
stevep 0:04a9f72bbca7 1940 if( index == 2 )
stevep 0:04a9f72bbca7 1941 {
stevep 0:04a9f72bbca7 1942 x = Sqrt(x, err);
stevep 0:04a9f72bbca7 1943
stevep 0:04a9f72bbca7 1944 return true;
stevep 0:04a9f72bbca7 1945 }
stevep 0:04a9f72bbca7 1946
stevep 0:04a9f72bbca7 1947 return false;
stevep 0:04a9f72bbca7 1948 }
stevep 0:04a9f72bbca7 1949
stevep 0:04a9f72bbca7 1950
stevep 0:04a9f72bbca7 1951 template<class ValueType>
stevep 0:04a9f72bbca7 1952 bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
stevep 0:04a9f72bbca7 1953 {
stevep 0:04a9f72bbca7 1954 if( !index.IsInteger() )
stevep 0:04a9f72bbca7 1955 {
stevep 0:04a9f72bbca7 1956 // index must be integer
stevep 0:04a9f72bbca7 1957 if( err )
stevep 0:04a9f72bbca7 1958 *err = err_improper_argument;
stevep 0:04a9f72bbca7 1959
stevep 0:04a9f72bbca7 1960 x.SetNan();
stevep 0:04a9f72bbca7 1961
stevep 0:04a9f72bbca7 1962 return true;
stevep 0:04a9f72bbca7 1963 }
stevep 0:04a9f72bbca7 1964
stevep 0:04a9f72bbca7 1965 return false;
stevep 0:04a9f72bbca7 1966 }
stevep 0:04a9f72bbca7 1967
stevep 0:04a9f72bbca7 1968
stevep 0:04a9f72bbca7 1969 template<class ValueType>
stevep 0:04a9f72bbca7 1970 bool RootCheckXZero(ValueType & x, ErrorCode * err)
stevep 0:04a9f72bbca7 1971 {
stevep 0:04a9f72bbca7 1972 if( x.IsZero() )
stevep 0:04a9f72bbca7 1973 {
stevep 0:04a9f72bbca7 1974 // root(0;index) is zero (if index!=0)
stevep 0:04a9f72bbca7 1975 // RootCheckIndexZero() must be called beforehand
stevep 0:04a9f72bbca7 1976 x.SetZero();
stevep 0:04a9f72bbca7 1977
stevep 0:04a9f72bbca7 1978 if( err )
stevep 0:04a9f72bbca7 1979 *err = err_ok;
stevep 0:04a9f72bbca7 1980
stevep 0:04a9f72bbca7 1981 return true;
stevep 0:04a9f72bbca7 1982 }
stevep 0:04a9f72bbca7 1983
stevep 0:04a9f72bbca7 1984 return false;
stevep 0:04a9f72bbca7 1985 }
stevep 0:04a9f72bbca7 1986
stevep 0:04a9f72bbca7 1987
stevep 0:04a9f72bbca7 1988 template<class ValueType>
stevep 0:04a9f72bbca7 1989 bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign)
stevep 0:04a9f72bbca7 1990 {
stevep 0:04a9f72bbca7 1991 *change_sign = false;
stevep 0:04a9f72bbca7 1992
stevep 0:04a9f72bbca7 1993 if( index.Mod2() )
stevep 0:04a9f72bbca7 1994 {
stevep 0:04a9f72bbca7 1995 // index is odd (1,3,5...)
stevep 0:04a9f72bbca7 1996 if( x.IsSign() )
stevep 0:04a9f72bbca7 1997 {
stevep 0:04a9f72bbca7 1998 *change_sign = true;
stevep 0:04a9f72bbca7 1999 x.Abs();
stevep 0:04a9f72bbca7 2000 }
stevep 0:04a9f72bbca7 2001 }
stevep 0:04a9f72bbca7 2002 else
stevep 0:04a9f72bbca7 2003 {
stevep 0:04a9f72bbca7 2004 // index is even
stevep 0:04a9f72bbca7 2005 // x cannot be negative
stevep 0:04a9f72bbca7 2006 if( x.IsSign() )
stevep 0:04a9f72bbca7 2007 {
stevep 0:04a9f72bbca7 2008 if( err )
stevep 0:04a9f72bbca7 2009 *err = err_improper_argument;
stevep 0:04a9f72bbca7 2010
stevep 0:04a9f72bbca7 2011 x.SetNan();
stevep 0:04a9f72bbca7 2012
stevep 0:04a9f72bbca7 2013 return true;
stevep 0:04a9f72bbca7 2014 }
stevep 0:04a9f72bbca7 2015 }
stevep 0:04a9f72bbca7 2016
stevep 0:04a9f72bbca7 2017 return false;
stevep 0:04a9f72bbca7 2018 }
stevep 0:04a9f72bbca7 2019
stevep 0:04a9f72bbca7 2020
stevep 0:04a9f72bbca7 2021 template<class ValueType>
stevep 0:04a9f72bbca7 2022 uint RootCorrectInteger(ValueType & old_x, ValueType & x, const ValueType & index)
stevep 0:04a9f72bbca7 2023 {
stevep 0:04a9f72bbca7 2024 if( !old_x.IsInteger() || x.IsInteger() || !index.exponent.IsSign() )
stevep 0:04a9f72bbca7 2025 return 0;
stevep 0:04a9f72bbca7 2026
stevep 0:04a9f72bbca7 2027 // old_x is integer,
stevep 0:04a9f72bbca7 2028 // x is not integer,
stevep 0:04a9f72bbca7 2029 // index is relatively small (index.exponent<0 or index.exponent<=0)
stevep 0:04a9f72bbca7 2030 // (because we're using a special powering algorithm Big::PowUInt())
stevep 0:04a9f72bbca7 2031
stevep 0:04a9f72bbca7 2032 uint c = 0;
stevep 0:04a9f72bbca7 2033
stevep 0:04a9f72bbca7 2034 ValueType temp(x);
stevep 0:04a9f72bbca7 2035 c += temp.Round();
stevep 0:04a9f72bbca7 2036
stevep 0:04a9f72bbca7 2037 ValueType temp_round(temp);
stevep 0:04a9f72bbca7 2038 c += temp.PowUInt(index);
stevep 0:04a9f72bbca7 2039
stevep 0:04a9f72bbca7 2040 if( temp == old_x )
stevep 0:04a9f72bbca7 2041 x = temp_round;
stevep 0:04a9f72bbca7 2042
stevep 0:04a9f72bbca7 2043 return (c==0)? 0 : 1;
stevep 0:04a9f72bbca7 2044 }
stevep 0:04a9f72bbca7 2045
stevep 0:04a9f72bbca7 2046
stevep 0:04a9f72bbca7 2047
stevep 0:04a9f72bbca7 2048 } // namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 2049
stevep 0:04a9f72bbca7 2050
stevep 0:04a9f72bbca7 2051
stevep 0:04a9f72bbca7 2052 /*!
stevep 0:04a9f72bbca7 2053 indexth Root of x
stevep 0:04a9f72bbca7 2054 index must be integer and not negative <0;1;2;3....)
stevep 0:04a9f72bbca7 2055
stevep 0:04a9f72bbca7 2056 if index==0 the result is one
stevep 0:04a9f72bbca7 2057 if x==0 the result is zero and we assume root(0;0) is not defined
stevep 0:04a9f72bbca7 2058
stevep 0:04a9f72bbca7 2059 if index is even (2;4;6...) the result is x^(1/index) and x>0
stevep 0:04a9f72bbca7 2060 if index is odd (1;2;3;...) the result is either
stevep 0:04a9f72bbca7 2061 -(abs(x)^(1/index)) if x<0 or
stevep 0:04a9f72bbca7 2062 x^(1/index)) if x>0
stevep 0:04a9f72bbca7 2063
stevep 0:04a9f72bbca7 2064 (for index==1 the result is equal x)
stevep 0:04a9f72bbca7 2065 */
stevep 0:04a9f72bbca7 2066 template<class ValueType>
stevep 0:04a9f72bbca7 2067 ValueType Root(ValueType x, const ValueType & index, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 2068 {
stevep 0:04a9f72bbca7 2069 using namespace auxiliaryfunctions;
stevep 0:04a9f72bbca7 2070
stevep 0:04a9f72bbca7 2071 if( x.IsNan() || index.IsNan() )
stevep 0:04a9f72bbca7 2072 {
stevep 0:04a9f72bbca7 2073 if( err )
stevep 0:04a9f72bbca7 2074 *err = err_improper_argument;
stevep 0:04a9f72bbca7 2075
stevep 0:04a9f72bbca7 2076 x.SetNan();
stevep 0:04a9f72bbca7 2077
stevep 0:04a9f72bbca7 2078 return x;
stevep 0:04a9f72bbca7 2079 }
stevep 0:04a9f72bbca7 2080
stevep 0:04a9f72bbca7 2081 if( RootCheckIndexSign(x, index, err) ) return x;
stevep 0:04a9f72bbca7 2082 if( RootCheckIndexZero(x, index, err) ) return x;
stevep 0:04a9f72bbca7 2083 if( RootCheckIndexOne ( index, err) ) return x;
stevep 0:04a9f72bbca7 2084 if( RootCheckIndexTwo (x, index, err) ) return x;
stevep 0:04a9f72bbca7 2085 if( RootCheckIndexFrac(x, index, err) ) return x;
stevep 0:04a9f72bbca7 2086 if( RootCheckXZero (x, err) ) return x;
stevep 0:04a9f72bbca7 2087
stevep 0:04a9f72bbca7 2088 // index integer and index!=0
stevep 0:04a9f72bbca7 2089 // x!=0
stevep 0:04a9f72bbca7 2090
stevep 0:04a9f72bbca7 2091 ValueType old_x(x);
stevep 0:04a9f72bbca7 2092 bool change_sign;
stevep 0:04a9f72bbca7 2093
stevep 0:04a9f72bbca7 2094 if( RootCheckIndex(x, index, err, &change_sign ) ) return x;
stevep 0:04a9f72bbca7 2095
stevep 0:04a9f72bbca7 2096 ValueType temp;
stevep 0:04a9f72bbca7 2097 uint c = 0;
stevep 0:04a9f72bbca7 2098
stevep 0:04a9f72bbca7 2099 // we're using the formula: root(x ; n) = exp( ln(x) / n )
stevep 0:04a9f72bbca7 2100 c += temp.Ln(x);
stevep 0:04a9f72bbca7 2101 c += temp.Div(index);
stevep 0:04a9f72bbca7 2102 c += x.Exp(temp);
stevep 0:04a9f72bbca7 2103
stevep 0:04a9f72bbca7 2104 if( change_sign )
stevep 0:04a9f72bbca7 2105 {
stevep 0:04a9f72bbca7 2106 // x is different from zero
stevep 0:04a9f72bbca7 2107 x.SetSign();
stevep 0:04a9f72bbca7 2108 }
stevep 0:04a9f72bbca7 2109
stevep 0:04a9f72bbca7 2110 c += RootCorrectInteger(old_x, x, index);
stevep 0:04a9f72bbca7 2111
stevep 0:04a9f72bbca7 2112 if( err )
stevep 0:04a9f72bbca7 2113 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 2114
stevep 0:04a9f72bbca7 2115 return x;
stevep 0:04a9f72bbca7 2116 }
stevep 0:04a9f72bbca7 2117
stevep 0:04a9f72bbca7 2118
stevep 0:04a9f72bbca7 2119
stevep 0:04a9f72bbca7 2120 /*!
stevep 0:04a9f72bbca7 2121 absolute value of x
stevep 0:04a9f72bbca7 2122 e.g. -2 = 2
stevep 0:04a9f72bbca7 2123 2 = 2
stevep 0:04a9f72bbca7 2124 */
stevep 0:04a9f72bbca7 2125 template<class ValueType>
stevep 0:04a9f72bbca7 2126 ValueType Abs(const ValueType & x)
stevep 0:04a9f72bbca7 2127 {
stevep 0:04a9f72bbca7 2128 ValueType result( x );
stevep 0:04a9f72bbca7 2129 result.Abs();
stevep 0:04a9f72bbca7 2130
stevep 0:04a9f72bbca7 2131 return result;
stevep 0:04a9f72bbca7 2132 }
stevep 0:04a9f72bbca7 2133
stevep 0:04a9f72bbca7 2134
stevep 0:04a9f72bbca7 2135 /*!
stevep 0:04a9f72bbca7 2136 it returns the sign of the value
stevep 0:04a9f72bbca7 2137 e.g. -2 = -1
stevep 0:04a9f72bbca7 2138 0 = 0
stevep 0:04a9f72bbca7 2139 10 = 1
stevep 0:04a9f72bbca7 2140 */
stevep 0:04a9f72bbca7 2141 template<class ValueType>
stevep 0:04a9f72bbca7 2142 ValueType Sgn(ValueType x)
stevep 0:04a9f72bbca7 2143 {
stevep 0:04a9f72bbca7 2144 x.Sgn();
stevep 0:04a9f72bbca7 2145
stevep 0:04a9f72bbca7 2146 return x;
stevep 0:04a9f72bbca7 2147 }
stevep 0:04a9f72bbca7 2148
stevep 0:04a9f72bbca7 2149
stevep 0:04a9f72bbca7 2150 /*!
stevep 0:04a9f72bbca7 2151 the remainder from a division
stevep 0:04a9f72bbca7 2152
stevep 0:04a9f72bbca7 2153 e.g.
stevep 0:04a9f72bbca7 2154 mod( 12.6 ; 3) = 0.6 because 12.6 = 3*4 + 0.6
stevep 0:04a9f72bbca7 2155 mod(-12.6 ; 3) = -0.6 bacause -12.6 = 3*(-4) + (-0.6)
stevep 0:04a9f72bbca7 2156 mod( 12.6 ; -3) = 0.6
stevep 0:04a9f72bbca7 2157 mod(-12.6 ; -3) = -0.6
stevep 0:04a9f72bbca7 2158 */
stevep 0:04a9f72bbca7 2159 template<class ValueType>
stevep 0:04a9f72bbca7 2160 ValueType Mod(ValueType a, const ValueType & b, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 2161 {
stevep 0:04a9f72bbca7 2162 if( a.IsNan() || b.IsNan() )
stevep 0:04a9f72bbca7 2163 {
stevep 0:04a9f72bbca7 2164 if( err )
stevep 0:04a9f72bbca7 2165 *err = err_improper_argument;
stevep 0:04a9f72bbca7 2166
stevep 0:04a9f72bbca7 2167 a.SetNan();
stevep 0:04a9f72bbca7 2168
stevep 0:04a9f72bbca7 2169 return a;
stevep 0:04a9f72bbca7 2170 }
stevep 0:04a9f72bbca7 2171
stevep 0:04a9f72bbca7 2172 uint c = a.Mod(b);
stevep 0:04a9f72bbca7 2173
stevep 0:04a9f72bbca7 2174 if( err )
stevep 0:04a9f72bbca7 2175 *err = c ? err_overflow : err_ok;
stevep 0:04a9f72bbca7 2176
stevep 0:04a9f72bbca7 2177 return a;
stevep 0:04a9f72bbca7 2178 }
stevep 0:04a9f72bbca7 2179
stevep 0:04a9f72bbca7 2180
stevep 0:04a9f72bbca7 2181
stevep 0:04a9f72bbca7 2182 namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 2183 {
stevep 0:04a9f72bbca7 2184
stevep 0:04a9f72bbca7 2185 /*!
stevep 0:04a9f72bbca7 2186 this function is used to store factorials in a given container
stevep 0:04a9f72bbca7 2187 'more' means how many values should be added at the end
stevep 0:04a9f72bbca7 2188
stevep 0:04a9f72bbca7 2189 e.g.
stevep 0:04a9f72bbca7 2190 std::vector<ValueType> fact;
stevep 0:04a9f72bbca7 2191 SetFactorialSequence(fact, 3);
stevep 0:04a9f72bbca7 2192 // now the container has three values: 1 1 2
stevep 0:04a9f72bbca7 2193
stevep 0:04a9f72bbca7 2194 SetFactorialSequence(fact, 2);
stevep 0:04a9f72bbca7 2195 // now the container has five values: 1 1 2 6 24
stevep 0:04a9f72bbca7 2196 */
stevep 0:04a9f72bbca7 2197 template<class ValueType>
stevep 0:04a9f72bbca7 2198 void SetFactorialSequence(std::vector<ValueType> & fact, uint more = 20)
stevep 0:04a9f72bbca7 2199 {
stevep 0:04a9f72bbca7 2200 if( more == 0 )
stevep 0:04a9f72bbca7 2201 more = 1;
stevep 0:04a9f72bbca7 2202
stevep 0:04a9f72bbca7 2203 uint start = static_cast<uint>(fact.size());
stevep 0:04a9f72bbca7 2204 fact.resize(fact.size() + more);
stevep 0:04a9f72bbca7 2205
stevep 0:04a9f72bbca7 2206 if( start == 0 )
stevep 0:04a9f72bbca7 2207 {
stevep 0:04a9f72bbca7 2208 fact[0] = 1;
stevep 0:04a9f72bbca7 2209 ++start;
stevep 0:04a9f72bbca7 2210 }
stevep 0:04a9f72bbca7 2211
stevep 0:04a9f72bbca7 2212 for(uint i=start ; i<fact.size() ; ++i)
stevep 0:04a9f72bbca7 2213 {
stevep 0:04a9f72bbca7 2214 fact[i] = fact[i-1];
stevep 0:04a9f72bbca7 2215 fact[i].MulInt(i);
stevep 0:04a9f72bbca7 2216 }
stevep 0:04a9f72bbca7 2217 }
stevep 0:04a9f72bbca7 2218
stevep 0:04a9f72bbca7 2219
stevep 0:04a9f72bbca7 2220 /*!
stevep 0:04a9f72bbca7 2221 an auxiliary function used to calculate Bernoulli numbers
stevep 0:04a9f72bbca7 2222
stevep 0:04a9f72bbca7 2223 this function returns a sum:
stevep 0:04a9f72bbca7 2224 sum(m) = sum_{k=0}^{m-1} {2^k * (m k) * B(k)} k in [0, m-1] (m k) means binomial coefficient = (m! / (k! * (m-k)!))
stevep 0:04a9f72bbca7 2225
stevep 0:04a9f72bbca7 2226 you should have sufficient factorials in cgamma.fact
stevep 0:04a9f72bbca7 2227 (cgamma.fact should have at least m items)
stevep 0:04a9f72bbca7 2228
stevep 0:04a9f72bbca7 2229 n_ should be equal 2
stevep 0:04a9f72bbca7 2230 */
stevep 0:04a9f72bbca7 2231 template<class ValueType>
stevep 0:04a9f72bbca7 2232 ValueType SetBernoulliNumbersSum(CGamma<ValueType> & cgamma, const ValueType & n_, uint m,
stevep 0:04a9f72bbca7 2233 const volatile StopCalculating * stop = 0)
stevep 0:04a9f72bbca7 2234 {
stevep 0:04a9f72bbca7 2235 ValueType k_, temp, temp2, temp3, sum;
stevep 0:04a9f72bbca7 2236
stevep 0:04a9f72bbca7 2237 sum.SetZero();
stevep 0:04a9f72bbca7 2238
stevep 0:04a9f72bbca7 2239 for(uint k=0 ; k<m ; ++k) // k<m means k<=m-1
stevep 0:04a9f72bbca7 2240 {
stevep 0:04a9f72bbca7 2241 if( stop && (k & 15)==0 ) // means: k % 16 == 0
stevep 0:04a9f72bbca7 2242 if( stop->WasStopSignal() )
stevep 0:04a9f72bbca7 2243 return ValueType(); // NaN
stevep 0:04a9f72bbca7 2244
stevep 0:04a9f72bbca7 2245 if( k>1 && (k & 1) == 1 ) // for that k the Bernoulli number is zero
stevep 0:04a9f72bbca7 2246 continue;
stevep 0:04a9f72bbca7 2247
stevep 0:04a9f72bbca7 2248 k_ = k;
stevep 0:04a9f72bbca7 2249
stevep 0:04a9f72bbca7 2250 temp = n_; // n_ is equal 2
stevep 0:04a9f72bbca7 2251 temp.Pow(k_);
stevep 0:04a9f72bbca7 2252 // temp = 2^k
stevep 0:04a9f72bbca7 2253
stevep 0:04a9f72bbca7 2254 temp2 = cgamma.fact[m];
stevep 0:04a9f72bbca7 2255 temp3 = cgamma.fact[k];
stevep 0:04a9f72bbca7 2256 temp3.Mul(cgamma.fact[m-k]);
stevep 0:04a9f72bbca7 2257 temp2.Div(temp3);
stevep 0:04a9f72bbca7 2258 // temp2 = (m k) = m! / ( k! * (m-k)! )
stevep 0:04a9f72bbca7 2259
stevep 0:04a9f72bbca7 2260 temp.Mul(temp2);
stevep 0:04a9f72bbca7 2261 temp.Mul(cgamma.bern[k]);
stevep 0:04a9f72bbca7 2262
stevep 0:04a9f72bbca7 2263 sum.Add(temp);
stevep 0:04a9f72bbca7 2264 // sum += 2^k * (m k) * B(k)
stevep 0:04a9f72bbca7 2265
stevep 0:04a9f72bbca7 2266 if( sum.IsNan() )
stevep 0:04a9f72bbca7 2267 break;
stevep 0:04a9f72bbca7 2268 }
stevep 0:04a9f72bbca7 2269
stevep 0:04a9f72bbca7 2270 return sum;
stevep 0:04a9f72bbca7 2271 }
stevep 0:04a9f72bbca7 2272
stevep 0:04a9f72bbca7 2273
stevep 0:04a9f72bbca7 2274 /*!
stevep 0:04a9f72bbca7 2275 an auxiliary function used to calculate Bernoulli numbers
stevep 0:04a9f72bbca7 2276 start is >= 2
stevep 0:04a9f72bbca7 2277
stevep 0:04a9f72bbca7 2278 we use the recurrence formula:
stevep 0:04a9f72bbca7 2279 B(m) = 1 / (2*(1 - 2^m)) * sum(m)
stevep 0:04a9f72bbca7 2280 where sum(m) is calculated by SetBernoulliNumbersSum()
stevep 0:04a9f72bbca7 2281 */
stevep 0:04a9f72bbca7 2282 template<class ValueType>
stevep 0:04a9f72bbca7 2283 bool SetBernoulliNumbersMore(CGamma<ValueType> & cgamma, uint start, const volatile StopCalculating * stop = 0)
stevep 0:04a9f72bbca7 2284 {
stevep 0:04a9f72bbca7 2285 ValueType denominator, temp, temp2, temp3, m_, sum, sum2, n_, k_;
stevep 0:04a9f72bbca7 2286
stevep 0:04a9f72bbca7 2287 const uint n = 2;
stevep 0:04a9f72bbca7 2288 n_ = n;
stevep 0:04a9f72bbca7 2289
stevep 0:04a9f72bbca7 2290 // start is >= 2
stevep 0:04a9f72bbca7 2291 for(uint m=start ; m<cgamma.bern.size() ; ++m)
stevep 0:04a9f72bbca7 2292 {
stevep 0:04a9f72bbca7 2293 if( (m & 1) == 1 )
stevep 0:04a9f72bbca7 2294 {
stevep 0:04a9f72bbca7 2295 cgamma.bern[m].SetZero();
stevep 0:04a9f72bbca7 2296 }
stevep 0:04a9f72bbca7 2297 else
stevep 0:04a9f72bbca7 2298 {
stevep 0:04a9f72bbca7 2299 m_ = m;
stevep 0:04a9f72bbca7 2300
stevep 0:04a9f72bbca7 2301 temp = n_; // n_ = 2
stevep 0:04a9f72bbca7 2302 temp.Pow(m_);
stevep 0:04a9f72bbca7 2303 // temp = 2^m
stevep 0:04a9f72bbca7 2304
stevep 0:04a9f72bbca7 2305 denominator.SetOne();
stevep 0:04a9f72bbca7 2306 denominator.Sub(temp);
stevep 0:04a9f72bbca7 2307 if( denominator.exponent.AddOne() ) // it means: denominator.MulInt(2)
stevep 0:04a9f72bbca7 2308 denominator.SetNan();
stevep 0:04a9f72bbca7 2309
stevep 0:04a9f72bbca7 2310 // denominator = 2 * (1 - 2^m)
stevep 0:04a9f72bbca7 2311
stevep 0:04a9f72bbca7 2312 cgamma.bern[m] = SetBernoulliNumbersSum(cgamma, n_, m, stop);
stevep 0:04a9f72bbca7 2313
stevep 0:04a9f72bbca7 2314 if( stop && stop->WasStopSignal() )
stevep 0:04a9f72bbca7 2315 {
stevep 0:04a9f72bbca7 2316 cgamma.bern.resize(m); // valid numbers are in [0, m-1]
stevep 0:04a9f72bbca7 2317 return false;
stevep 0:04a9f72bbca7 2318 }
stevep 0:04a9f72bbca7 2319
stevep 0:04a9f72bbca7 2320 cgamma.bern[m].Div(denominator);
stevep 0:04a9f72bbca7 2321 }
stevep 0:04a9f72bbca7 2322 }
stevep 0:04a9f72bbca7 2323
stevep 0:04a9f72bbca7 2324 return true;
stevep 0:04a9f72bbca7 2325 }
stevep 0:04a9f72bbca7 2326
stevep 0:04a9f72bbca7 2327
stevep 0:04a9f72bbca7 2328 /*!
stevep 0:04a9f72bbca7 2329 this function is used to calculate Bernoulli numbers,
stevep 0:04a9f72bbca7 2330 returns false if there was a stop signal,
stevep 0:04a9f72bbca7 2331 'more' means how many values should be added at the end
stevep 0:04a9f72bbca7 2332
stevep 0:04a9f72bbca7 2333 e.g.
stevep 0:04a9f72bbca7 2334 typedef Big<1,2> MyBig;
stevep 0:04a9f72bbca7 2335 CGamma<MyBig> cgamma;
stevep 0:04a9f72bbca7 2336 SetBernoulliNumbers(cgamma, 3);
stevep 0:04a9f72bbca7 2337 // now we have three first Bernoulli numbers: 1 -0.5 0.16667
stevep 0:04a9f72bbca7 2338
stevep 0:04a9f72bbca7 2339 SetBernoulliNumbers(cgamma, 4);
stevep 0:04a9f72bbca7 2340 // now we have 7 Bernoulli numbers: 1 -0.5 0.16667 0 -0.0333 0 0.0238
stevep 0:04a9f72bbca7 2341 */
stevep 0:04a9f72bbca7 2342 template<class ValueType>
stevep 0:04a9f72bbca7 2343 bool SetBernoulliNumbers(CGamma<ValueType> & cgamma, uint more = 20, const volatile StopCalculating * stop = 0)
stevep 0:04a9f72bbca7 2344 {
stevep 0:04a9f72bbca7 2345 if( more == 0 )
stevep 0:04a9f72bbca7 2346 more = 1;
stevep 0:04a9f72bbca7 2347
stevep 0:04a9f72bbca7 2348 uint start = static_cast<uint>(cgamma.bern.size());
stevep 0:04a9f72bbca7 2349 cgamma.bern.resize(cgamma.bern.size() + more);
stevep 0:04a9f72bbca7 2350
stevep 0:04a9f72bbca7 2351 if( start == 0 )
stevep 0:04a9f72bbca7 2352 {
stevep 0:04a9f72bbca7 2353 cgamma.bern[0].SetOne();
stevep 0:04a9f72bbca7 2354 ++start;
stevep 0:04a9f72bbca7 2355 }
stevep 0:04a9f72bbca7 2356
stevep 0:04a9f72bbca7 2357 if( cgamma.bern.size() == 1 )
stevep 0:04a9f72bbca7 2358 return true;
stevep 0:04a9f72bbca7 2359
stevep 0:04a9f72bbca7 2360 if( start == 1 )
stevep 0:04a9f72bbca7 2361 {
stevep 0:04a9f72bbca7 2362 cgamma.bern[1].Set05();
stevep 0:04a9f72bbca7 2363 cgamma.bern[1].ChangeSign();
stevep 0:04a9f72bbca7 2364 ++start;
stevep 0:04a9f72bbca7 2365 }
stevep 0:04a9f72bbca7 2366
stevep 0:04a9f72bbca7 2367 // we should have sufficient factorials in cgamma.fact
stevep 0:04a9f72bbca7 2368 if( cgamma.fact.size() < cgamma.bern.size() )
stevep 0:04a9f72bbca7 2369 SetFactorialSequence(cgamma.fact, static_cast<uint>(cgamma.bern.size() - cgamma.fact.size()));
stevep 0:04a9f72bbca7 2370
stevep 0:04a9f72bbca7 2371
stevep 0:04a9f72bbca7 2372 return SetBernoulliNumbersMore(cgamma, start, stop);
stevep 0:04a9f72bbca7 2373 }
stevep 0:04a9f72bbca7 2374
stevep 0:04a9f72bbca7 2375
stevep 0:04a9f72bbca7 2376 /*!
stevep 0:04a9f72bbca7 2377 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2378
stevep 0:04a9f72bbca7 2379 we calculate a sum:
stevep 0:04a9f72bbca7 2380 sum(n) = sum_{m=2} { B(m) / ( (m^2 - m) * n^(m-1) ) } = 1/(12*n) - 1/(360*n^3) + 1/(1260*n^5) + ...
stevep 0:04a9f72bbca7 2381 B(m) means a mth Bernoulli number
stevep 0:04a9f72bbca7 2382 the sum starts from m=2, we calculate as long as the value will not change after adding a next part
stevep 0:04a9f72bbca7 2383 */
stevep 0:04a9f72bbca7 2384 template<class ValueType>
stevep 0:04a9f72bbca7 2385 ValueType GammaFactorialHighSum(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
stevep 0:04a9f72bbca7 2386 const volatile StopCalculating * stop)
stevep 0:04a9f72bbca7 2387 {
stevep 0:04a9f72bbca7 2388 ValueType temp, temp2, denominator, sum, oldsum;
stevep 0:04a9f72bbca7 2389
stevep 0:04a9f72bbca7 2390 sum.SetZero();
stevep 0:04a9f72bbca7 2391
stevep 0:04a9f72bbca7 2392 for(uint m=2 ; m<TTMATH_ARITHMETIC_MAX_LOOP ; m+=2)
stevep 0:04a9f72bbca7 2393 {
stevep 0:04a9f72bbca7 2394 if( stop && (m & 3)==0 ) // (m & 3)==0 means: (m % 4)==0
stevep 0:04a9f72bbca7 2395 if( stop->WasStopSignal() )
stevep 0:04a9f72bbca7 2396 {
stevep 0:04a9f72bbca7 2397 err = err_interrupt;
stevep 0:04a9f72bbca7 2398 return ValueType(); // NaN
stevep 0:04a9f72bbca7 2399 }
stevep 0:04a9f72bbca7 2400
stevep 0:04a9f72bbca7 2401 temp = (m-1);
stevep 0:04a9f72bbca7 2402 denominator = n;
stevep 0:04a9f72bbca7 2403 denominator.Pow(temp);
stevep 0:04a9f72bbca7 2404 // denominator = n ^ (m-1)
stevep 0:04a9f72bbca7 2405
stevep 0:04a9f72bbca7 2406 temp = m;
stevep 0:04a9f72bbca7 2407 temp2 = temp;
stevep 0:04a9f72bbca7 2408 temp.Mul(temp2);
stevep 0:04a9f72bbca7 2409 temp.Sub(temp2);
stevep 0:04a9f72bbca7 2410 // temp = m^2 - m
stevep 0:04a9f72bbca7 2411
stevep 0:04a9f72bbca7 2412 denominator.Mul(temp);
stevep 0:04a9f72bbca7 2413 // denominator = (m^2 - m) * n ^ (m-1)
stevep 0:04a9f72bbca7 2414
stevep 0:04a9f72bbca7 2415 if( m >= cgamma.bern.size() )
stevep 0:04a9f72bbca7 2416 {
stevep 0:04a9f72bbca7 2417 if( !SetBernoulliNumbers(cgamma, m - cgamma.bern.size() + 1 + 3, stop) ) // 3 more than needed
stevep 0:04a9f72bbca7 2418 {
stevep 0:04a9f72bbca7 2419 // there was the stop signal
stevep 0:04a9f72bbca7 2420 err = err_interrupt;
stevep 0:04a9f72bbca7 2421 return ValueType(); // NaN
stevep 0:04a9f72bbca7 2422 }
stevep 0:04a9f72bbca7 2423 }
stevep 0:04a9f72bbca7 2424
stevep 0:04a9f72bbca7 2425 temp = cgamma.bern[m];
stevep 0:04a9f72bbca7 2426 temp.Div(denominator);
stevep 0:04a9f72bbca7 2427
stevep 0:04a9f72bbca7 2428 oldsum = sum;
stevep 0:04a9f72bbca7 2429 sum.Add(temp);
stevep 0:04a9f72bbca7 2430
stevep 0:04a9f72bbca7 2431 if( sum.IsNan() || oldsum==sum )
stevep 0:04a9f72bbca7 2432 break;
stevep 0:04a9f72bbca7 2433 }
stevep 0:04a9f72bbca7 2434
stevep 0:04a9f72bbca7 2435 return sum;
stevep 0:04a9f72bbca7 2436 }
stevep 0:04a9f72bbca7 2437
stevep 0:04a9f72bbca7 2438
stevep 0:04a9f72bbca7 2439 /*!
stevep 0:04a9f72bbca7 2440 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2441
stevep 0:04a9f72bbca7 2442 we calculate a helper function GammaFactorialHigh() by using Stirling's series:
stevep 0:04a9f72bbca7 2443 n! = (n/e)^n * sqrt(2*pi*n) * exp( sum(n) )
stevep 0:04a9f72bbca7 2444 where n is a real number (not only an integer) and is sufficient large (greater than TTMATH_GAMMA_BOUNDARY)
stevep 0:04a9f72bbca7 2445 and sum(n) is calculated by GammaFactorialHighSum()
stevep 0:04a9f72bbca7 2446 */
stevep 0:04a9f72bbca7 2447 template<class ValueType>
stevep 0:04a9f72bbca7 2448 ValueType GammaFactorialHigh(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
stevep 0:04a9f72bbca7 2449 const volatile StopCalculating * stop)
stevep 0:04a9f72bbca7 2450 {
stevep 0:04a9f72bbca7 2451 ValueType temp, temp2, temp3, denominator, sum;
stevep 0:04a9f72bbca7 2452
stevep 0:04a9f72bbca7 2453 temp.Set2Pi();
stevep 0:04a9f72bbca7 2454 temp.Mul(n);
stevep 0:04a9f72bbca7 2455 temp2 = Sqrt(temp);
stevep 0:04a9f72bbca7 2456 // temp2 = sqrt(2*pi*n)
stevep 0:04a9f72bbca7 2457
stevep 0:04a9f72bbca7 2458 temp = n;
stevep 0:04a9f72bbca7 2459 temp3.SetE();
stevep 0:04a9f72bbca7 2460 temp.Div(temp3);
stevep 0:04a9f72bbca7 2461 temp.Pow(n);
stevep 0:04a9f72bbca7 2462 // temp = (n/e)^n
stevep 0:04a9f72bbca7 2463
stevep 0:04a9f72bbca7 2464 sum = GammaFactorialHighSum(n, cgamma, err, stop);
stevep 0:04a9f72bbca7 2465 temp3.Exp(sum);
stevep 0:04a9f72bbca7 2466 // temp3 = exp(sum)
stevep 0:04a9f72bbca7 2467
stevep 0:04a9f72bbca7 2468 temp.Mul(temp2);
stevep 0:04a9f72bbca7 2469 temp.Mul(temp3);
stevep 0:04a9f72bbca7 2470
stevep 0:04a9f72bbca7 2471 return temp;
stevep 0:04a9f72bbca7 2472 }
stevep 0:04a9f72bbca7 2473
stevep 0:04a9f72bbca7 2474
stevep 0:04a9f72bbca7 2475 /*!
stevep 0:04a9f72bbca7 2476 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2477
stevep 0:04a9f72bbca7 2478 Gamma(x) = GammaFactorialHigh(x-1)
stevep 0:04a9f72bbca7 2479 */
stevep 0:04a9f72bbca7 2480 template<class ValueType>
stevep 0:04a9f72bbca7 2481 ValueType GammaPlusHigh(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
stevep 0:04a9f72bbca7 2482 {
stevep 0:04a9f72bbca7 2483 ValueType one;
stevep 0:04a9f72bbca7 2484
stevep 0:04a9f72bbca7 2485 one.SetOne();
stevep 0:04a9f72bbca7 2486 n.Sub(one);
stevep 0:04a9f72bbca7 2487
stevep 0:04a9f72bbca7 2488 return GammaFactorialHigh(n, cgamma, err, stop);
stevep 0:04a9f72bbca7 2489 }
stevep 0:04a9f72bbca7 2490
stevep 0:04a9f72bbca7 2491
stevep 0:04a9f72bbca7 2492 /*!
stevep 0:04a9f72bbca7 2493 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2494
stevep 0:04a9f72bbca7 2495 we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
stevep 0:04a9f72bbca7 2496 we use the formula:
stevep 0:04a9f72bbca7 2497 gamma(n) = (n-1)! = 1 * 2 * 3 * ... * (n-1)
stevep 0:04a9f72bbca7 2498 */
stevep 0:04a9f72bbca7 2499 template<class ValueType>
stevep 0:04a9f72bbca7 2500 ValueType GammaPlusLowIntegerInt(uint n, CGamma<ValueType> & cgamma)
stevep 0:04a9f72bbca7 2501 {
stevep 0:04a9f72bbca7 2502 TTMATH_ASSERT( n > 0 )
stevep 0:04a9f72bbca7 2503
stevep 0:04a9f72bbca7 2504 if( n - 1 < static_cast<uint>(cgamma.fact.size()) )
stevep 0:04a9f72bbca7 2505 return cgamma.fact[n - 1];
stevep 0:04a9f72bbca7 2506
stevep 0:04a9f72bbca7 2507 ValueType res;
stevep 0:04a9f72bbca7 2508 uint start = 2;
stevep 0:04a9f72bbca7 2509
stevep 0:04a9f72bbca7 2510 if( cgamma.fact.size() < 2 )
stevep 0:04a9f72bbca7 2511 {
stevep 0:04a9f72bbca7 2512 res.SetOne();
stevep 0:04a9f72bbca7 2513 }
stevep 0:04a9f72bbca7 2514 else
stevep 0:04a9f72bbca7 2515 {
stevep 0:04a9f72bbca7 2516 start = static_cast<uint>(cgamma.fact.size());
stevep 0:04a9f72bbca7 2517 res = cgamma.fact[start-1];
stevep 0:04a9f72bbca7 2518 }
stevep 0:04a9f72bbca7 2519
stevep 0:04a9f72bbca7 2520 for(uint i=start ; i<n ; ++i)
stevep 0:04a9f72bbca7 2521 res.MulInt(i);
stevep 0:04a9f72bbca7 2522
stevep 0:04a9f72bbca7 2523 return res;
stevep 0:04a9f72bbca7 2524 }
stevep 0:04a9f72bbca7 2525
stevep 0:04a9f72bbca7 2526
stevep 0:04a9f72bbca7 2527 /*!
stevep 0:04a9f72bbca7 2528 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2529
stevep 0:04a9f72bbca7 2530 we use this function when n is integer and a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
stevep 0:04a9f72bbca7 2531 */
stevep 0:04a9f72bbca7 2532 template<class ValueType>
stevep 0:04a9f72bbca7 2533 ValueType GammaPlusLowInteger(const ValueType & n, CGamma<ValueType> & cgamma)
stevep 0:04a9f72bbca7 2534 {
stevep 0:04a9f72bbca7 2535 sint n_;
stevep 0:04a9f72bbca7 2536
stevep 0:04a9f72bbca7 2537 n.ToInt(n_);
stevep 0:04a9f72bbca7 2538
stevep 0:04a9f72bbca7 2539 return GammaPlusLowIntegerInt(n_, cgamma);
stevep 0:04a9f72bbca7 2540 }
stevep 0:04a9f72bbca7 2541
stevep 0:04a9f72bbca7 2542
stevep 0:04a9f72bbca7 2543 /*!
stevep 0:04a9f72bbca7 2544 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2545
stevep 0:04a9f72bbca7 2546 we use this function when n is a small value (from 0 to TTMATH_GAMMA_BOUNDARY]
stevep 0:04a9f72bbca7 2547 we use a recurrence formula:
stevep 0:04a9f72bbca7 2548 gamma(z+1) = z * gamma(z)
stevep 0:04a9f72bbca7 2549 then: gamma(z) = gamma(z+1) / z
stevep 0:04a9f72bbca7 2550
stevep 0:04a9f72bbca7 2551 e.g.
stevep 0:04a9f72bbca7 2552 gamma(3.89) = gamma(2001.89) / ( 3.89 * 4.89 * 5.89 * ... * 1999.89 * 2000.89 )
stevep 0:04a9f72bbca7 2553 */
stevep 0:04a9f72bbca7 2554 template<class ValueType>
stevep 0:04a9f72bbca7 2555 ValueType GammaPlusLow(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
stevep 0:04a9f72bbca7 2556 {
stevep 0:04a9f72bbca7 2557 ValueType one, denominator, temp, boundary;
stevep 0:04a9f72bbca7 2558
stevep 0:04a9f72bbca7 2559 if( n.IsInteger() )
stevep 0:04a9f72bbca7 2560 return GammaPlusLowInteger(n, cgamma);
stevep 0:04a9f72bbca7 2561
stevep 0:04a9f72bbca7 2562 one.SetOne();
stevep 0:04a9f72bbca7 2563 denominator = n;
stevep 0:04a9f72bbca7 2564 boundary = TTMATH_GAMMA_BOUNDARY;
stevep 0:04a9f72bbca7 2565
stevep 0:04a9f72bbca7 2566 while( n < boundary )
stevep 0:04a9f72bbca7 2567 {
stevep 0:04a9f72bbca7 2568 n.Add(one);
stevep 0:04a9f72bbca7 2569 denominator.Mul(n);
stevep 0:04a9f72bbca7 2570 }
stevep 0:04a9f72bbca7 2571
stevep 0:04a9f72bbca7 2572 n.Add(one);
stevep 0:04a9f72bbca7 2573
stevep 0:04a9f72bbca7 2574 // now n is sufficient big
stevep 0:04a9f72bbca7 2575 temp = GammaPlusHigh(n, cgamma, err, stop);
stevep 0:04a9f72bbca7 2576 temp.Div(denominator);
stevep 0:04a9f72bbca7 2577
stevep 0:04a9f72bbca7 2578 return temp;
stevep 0:04a9f72bbca7 2579 }
stevep 0:04a9f72bbca7 2580
stevep 0:04a9f72bbca7 2581
stevep 0:04a9f72bbca7 2582 /*!
stevep 0:04a9f72bbca7 2583 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2584 */
stevep 0:04a9f72bbca7 2585 template<class ValueType>
stevep 0:04a9f72bbca7 2586 ValueType GammaPlus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
stevep 0:04a9f72bbca7 2587 {
stevep 0:04a9f72bbca7 2588 if( n > TTMATH_GAMMA_BOUNDARY )
stevep 0:04a9f72bbca7 2589 return GammaPlusHigh(n, cgamma, err, stop);
stevep 0:04a9f72bbca7 2590
stevep 0:04a9f72bbca7 2591 return GammaPlusLow(n, cgamma, err, stop);
stevep 0:04a9f72bbca7 2592 }
stevep 0:04a9f72bbca7 2593
stevep 0:04a9f72bbca7 2594
stevep 0:04a9f72bbca7 2595 /*!
stevep 0:04a9f72bbca7 2596 an auxiliary function used to calculate the Gamma() function
stevep 0:04a9f72bbca7 2597
stevep 0:04a9f72bbca7 2598 this function is used when n is negative
stevep 0:04a9f72bbca7 2599 we use the reflection formula:
stevep 0:04a9f72bbca7 2600 gamma(1-z) * gamma(z) = pi / sin(pi*z)
stevep 0:04a9f72bbca7 2601 then: gamma(z) = pi / (sin(pi*z) * gamma(1-z))
stevep 0:04a9f72bbca7 2602
stevep 0:04a9f72bbca7 2603 */
stevep 0:04a9f72bbca7 2604 template<class ValueType>
stevep 0:04a9f72bbca7 2605 ValueType GammaMinus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
stevep 0:04a9f72bbca7 2606 {
stevep 0:04a9f72bbca7 2607 ValueType pi, denominator, temp, temp2;
stevep 0:04a9f72bbca7 2608
stevep 0:04a9f72bbca7 2609 if( n.IsInteger() )
stevep 0:04a9f72bbca7 2610 {
stevep 0:04a9f72bbca7 2611 // gamma function is not defined when n is negative and integer
stevep 0:04a9f72bbca7 2612 err = err_improper_argument;
stevep 0:04a9f72bbca7 2613 return temp; // NaN
stevep 0:04a9f72bbca7 2614 }
stevep 0:04a9f72bbca7 2615
stevep 0:04a9f72bbca7 2616 pi.SetPi();
stevep 0:04a9f72bbca7 2617
stevep 0:04a9f72bbca7 2618 temp = pi;
stevep 0:04a9f72bbca7 2619 temp.Mul(n);
stevep 0:04a9f72bbca7 2620 temp2 = Sin(temp);
stevep 0:04a9f72bbca7 2621 // temp2 = sin(pi * n)
stevep 0:04a9f72bbca7 2622
stevep 0:04a9f72bbca7 2623 temp.SetOne();
stevep 0:04a9f72bbca7 2624 temp.Sub(n);
stevep 0:04a9f72bbca7 2625 temp = GammaPlus(temp, cgamma, err, stop);
stevep 0:04a9f72bbca7 2626 // temp = gamma(1 - n)
stevep 0:04a9f72bbca7 2627
stevep 0:04a9f72bbca7 2628 temp.Mul(temp2);
stevep 0:04a9f72bbca7 2629 pi.Div(temp);
stevep 0:04a9f72bbca7 2630
stevep 0:04a9f72bbca7 2631 return pi;
stevep 0:04a9f72bbca7 2632 }
stevep 0:04a9f72bbca7 2633
stevep 0:04a9f72bbca7 2634 } // namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 2635
stevep 0:04a9f72bbca7 2636
stevep 0:04a9f72bbca7 2637
stevep 0:04a9f72bbca7 2638 /*!
stevep 0:04a9f72bbca7 2639 this function calculates the Gamma function
stevep 0:04a9f72bbca7 2640
stevep 0:04a9f72bbca7 2641 it's multithread safe, you should create a CGamma<> object and use it whenever you call the Gamma()
stevep 0:04a9f72bbca7 2642 e.g.
stevep 0:04a9f72bbca7 2643 typedef Big<1,2> MyBig;
stevep 0:04a9f72bbca7 2644 MyBig x=234, y=345.53;
stevep 0:04a9f72bbca7 2645 CGamma<MyBig> cgamma;
stevep 0:04a9f72bbca7 2646 std::cout << Gamma(x, cgamma) << std::endl;
stevep 0:04a9f72bbca7 2647 std::cout << Gamma(y, cgamma) << std::endl;
stevep 0:04a9f72bbca7 2648 in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers),
stevep 0:04a9f72bbca7 2649 and they will be reused in next calls to the function
stevep 0:04a9f72bbca7 2650
stevep 0:04a9f72bbca7 2651 each thread should have its own CGamma<> object, and you can use these objects with Factorial() function too
stevep 0:04a9f72bbca7 2652 */
stevep 0:04a9f72bbca7 2653 template<class ValueType>
stevep 0:04a9f72bbca7 2654 ValueType Gamma(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
stevep 0:04a9f72bbca7 2655 const volatile StopCalculating * stop = 0)
stevep 0:04a9f72bbca7 2656 {
stevep 0:04a9f72bbca7 2657 using namespace auxiliaryfunctions;
stevep 0:04a9f72bbca7 2658
stevep 0:04a9f72bbca7 2659 ValueType result;
stevep 0:04a9f72bbca7 2660 ErrorCode err_tmp;
stevep 0:04a9f72bbca7 2661
stevep 0:04a9f72bbca7 2662 if( n.IsNan() )
stevep 0:04a9f72bbca7 2663 {
stevep 0:04a9f72bbca7 2664 if( err )
stevep 0:04a9f72bbca7 2665 *err = err_improper_argument;
stevep 0:04a9f72bbca7 2666
stevep 0:04a9f72bbca7 2667 return n;
stevep 0:04a9f72bbca7 2668 }
stevep 0:04a9f72bbca7 2669
stevep 0:04a9f72bbca7 2670 if( cgamma.history.Get(n, result, err_tmp) )
stevep 0:04a9f72bbca7 2671 {
stevep 0:04a9f72bbca7 2672 if( err )
stevep 0:04a9f72bbca7 2673 *err = err_tmp;
stevep 0:04a9f72bbca7 2674
stevep 0:04a9f72bbca7 2675 return result;
stevep 0:04a9f72bbca7 2676 }
stevep 0:04a9f72bbca7 2677
stevep 0:04a9f72bbca7 2678 err_tmp = err_ok;
stevep 0:04a9f72bbca7 2679
stevep 0:04a9f72bbca7 2680 if( n.IsSign() )
stevep 0:04a9f72bbca7 2681 {
stevep 0:04a9f72bbca7 2682 result = GammaMinus(n, cgamma, err_tmp, stop);
stevep 0:04a9f72bbca7 2683 }
stevep 0:04a9f72bbca7 2684 else
stevep 0:04a9f72bbca7 2685 if( n.IsZero() )
stevep 0:04a9f72bbca7 2686 {
stevep 0:04a9f72bbca7 2687 err_tmp = err_improper_argument;
stevep 0:04a9f72bbca7 2688 result.SetNan();
stevep 0:04a9f72bbca7 2689 }
stevep 0:04a9f72bbca7 2690 else
stevep 0:04a9f72bbca7 2691 {
stevep 0:04a9f72bbca7 2692 result = GammaPlus(n, cgamma, err_tmp, stop);
stevep 0:04a9f72bbca7 2693 }
stevep 0:04a9f72bbca7 2694
stevep 0:04a9f72bbca7 2695 if( result.IsNan() && err_tmp==err_ok )
stevep 0:04a9f72bbca7 2696 err_tmp = err_overflow;
stevep 0:04a9f72bbca7 2697
stevep 0:04a9f72bbca7 2698 if( err )
stevep 0:04a9f72bbca7 2699 *err = err_tmp;
stevep 0:04a9f72bbca7 2700
stevep 0:04a9f72bbca7 2701 if( stop && !stop->WasStopSignal() )
stevep 0:04a9f72bbca7 2702 cgamma.history.Add(n, result, err_tmp);
stevep 0:04a9f72bbca7 2703
stevep 0:04a9f72bbca7 2704 return result;
stevep 0:04a9f72bbca7 2705 }
stevep 0:04a9f72bbca7 2706
stevep 0:04a9f72bbca7 2707
stevep 0:04a9f72bbca7 2708 /*!
stevep 0:04a9f72bbca7 2709 this function calculates the Gamma function
stevep 0:04a9f72bbca7 2710
stevep 0:04a9f72bbca7 2711 note: this function should be used only in a single-thread environment
stevep 0:04a9f72bbca7 2712 */
stevep 0:04a9f72bbca7 2713 template<class ValueType>
stevep 0:04a9f72bbca7 2714 ValueType Gamma(const ValueType & n, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 2715 {
stevep 0:04a9f72bbca7 2716 // warning: this static object is not thread safe
stevep 0:04a9f72bbca7 2717 static CGamma<ValueType> cgamma;
stevep 0:04a9f72bbca7 2718
stevep 0:04a9f72bbca7 2719 return Gamma(n, cgamma, err);
stevep 0:04a9f72bbca7 2720 }
stevep 0:04a9f72bbca7 2721
stevep 0:04a9f72bbca7 2722
stevep 0:04a9f72bbca7 2723
stevep 0:04a9f72bbca7 2724 namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 2725 {
stevep 0:04a9f72bbca7 2726
stevep 0:04a9f72bbca7 2727 /*!
stevep 0:04a9f72bbca7 2728 an auxiliary function for calculating the factorial function
stevep 0:04a9f72bbca7 2729
stevep 0:04a9f72bbca7 2730 we use the formula:
stevep 0:04a9f72bbca7 2731 x! = gamma(x+1)
stevep 0:04a9f72bbca7 2732 */
stevep 0:04a9f72bbca7 2733 template<class ValueType>
stevep 0:04a9f72bbca7 2734 ValueType Factorial2(ValueType x,
stevep 0:04a9f72bbca7 2735 CGamma<ValueType> * cgamma = 0,
stevep 0:04a9f72bbca7 2736 ErrorCode * err = 0,
stevep 0:04a9f72bbca7 2737 const volatile StopCalculating * stop = 0)
stevep 0:04a9f72bbca7 2738 {
stevep 0:04a9f72bbca7 2739 ValueType result, one;
stevep 0:04a9f72bbca7 2740
stevep 0:04a9f72bbca7 2741 if( x.IsNan() || x.IsSign() || !x.IsInteger() )
stevep 0:04a9f72bbca7 2742 {
stevep 0:04a9f72bbca7 2743 if( err )
stevep 0:04a9f72bbca7 2744 *err = err_improper_argument;
stevep 0:04a9f72bbca7 2745
stevep 0:04a9f72bbca7 2746 x.SetNan();
stevep 0:04a9f72bbca7 2747
stevep 0:04a9f72bbca7 2748 return x;
stevep 0:04a9f72bbca7 2749 }
stevep 0:04a9f72bbca7 2750
stevep 0:04a9f72bbca7 2751 one.SetOne();
stevep 0:04a9f72bbca7 2752 x.Add(one);
stevep 0:04a9f72bbca7 2753
stevep 0:04a9f72bbca7 2754 if( cgamma )
stevep 0:04a9f72bbca7 2755 return Gamma(x, *cgamma, err, stop);
stevep 0:04a9f72bbca7 2756
stevep 0:04a9f72bbca7 2757 return Gamma(x, err);
stevep 0:04a9f72bbca7 2758 }
stevep 0:04a9f72bbca7 2759
stevep 0:04a9f72bbca7 2760 } // namespace auxiliaryfunctions
stevep 0:04a9f72bbca7 2761
stevep 0:04a9f72bbca7 2762
stevep 0:04a9f72bbca7 2763
stevep 0:04a9f72bbca7 2764 /*!
stevep 0:04a9f72bbca7 2765 the factorial from given 'x'
stevep 0:04a9f72bbca7 2766 e.g.
stevep 0:04a9f72bbca7 2767 Factorial(4) = 4! = 1*2*3*4
stevep 0:04a9f72bbca7 2768
stevep 0:04a9f72bbca7 2769 it's multithread safe, you should create a CGamma<> object and use it whenever you call the Factorial()
stevep 0:04a9f72bbca7 2770 e.g.
stevep 0:04a9f72bbca7 2771 typedef Big<1,2> MyBig;
stevep 0:04a9f72bbca7 2772 MyBig x=234, y=54345;
stevep 0:04a9f72bbca7 2773 CGamma<MyBig> cgamma;
stevep 0:04a9f72bbca7 2774 std::cout << Factorial(x, cgamma) << std::endl;
stevep 0:04a9f72bbca7 2775 std::cout << Factorial(y, cgamma) << std::endl;
stevep 0:04a9f72bbca7 2776 in the CGamma<> object the function stores some coefficients (factorials, Bernoulli numbers),
stevep 0:04a9f72bbca7 2777 and they will be reused in next calls to the function
stevep 0:04a9f72bbca7 2778
stevep 0:04a9f72bbca7 2779 each thread should have its own CGamma<> object, and you can use these objects with Gamma() function too
stevep 0:04a9f72bbca7 2780 */
stevep 0:04a9f72bbca7 2781 template<class ValueType>
stevep 0:04a9f72bbca7 2782 ValueType Factorial(const ValueType & x, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
stevep 0:04a9f72bbca7 2783 const volatile StopCalculating * stop = 0)
stevep 0:04a9f72bbca7 2784 {
stevep 0:04a9f72bbca7 2785 return auxiliaryfunctions::Factorial2(x, &cgamma, err, stop);
stevep 0:04a9f72bbca7 2786 }
stevep 0:04a9f72bbca7 2787
stevep 0:04a9f72bbca7 2788
stevep 0:04a9f72bbca7 2789 /*!
stevep 0:04a9f72bbca7 2790 the factorial from given 'x'
stevep 0:04a9f72bbca7 2791 e.g.
stevep 0:04a9f72bbca7 2792 Factorial(4) = 4! = 1*2*3*4
stevep 0:04a9f72bbca7 2793
stevep 0:04a9f72bbca7 2794 note: this function should be used only in a single-thread environment
stevep 0:04a9f72bbca7 2795 */
stevep 0:04a9f72bbca7 2796 template<class ValueType>
stevep 0:04a9f72bbca7 2797 ValueType Factorial(const ValueType & x, ErrorCode * err = 0)
stevep 0:04a9f72bbca7 2798 {
stevep 0:04a9f72bbca7 2799 return auxiliaryfunctions::Factorial2(x, (CGamma<ValueType>*)0, err, 0);
stevep 0:04a9f72bbca7 2800 }
stevep 0:04a9f72bbca7 2801
stevep 0:04a9f72bbca7 2802
stevep 0:04a9f72bbca7 2803 /*!
stevep 0:04a9f72bbca7 2804 this method prepares some coefficients: factorials and Bernoulli numbers
stevep 0:04a9f72bbca7 2805 stored in 'fact' and 'bern' objects
stevep 0:04a9f72bbca7 2806
stevep 0:04a9f72bbca7 2807 we're defining the method here because we're using Gamma() function which
stevep 0:04a9f72bbca7 2808 is not available in ttmathobjects.h
stevep 0:04a9f72bbca7 2809
stevep 0:04a9f72bbca7 2810 read the doc info in ttmathobjects.h file where CGamma<> struct is declared
stevep 0:04a9f72bbca7 2811 */
stevep 0:04a9f72bbca7 2812 template<class ValueType>
stevep 0:04a9f72bbca7 2813 void CGamma<ValueType>::InitAll()
stevep 0:04a9f72bbca7 2814 {
stevep 0:04a9f72bbca7 2815 ValueType x = TTMATH_GAMMA_BOUNDARY + 1;
stevep 0:04a9f72bbca7 2816
stevep 0:04a9f72bbca7 2817 // history.Remove(x) removes only one object
stevep 0:04a9f72bbca7 2818 // we must be sure that there are not others objects with the key 'x'
stevep 0:04a9f72bbca7 2819 while( history.Remove(x) )
stevep 0:04a9f72bbca7 2820 {
stevep 0:04a9f72bbca7 2821 }
stevep 0:04a9f72bbca7 2822
stevep 0:04a9f72bbca7 2823 // the simplest way to initialize is to call the Gamma function with (TTMATH_GAMMA_BOUNDARY + 1)
stevep 0:04a9f72bbca7 2824 // when x is larger then fewer coefficients we need
stevep 0:04a9f72bbca7 2825 Gamma(x, *this);
stevep 0:04a9f72bbca7 2826 }
stevep 0:04a9f72bbca7 2827
stevep 0:04a9f72bbca7 2828
stevep 0:04a9f72bbca7 2829
stevep 0:04a9f72bbca7 2830 } // namespace
stevep 0:04a9f72bbca7 2831
stevep 0:04a9f72bbca7 2832
stevep 0:04a9f72bbca7 2833 /*!
stevep 0:04a9f72bbca7 2834 this is for convenience for the user
stevep 0:04a9f72bbca7 2835 he can only use '#include <ttmath/ttmath.h>'
stevep 0:04a9f72bbca7 2836 */
stevep 0:04a9f72bbca7 2837 #include "ttmathparser.h"
stevep 0:04a9f72bbca7 2838
stevep 0:04a9f72bbca7 2839 // Dec is not finished yet
stevep 0:04a9f72bbca7 2840 //#include "ttmathdec.h"
stevep 0:04a9f72bbca7 2841
stevep 0:04a9f72bbca7 2842
stevep 0:04a9f72bbca7 2843
stevep 0:04a9f72bbca7 2844 #ifdef _MSC_VER
stevep 0:04a9f72bbca7 2845 //warning C4127: conditional expression is constant
stevep 0:04a9f72bbca7 2846 #pragma warning( default: 4127 )
stevep 0:04a9f72bbca7 2847 //warning C4702: unreachable code
stevep 0:04a9f72bbca7 2848 #pragma warning( default: 4702 )
stevep 0:04a9f72bbca7 2849 //warning C4800: forcing value to bool 'true' or 'false' (performance warning)
stevep 0:04a9f72bbca7 2850 #pragma warning( default: 4800 )
stevep 0:04a9f72bbca7 2851 #endif
stevep 0:04a9f72bbca7 2852
stevep 0:04a9f72bbca7 2853 #endif
stevep 0:04a9f72bbca7 2854