← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty-proposed

« back to all changes in this revision

Viewing changes to util/bigfloat/bigfloat.h

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
util/bigfloat/bigfloat.h
 
1
#ifndef _BIGFLOAT_H
 
2
#define _BIGFLOAT_H
 
3
///
 
4
/// This file is part of Rheolef.
 
5
///
 
6
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
7
///
 
8
/// Rheolef is free software; you can redistribute it and/or modify
 
9
/// it under the terms of the GNU General Public License as published by
 
10
/// the Free Software Foundation; either version 2 of the License, or
 
11
/// (at your option) any later version.
 
12
///
 
13
/// Rheolef is distributed in the hope that it will be useful,
 
14
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
15
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
16
/// GNU General Public License for more details.
 
17
///
 
18
/// You should have received a copy of the GNU General Public License
 
19
/// along with Rheolef; if not, write to the Free Software
 
20
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
21
/// 
 
22
/// =========================================================================
 
23
//
 
24
// bigfloat<N> class with N digits precision
 
25
// as a remplacement for the double
 
26
//
 
27
//    class bigfloat<N>
 
28
//    class numeric_limits <bigfloat<N> >
 
29
//    class numeric_flags  <bigfloat<N> >
 
30
//
 
31
// author: Pierre.Saramito@imag.fr
 
32
//
 
33
// date: 25 january 2000
 
34
//
 
35
// IMPLEMENTATION NOTE:
 
36
// this is a wrapper class for cln_F:
 
37
// provides an interface comparable to double.
 
38
//
 
39
// TODO: ostream
 
40
//   knows setprecision(n) but 
 
41
// does not respond to 
 
42
//   scientific & fixed & uppercase
 
43
//
 
44
// TODO: complete the math lib, as IEEE
 
45
// -> see the blitz test suite.
 
46
//
 
47
 
 
48
#ifdef HAVE_CONFIG_H
 
49
#include "config.h"
 
50
#endif
 
51
 
 
52
#include <cmath>
 
53
#include <cln/cln.h>
 
54
#ifdef _RHEOLEF_HAVE_GINAC
 
55
#include <ginac/ginac.h>
 
56
#endif // _RHEOLEF_HAVE_GINAC
 
57
#include "rheolef/numeric_limits.h"
 
58
#include "rheolef/numeric_flags.h"
 
59
#include "all_float_io.h"
 
60
#include <sstream>
 
61
 
 
62
template <int N>
 
