-
Notifications
You must be signed in to change notification settings - Fork 2
/
ladiv.h
72 lines (64 loc) · 1.67 KB
/
ladiv.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
//
// ladiv.h
// Linear Algebra Template Library
//
// Created by Rodney James on 3/12/12.
// Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//
#ifndef _ladiv_h
#define _ladiv_h
/// @file ladiv.h Performs complex division of two scalars taking care to avoid underflow.
#include <cmath>
#include "latl.h"
namespace LATL
{
/// @brief Performs complex division in real arithmetic.
///
/// p + iq = (a +ib) / (c + id)
///
/// @tparam real_t Floating point type.
/// @param a Real part of numerator.
/// @param b Imaginary part of numerator.
/// @param c Real part of denominator.
/// @param d Imaginary part of denominator.
/// @param p Real part of quotient.
/// @param q Imaginary part of quotient.
/// @ingroup AUX
template<typename real_t>
void LADIV(real_t a, real_t b, real_t c, real_t d, real_t &p, real_t &q)
{
using std::abs;
real_t e,f;
if(abs(d)<abs(c))
{
e=d/c;
f=c+d*e;
p=(a+b*e)/f;
q=(b-a*e)/f;
}
else
{
e=c/d;
f=c+d*e;
p=(b+a*e)/f;
q=(-a+b*e)/f;
}
}
/// @brief Performs complex division in real arithmetic with complex arguments.
/// @return x/y
/// @tparam real_t Floating point type.
/// @param x Complex numerator.
/// @param y Complex denominator.
/// @ingroup AUX
template<typename real_t>
complex<real_t> LADIV(complex<real_t> x, complex<real_t> y)
{
using std::abs;
using std::real;
using std::imag;
real_t zr,zi;
LADIV<real_t>(real(x),imag(x),real(y),imag(y),zr,zi);
return complex<real_t>(zr,zi);
}
}
#endif