← Back to branch summary

~ubuntu-branches/debian/sid/pgadmin3/sid

« back to all changes in this revision

Viewing changes to pgadmin/pgscript/utilities/m_apm/mapm_exp.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Gerfried Fuchs
  • Date: 2009-07-30 12:27:16 UTC
  • mfrom: (1.1.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20090730122716-fddbh42on721bbs2
Tags: 1.10.0-1
* New upstream release.
* Adjusted watch file to match release candidates.
* Updated to Standards-Version 3.8.2:
  - Moved to Section: database.
  - Add DEB_BUILD_OPTIONS support for parallel building.
  - Move from findstring to filter suggestion for DEB_BUILD_OPTIONS parsing.
* pgagent got split into its own separate source package by upstream.
* Exclude Docs.vcproj from installation.
* Move doc-base.enus from pgadmin3 to pgadmin3-data package, the files are
  in there too.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pgadmin/pgscript/utilities/m_apm/mapm_exp.cpp
 
1
 
 
2
/* 
 
3
 *  M_APM  -  mapm_exp.c
 
4
 *
 
5
 *  Copyright (C) 1999 - 2007   Michael C. Ring
 
6
 *
 
7
 *  Permission to use, copy, and distribute this software and its
 
8
 *  documentation for any purpose with or without fee is hereby granted,
 
9
 *  provided that the above copyright notice appear in all copies and
 
10
 *  that both that copyright notice and this permission notice appear
 
11
 *  in supporting documentation.
 
12
 *
 
13
 *  Permission to modify the software is granted. Permission to distribute
 
14
 *  the modified code is granted. Modifications are to be distributed by
 
15
 *  using the file 'license.txt' as a template to modify the file header.
 
16
 *  'license.txt' is available in the official MAPM distribution.
 
17
 *
 
18
 *  This software is provided "as is" without express or implied warranty.
 
19
 */
 
20
 
 
21
/*
 
22
 *
 
23
 *      This file contains the EXP function.
 
24
 *
 
25
 */
 
26
 
 
27
#include "pgAdmin3.h"
 
28
#include "pgscript/utilities/mapm-lib/m_apm_lc.h"
 
29
 
 
30
static  M_APM  MM_exp_log2R;
 
31
static  M_APM  MM_exp_512R;
 
32
static  int    MM_firsttime1 = TRUE;
 
33
 
 
34
/****************************************************************************/
 
35
void    M_free_all_exp()
 
36
{
 
37
if (MM_firsttime1 == FALSE)
 
38
  {
 
39
   m_apm_free(MM_exp_log2R);
 
40
   m_apm_free(MM_exp_512R);
 
41
 
 
42
   MM_firsttime1 = TRUE;
 
43
  }
 
44
}
 
45
/****************************************************************************/
 
46
void    m_apm_exp(M_APM r, int places, M_APM x)
 
47
{
 
48
M_APM   tmp7, tmp8, tmp9;
 
49
int     dplaces, nn, ii;
 
50
 
 
51
if (MM_firsttime1)
 
52
  {
 
53
   MM_firsttime1 = FALSE;
 
54
 
 
55
   MM_exp_log2R = m_apm_init();
 
56
   MM_exp_512R  = m_apm_init();
 
57
 
 
58
   m_apm_set_string(MM_exp_log2R, "1.44269504089");   /* ~ 1 / log(2) */
 
59
   m_apm_set_string(MM_exp_512R,  "1.953125E-3");     /*   1 / 512    */
 
60
  }
 
61
 
 
62
tmp7 = M_get_stack_var();
 
63
tmp8 = M_get_stack_var();
 
64
tmp9 = M_get_stack_var();
 
65
 
 
66
if (x->m_apm_sign == 0)         /* if input == 0, return '1' */
 
67
  {
 
68
   m_apm_copy(r, MM_One);
 
69
   M_restore_stack(3);
 
70
   return;
 
71
  }
 
72
 
 
73
if (x->m_apm_exponent <= -3)  /* already small enough so call _raw directly */
 
74
  {
 
75
   M_raw_exp(tmp9, (places + 6), x);
 
76
   m_apm_round(r, places, tmp9);
 
77
   M_restore_stack(3);
 
78
   return;
 
79
  }
 
80
 
 
81
/*
 
82
    From David H. Bailey's MPFUN Fortran package :
 
83
 
 
84
    exp (t) =  (1 + r + r^2 / 2! + r^3 / 3! + r^4 / 4! ...) ^ q * 2 ^ n
 
85
 
 
86
    where q = 256, r = t' / q, t' = t - n Log(2) and where n is chosen so
 
87
    that -0.5 Log(2) < t' <= 0.5 Log(2).  Reducing t mod Log(2) and
 
88
    dividing by 256 insures that -0.001 < r <= 0.001, which accelerates
 
89
    convergence in the above series.
 
90
 
 
91
    I use q = 512 and also limit how small 'r' can become. The 'r' used
 
92
    here is limited in magnitude from 1.95E-4 < |r| < 1.35E-3. Forcing
 
93
    'r' into a narrow range keeps the algorithm 'well behaved'.
 
94
 
 
95
    ( the range is [0.1 / 512] to [log(2) / 512] )
 
96
*/
 
97
 
 
98
if (M_exp_compute_nn(&nn, tmp7, x) != 0)
 
99
  {
 
100
   M_apm_log_error_msg(M_APM_RETURN, 
 
101
      "\'m_apm_exp\', Input too large, Overflow");
 
102
 
 
103
   M_set_to_zero(r);
 
104
   M_restore_stack(3);
 
105
   return;
 
106
  }
 
107
 
 
108
dplaces = places + 8;
 
109
 
 
110
/* check to make sure our log(2) is accurate enough */
 
111
 
 
112
M_check_log_places(dplaces);
 
113
 
 
114
m_apm_multiply(tmp8, tmp7, MM_lc_log2);
 
115
m_apm_subtract(tmp7, x, tmp8);
 
116
 
 
117
/*
 
118
 *     guarantee that |tmp7| is between 0.1 and 0.9999999....
 
119
 *     (in practice, the upper limit only reaches log(2), 0.693... )
 
120
 */
 
121
 
 
122
while (TRUE)
 
123
  {
 
124
   if (tmp7->m_apm_sign != 0)
 
125
     {
 
126
      if (tmp7->m_apm_exponent == 0)
 
127
        break;
 
128
     }
 
129
     
 
130
   if (tmp7->m_apm_sign >= 0)
 
131
     {
 
132
      nn++;
 
133
      m_apm_subtract(tmp8, tmp7, MM_lc_log2);
 
134
      m_apm_copy(tmp7, tmp8);
 
135
     }
 
136
   else
 
137
     {
 
138
      nn--;
 
139
      m_apm_add(tmp8, tmp7, MM_lc_log2);
 
140
      m_apm_copy(tmp7, tmp8);
 
141
     }
 
142
  }
 
143
 
 
144
m_apm_multiply(tmp9, tmp7, MM_exp_512R);
 
145
 
 
146
/* perform the series expansion ... */
 
147
 
 
148
M_raw_exp(tmp8, dplaces, tmp9);
 
149
 
 
150
/*
 
151
 *   raise result to the 512 power
 
152
 *
 
153
 *   note : x ^ 512  =  (((x ^ 2) ^ 2) ^ 2) ... 9 times
 
154
 */
 
155
 
 
156
ii = 9;
 
157
 
 
158
while (TRUE)
 
159
  {
 
160
   m_apm_multiply(tmp9, tmp8, tmp8);
 
161
   m_apm_round(tmp8, dplaces, tmp9);
 
162
 
 
163
   if (--ii == 0)
 
164
     break;
 
165
  }
 
166
 
 
167
/* now compute 2 ^ N */
 
168
 
 
169
m_apm_integer_pow(tmp7, dplaces, MM_Two, nn);
 
170
 
 
171
m_apm_multiply(tmp9, tmp7, tmp8);
 
172
m_apm_round(r, places, tmp9);
 
173
 
 
174
M_restore_stack(3);                    /* restore the 3 locals we used here */
 
175
}
 
176
/****************************************************************************/
 
177
/*
 
178
        compute  int *n  = round_to_nearest_int(a / log(2))
 
179
                 M_APM b = MAPM version of *n
 
180
 
 
181
        returns      0: OK
 
182
                 -1, 1: failure
 
183
*/
 
184
int     M_exp_compute_nn(int *n, M_APM b, M_APM a)
 
185
{
 
186
M_APM   tmp0, tmp1;
 
187
void    *vp;
 
188
char    *cp, sbuf[48];
 
189
int     kk;
 
190
 
 
191
*n   = 0;
 
192
vp   = NULL;
 
193
cp   = sbuf;
 
194
tmp0 = M_get_stack_var();
 
195
tmp1 = M_get_stack_var();
 
196
 
 
197
/* find 'n' and convert it to a normal C int            */
 
198
/* we just need an approx 1/log(2) for this calculation */
 
199
 
 
200
m_apm_multiply(tmp1, a, MM_exp_log2R);
 
201
 
 
202
/* round to the nearest int */
 
203
 
 
204
if (tmp1->m_apm_sign >= 0)
 
205
  {
 
206
   m_apm_add(tmp0, tmp1, MM_0_5);
 
207
   m_apm_floor(tmp1, tmp0);
 
208
  }
 
209
else
 
210
  {
 
211
   m_apm_subtract(tmp0, tmp1, MM_0_5);
 
212
   m_apm_ceil(tmp1, tmp0);
 
213
  }
 
214
 
 
215
kk = tmp1->m_apm_exponent;
 
216
if (kk >= 42)
 
217
  {
 
218
   if ((vp = (void *)MAPM_MALLOC((kk + 16) * sizeof(char))) == NULL)
 
219
     {
 
220
      /* fatal, this does not return */
 
221
 
 
222
      M_apm_log_error_msg(M_APM_FATAL, "\'M_exp_compute_nn\', Out of memory");
 
223
     }
 
224
 
 
225
   cp = (char *)vp;
 
226
  }
 
227
 
 
228
m_apm_to_integer_string(cp, tmp1);
 
229
*n = atoi(cp);
 
230
 
 
231
m_apm_set_long(b, (long)(*n));
 
232
 
 
233
kk = m_apm_compare(b, tmp1);
 
234
 
 
235
if (vp != NULL)
 
236
  MAPM_FREE(vp);
 
237
 
 
238
M_restore_stack(2);
 
239
return(kk);
 
240
}
 
241
/****************************************************************************/
 
242
/*
 
243
        calculate the exponential function using the following series :
 
244
 
 
245
                              x^2     x^3     x^4     x^5
 
246
        exp(x) == 1  +  x  +  ---  +  ---  +  ---  +  ---  ...
 
247
                               2!      3!      4!      5!
 
248
 
 
249
*/
 
250
void    M_raw_exp(M_APM rr, int places, M_APM xx)
 
251
{
 
252
M_APM   tmp0, digit, term;
 
253
int     tolerance,  local_precision, prev_exp;
 
254
long    m1;
 
255
 
 
256
tmp0  = M_get_stack_var();
 
257
term  = M_get_stack_var();
 
258
digit = M_get_stack_var();
 
259
 
 
260
local_precision = places + 8;
 
261
tolerance       = -(places + 4);
 
262
prev_exp        = 0;
 
263
 
 
264
m_apm_add(rr, MM_One, xx);
 
265
m_apm_copy(term, xx);
 
266
 
 
267
m1 = 2L;
 
268
 
 
269
while (TRUE)
 
270
  {
 
271
   m_apm_set_long(digit, m1);
 
272
   m_apm_multiply(tmp0, term, xx);
 
273
   m_apm_divide(term, local_precision, tmp0, digit);
 
274
   m_apm_add(tmp0, rr, term);
 
275
   m_apm_copy(rr, tmp0);
 
276
 
 
277
   if ((term->m_apm_exponent < tolerance) || (term->m_apm_sign == 0))
 
278
     break;
 
279
 
 
280
   if (m1 != 2L)
 
281
     {
 
282
      local_precision = local_precision + term->m_apm_exponent - prev_exp;
 
283
 
 
284
      if (local_precision < 20)
 
285
        local_precision = 20;
 
286
     }
 
287
 
 
288
   prev_exp = term->m_apm_exponent;
 
289
   m1++;
 
290
  }
 
291
 
 
292
M_restore_stack(3);                    /* restore the 3 locals we used here */
 
293
}
 
294
/****************************************************************************/