63
class bigfloat {
 
64
public:
 
65
 
 
66
// allocators/deallocators
 
67
 
 
68
  bigfloat(double x = 0);
 
69
  bigfloat(const char* s);
 
70
  bigfloat(int x);
 
71
  bigfloat(long x);
 
72
  bigfloat(unsigned int x);
 
73
  bigfloat(unsigned long x);
 
74
 
 
75
// type conversions
 
76
 
 
77
  operator double() const; // compiler do not send any warning, such as
 
78
  operator int() const;    // warning: initialization to `int' from `bigfloat' !!
 
79
#ifdef _RHEOLEF_HAVE_GINAC
 
80
  operator GiNaC::ex() const;
 
81
#endif // _RHEOLEF_HAVE_GINAC
 
82
 
 
83
// assignments
 
84
 
 
85
  bigfloat<N>& operator = (const bigfloat<N>&);
 
86
  bigfloat<N>& operator = (const char*&);
 
87
  bigfloat<N>& operator = (const double&);
 
88
  bigfloat<N>& operator = (const int&);
 
89
  bigfloat<N>& operator = (const long&);
 
90
  bigfloat<N>& operator = (const unsigned int&);
 
91
  bigfloat<N>& operator = (const unsigned long&);
 
92
 
 
93
// operators
 
94
  template <int M>
 
95
  friend bigfloat<M> operator - (const bigfloat<M>&);
 
96
 
 
97
# define bigfloat_binary_operator_basic(op,T1,T2)               \
 
98
  template <int M>                                              \
 
99
  friend bigfloat<M> operator op (const T1&, const T2&);
 
100
 
 
101
# define bigfloat_binary_operator(op)                           \
 
102
bigfloat_binary_operator_basic(op,bigfloat<M>,bigfloat<M>)      \
 
103
bigfloat_binary_operator_basic(op,bigfloat<M>,double)           \
 
104
bigfloat_binary_operator_basic(op,bigfloat<M>,int)              \
 
105
bigfloat_binary_operator_basic(op,bigfloat<M>,long)             \
 
106
bigfloat_binary_operator_basic(op,bigfloat<M>,unsigned int)     \
 
107
bigfloat_binary_operator_basic(op,bigfloat<M>,unsigned long)    \
 
108
bigfloat_binary_operator_basic(op,double,bigfloat<M>)           \
 
109
bigfloat_binary_operator_basic(op,int,bigfloat<M>)              \
 
110
bigfloat_binary_operator_basic(op,long,bigfloat<M>)             \
 
111
bigfloat_binary_operator_basic(op,unsigned int,bigfloat<M>)     \
 
112
bigfloat_binary_operator_basic(op,unsigned long,bigfloat<M>)
 
113
 
 
114
bigfloat_binary_operator(+)
 
115
bigfloat_binary_operator(-)
 
116
bigfloat_binary_operator(*)
 
117
bigfloat_binary_operator(/)
 
118
 
 
119
# undef bigfloat_binary_operator_basic
 
120
# undef bigfloat_binary_operator
 
121
 
 
122
// computed assignments
 
123
 
 
124
# define bigfloat_computed_assignment(op)                       \
 
125
  bigfloat<N>& operator op(const bigfloat<N>&);                 \
 
126
  bigfloat<N>& operator op(const double&);                      \
 
127
  bigfloat<N>& operator op(const int&);                         \
 
128
  bigfloat<N>& operator op(const long&);                        \
 
129
  bigfloat<N>& operator op(const unsigned int&);                \
 
130
  bigfloat<N>& operator op(const unsigned long&);
 
131
 
 
132
bigfloat_computed_assignment(+=)
 
133
bigfloat_computed_assignment(-=)
 
134
bigfloat_computed_assignment(*=)
 
135
bigfloat_computed_assignment(/=)
 
136
 
 
137
#undef bigfloat_computed_assignment
 
138
 
 
139
// comparators
 
140
 
 
141
#define bigfloat_comparator_basic(op)                                   \
 
142
  template <int M>                                                      \
 
143
  friend bool operator op (const bigfloat<M>&, const bigfloat<M>&);
 
144
#define bigfloat_comparator_both(op,T)                                  \
 
145
  template <int M>                                                      \
 
146
  friend bool operator op (const bigfloat<M>&, const T&);               \
 
147
  template <int M>                                                      \
 
148
  friend bool operator op (const T&, const bigfloat<M>&);
 
149
#define bigfloat_comparator(op)                                         \
 
150
  bigfloat_comparator_basic(op)                                         \
 
151
  bigfloat_comparator_both(op,double)                                   \
 
152
  bigfloat_comparator_both(op,int)                                      \
 
153
  bigfloat_comparator_both(op,long)                                     \
 
154
  bigfloat_comparator_both(op,unsigned int)                             \
 
155
  bigfloat_comparator_both(op,unsigned long)
 
156
 
 
157
bigfloat_comparator(==)
 
158
bigfloat_comparator(!=)
 
159
bigfloat_comparator(<)
 
160
bigfloat_comparator(>)
 
161
bigfloat_comparator(<=)
 
162
bigfloat_comparator(>=)
 
163
 
 
164
#undef bigfloat_comparator
 
165
#undef bigfloat_comparator_basic
 
166
#undef bigfloat_comparator_both
 
167
 
 
168
// math functions
 
169
 
 
170
  template <int M> friend bigfloat<M> fabs(const bigfloat<M>& x);
 
171
  template <int M> friend bigfloat<M> sqrt(const bigfloat<M>& x);
 
172
  template <int M> friend bigfloat<M> sin(const bigfloat<M>& x);
 
173
  template <int M> friend bigfloat<M> cos(const bigfloat<M>& x);
 
174
  template <int M> friend bigfloat<M> tan(const bigfloat<M>& x);
 
175
  template <int M> friend bigfloat<M> log(const bigfloat<M>& x);
 
176
  template <int M> friend bigfloat<M> exp(const bigfloat<M>& x);
 
177
  template <int M> friend bigfloat<M> sinh(const bigfloat<M>& x);
 
178
  template <int M> friend bigfloat<M> cosh(const bigfloat<M>& x);
 
179
  template <int M> friend bigfloat<M> tanh(const bigfloat<M>& x);
 
180
  template <int M> friend bigfloat<M> atan(const bigfloat<M>& x);
 
181
  template <int M> friend bigfloat<M> acos(const bigfloat<M>& x);
 
182
  template <int M> friend bigfloat<M> sqr(const bigfloat<M>& x);
 
183
  template <int M> friend bigfloat<M> floor(const bigfloat<M>& x);
 
184
  template <int M> friend bigfloat<M> ceil(const bigfloat<M>& x);
 
185
  template <int M> friend bigfloat<M> trunc(const bigfloat<M>& x);
 
186
  template <int M> friend bigfloat<M> round(const bigfloat<M>& x);
 
187
 
 
188
  template <int M> friend bigfloat<M> log10(const bigfloat<M>& x);
 
189
  template <int M> friend bigfloat<M> pow(const bigfloat<M>& x, const bigfloat<M>& y);
 
190
  template <int M> friend bigfloat<M> pow(const bigfloat<M>& x, int n);
 
191
 
 
192
// inputs/outputs
 
193
 
 
194
  template <int M> friend std::ostream& operator <<  (std::ostream& s, const bigfloat<M>& x);
 
195
  template <int M> friend std::istream& operator >>  (std::istream& s, bigfloat<M>& x);
 
196
 
 
197
// constants (may be const! -> core dump..)
 
198
 
 
199
#ifdef TODO
 
200
  static bigfloat<N> Pi;
 
201
  static bigfloat<N> one_on_log10;
 
202
#endif // TODO
 
203
 
 
204
// numeric limits
 
205
 
 
206
// Returns the smallest float format which guarantees at least
 
207
// n decimal digits in the mantissa (after the decimal point).
 
208
// Methode from: cln-1.0.2/src/float/misc/cl_float_format.cc
 
209
// Mindestens 1+n Dezimalstellen (inklusive Vorkommastelle)
 
210
// bedeutet mindestens ceiling((1+n)*ln(10)/ln(2)) Bin�rstellen.
 
211
// ln(10)/ln(2) = 3.321928095 = (bin�r) 11.01010010011010011110000100101111...
 
212
//                       = (bin�r) 100 - 0.10101101100101100001111011010001
 
213
// Durch diese Berechnungsmethode wird das Ergebnis sicher >= (1+n)*ln(10)/ln(2)
 
214
// sein, evtl. um ein paar Bit zu gro�, aber nicht zu klein.
 
215
 
 
216
  static const int digits
 
217
        = ((1+N) << 2) 
 
218
        - ((1+N) >> 1)  - ((1+N) >> 3)  - ((1+N) >> 5)
 
219
        - ((1+N) >> 6)  - ((1+N) >> 8)  - ((1+N) >> 9)
 
220
        - ((1+N) >> 12) - ((1+N) >> 14) - ((1+N) >> 15);
 
221
 
 
222
  static const int digits10             = N;
 
223
 
 
224
  static bigfloat<N> epsilon();
 
225
  static bigfloat<N> min();
 
226
  static bigfloat<N> max();
 
227
 
 
228
// computation of some constants
 
229
 
 
230
  static bigfloat<N> compute_pi();
 
231
 
 
232
// data
 
233
 
 
234
protected:
 
235
 
 
236
  cln::cl_F _f;
 
237
 
 
238
  // internal constructor
 
239
  bigfloat(const cln::cl_F& y) : _f(y) {}
 
240
 
 
241
public:
 
242
  // for debug
 
243
  std::ostream& dump(std::ostream&);
 
244
};
 
