← 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_lg3.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_lg3.cpp
 
1
 
 
2
/* 
 
3
 *  M_APM  -  mapm_lg3.c
 
4
 *
 
5
 *  Copyright (C) 2003 - 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 function to compute log(2), log(10),
 
24
 *      and 1/log(10) to the desired precision using an AGM algorithm.
 
25
 *
 
26
 */
 
27
 
 
28
#include "pgAdmin3.h"
 
29
#include "pgscript/utilities/mapm-lib/m_apm_lc.h"
 
30
 
 
31
/*
 
32
 *  using the 'R' function (defined below) for 'N' decimal places :
 
33
 *
 
34
 *
 
35
 *                          -N             -N
 
36
 *  log(2)  =  R(1, 0.5 * 10  )  -  R(1, 10  ) 
 
37
 *
 
38
 *
 
39
 *                          -N             -N
 
40
 *  log(10) =  R(1, 0.1 * 10  )  -  R(1, 10  ) 
 
41
 *
 
42
 *
 
43
 *  In general:
 
44
 *
 
45
 *                    -N                -N
 
46
 *  log(x)  =  R(1, 10  / x)  -  R(1, 10  ) 
 
47
 *
 
48
 *
 
49
 *  I found this on a web site which went into considerable detail
 
50
 *  on the history of log(2). This formula is algebraically identical
 
51
 *  to the formula specified in J. Borwein and P. Borwein's book 
 
52
 *  "PI and the AGM". (reference algorithm 7.2)
 
53
 */
 
54
 
 
55
/****************************************************************************/
 
56
/*
 
57
 *      check if our local copy of log(2) & log(10) is precise
 
58
 *      enough for our purpose. if not, calculate them so it's
 
59
 *      as precise as desired, accurate to at least 'places'.
 
60
 */
 
61
void    M_check_log_places(int places)
 
62
{
 
63
M_APM   tmp6, tmp7, tmp8, tmp9;
 
64
int     dplaces;
 
65
 
 
66
dplaces = places + 4;
 
67
 
 
68
if (dplaces > MM_lc_log_digits)
 
69
  {
 
70
   MM_lc_log_digits = dplaces + 4;
 
71
   
 
72
   tmp6 = M_get_stack_var();
 
73
   tmp7 = M_get_stack_var();
 
74
   tmp8 = M_get_stack_var();
 
75
   tmp9 = M_get_stack_var();
 
76
   
 
77
   dplaces += 6 + (int)log10((double)places);
 
78
   
 
79
   m_apm_copy(tmp7, MM_One);
 
80
   tmp7->m_apm_exponent = -places;
 
81
   
 
82
   M_log_AGM_R_func(tmp8, dplaces, MM_One, tmp7);
 
83
   
 
84
   m_apm_multiply(tmp6, tmp7, MM_0_5);
 
85
   
 
86
   M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp6);
 
87
   
 
88
   m_apm_subtract(MM_lc_log2, tmp9, tmp8);               /* log(2) */
 
89
 
 
90
   tmp7->m_apm_exponent -= 1;                            /* divide by 10 */
 
91
 
 
92
   M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp7);
 
93
 
 
94
   m_apm_subtract(MM_lc_log10, tmp9, tmp8);              /* log(10) */
 
95
   m_apm_reciprocal(MM_lc_log10R, dplaces, MM_lc_log10); /* 1 / log(10) */
 
96
 
 
97
   M_restore_stack(4);
 
98
  }
 
99
}
 
100
/****************************************************************************/
 
