Skip to content

Commit

Permalink
Improve accuracy of complex tan/tanh.
Browse files Browse the repository at this point in the history 
Test self assignment when using complex_adaptor and fix resulting bugs.
Fixes #262.
  • Loading branch information
@jzmaddock
jzmaddock committed Aug 6, 2020
1 parent 34881d5 commit ed66aa5
Show file tree
Hide file tree
Showing 2 changed files with 128 additions and 48 deletions.
98 changes: 76 additions & 22 deletions include/boost/multiprecision/complex_adaptor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -550,14 +550,15 @@ inline void eval_sin(complex_adaptor<Backend>& result, const complex_adaptor<Bac
using default_ops::eval_sin;
using default_ops::eval_sinh;

Backend t1, t2;
Backend t1, t2, t3;
eval_sin(t1, arg.real_data());
eval_cosh(t2, arg.imag_data());
eval_multiply(result.real_data(), t1, t2);
eval_multiply(t3, t1, t2);

eval_cos(t1, arg.real_data());
eval_sinh(t2, arg.imag_data());
eval_multiply(result.imag_data(), t1, t2);
result.real_data() = t3;
}

template <class Backend>
Expand All @@ -568,24 +569,85 @@ inline void eval_cos(complex_adaptor<Backend>& result, const complex_adaptor<Bac
using default_ops::eval_sin;
using default_ops::eval_sinh;

Backend t1, t2;
Backend t1, t2, t3;
eval_cos(t1, arg.real_data());
eval_cosh(t2, arg.imag_data());
eval_multiply(result.real_data(), t1, t2);
eval_multiply(t3, t1, t2);

eval_sin(t1, arg.real_data());
eval_sinh(t2, arg.imag_data());
eval_multiply(result.imag_data(), t1, t2);
result.imag_data().negate();
result.real_data() = t3;
}

template <class T>
void tanh_imp(const T& r, const T& i, T& r_result, T& i_result)
{
using default_ops::eval_tan;
using default_ops::eval_sinh;
using default_ops::eval_add;
using default_ops::eval_fpclassify;
using default_ops::eval_get_sign;

typedef typename boost::mpl::front<typename T::unsigned_types>::type ui_type;
ui_type one(1);
//
// Set:
// t = tan(i);
// s = sinh(r);
// b = s * (1 + t^2);
// d = 1 + b * s;
//
T t, s, b, d;
eval_tan(t, i);
eval_sinh(s, r);
eval_multiply(d, t, t);
eval_add(d, one);
eval_multiply(b, d, s);
eval_multiply(d, b, s);
eval_add(d, one);

if (eval_fpclassify(d) == FP_INFINITE)
{
r_result = one;
if (eval_get_sign(s) < 0)
r_result.negate();
//
// Imaginary part is a signed zero:
//
ui_type zero(0);
i_result = zero;
if (eval_get_sign(t) < 0)
i_result.negate();
}
//
// Real part is sqrt(1 + s^2) * b / d;
// Imaginary part is t / d;
//
eval_divide(i_result, t, d);
//
// variable t is now spare, as is r_result.
//
eval_multiply(t, s, s);
eval_add(t, one);
eval_sqrt(r_result, t);
eval_multiply(t, r_result, b);
eval_divide(r_result, t, d);
}

template <class Backend>
inline void eval_tanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
{
tanh_imp(arg.real_data(), arg.imag_data(), result.real_data(), result.imag_data());
}
template <class Backend>
inline void eval_tan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
{
complex_adaptor<Backend> c;
eval_cos(c, arg);
eval_sin(result, arg);
eval_divide(result, c);
Backend t(arg.imag_data());
t.negate();
tanh_imp(t, arg.real_data(), result.imag_data(), result.real_data());
result.imag_data().negate();
}

template <class Backend>
Expand Down Expand Up @@ -659,14 +721,15 @@ inline void eval_sinh(complex_adaptor<Backend>& result, const complex_adaptor<Ba
using default_ops::eval_sin;
using default_ops::eval_sinh;

Backend t1, t2;
Backend t1, t2, t3;
eval_cos(t1, arg.imag_data());
eval_sinh(t2, arg.real_data());
eval_multiply(result.real_data(), t1, t2);
eval_multiply(t3, t1, t2);

eval_cosh(t1, arg.real_data());
eval_sin(t2, arg.imag_data());
eval_multiply(result.imag_data(), t1, t2);
result.real_data() = t3;
}