245
//
 
246
// TODO: complete the fields...
 
247
//
 
248
namespace std {
 
249
template<int N>
 
250
class numeric_limits <bigfloat<N> > {
 
251
public:
 
252
    static const int  digits         = bigfloat<N>::digits;
 
253
    static const int  digits10       = bigfloat<N>::digits10;
 
254
 
 
255
    // have no representation...
 
256
    static bigfloat<N> min ()          { return bigfloat<N>::min(); }
 
257
    static bigfloat<N> max ()          { return bigfloat<N>::max(); }
 
258
 
 
259
    static bigfloat<N> epsilon ()      { return bigfloat<N>::epsilon(); }
 
260
    static bigfloat<N> denorm_min ()   { return bigfloat<N>::min(); }
 
261
 
 
262
    // CLN does not implement features like NaNs, denormalized numbers
 
263
    // and gradual underflow.
 
264
    // If the exponent range of some floating-point type is too limited 
 
265
    // for your application, choose another floating-point type with
 
266
    // larger exponent range. 
 
267
 
 
268
#ifdef TODO
 
269
    static const int  min_exponent   = INT_MIN;
 
270
    static const int  max_exponent   = INT_MAX;
 
271
    static const int  min_exponent10 = INT_MIN;
 
272
    static const int  max_exponent10 = INT_MAX;
 
273
 
 
274
    __NUMERIC_LIMITS_FLOAT(<bigfloat<N>)
 
275
 
 
276
#endif // TODO
 
277
};
 
278
} // namespace std
 
