1 /* 2 3 PETSc mathematics include file. Defines certain basic mathematical 4 constants and functions for working with single, double, and quad precision 5 floating point numbers as well as complex single and double. 6 7 This file is included by petscsys.h and should not be used directly. 8 9 */ 10 11 #if !defined(__PETSCMATH_H) 12 #define __PETSCMATH_H 13 #include <math.h> 14 15 /* 16 17 Defines operations that are different for complex and real numbers; 18 note that one cannot mix the use of complex and real in the same 19 PETSc program. All PETSc objects in one program are built around the object 20 PetscScalar which is either always a real or a complex. 21 22 */ 23 24 #define PetscExpPassiveScalar(a) PetscExpScalar() 25 #if defined(PETSC_USE_REAL_SINGLE) 26 #define MPIU_REAL MPI_FLOAT 27 typedef float PetscReal; 28 #define PetscSqrtReal(a) sqrt(a) 29 #define PetscExpReal(a) exp(a) 30 #define PetscLogReal(a) log(a) 31 #define PetscLog10Real(a) log10(a) 32 #define PetscSinReal(a) sin(a) 33 #define PetscCosReal(a) cos(a) 34 #define PetscTanReal(a) tan(a) 35 #define PetscAsinReal(a) asin(a) 36 #define PetscAcosReal(a) acos(a) 37 #define PetscAtanReal(a) atan(a) 38 #define PetscAtan2Real(a,b) atan2(a,b) 39 #define PetscSinhReal(a) sinh(a) 40 #define PetscCoshReal(a) cosh(a) 41 #define PetscTanhReal(a) tanh(a) 42 #define PetscPowReal(a,b) pow(a,b) 43 #define PetscCeilReal(a) ceil(a) 44 #define PetscFloorReal(a) floor(a) 45 #define PetscFmodReal(a,b) fmod(a,b) 46 #define PetscTGamma(a) tgammaf(a) 47 #elif defined(PETSC_USE_REAL_DOUBLE) 48 #define MPIU_REAL MPI_DOUBLE 49 typedef double PetscReal; 50 #define PetscSqrtReal(a) sqrt(a) 51 #define PetscExpReal(a) exp(a) 52 #define PetscLogReal(a) log(a) 53 #define PetscLog10Real(a) log10(a) 54 #define PetscSinReal(a) sin(a) 55 #define PetscCosReal(a) cos(a) 56 #define PetscTanReal(a) tan(a) 57 #define PetscAsinReal(a) asin(a) 58 #define PetscAcosReal(a) acos(a) 59 #define PetscAtanReal(a) atan(a) 60 #define PetscAtan2Real(a,b) atan2(a,b) 61 #define PetscSinhReal(a) sinh(a) 62 #define PetscCoshReal(a) cosh(a) 63 #define PetscTanhReal(a) tanh(a) 64 #define PetscPowReal(a,b) pow(a,b) 65 #define PetscCeilReal(a) ceil(a) 66 #define PetscFloorReal(a) floor(a) 67 #define PetscFmodReal(a,b) fmod(a,b) 68 #define PetscTGamma(a) tgamma(a) 69 #elif defined(PETSC_USE_REAL___FLOAT128) 70 #if defined(__cplusplus) 71 extern "C" { 72 #endif 73 #include <quadmath.h> 74 #if defined(__cplusplus) 75 } 76 #endif 77 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128 PetscAttrMPITypeTag(__float128); 78 #define MPIU_REAL MPIU___FLOAT128 79 typedef __float128 PetscReal; 80 #define PetscSqrtReal(a) sqrtq(a) 81 #define PetscExpReal(a) expq(a) 82 #define PetscLogReal(a) logq(a) 83 #define PetscLog10Real(a) log10q(a) 84 #define PetscSinReal(a) sinq(a) 85 #define PetscCosReal(a) cosq(a) 86 #define PetscTanReal(a) tanq(a) 87 #define PetscAsinReal(a) asinq(a) 88 #define PetscAcosReal(a) acosq(a) 89 #define PetscAtanReal(a) atanq(a) 90 #define PetscAtan2Real(a,b) atan2q(a,b) 91 #define PetscSinhReal(a) sinhq(a) 92 #define PetscCoshReal(a) coshq(a) 93 #define PetscTanhReal(a) tanhq(a) 94 #define PetscPowReal(a,b) powq(a,b) 95 #define PetscCeilReal(a) ceilq(a) 96 #define PetscFloorReal(a) floorq(a) 97 #define PetscFmodReal(a,b) fmodq(a,b) 98 #define PetscTGamma(a) tgammaq(a) 99 #endif /* PETSC_USE_REAL_* */ 100 101 /* 102 Complex number definitions 103 */ 104 #if defined(__cplusplus) && defined(PETSC_HAVE_CXX_COMPLEX) && !defined(PETSC_USE_REAL___FLOAT128) 105 #if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX)) || defined(PETSC_DESIRE_COMPLEX) 106 #define PETSC_HAVE_COMPLEX 1 107 /* C++ support of complex number */ 108 #if defined(PETSC_HAVE_CUSP) 109 #define complexlib cusp 110 #include <cusp/complex.h> 111 #else 112 #define complexlib std 113 #include <complex> 114 #endif 115 116 #define PetscRealPartComplex(a) (a).real() 117 #define PetscImaginaryPartComplex(a) (a).imag() 118 #define PetscAbsComplex(a) complexlib::abs(a) 119 #define PetscConjComplex(a) complexlib::conj(a) 120 #define PetscSqrtComplex(a) complexlib::sqrt(a) 121 #define PetscPowComplex(a,b) complexlib::pow(a,b) 122 #define PetscExpComplex(a) complexlib::exp(a) 123 #define PetscLogComplex(a) complexlib::log(a) 124 #define PetscSinComplex(a) complexlib::sin(a) 125 #define PetscCosComplex(a) complexlib::cos(a) 126 #define PetscAsinComplex(a) complexlib::asin(a) 127 #define PetscAcosComplex(a) complexlib::acos(a) 128 #define PetscTanComplex(a) complexlib::tan(a) 129 #define PetscSinhComplex(a) complexlib::sinh(a) 130 #define PetscCoshComplex(a) complexlib::cosh(a) 131 #define PetscTanhComplex(a) complexlib::tanh(a) 132 133 #if defined(PETSC_USE_REAL_SINGLE) 134 typedef complexlib::complex<float> PetscComplex; 135 #elif defined(PETSC_USE_REAL_DOUBLE) 136 typedef complexlib::complex<double> PetscComplex; 137 #elif defined(PETSC_USE_REAL___FLOAT128) 138 typedef complexlib::complex<__float128> PetscComplex; /* Notstandard and not expected to work, use __complex128 */ 139 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128; 140 #endif /* PETSC_USE_REAL_ */ 141 #endif /* PETSC_USE_COMPLEX || PETSC_DESIRE_COMPLEX */ 142 143 #elif !defined(__cplusplus) && defined(PETSC_HAVE_C99_COMPLEX) 144 /* Use C99 _Complex for the type. Do not include complex.h by default to define "complex" because of symbol conflicts in Hypre. */ 145 /* Compilation units that can safely use complex should define PETSC_DESIRE_COMPLEX before including any headers */ 146 #if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX)) || defined(PETSC_DESIRE_COMPLEX) 147 #define PETSC_HAVE_COMPLEX 1 148 #include <complex.h> 149 150 #if defined(PETSC_USE_REAL_SINGLE) 151 typedef float _Complex PetscComplex; 152 153 #define PetscRealPartComplex(a) crealf(a) 154 #define PetscImaginaryPartComplex(a) cimagf(a) 155 #define PetscAbsComplex(a) cabsf(a) 156 #define PetscConjComplex(a) conjf(a) 157 #define PetscSqrtComplex(a) csqrtf(a) 158 #define PetscPowComplex(a,b) cpowf(a,b) 159 #define PetscExpComplex(a) cexpf(a) 160 #define PetscLogComplex(a) clogf(a) 161 #define PetscSinComplex(a) csinf(a) 162 #define PetscCosComplex(a) ccosf(a) 163 #define PetscAsinComplex(a) casinf(a) 164 #define PetscAcosComplex(a) cacosf(a) 165 #define PetscTanComplex(a) ctanf(a) 166 #define PetscSinhComplex(a) csinhf(a) 167 #define PetscCoshComplex(a) ccoshf(a) 168 #define PetscTanhComplex(a) ctanhf(a) 169 170 #elif defined(PETSC_USE_REAL_DOUBLE) 171 typedef double _Complex PetscComplex; 172 173 #define PetscRealPartComplex(a) creal(a) 174 #define PetscImaginaryPartComplex(a) cimag(a) 175 #define PetscAbsComplex(a) cabs(a) 176 #define PetscConjComplex(a) conj(a) 177 #define PetscSqrtComplex(a) csqrt(a) 178 #define PetscPowComplex(a,b) cpow(a,b) 179 #define PetscExpComplex(a) cexp(a) 180 #define PetscLogComplex(a) clog(a) 181 #define PetscSinComplex(a) csin(a) 182 #define PetscCosComplex(a) ccos(a) 183 #define PetscAsinComplex(a) casin(a) 184 #define PetscAcosComplex(a) cacos(a) 185 #define PetscTanComplex(a) ctan(a) 186 #define PetscSinhComplex(a) csinh(a) 187 #define PetscCoshComplex(a) ccosh(a) 188 #define PetscTanhComplex(a) ctanh(a) 189 190 #elif defined(PETSC_USE_REAL___FLOAT128) 191 typedef __complex128 PetscComplex; 192 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128 PetscAttrMPITypeTag(__complex128); 193 194 #define PetscRealPartComplex(a) crealq(a) 195 #define PetscImaginaryPartComplex(a) cimagq(a) 196 #define PetscAbsComplex(a) cabsq(a) 197 #define PetscConjComplex(a) conjq(a) 198 #define PetscSqrtComplex(a) csqrtq(a) 199 #define PetscPowComplex(a,b) cpowq(a,b) 200 #define PetscExpComplex(a) cexpq(a) 201 #define PetscLogComplex(a) clogq(a) 202 #define PetscSinComplex(a) csinq(a) 203 #define PetscCosComplex(a) ccosq(a) 204 #define PetscAsinComplex(a) casinq(a) 205 #define PetscAcosComplex(a) cacosq(a) 206 #define PetscTanComplex(a) ctanq(a) 207 #define PetscSinhComplex(a) csinhq(a) 208 #define PetscCoshComplex(a) ccoshq(a) 209 #define PetscTanhComplex(a) ctanhq(a) 210 211 #endif /* PETSC_USE_REAL_* */ 212 #elif (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX)) 213 #error "PETSc was configured --with-scalar-type=complex, but a language-appropriate complex library is not available" 214 #endif /* PETSC_USE_COMPLEX || PETSC_DESIRE_COMPLEX */ 215 #endif /* (__cplusplus && PETSC_HAVE_CXX_COMPLEX) else-if (!__cplusplus && PETSC_HAVE_C99_COMPLEX) */ 216 217 #if defined(PETSC_HAVE_COMPLEX) 218 #if defined(PETSC_HAVE_MPI_C_DOUBLE_COMPLEX) 219 #define MPIU_C_DOUBLE_COMPLEX MPI_C_DOUBLE_COMPLEX 220 #define MPIU_C_COMPLEX MPI_C_COMPLEX 221 #else 222 # if defined(__cplusplus) && defined(PETSC_HAVE_CXX_COMPLEX) 223 typedef complexlib::complex<double> petsc_mpiu_c_double_complex; 224 typedef complexlib::complex<float> petsc_mpiu_c_complex; 225 # elif !defined(__cplusplus) && defined(PETSC_HAVE_C99_COMPLEX) 226 typedef double _Complex petsc_mpiu_c_double_complex; 227 typedef float _Complex petsc_mpiu_c_complex; 228 # else 229 typedef struct {double real,imag;} petsc_mpiu_c_double_complex; 230 typedef struct {float real,imag;} petsc_mpiu_c_complex; 231 # endif 232 PETSC_EXTERN MPI_Datatype MPIU_C_DOUBLE_COMPLEX PetscAttrMPITypeTagLayoutCompatible(petsc_mpiu_c_double_complex); 233 PETSC_EXTERN MPI_Datatype MPIU_C_COMPLEX PetscAttrMPITypeTagLayoutCompatible(petsc_mpiu_c_complex); 234 #endif /* PETSC_HAVE_MPI_C_DOUBLE_COMPLEX */ 235 #endif /* PETSC_HAVE_COMPLEX */ 236 237 #if defined(PETSC_HAVE_COMPLEX) 238 # if defined(PETSC_USE_REAL_SINGLE) 239 # define MPIU_COMPLEX MPIU_C_COMPLEX 240 # elif defined(PETSC_USE_REAL_DOUBLE) 241 # define MPIU_COMPLEX MPIU_C_DOUBLE_COMPLEX 242 # elif defined(PETSC_USE_REAL___FLOAT128) 243 # define MPIU_COMPLEX MPIU___COMPLEX128 244 # endif /* PETSC_USE_REAL_* */ 245 #endif 246 247 #if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX)) 248 typedef PetscComplex PetscScalar; 249 #define PetscRealPart(a) PetscRealPartComplex(a) 250 #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a) 251 #define PetscAbsScalar(a) PetscAbsComplex(a) 252 #define PetscConj(a) PetscConjComplex(a) 253 #define PetscSqrtScalar(a) PetscSqrtComplex(a) 254 #define PetscPowScalar(a,b) PetscPowComplex(a,b) 255 #define PetscExpScalar(a) PetscExpComplex(a) 256 #define PetscLogScalar(a) PetscLogComplex(a) 257 #define PetscSinScalar(a) PetscSinComplex(a) 258 #define PetscCosScalar(a) PetscCosComplex(a) 259 #define PetscAsinScalar(a) PetscAsinComplex(a) 260 #define PetscAcosScalar(a) PetscAcosComplex(a) 261 #define PetscTanScalar(a) PetscTanComplex(a) 262 #define PetscSinhScalar(a) PetscSinhComplex(a) 263 #define PetscCoshScalar(a) PetscCoshComplex(a) 264 #define PetscTanhScalar(a) PetscTanhComplex(a) 265 #define MPIU_SCALAR MPIU_COMPLEX 266 267 /* 268 real number definitions 269 */ 270 #else /* PETSC_USE_COMPLEX */ 271 typedef PetscReal PetscScalar; 272 #define MPIU_SCALAR MPIU_REAL 273 274 #define PetscRealPart(a) (a) 275 #define PetscImaginaryPart(a) ((PetscReal)0.) 276 PETSC_STATIC_INLINE PetscReal PetscAbsScalar(PetscScalar a) {return a < 0.0 ? -a : a;} 277 #define PetscConj(a) (a) 278 #if !defined(PETSC_USE_REAL___FLOAT128) 279 #define PetscSqrtScalar(a) sqrt(a) 280 #define PetscPowScalar(a,b) pow(a,b) 281 #define PetscExpScalar(a) exp(a) 282 #define PetscLogScalar(a) log(a) 283 #define PetscSinScalar(a) sin(a) 284 #define PetscCosScalar(a) cos(a) 285 #define PetscAsinScalar(a) asin(a) 286 #define PetscAcosScalar(a) acos(a) 287 #define PetscTanScalar(a) tan(a) 288 #define PetscSinhScalar(a) sinh(a) 289 #define PetscCoshScalar(a) cosh(a) 290 #define PetscTanhScalar(a) tanh(a) 291 #else /* PETSC_USE_REAL___FLOAT128 */ 292 #define PetscSqrtScalar(a) sqrtq(a) 293 #define PetscPowScalar(a,b) powq(a,b) 294 #define PetscExpScalar(a) expq(a) 295 #define PetscLogScalar(a) logq(a) 296 #define PetscSinScalar(a) sinq(a) 297 #define PetscCosScalar(a) cosq(a) 298 #define PetscAsinScalar(a) asinq(a) 299 #define PetscAcosScalar(a) acosq(a) 300 #define PetscTanScalar(a) tanq(a) 301 #define PetscSinhScalar(a) sinhq(a) 302 #define PetscCoshScalar(a) coshq(a) 303 #define PetscTanhScalar(a) tanhq(a) 304 #endif /* PETSC_USE_REAL___FLOAT128 */ 305 306 #endif /* PETSC_USE_COMPLEX */ 307 308 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1) 309 #define PetscAbs(a) (((a) >= 0) ? (a) : -(a)) 310 311 /* --------------------------------------------------------------------------*/ 312 313 /* 314 Certain objects may be created using either single or double precision. 315 This is currently not used. 316 */ 317 typedef enum { PETSC_SCALAR_DOUBLE,PETSC_SCALAR_SINGLE, PETSC_SCALAR_LONG_DOUBLE } PetscScalarPrecision; 318 319 #if defined(PETSC_HAVE_COMPLEX) 320 /* PETSC_i is the imaginary number, i */ 321 PETSC_EXTERN PetscComplex PETSC_i; 322 #endif 323 324 /*MC 325 PetscMin - Returns minimum of two numbers 326 327 Synopsis: 328 #include <petscmath.h> 329 type PetscMin(type v1,type v2) 330 331 Not Collective 332 333 Input Parameter: 334 + v1 - first value to find minimum of 335 - v2 - second value to find minimum of 336 337 Notes: type can be integer or floating point value 338 339 Level: beginner 340 341 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 342 343 M*/ 344 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 345 346 /*MC 347 PetscMax - Returns maxium of two numbers 348 349 Synopsis: 350 #include <petscmath.h> 351 type max PetscMax(type v1,type v2) 352 353 Not Collective 354 355 Input Parameter: 356 + v1 - first value to find maximum of 357 - v2 - second value to find maximum of 358 359 Notes: type can be integer or floating point value 360 361 Level: beginner 362 363 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 364 365 M*/ 366 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 367 368 /*MC 369 PetscClipInterval - Returns a number clipped to be within an interval 370 371 Synopsis: 372 #include <petscmath.h> 373 type clip PetscClipInterval(type x,type a,type b) 374 375 Not Collective 376 377 Input Parameter: 378 + x - value to use if within interval (a,b) 379 . a - lower end of interval 380 - b - upper end of interval 381 382 Notes: type can be integer or floating point value 383 384 Level: beginner 385 386 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 387 388 M*/ 389 #define PetscClipInterval(x,a,b) (PetscMax((a),PetscMin((x),(b)))) 390 391 /*MC 392 PetscAbsInt - Returns the absolute value of an integer 393 394 Synopsis: 395 #include <petscmath.h> 396 int abs PetscAbsInt(int v1) 397 398 Not Collective 399 400 Input Parameter: 401 . v1 - the integer 402 403 Level: beginner 404 405 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 406 407 M*/ 408 #define PetscAbsInt(a) (((a)<0) ? -(a) : (a)) 409 410 /*MC 411 PetscAbsReal - Returns the absolute value of an real number 412 413 Synopsis: 414 #include <petscmath.h> 415 Real abs PetscAbsReal(PetscReal v1) 416 417 Not Collective 418 419 Input Parameter: 420 . v1 - the double 421 422 423 Level: beginner 424 425 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 426 427 M*/ 428 #define PetscAbsReal(a) (((a)<0) ? -(a) : (a)) 429 430 /*MC 431 PetscSqr - Returns the square of a number 432 433 Synopsis: 434 #include <petscmath.h> 435 type sqr PetscSqr(type v1) 436 437 Not Collective 438 439 Input Parameter: 440 . v1 - the value 441 442 Notes: type can be integer or floating point value 443 444 Level: beginner 445 446 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 447 448 M*/ 449 #define PetscSqr(a) ((a)*(a)) 450 451 /* ----------------------------------------------------------------------------*/ 452 /* 453 Basic constants 454 */ 455 #if defined(PETSC_USE_REAL___FLOAT128) 456 #define PETSC_PI M_PIq 457 #elif defined(M_PI) 458 #define PETSC_PI M_PI 459 #else 460 #define PETSC_PI 3.14159265358979323846264338327950288419716939937510582 461 #endif 462 463 #if !defined(PETSC_USE_64BIT_INDICES) 464 #define PETSC_MAX_INT 2147483647 465 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 466 #else 467 #define PETSC_MAX_INT 9223372036854775807L 468 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 469 #endif 470 471 #if defined(PETSC_USE_REAL_SINGLE) 472 # define PETSC_MAX_REAL 3.40282346638528860e+38F 473 # define PETSC_MIN_REAL -PETSC_MAX_REAL 474 # define PETSC_MACHINE_EPSILON 1.19209290e-07F 475 # define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F 476 # define PETSC_SMALL 1.e-5 477 #elif defined(PETSC_USE_REAL_DOUBLE) 478 # define PETSC_MAX_REAL 1.7976931348623157e+308 479 # define PETSC_MIN_REAL -PETSC_MAX_REAL 480 # define PETSC_MACHINE_EPSILON 2.2204460492503131e-16 481 # define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08 482 # define PETSC_SMALL 1.e-10 483 #elif defined(PETSC_USE_REAL___FLOAT128) 484 # define PETSC_MAX_REAL FLT128_MAX 485 # define PETSC_MIN_REAL -FLT128_MAX 486 # define PETSC_MACHINE_EPSILON FLT128_EPSILON 487 # define PETSC_SQRT_MACHINE_EPSILON 1.38777878078e-17 488 # define PETSC_SMALL 1.e-20 489 #endif 490 491 #define PETSC_INFINITY PETSC_MAX_REAL/4.0 492 #define PETSC_NINFINITY -PETSC_INFINITY 493 494 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanScalar(PetscScalar); 495 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanReal(PetscReal); 496 PETSC_EXTERN PetscBool PetscIsNormalScalar(PetscScalar); 497 PETSC_EXTERN PetscBool PetscIsNormalReal(PetscReal); 498 499 /* ----------------------------------------------------------------------------*/ 500 #define PassiveReal PetscReal 501 #define PassiveScalar PetscScalar 502 503 /* 504 These macros are currently hardwired to match the regular data types, so there is no support for a different 505 MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again. 506 */ 507 #define MPIU_MATSCALAR MPIU_SCALAR 508 typedef PetscScalar MatScalar; 509 typedef PetscReal MatReal; 510 511 struct petsc_mpiu_2scalar {PetscScalar a,b;}; 512 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2scalar); 513 #if defined(PETSC_USE_64BIT_INDICES) || !defined(MPI_2INT) 514 struct petsc_mpiu_2int {PetscInt a,b;}; 515 PETSC_EXTERN MPI_Datatype MPIU_2INT PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2int); 516 #else 517 #define MPIU_2INT MPI_2INT 518 #endif 519 520 PETSC_STATIC_INLINE PetscInt PetscPowInt(PetscInt base,PetscInt power) 521 { 522 PetscInt result = 1; 523 while (power) { 524 if (power & 1) result *= base; 525 power >>= 1; 526 base *= base; 527 } 528 return result; 529 } 530 531 PETSC_STATIC_INLINE PetscReal PetscPowRealInt(PetscReal base,PetscInt power) 532 { 533 PetscReal result = 1; 534 if (power < 0) { 535 power = -power; 536 if (base != 0.0) base = 1./base; 537 } 538 while (power) { 539 if (power & 1) result *= base; 540 power >>= 1; 541 base *= base; 542 } 543 return result; 544 } 545 546 PETSC_STATIC_INLINE PetscScalar PetscPowScalarInt(PetscScalar base,PetscInt power) 547 { 548 PetscScalar result = 1; 549 if (power < 0) { 550 power = -power; 551 if (base != 0.0) base = 1./base; 552 } 553 while (power) { 554 if (power & 1) result *= base; 555 power >>= 1; 556 base *= base; 557 } 558 return result; 559 } 560 561 PETSC_STATIC_INLINE PetscScalar PetscPowScalarReal(PetscScalar base,PetscReal power) 562 { 563 PetscScalar cpower = power; 564 return PetscPowScalar(base,cpower); 565 } 566 #endif 567