template <class Backend>
Expand All @@ -677,24 +740,15 @@ inline void eval_cosh(complex_adaptor<Backend>& result, const complex_adaptor<Ba
using default_ops::eval_sin;
using default_ops::eval_sinh;

Backend t1, t2;
Backend t1, t2, t3;
eval_cos(t1, arg.imag_data());
eval_cosh(t2, arg.real_data());
eval_multiply(result.real_data(), t1, t2);
eval_multiply(t3, t1, t2);

eval_sin(t1, arg.imag_data());
eval_sinh(t2, arg.real_data());
eval_multiply(result.imag_data(), t1, t2);
}

template <class Backend>
inline void eval_tanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
{
using default_ops::eval_divide;
complex_adaptor<Backend> s, c;
eval_sinh(s, arg);
eval_cosh(c, arg);
eval_divide(result, s, c);
result.real_data() = t3;
}

template <class Backend>
Expand Down
78 changes: 52 additions & 26 deletions test/test_arithmetic.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2386,83 +2386,109 @@ typename boost::enable_if_c<boost::multiprecision::number_category<Real>::value
BOOST_CHECK_CLOSE_FRACTION(real_type("-8.8931513442797186948734782808862447235385767991868219480917324534839621090167050538805196124711247247992169338"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-1.3826999557878897572499699021550296885662132089951379549068064961882821777067532977546360861176011175070188118"), imag(a), tol * 3);

a = sqrt(c);
a = c;
a = sqrt(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.674149228035540040448039300849051821674708677883920366727287836003399240343274891876712629708287692163156802065"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.8959774761298381247157337552900434410433241995549314932449006989874470582160955817053273057885402621549320588976"), imag(a), tol);
a = sin(c);
a = c;
a = sin(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("9.154499146911429573467299544609832559158860568765182977899828142590020335321896403936690014669532606510294425039"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-4.168906959966564350754813058853754843573565604758055889965478710592666260138453299795649308385497563475115931624"), imag(a), tol);
a = cos(c);
a = c;
a = cos(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("-4.1896256909688072301325550196159737286219454041279210357407905058369727912162626993926269783331491034500484583"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-9.1092278937553365979791972627788621213326202389201695649104967309554222940748568716960841549279996556547993373"), imag(a), tol);
a = tan(c);
a = c;
a = tan(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("-0.0037640256415042482927512211303226908396306202016580864328644932511249097100916559688254811519914564480500042311"), real(a), tol * 5);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.0032386273536098014463585978219272598077897241071003399272426939850671219193120708438426543945017427085738411"), imag(a), tol);
a = asin(c);
a = c;
a = asin(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.5706527843210994007102838796856696501828032450960401365302732598209740064262509342420347149436326252483895113827"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.983387029916535432347076902894039565014248302909345356125267430944752731616095111727103650117987412058949254132"), imag(a), tol);
a = acos(c);
a = c;
a = acos(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.000143542473797218521037811954081791915781454591512773957199036332934196716853565071982697727425908742684531873"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-1.983387029916535432347076902894039565014248302909345356125267430944752731616095111727103650117987412058949254132"), imag(a), tol);
a = atan(c);
a = c;
a = atan(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.409921049596575522530619384460420782588207051908724814771070766475530084440199227135813201495737846771570458568"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.2290726829685387662958818029420027678625253049770656169479919704951963414344907622560676377741902308144912055002"), imag(a), tol);
a = sqrt(c + 0);
a = c;
a = sqrt(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.674149228035540040448039300849051821674708677883920366727287836003399240343274891876712629708287692163156802065"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.8959774761298381247157337552900434410433241995549314932449006989874470582160955817053273057885402621549320588976"), imag(a), tol);
a = sin(c + 0);
a = c;
a = sin(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("9.154499146911429573467299544609832559158860568765182977899828142590020335321896403936690014669532606510294425039"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-4.168906959966564350754813058853754843573565604758055889965478710592666260138453299795649308385497563475115931624"), imag(a), tol);
a = cos(c + 0);
a = c;
a = cos(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("-4.1896256909688072301325550196159737286219454041279210357407905058369727912162626993926269783331491034500484583"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-9.1092278937553365979791972627788621213326202389201695649104967309554222940748568716960841549279996556547993373"), imag(a), tol);
a = tan(c + 0);
a = c;
a = tan(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("-0.0037640256415042482927512211303226908396306202016580864328644932511249097100916559688254811519914564480500042311"), real(a), tol * 5);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.0032386273536098014463585978219272598077897241071003399272426939850671219193120708438426543945017427085738411"), imag(a), tol);
a = asin(c + 0);
a = c;
a = asin(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.5706527843210994007102838796856696501828032450960401365302732598209740064262509342420347149436326252483895113827"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.983387029916535432347076902894039565014248302909345356125267430944752731616095111727103650117987412058949254132"), imag(a), tol);
a = acos(c + 0);
a = c;
a = acos(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.000143542473797218521037811954081791915781454591512773957199036332934196716853565071982697727425908742684531873"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-1.983387029916535432347076902894039565014248302909345356125267430944752731616095111727103650117987412058949254132"), imag(a), tol);
a = atan(c + 0);
a = c;
a = atan(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.409921049596575522530619384460420782588207051908724814771070766475530084440199227135813201495737846771570458568"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.2290726829685387662958818029420027678625253049770656169479919704951963414344907622560676377741902308144912055002"), imag(a), tol);

a = sinh(c);
a = c;
a = sinh(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("-3.5905645899857799520125654477948167931949136757293015099986213974178826801534614215227593814301490087307920223"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.53092108624851980526704009066067655967277345095149103008706855371803528753067068552935673000832252607835087747"), imag(a), tol);
a = cosh(c);
a = c;
a = cosh(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("-3.7245455049153225654739707032559725286749657732153307267858945686649501059065292889110148294141744084833329553"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.51182256998738460883446384980187563424555660949074386745538379123585339045741119409984041226187262097496424111"), imag(a), tol);
a = tanh(c);
a = c;
a = tanh(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.965385879022133124278480269394560685879729650005757773636908240066639772853967550095754361348005358178253777920"), real(a), tol * 5);
BOOST_CHECK_CLOSE_FRACTION(real_type("-0.00988437503832249372031403430350121097961813353467039031861010606115560355679254344335582852193041894874685555114"), imag(a), tol);
a = asinh(c);
a = c;
a = asinh(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076230641308014923"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.9646585044076027920454110594995323555197773725073316527132580297155508786089335572049608301897631767195194427315"), imag(a), tol);
a = acosh(c);
a = c;
a = acosh(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.983387029916535432347076902894039565014248302909345356125267430944752731616095111727103650117987412058949254132"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.000143542473797218521037811954081791915781454591512773957199036332934196716853565071982697727425908742684531873"), imag(a), tol);
a = atanh(c);
a = c;
a = atanh(a);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.1469466662255297520474327851547159424423449403442452953891851939502023996823900422792744078835711416939934387775"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.338972522294493561124193575909144241084316172544492778582005751793809271060233646663717270678614587712809117131"), imag(a), tol);
a = sinh(c + 0);
a = c;
a = sinh(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("-3.5905645899857799520125654477948167931949136757293015099986213974178826801534614215227593814301490087307920223"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.53092108624851980526704009066067655967277345095149103008706855371803528753067068552935673000832252607835087747"), imag(a), tol);
a = cosh(c + 0);
a = c;
a = cosh(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("-3.7245455049153225654739707032559725286749657732153307267858945686649501059065292889110148294141744084833329553"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.51182256998738460883446384980187563424555660949074386745538379123585339045741119409984041226187262097496424111"), imag(a), tol);
a = tanh(c + 0);
a = c;
a = tanh(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.965385879022133124278480269394560685879729650005757773636908240066639772853967550095754361348005358178253777920"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("-0.00988437503832249372031403430350121097961813353467039031861010606115560355679254344335582852193041894874685555114"), imag(a), tol);
a = asinh(c + 0);
a = c;
a = asinh(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076230641308014923"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.9646585044076027920454110594995323555197773725073316527132580297155508786089335572049608301897631767195194427315"), imag(a), tol);
a = acosh(c + 0);
a = c;
a = acosh(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.983387029916535432347076902894039565014248302909345356125267430944752731616095111727103650117987412058949254132"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.000143542473797218521037811954081791915781454591512773957199036332934196716853565071982697727425908742684531873"), imag(a), tol);
a = atanh(c + 0);
a = c;
a = atanh(a + 0);
BOOST_CHECK_CLOSE_FRACTION(real_type("0.1469466662255297520474327851547159424423449403442452953891851939502023996823900422792744078835711416939934387775"), real(a), tol);
BOOST_CHECK_CLOSE_FRACTION(real_type("1.338972522294493561124193575909144241084316172544492778582005751793809271060233646663717270678614587712809117131"), imag(a), tol);
}
Expand Down

0 comments on commit ed66aa5

Please sign in to comment.