279
template <int N>
 
280
class numeric_flags<bigfloat<N> > {
 
281
public:
 
282
    static bool output_as_double()       { return _output_as_double; }
 
283
    static void output_as_double(bool f) { _output_as_double = f; }
 
284
protected:
 
285
    static bool _output_as_double;
 
286
};
 
287
 
 
288
// flag initialisation
 
289
 
 
290
template <int N>
 
291
bool numeric_flags<bigfloat<N> >::_output_as_double = false;
 
292
 
 
293
 
 
294
// numeric limits
 
295
 
 
296
template <int N>
 
297
inline
 
298
bigfloat<N>
 
299
bigfloat<N>::min() 
 
300
{
 
301
    return bigfloat<N>(most_negative_float(cln::float_format_t(digits)));
 
302
}
 
303
template <int N>
 
304
inline
 
305
bigfloat<N>
 
306
bigfloat<N>::max() 
 
307
{
 
308
    return bigfloat<N>(most_positive_float(cln::float_format_t(digits)));
 
309
}
 
310
template <int N>
 
311
inline
 
312
bigfloat<N>
 
313
bigfloat<N>::epsilon() 
 
314
{
 
315
    return bigfloat<N>(float_epsilon(cln::float_format_t(digits)));
 
316
}
 
317
 
 
318
// constants
 
319
 
 
320
template <int N>
 
321
inline
 
322
bigfloat<N> 
 
323
bigfloat<N>::compute_pi()
 
324
{
 
325
        return bigfloat<N>(cln::pi(cln::float_format_t(digits)));
 
326
}
 
327
 
 
328
// allocators/deallocators
 
329
 
 
330
#define bigfloat_cstor(T)                               \
 
331
template <int N>                                        \
 
332
inline                                                  \
 
333
bigfloat<N>::bigfloat(T x)                              \
 
334
        : _f(cln::cl_float(x, cln::float_format_t(digits)))     \
 
335
{                                                       \
 
336
}
 
337
bigfloat_cstor(double)
 
338
bigfloat_cstor(int)
 
339
bigfloat_cstor(unsigned int)
 
340
#undef  bigfloat_cstor
 
341
 
 
342
template <int N>
 
343
inline
 
344
bigfloat<N>::bigfloat(const char* s)
 
345
        : _f(cln::cl_float(cln::cl_F(s), cln::float_format_t(digits)))
 
346
{
 
347
}
 
348
 
 
349
// type conversions
 
350
 
 
351
template <int N>
 
352
inline
 
353
bigfloat<N>::operator double() const
 
354
{
 
355
        return cln::double_approx(_f);
 
356
}
 
357
template <int N>
 
358
inline
 
359
bigfloat<N>::operator int() const
 
360
{
 
361
        return int(double(*this));
 
362
}
 
