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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Defining a Special Function.</title>
<link rel="stylesheet" href="../../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.77.1">
<link rel="home" href="../../../../index.html" title="Chapter 1. Boost.Multiprecision">
<link rel="up" href="../fp_eg.html" title="Examples">
<link rel="prev" href="aos.html" title="Area of Circle">
<link rel="next" href="nd.html" title="Calculating a Derivative">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="aos.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../fp_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="nd.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h5 class="title">
<a name="boost_multiprecision.tut.floats.fp_eg.jel"></a><a class="link" href="jel.html" title="Defining a Special Function.">Defining
a Special Function.</a>
</h5></div></div></div>
<p>
In this example we'll show several implementations of the <a href="http://mathworld.wolfram.com/LambdaFunction.html" target="_top">Jahnke
and Emden Lambda function</a>, each implementation a little more
sophisticated than the last.
</p>
<p>
The Jahnke-Emden Lambda function is defined by the equation:
</p>
<div class="blockquote"><blockquote class="blockquote"><p>
<span class="emphasis"><em>JahnkeEmden(v, z) = Γ(v+1) * J<sub>v</sub>(z) / (z / 2)<sup>v</sup></em></span>
</p></blockquote></div>
<p>
If we were to implement this at double precision using Boost.Math's facilities
for the Gamma and Bessel function calls it would look like this:
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">JEL1</span><span class="special">(</span><span class="keyword">double</span> <span class="identifier">v</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tgamma</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span><span class="special">(</span><span class="identifier">z</span> <span class="special">/</span> <span class="number">2</span><span class="special">,</span> <span class="identifier">v</span><span class="special">);</span>
<span class="special">}</span>
</pre>
<p>
Calling this function as:
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">scientific</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">JEL1</span><span class="special">(</span><span class="number">2.5</span><span class="special">,</span> <span class="number">0.5</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
Yields the output:
</p>
<pre class="programlisting">9.822663964796047e-001</pre>
<p>
Now let's implement the function again, but this time using the multiprecision
type <code class="computeroutput"><span class="identifier">cpp_dec_float_50</span></code>
as the argument type:
</p>
<pre class="programlisting"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span>
<span class="identifier">JEL2</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tgamma</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">pow</span><span class="special">(</span><span class="identifier">z</span> <span class="special">/</span> <span class="number">2</span><span class="special">,</span> <span class="identifier">v</span><span class="special">);</span>
<span class="special">}</span>
</pre>
<p>
The implementation is almost the same as before, but with one key difference
- we can no longer call <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span></code>,
instead we must call the version inside the <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span></code>
namespace. In point of fact, we could have omitted the namespace prefix
on the call to <code class="computeroutput"><span class="identifier">pow</span></code> since
the right overload would have been found via <a href="http://en.wikipedia.org/wiki/Argument-dependent_name_lookup" target="_top">argument
dependent lookup</a> in any case.
</p>
<p>
Note also that the first argument to <code class="computeroutput"><span class="identifier">pow</span></code>
along with the argument to <code class="computeroutput"><span class="identifier">tgamma</span></code>
in the above code are actually expression templates. The <code class="computeroutput"><span class="identifier">pow</span></code> and <code class="computeroutput"><span class="identifier">tgamma</span></code>
functions will handle these arguments just fine.
</p>
<p>
Here's an example of how the function may be called:
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">scientific</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_dec_float_50</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">JEL2</span><span class="special">(</span><span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">2.5</span><span class="special">),</span> <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">0.5</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
Which outputs:
</p>
<pre class="programlisting">9.82266396479604757017335009796882833995903762577173e-01</pre>
<p>
Now that we've seen some non-template examples, lets repeat the code
again, but this time as a template that can be called either with a builtin
type (<code class="computeroutput"><span class="keyword">float</span></code>, <code class="computeroutput"><span class="keyword">double</span></code> etc), or with a multiprecision
type:
</p>
<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Float</span><span class="special">></span>
<span class="identifier">Float</span> <span class="identifier">JEL3</span><span class="special">(</span><span class="identifier">Float</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">Float</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span><span class="special">;</span>
<span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tgamma</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">z</span> <span class="special">/</span> <span class="number">2</span><span class="special">,</span> <span class="identifier">v</span><span class="special">);</span>
<span class="special">}</span>
</pre>
<p>
Once again the code is almost the same as before, but the call to <code class="computeroutput"><span class="identifier">pow</span></code> has changed yet again. We need
the call to resolve to either <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span></code>
(when the argument is a builtin type), or to <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">pow</span></code>
(when the argument is a multiprecision type). We do that by making the
call unqualified so that versions of <code class="computeroutput"><span class="identifier">pow</span></code>
defined in the same namespace as type <code class="computeroutput"><span class="identifier">Float</span></code>
are found via argument dependent lookup, while the <code class="computeroutput"><span class="keyword">using</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span></code> directive makes the standard library
versions visible for builtin floating point types.
</p>
<p>
Let's call the function with both <code class="computeroutput"><span class="keyword">double</span></code>
and multiprecision arguments:
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">scientific</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">JEL3</span><span class="special">(</span><span class="number">2.5</span><span class="special">,</span> <span class="number">0.5</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">scientific</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_dec_float_50</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">JEL3</span><span class="special">(</span><span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">2.5</span><span class="special">),</span> <span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="number">0.5</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
Which outputs:
</p>
<pre class="programlisting">9.822663964796047e-001
9.82266396479604757017335009796882833995903762577173e-01
</pre>
<p>
Unfortunately there is a problem with this version: if we were to call
it like this:
</p>
<pre class="programlisting"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span> <span class="identifier">v</span><span class="special">(</span><span class="number">2</span><span class="special">),</span> <span class="identifier">z</span><span class="special">(</span><span class="number">0.5</span><span class="special">);</span>
<span class="identifier">JEL3</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">0.5</span><span class="special">,</span> <span class="identifier">z</span><span class="special">);</span>
</pre>
<p>
Then we would get a long and inscrutable error message from the compiler:
the problem here is that the first argument to <code class="computeroutput"><span class="identifier">JEL3</span></code>
is not a number type, but an expression template. We could obviously
add a typecast to fix the issue:
</p>
<pre class="programlisting"><span class="identifier">JEL</span><span class="special">(</span><span class="identifier">cpp_dec_float_50</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">0.5</span><span class="special">),</span> <span class="identifier">z</span><span class="special">);</span>
</pre>
<p>
However, if we want the function JEL to be truly reusable, then a better
solution might be preferred. To achieve this we can borrow some code
from Boost.Math which calculates the return type of mixed-argument functions,
here's how the new code looks now:
</p>
<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Float1</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Float2</span><span class="special">></span>
<span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">promote_args</span><span class="special"><</span><span class="identifier">Float1</span><span class="special">,</span> <span class="identifier">Float2</span><span class="special">>::</span><span class="identifier">type</span>
<span class="identifier">JEL4</span><span class="special">(</span><span class="identifier">Float1</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">Float2</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span><span class="special">;</span>
<span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tgamma</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">z</span> <span class="special">/</span> <span class="number">2</span><span class="special">,</span> <span class="identifier">v</span><span class="special">);</span>
<span class="special">}</span>
</pre>
<p>
As you can see the two arguments to the function are now separate template
types, and the return type is computed using the <code class="computeroutput"><span class="identifier">promote_args</span></code>
metafunction from Boost.Math.
</p>
<p>
Now we can call:
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">scientific</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">cpp_dec_float_100</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">JEL4</span><span class="special">(</span><span class="identifier">cpp_dec_float_100</span><span class="special">(</span><span class="number">2</span><span class="special">)</span> <span class="special">+</span> <span class="number">0.5</span><span class="special">,</span> <span class="identifier">cpp_dec_float_100</span><span class="special">(</span><span class="number">0.5</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
And get 100 digits of output:
</p>
<pre class="programlisting">9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01</pre>
<p>
As a bonus, we can now call the function not just with expression templates,
but with other mixed types as well: for example <code class="computeroutput"><span class="keyword">float</span></code>
and <code class="computeroutput"><span class="keyword">double</span></code> or <code class="computeroutput"><span class="keyword">int</span></code> and <code class="computeroutput"><span class="keyword">double</span></code>,
and the correct return type will be computed in each case.
</p>
<p>
Note that while in this case we didn't have to change the body of the
function, in the general case any function like this which creates local
variables internally would have to use <code class="computeroutput"><span class="identifier">promote_args</span></code>
to work out what type those variables should be, for example:
</p>
<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Float1</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Float2</span><span class="special">></span>
<span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">promote_args</span><span class="special"><</span><span class="identifier">Float1</span><span class="special">,</span> <span class="identifier">Float2</span><span class="special">>::</span><span class="identifier">type</span>
<span class="identifier">JEL5</span><span class="special">(</span><span class="identifier">Float1</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">Float2</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">pow</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">promote_args</span><span class="special"><</span><span class="identifier">Float1</span><span class="special">,</span> <span class="identifier">Float2</span><span class="special">>::</span><span class="identifier">type</span> <span class="identifier">variable_type</span><span class="special">;</span>
<span class="identifier">variable_type</span> <span class="identifier">t</span> <span class="special">=</span> <span class="identifier">pow</span><span class="special">(</span><span class="identifier">z</span> <span class="special">/</span> <span class="number">2</span><span class="special">,</span> <span class="identifier">v</span><span class="special">);</span>
<span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tgamma</span><span class="special">(</span><span class="identifier">v</span> <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">t</span><span class="special">;</span>
<span class="special">}</span>
</pre>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2002-2013 John Maddock and Christopher Kormanyos<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="aos.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../fp_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="nd.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
