Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64044828
polynomial.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, May 24, 04:44
Size
11 KB
Mime Type
text/x-c++
Expires
Sun, May 26, 04:44 (2 d)
Engine
blob
Format
Raw Data
Handle
17848578
Attached To
rTAMAAS tamaas
polynomial.hh
View Options
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Expolit is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Expolit is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
* more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Expolit. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef POLYNOMIALS_HH
#define POLYNOMIALS_HH
#include "operations.hh"
#include "types.hh"
#include <algorithm>
#include <array>
#include <numeric>
#include <ostream>
namespace
expolit
{
template
<
typename
From
,
typename
To
>
using
compatible
=
std
::
enable_if_t
<
std
::
is_convertible
<
From
,
To
>::
value
>
;
/**
* @brief Polynomial type
*
* Type for polynomial expressions of arbitrary order, variable type and
* coefficient type
*/
template
<
typename
CoeffType
,
UInt
order
>
struct
Polynomial
:
Expression
<
Polynomial
<
CoeffType
,
order
>>
{
/// Number of terms
constexpr
static
std
::
size_t
terms
=
order
+
1
;
constexpr
Polynomial
()
=
default
;
constexpr
explicit
Polynomial
(
const
std
::
array
<
CoeffType
,
terms
>&
coeffs
)
:
coeffs
(
coeffs
)
{}
/// Evaluating polynomial
template
<
typename
Val
,
typename
=
compatible
<
Val
,
CoeffType
>>
constexpr
auto
operator
()(
Val
z
)
const
{
using
Res
=
decltype
(
std
::
declval
<
CoeffType
>
()
*
z
);
Res
res
(
coeffs
.
front
());
Val
power
(
z
);
for
(
UInt
i
=
1
;
i
<
terms
;
++
i
)
{
res
+=
coeffs
[
i
]
*
power
;
power
*=
z
;
}
return
res
;
}
/// Storing coefficients
std
::
array
<
CoeffType
,
terms
>
coeffs
=
{
0
};
};
template
<
typename
CoeffType
>
using
Constant
=
Polynomial
<
CoeffType
,
0
>
;
namespace
detail
{
template
<
typename
Poly
>
struct
get_terms
:
std
::
integral_constant
<
std
::
size_t
,
std
::
remove_reference_t
<
std
::
remove_cv_t
<
Poly
>>::
terms
>
{
};
}
// namespace detail
/* -------------------------------------------------------------------------- */
/* Algebraic rules */
/* -------------------------------------------------------------------------- */
/// Sum with scalar
template
<
typename
CT
,
UInt
order
,
typename
DT
,
typename
=
compatible
<
DT
,
CT
>>
constexpr
auto
operator
+
(
DT
x
,
const
Polynomial
<
CT
,
order
>&
p
)
{
Polynomial
<
decltype
(
std
::
declval
<
CT
>
()
+
x
),
order
>
res
(
p
);
std
::
get
<
0
>
(
res
.
coeffs
)
+=
x
;
return
res
;
}
/// Sum with scalar
template
<
typename
CT
,
UInt
order
,
typename
DT
,
typename
=
compatible
<
DT
,
CT
>>
constexpr
auto
operator
+
(
const
Polynomial
<
CT
,
order
>&
p
,
DT
x
)
{
return
x
+
p
;
}
namespace
detail
{
template
<
std
::
size_t
i
>
struct
for_scalar
{
template
<
typename
CT
,
std
::
size_t
n
,
typename
DT
>
constexpr
static
void
prod
(
std
::
array
<
CT
,
n
>&
res
,
DT
x
)
{
std
::
get
<
i
>
(
res
)
*=
x
;
for_scalar
<
i
-
1
>::
prod
(
res
,
x
);
}
template
<
typename
CT
,
std
::
size_t
n
,
typename
DT
>
constexpr
static
void
div
(
std
::
array
<
CT
,
n
>&
res
,
DT
x
)
{
std
::
get
<
i
>
(
res
)
/=
x
;
for_scalar
<
i
-
1
>::
prod
(
res
,
x
);
}
};
template
<>
struct
for_scalar
<
0
>
{
template
<
typename
CT
,
std
::
size_t
n
,
typename
DT
>
constexpr
static
void
prod
(
std
::
array
<
CT
,
n
>&
res
,
DT
x
)
{
std
::
get
<
0
>
(
res
)
*=
x
;
}
template
<
typename
CT
,
std
::
size_t
n
,
typename
DT
>
constexpr
static
void
div
(
std
::
array
<
CT
,
n
>&
res
,
DT
x
)
{
std
::
get
<
0
>
(
res
)
/=
x
;
}
};
}
// namespace detail
/// Product with scalar
template
<
typename
CT
,
UInt
order
,
typename
DT
,
typename
=
compatible
<
DT
,
CT
>>
constexpr
auto
operator
*
(
DT
x
,
const
Polynomial
<
CT
,
order
>&
p
)
{
Polynomial
<
decltype
(
std
::
declval
<
CT
>
()
*
x
),
order
>
res
(
p
);
detail
::
for_scalar
<
res
.
terms
-
1
>::
prod
(
res
.
coeffs
,
x
);
return
res
;
}
/// Product with scalar
template
<
typename
CT
,
UInt
order
,
typename
DT
,
typename
=
compatible
<
DT
,
CT
>>
constexpr
auto
operator
*
(
const
Polynomial
<
CT
,
order
>&
p
,
DT
x
)
{
return
x
*
p
;
}
/// Division with scalar
template
<
typename
CT
,
UInt
order
,
typename
DT
,
typename
=
compatible
<
DT
,
CT
>>
constexpr
auto
operator
/
(
const
Polynomial
<
CT
,
order
>&
p
,
DT
x
)
{
Polynomial
<
decltype
(
std
::
declval
<
CT
>
()
/
x
),
order
>
res
(
p
);
detail
::
for_scalar
<
res
.
terms
-
1
>::
div
(
res
.
coeffs
,
x
);
return
res
;
}
/// Division with constant
template
<
typename
CT
,
UInt
order
,
typename
CT2
,
typename
=
compatible
<
CT2
,
CT
>>
constexpr
auto
operator
/
(
const
Polynomial
<
CT
,
order
>&
p
,
const
Constant
<
CT2
>&
x
)
{
Polynomial
<
decltype
(
std
::
declval
<
CT
>
()
/
x
.
coeffs
.
front
()),
order
>
res
(
p
);
detail
::
for_scalar
<
res
.
terms
-
1
>::
div
(
res
.
coeffs
,
std
::
get
<
0
>
(
x
.
coeffs
));
return
res
;
}
namespace
detail
{
template
<
std
::
size_t
i
,
std
::
size_t
j
>
struct
coeff_product
{
template
<
typename
CT1
,
typename
CT2
,
typename
CT3
,
std
::
size_t
n
,
std
::
size_t
m
,
std
::
size_t
l
>
static
constexpr
void
prod
(
std
::
array
<
CT1
,
n
>&
res
,
const
std
::
array
<
CT2
,
m
>&
a
,
const
std
::
array
<
CT3
,
l
>&
b
)
{
std
::
get
<
i
+
j
>
(
res
)
+=
a
[
i
]
*
b
[
j
];
coeff_product
<
i
,
j
-
1
>::
prod
(
res
,
a
,
b
);
}
};
template
<
std
::
size_t
i
>
struct
coeff_product
<
i
,
0
>
{
template
<
typename
CT1
,
typename
CT2
,
typename
CT3
,
std
::
size_t
n
,
std
::
size_t
m
,
std
::
size_t
l
>
static
constexpr
void
prod
(
std
::
array
<
CT1
,
n
>&
res
,
const
std
::
array
<
CT2
,
m
>&
a
,
const
std
::
array
<
CT3
,
l
>&
b
)
{
std
::
get
<
i
>
(
res
)
+=
a
[
i
]
*
b
[
0
];
coeff_product
<
i
-
1
,
l
-
1
>::
prod
(
res
,
a
,
b
);
}
};
template
<>
struct
coeff_product
<
0
,
0
>
{
template
<
typename
CT1
,
typename
CT2
,
typename
CT3
,
std
::
size_t
n
,
std
::
size_t
m
,
std
::
size_t
l
>
static
constexpr
void
prod
(
std
::
array
<
CT1
,
n
>&
res
,
const
std
::
array
<
CT2
,
m
>&
a
,
const
std
::
array
<
CT3
,
l
>&
b
)
{
std
::
get
<
0
>
(
res
)
+=
a
[
0
]
*
b
[
0
];
}
};
}
// namespace detail
/// Product of polynomials
template
<
typename
CT1
,
typename
CT2
,
UInt
o1
,
UInt
o2
>
constexpr
auto
operator
*
(
const
Polynomial
<
CT1
,
o1
>&
p1
,
const
Polynomial
<
CT2
,
o2
>&
p2
)
{
Polynomial
<
decltype
(
std
::
declval
<
CT1
>
()
*
std
::
declval
<
CT2
>
()),
o1
+
o2
>
res
;
detail
::
coeff_product
<
detail
::
get_terms
<
decltype
(
p1
)
>::
value
-
1
,
detail
::
get_terms
<
decltype
(
p2
)
>::
value
-
1
>::
prod
(
res
.
coeffs
,
p1
.
coeffs
,
p2
.
coeffs
);
return
res
;
}
namespace
detail
{
template
<
std
::
size_t
i
>
struct
coeff_sum
{
template
<
typename
CT1
,
typename
CT2
,
typename
CT3
,
std
::
size_t
n
,
std
::
size_t
m
,
std
::
size_t
l
>
static
constexpr
void
sum
(
std
::
array
<
CT1
,
n
>&
res
,
const
std
::
array
<
CT2
,
m
>&
a
,
const
std
::
array
<
CT3
,
l
>&
b
)
{
std
::
get
<
i
>
(
res
)
+=
((
i
<
m
)
?
a
[
i
]
:
0
)
+
((
i
<
l
)
?
b
[
i
]
:
0
);
coeff_sum
<
i
-
1
>::
sum
(
res
,
a
,
b
);
}
};
template
<>
struct
coeff_sum
<
0
>
{
template
<
typename
CT1
,
typename
CT2
,
typename
CT3
,
std
::
size_t
n
,
std
::
size_t
m
,
std
::
size_t
l
>
static
constexpr
void
sum
(
std
::
array
<
CT1
,
n
>&
res
,
const
std
::
array
<
CT2
,
m
>&
a
,
const
std
::
array
<
CT3
,
l
>&
b
)
{
std
::
get
<
0
>
(
res
)
+=
a
[
0
]
+
b
[
0
];
}
};
}
// namespace detail
/// Sum of polynomials
template
<
typename
CT1
,
typename
CT2
,
UInt
o1
,
UInt
o2
>
constexpr
auto
operator
+
(
const
Polynomial
<
CT1
,
o1
>&
p1
,
const
Polynomial
<
CT2
,
o2
>&
p2
)
{
Polynomial
<
decltype
(
std
::
declval
<
CT1
>
()
+
std
::
declval
<
CT2
>
()),
std
::
max
(
o1
,
o2
)
>
res
;
detail
::
coeff_sum
<
res
.
terms
-
1
>::
sum
(
res
.
coeffs
,
p1
.
coeffs
,
p2
.
coeffs
);
return
res
;
}
/* -------------------------------------------------------------------------- */
/* Simplification rules */
/* -------------------------------------------------------------------------- */
#define POLY_SIMPLIFICATION(opsign, std_op) \
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const Polynomial<CT1, o1>& p1, \
const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) { \
return p2.op(p2.op(p1, p2.operands.first), p2.operands.second); \
} \
\
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2, \
const Polynomial<CT1, o1>& p1) { \
return p2.op(p1, p2); \
} \
\
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const Polynomial<CT1, o1>& p1, \
const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) { \
return p2.op(p1, p2.commute()); \
} \
\
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2, \
const Polynomial<CT1, o1>& p1) { \
return p2.op(p1, p2.commute()); \
}
POLY_SIMPLIFICATION
(
+
,
std
::
plus
<>
)
POLY_SIMPLIFICATION
(
*
,
std
::
multiplies
<>
)
#undef POLY_SIMPLIFICATION
template
<
typename
CT1
,
UInt
o
,
typename
CT2
,
typename
Derived
>
constexpr
auto
operator
/
(
const
Product
<
Polynomial
<
CT1
,
o
>
,
Derived
>&
p
,
const
Constant
<
CT2
>&
c
)
{
return
p
.
operands
.
second
*
(
p
.
operands
.
first
/
c
);
}
template
<
typename
CT1
,
UInt
o
,
typename
CT2
,
typename
Derived
>
constexpr
auto
operator
/
(
const
Product
<
Derived
,
Polynomial
<
CT1
,
o
>>&
p
,
const
Constant
<
CT2
>&
c
)
{
return
p
.
operands
.
second
*
(
p
.
operands
.
first
/
c
);
}
template
<
typename
CT1
,
UInt
o
,
typename
CT2
,
typename
Derived
>
constexpr
auto
operator
*
(
const
Polynomial
<
CT1
,
o
>&
p
,
const
Division
<
Derived
,
Constant
<
CT2
>>&
d
)
{
return
(
p
/
d
.
operands
.
second
)
*
d
.
operands
.
first
;
}
template
<
typename
CT1
,
UInt
o
,
typename
CT2
,
typename
Derived
>
constexpr
auto
operator
*
(
const
Division
<
Derived
,
Constant
<
CT2
>>&
d
,
const
Polynomial
<
CT1
,
o
>&
p
)
{
return
(
p
/
d
.
operands
.
second
)
*
d
.
operands
.
first
;
}
/* -------------------------------------------------------------------------- */
/* Output */
/* -------------------------------------------------------------------------- */
template
<
typename
CT
,
UInt
order
>
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
Polynomial
<
CT
,
order
>&
p
)
{
os
<<
p
.
coeffs
[
0
];
if
(
p
.
terms
>=
2
)
{
os
<<
" + "
<<
p
.
coeffs
[
1
]
<<
"*z"
;
for
(
UInt
i
=
2
;
i
<
p
.
terms
;
++
i
)
os
<<
" + "
<<
p
.
coeffs
[
i
]
<<
"*z^"
<<
i
;
}
return
os
;
}
}
// namespace expolit
#endif
Event Timeline
Log In to Comment