363
#ifdef _RHEOLEF_HAVE_GINAC
 
364
template <int N>
 
365
inline
 
366
bigfloat<N>::operator GiNaC::ex() const
 
367
{
 
368
        std::stringstream s;
 
369
        s << std::setprecision(N) << *this;
 
370
        // std::string sx = s.str();
 
371
        GiNaC::ex x;
 
372
        s >> x;
 
373
        return x;
 
374
}
 
375
#endif // _RHEOLEF_HAVE_GINAC
 
376
 
 
377
// assignments
 
378
 
 
379
template <int N>
 
380
inline
 
381
bigfloat<N>&
 
382
bigfloat<N>::operator = (const bigfloat<N>& x) 
 
383
{
 
384
        _f = x._f;
 
385
        return *this;
 
386
}
 
387
#define bigfloat_assignment(T)                          \
 
388
template <int N>                                        \
 
389
inline                                                  \
 
390
bigfloat<N>&                                            \
 
391
bigfloat<N>::operator = (const T& x)                    \
 
392
{                                                       \
 
393
        _f = cln::cl_float(x, cln::float_format_t(digits));     \
 
394
        return *this;                                   \
 
395
}
 
396
bigfloat_assignment(double)
 
397
bigfloat_assignment(int)
 
398
bigfloat_assignment(unsigned int)
 
399
#undef bigfloat_assignment
 
400
 
 
401
template <int N>
 
402
inline
 
403
bigfloat<N>&
 
404
bigfloat<N>::operator = (const char*& x)
 
405
{
 
406
        _f = cln::cl_float(cln::cl_F(x), cln::float_format_t(digits));
 
407
        return *this;
 
408
}
 
409
 
 
410
// operators
 
411
 
 
412
template <int N>
 
413
inline
 
414
bigfloat<N>
 
415
operator - (const bigfloat<N>& x)
 
416
{
 
417
        return bigfloat<N>(-x._f);
 
418
}
 
419
 
 
420
# define bigfloat_binary_operator_basic(op)                             \
 
421
template <int N>                                                        \
 
422
inline                                                                  \
 
423
bigfloat<N>                                                             \
 
424
operator op (const bigfloat<N>& x, const bigfloat<N>& y)                \
 
425
{                                                                       \
 
426
        return bigfloat<N>(x._f op y._f);                               \
 
427
}
 
428
# define bigfloat_binary_operator_left(op,T1)                           \
 
429
template <int N>                                                        \
 
430
inline                                                                  \
 
431
bigfloat<N>                                                             \
 
432
operator op (const T1& x, const bigfloat<N>& y)                         \
 
433
{                                                                       \
 
434
        return bigfloat<N>(x) op y;                                     \
 
435
}
 
436
# define bigfloat_binary_operator_right(op,T2)                          \
 
437
template <int N>                                                        \
 
438
inline                                                                  \
 
439
bigfloat<N>                                                             \
 
440
operator op (const bigfloat<N>& x, const T2& y)                         \
 
441
{                                                                       \
 
442
        return x op bigfloat<N>(y);                                     \
 
443
}
 
444
# define bigfloat_binary_operator_both(op,T)                            \
 
445
bigfloat_binary_operator_left(op,T)                                     \
 
446
bigfloat_binary_operator_right(op,T)
 
447
 
 
448
# define bigfloat_binary_operator(op)                                   \
 
449
bigfloat_binary_operator_basic(op)                                      \
 
450
bigfloat_binary_operator_both(op,double)                                \
 
451
bigfloat_binary_operator_both(op,int)                                   \
 
452
bigfloat_binary_operator_both(op,long)                                  \
 
453
bigfloat_binary_operator_both(op,unsigned int)                          \
 
454
bigfloat_binary_operator_both(op,unsigned long)
 
455
 
 
456
bigfloat_binary_operator(+)
 
457
bigfloat_binary_operator(-)
 
458
bigfloat_binary_operator(*)
 
459
bigfloat_binary_operator(/)
 
460
 
 
461
# undef bigfloat_binary_operator_basic
 
462
# undef bigfloat_binary_operator
 
463
 
 
464
// computed assignments
 