101
 
 
102
/*
 
103
 *      define a notation for a function 'R' :
 
104
 *
 
105
 *
 
106
 *
 
107
 *                                    1
 
108
 *      R (a0, b0)  =  ------------------------------
 
109
 *
 
110
 *                          ----
 
111
 *                           \ 
 
112
 *                            \     n-1      2    2
 
113
 *                      1  -   |   2    *  (a  - b )
 
114
 *                            /              n    n
 
115
 *                           /
 
116
 *                          ----
 
117
 *                         n >= 0
 
118
 *
 
119
 *
 
120
 *      where a, b are the classic AGM iteration :
 
121
 *
 
122
 *     
 
123
 *      a    =  0.5 * (a  + b )
 
124
 *       n+1            n    n
 
125
 *
 
126
 *
 
127
 *      b    =  sqrt(a  * b )
 
128
 *       n+1          n    n
 
129
 *
 
130
 *
 
131
 *
 
132
 *      define a variable 'c' for more efficient computation :
 
133
 *
 
134
 *                                      2     2     2
 
135
 *      c    =  0.5 * (a  - b )    ,   c  =  a  -  b
 
136
 *       n+1            n    n          n     n     n
 
137
 *
 
138
 */
 
139
 
 
140
/****************************************************************************/
 
141
void    M_log_AGM_R_func(M_APM rr, int places, M_APM aa, M_APM bb)
 
142
{
 
143
M_APM   tmp1, tmp2, tmp3, tmp4, tmpC2, sum, pow_2, tmpA0, tmpB0;
 
144
int     tolerance, dplaces;
 
145
 
 
146
tmpA0 = M_get_stack_var();
 
147
tmpB0 = M_get_stack_var();
 
148
tmpC2 = M_get_stack_var();
 
149
tmp1  = M_get_stack_var();
 
150
tmp2  = M_get_stack_var();
 
151
tmp3  = M_get_stack_var();
 
152
tmp4  = M_get_stack_var();
 
153
sum   = M_get_stack_var();
 
154
pow_2 = M_get_stack_var();
 
155
 
 
156
tolerance = places + 8;
 
157
dplaces   = places + 16;
 
158
 
 
159
m_apm_copy(tmpA0, aa);
 
160
m_apm_copy(tmpB0, bb);
 
161
m_apm_copy(pow_2, MM_0_5);
 
162
 
 
163
m_apm_multiply(tmp1, aa, aa);               /* 0.5 * [ a ^ 2 - b ^ 2 ] */
 
164
m_apm_multiply(tmp2, bb, bb);
 
165
m_apm_subtract(tmp3, tmp1, tmp2);
 
166
m_apm_multiply(sum, MM_0_5, tmp3);
 
167
 
 
168
while (TRUE)
 
169
  {
 
170
   m_apm_subtract(tmp1, tmpA0, tmpB0);      /* C n+1 = 0.5 * [ An - Bn ] */
 
171
   m_apm_multiply(tmp4, MM_0_5, tmp1);      /* C n+1 */
 
172
   m_apm_multiply(tmpC2, tmp4, tmp4);       /* C n+1 ^ 2 */
 
173
 
 
174
   /* do the AGM */
 
175
 
 
176
   m_apm_add(tmp1, tmpA0, tmpB0);
 
177
   m_apm_multiply(tmp3, MM_0_5, tmp1);
 
178
 
 
179
   m_apm_multiply(tmp2, tmpA0, tmpB0);
 
180
   m_apm_sqrt(tmpB0, dplaces, tmp2);
 
181
 
 
182
   m_apm_round(tmpA0, dplaces, tmp3);
 
183
 
 
184
   /* end AGM */
 
185
 
 
186
   m_apm_multiply(tmp2, MM_Two, pow_2);
 
187
   m_apm_copy(pow_2, tmp2);
 
188
 
 
189
   m_apm_multiply(tmp1, tmpC2, pow_2);
 
190
   m_apm_add(tmp3, sum, tmp1);
 
191
 
 
192
   if ((tmp1->m_apm_sign == 0) || 
 
193
      ((-2 * tmp1->m_apm_exponent) > tolerance))
 
194
     break;
 
195
 
 
196
   m_apm_round(sum, dplaces, tmp3);
 
197
  }
 
198
 
 
199
m_apm_subtract(tmp4, MM_One, tmp3);
 
200
m_apm_reciprocal(rr, places, tmp4);
 
201
 
 
202
M_restore_stack(9);
 
203
}
 
204
/****************************************************************************/