465
 
 
466
# define bigfloat_computed_assignment_basic(op,T)               \
 
467
template <int N>                                                \
 
468
inline                                                          \
 
469
bigfloat<N>&                                                    \
 
470
bigfloat<N>::operator op##= (const T& x)                        \
 
471
{                                                               \
 
472
    (*this) = (*this) op x;                                     \
 
473
    return *this;                                               \
 
474
}
 
475
# define bigfloat_computed_assignment(op)                       \
 
476
bigfloat_computed_assignment_basic(op,bigfloat<N>)              \
 
477
bigfloat_computed_assignment_basic(op,double)                   \
 
478
bigfloat_computed_assignment_basic(op,int)                      \
 
479
bigfloat_computed_assignment_basic(op,long)                     \
 
480
bigfloat_computed_assignment_basic(op,unsigned int)             \
 
481
bigfloat_computed_assignment_basic(op,unsigned long)
 
482
 
 
483
bigfloat_computed_assignment(+)
 
484
bigfloat_computed_assignment(-)
 
485
bigfloat_computed_assignment(*)
 
486
bigfloat_computed_assignment(/)
 
487
 
 
488
# undef bigfloat_computed_assignment_basic
 
489
# undef bigfloat_computed_assignment
 
490
 
 
491
// comparators
 
492
 
 
493
#define bigfloat_comparator_basic(op)                           \
 
494
template <int N>                                                \
 
495
inline                                                          \
 
496
bool                                                            \
 
497
operator op (const bigfloat<N>& x, const bigfloat<N>& y)        \
 
498
{                                                               \
 
499
        return (cln::compare(x._f,y._f) op 0);                  \
 
500
}
 
501
#define bigfloat_comparator_both(op,T)                          \
 
502
template <int N>                                                \
 
503
inline                                                          \
 
504
bool                                                            \
 
505
operator op (const bigfloat<N>& x, const T& y)                  \
 
506
{                                                               \
 
507
        return (x op bigfloat<N>(y));                           \
 
508
}                                                               \
 
509
template <int N>                                                \
 
510
inline                                                          \
 
511
bool                                                            \
 
512
operator op (const T& x, const bigfloat<N>& y)                  \
 
513
{                                                               \
 
514
        return (bigfloat<N>(x) op y);                           \
 
515
}
 
516
#define bigfloat_comparator(op)                                 \
 
517
  bigfloat_comparator_basic(op)                                 \
 
518
  bigfloat_comparator_both(op,double)                           \
 
519
  bigfloat_comparator_both(op,int)                              \
 
520
  bigfloat_comparator_both(op,long)                             \
 
521
  bigfloat_comparator_both(op,unsigned int)                     \
 
522
  bigfloat_comparator_both(op,unsigned long)
 
523
 
 
524
bigfloat_comparator(==)                                         \
 
525
bigfloat_comparator(!=)                                         \
 
526
bigfloat_comparator(<)                                          \
 
527
bigfloat_comparator(>)                                          \
 
528
bigfloat_comparator(<=)                                         \
 
529
bigfloat_comparator(>=)
 
530
 
 
531
#undef bigfloat_comparator_basic
 
532
#undef bigfloat_comparator_both
 
533
#undef bigfloat_comparator
 
534
 
 
535
// math functions
 
536
 
 
537
#ifdef TODO
 
538
doubledouble hypot(const doubledouble, const doubledouble);
 
539
doubledouble cub(const doubledouble&);
 
540
doubledouble rint(const doubledouble&);
 
541
doubledouble fmod(const doubledouble&, const int);
 
542
doubledouble modf(const doubledouble&, doubledouble *ip);
 
543
void sincos(const doubledouble x, doubledouble& sinx, doubledouble& cosx);
 
544
doubledouble atan2(const doubledouble&, const doubledouble&);
 
545
doubledouble asin(const doubledouble&);
 
546
doubledouble erf(const doubledouble);
 
547
doubledouble erfc(const doubledouble);
 
548
doubledouble gamma(const doubledouble);
 
549
int  digits(const doubledouble&,const doubledouble&);
 
550
doubledouble modr(const doubledouble a, const doubledouble b, int& n, doubledouble& rem);
 
551
#endif // TODO
 
552
 
 
553
#define bigfloat_math_function(fct, cln_fct)            \
 
554
template <int N>                                        \
 
555
inline                                                  \
 
556
bigfloat<N>                                             \
 
557
fct (const bigfloat<N>& x)                              \
 
558
{                                                       \
 
559
        return cln_fct (x._f);                          \
 
560
}
 
561
bigfloat_math_function(fabs, abs)
 
562
bigfloat_math_function(sqrt, sqrt)
 
563
bigfloat_math_function(sin, sin)
 
564
bigfloat_math_function(cos, cos)
 
565
bigfloat_math_function(tan, tan)
 
566
bigfloat_math_function(log, ln)
 
567
bigfloat_math_function(exp, exp)
 
568
bigfloat_math_function(sinh, sinh)
 
569
bigfloat_math_function(cosh, cosh)
 
570
bigfloat_math_function(tanh, tanh)
 
571
bigfloat_math_function(atan, atan)
 
572
bigfloat_math_function(sqr, square)
 
573
bigfloat_math_function(floor, ffloor)
 
574
bigfloat_math_function(ceil, fceiling)
 
575
bigfloat_math_function(trunc, ftruncate)
 
576
bigfloat_math_function(round, fround)
 
577
#undef bigfloat_math_function
 
578
 
 
579
template <int N>
 
580
inline
 
581
bigfloat<N>
 
582
log10 (const bigfloat<N>& x)
 
583
{
 
584
#ifdef TODO
 
585
        // How to store one_on_log10 one time for all ?
 
586
        return log(x) * bigfloat<N>::one_on_log10;
 
587
#else
 
588
        return log(x) / log(bigfloat<N>(10));
 
589
#endif // TODO
 
590
}
 
591
template <int N>
 
592
inline
 
593
bigfloat<N>
 
594
pow (const bigfloat<N>& x, const bigfloat<N>& y)
 
595
{
 
596
        if (y == trunc(y)) {
 
597
            return pow (x, int(y));
 
598
        }
 
599
        if (x == bigfloat<N>(0)) return bigfloat<N>(1);
 
600
        if (x <  bigfloat<N>(0)) {
 
601
            fatal_macro ("argument x="<<x<<" of pow(x,y) may be >= 0"); 
 
602
        }
 
603
        return exp(y*log(x));
 
604
}
 
605
template <int N>
 
606
inline
 
607
bigfloat<N>
 
608
pow (const bigfloat<N>& x, int n)
 
609
{
 
610
        if (x == bigfloat<N>(0)) {
 
611
            if (n == 0) return bigfloat<N>(1);
 
612
            if (n <  0) return bigfloat<N>(1) / bigfloat<N>(0);
 
613
            else        return bigfloat<N>(0);
 
614
        }
 
615
        if (x <  bigfloat<N>(0)) {
 
616
            if (n % 2 == 1) {
 
617
                return - exp(n*log(-x));
 
618
            } else {
 
619
                return exp(n*log(-x));
 
620
            }
 
621
        }
 
622
        return exp(n*log(x));
 
623
}
 
624
template <int N>
 
625
inline
 
626
bigfloat<N>
 
627
acos (const bigfloat<N>& x)
 
628
{
 
629
        return bigfloat<N>(cln::cl_float(cln::realpart(acos(x._f))));
 
630
}
 
631
 
 
632
// inputs/outputs
 
633
 
 
634
template<int N>
 
635
inline
 
636
std::istream&
 
637
operator >> (std::istream& s, bigfloat<N>& x)
 
638
{
 
639
  return all_float_read (s, x);
 
640
}
 
641
template <int N>
 
642
std::ostream&
 
643
operator << (std::ostream& s, const bigfloat<N>& x)
 
644
{
 
645
  return all_float_write (s, x);
 
646
}
 
647
// for debug
 
648
template <int N>
 
649
inline
 
650
std::ostream&
 
651
bigfloat<N>::dump(std::ostream& s)
 
652
{
 
653
    return s << _f;
 
654
}
 
655
#endif //_BIGFLOAT_H