1*59d799daSIngo Weinhold
2*59d799daSIngo Weinhold /*
3*59d799daSIngo Weinhold * M_APM - mapm_lg4.c
4*59d799daSIngo Weinhold *
5*59d799daSIngo Weinhold * Copyright (C) 2003 - 2007 Michael C. Ring
6*59d799daSIngo Weinhold *
7*59d799daSIngo Weinhold * Permission to use, copy, and distribute this software and its
8*59d799daSIngo Weinhold * documentation for any purpose with or without fee is hereby granted,
9*59d799daSIngo Weinhold * provided that the above copyright notice appear in all copies and
10*59d799daSIngo Weinhold * that both that copyright notice and this permission notice appear
11*59d799daSIngo Weinhold * in supporting documentation.
12*59d799daSIngo Weinhold *
13*59d799daSIngo Weinhold * Permission to modify the software is granted. Permission to distribute
14*59d799daSIngo Weinhold * the modified code is granted. Modifications are to be distributed by
15*59d799daSIngo Weinhold * using the file 'license.txt' as a template to modify the file header.
16*59d799daSIngo Weinhold * 'license.txt' is available in the official MAPM distribution.
17*59d799daSIngo Weinhold *
18*59d799daSIngo Weinhold * This software is provided "as is" without express or implied warranty.
19*59d799daSIngo Weinhold */
20*59d799daSIngo Weinhold
21*59d799daSIngo Weinhold /*
22*59d799daSIngo Weinhold * $Id: mapm_lg4.c,v 1.3 2007/12/03 01:43:32 mike Exp $
23*59d799daSIngo Weinhold *
24*59d799daSIngo Weinhold * This file contains the LOG_NEAR_1 function.
25*59d799daSIngo Weinhold *
26*59d799daSIngo Weinhold * $Log: mapm_lg4.c,v $
27*59d799daSIngo Weinhold * Revision 1.3 2007/12/03 01:43:32 mike
28*59d799daSIngo Weinhold * Update license
29*59d799daSIngo Weinhold *
30*59d799daSIngo Weinhold * Revision 1.2 2003/06/02 18:08:45 mike
31*59d799daSIngo Weinhold * tweak decimal places and add comments
32*59d799daSIngo Weinhold *
33*59d799daSIngo Weinhold * Revision 1.1 2003/06/02 17:27:26 mike
34*59d799daSIngo Weinhold * Initial revision
35*59d799daSIngo Weinhold */
36*59d799daSIngo Weinhold
37*59d799daSIngo Weinhold #include "m_apm_lc.h"
38*59d799daSIngo Weinhold
39*59d799daSIngo Weinhold /****************************************************************************/
40*59d799daSIngo Weinhold /*
41*59d799daSIngo Weinhold calculate log (1 + x) with the following series:
42*59d799daSIngo Weinhold
43*59d799daSIngo Weinhold x
44*59d799daSIngo Weinhold y = ----- ( |y| < 1 )
45*59d799daSIngo Weinhold x + 2
46*59d799daSIngo Weinhold
47*59d799daSIngo Weinhold
48*59d799daSIngo Weinhold [ 1 + y ] y^3 y^5 y^7
49*59d799daSIngo Weinhold log [-------] = 2 * [ y + --- + --- + --- ... ]
50*59d799daSIngo Weinhold [ 1 - y ] 3 5 7
51*59d799daSIngo Weinhold
52*59d799daSIngo Weinhold */
M_log_near_1(M_APM rr,int places,M_APM xx)53*59d799daSIngo Weinhold void M_log_near_1(M_APM rr, int places, M_APM xx)
54*59d799daSIngo Weinhold {
55*59d799daSIngo Weinhold M_APM tmp0, tmp1, tmp2, tmpS, term;
56*59d799daSIngo Weinhold int tolerance, dplaces, local_precision;
57*59d799daSIngo Weinhold long m1;
58*59d799daSIngo Weinhold
59*59d799daSIngo Weinhold tmp0 = M_get_stack_var();
60*59d799daSIngo Weinhold tmp1 = M_get_stack_var();
61*59d799daSIngo Weinhold tmp2 = M_get_stack_var();
62*59d799daSIngo Weinhold tmpS = M_get_stack_var();
63*59d799daSIngo Weinhold term = M_get_stack_var();
64*59d799daSIngo Weinhold
65*59d799daSIngo Weinhold tolerance = xx->m_apm_exponent - (places + 6);
66*59d799daSIngo Weinhold dplaces = (places + 12) - xx->m_apm_exponent;
67*59d799daSIngo Weinhold
68*59d799daSIngo Weinhold m_apm_add(tmp0, xx, MM_Two);
69*59d799daSIngo Weinhold m_apm_divide(tmpS, (dplaces + 6), xx, tmp0);
70*59d799daSIngo Weinhold
71*59d799daSIngo Weinhold m_apm_copy(term, tmpS);
72*59d799daSIngo Weinhold m_apm_multiply(tmp0, tmpS, tmpS);
73*59d799daSIngo Weinhold m_apm_round(tmp2, (dplaces + 6), tmp0);
74*59d799daSIngo Weinhold
75*59d799daSIngo Weinhold m1 = 3L;
76*59d799daSIngo Weinhold
77*59d799daSIngo Weinhold while (TRUE)
78*59d799daSIngo Weinhold {
79*59d799daSIngo Weinhold m_apm_multiply(tmp0, term, tmp2);
80*59d799daSIngo Weinhold
81*59d799daSIngo Weinhold if ((tmp0->m_apm_exponent < tolerance) || (tmp0->m_apm_sign == 0))
82*59d799daSIngo Weinhold break;
83*59d799daSIngo Weinhold
84*59d799daSIngo Weinhold local_precision = dplaces + tmp0->m_apm_exponent;
85*59d799daSIngo Weinhold
86*59d799daSIngo Weinhold if (local_precision < 20)
87*59d799daSIngo Weinhold local_precision = 20;
88*59d799daSIngo Weinhold
89*59d799daSIngo Weinhold m_apm_set_long(tmp1, m1);
90*59d799daSIngo Weinhold m_apm_round(term, local_precision, tmp0);
91*59d799daSIngo Weinhold m_apm_divide(tmp0, local_precision, term, tmp1);
92*59d799daSIngo Weinhold m_apm_add(tmp1, tmpS, tmp0);
93*59d799daSIngo Weinhold m_apm_copy(tmpS, tmp1);
94*59d799daSIngo Weinhold m1 += 2;
95*59d799daSIngo Weinhold }
96*59d799daSIngo Weinhold
97*59d799daSIngo Weinhold m_apm_multiply(tmp0, MM_Two, tmpS);
98*59d799daSIngo Weinhold m_apm_round(rr, places, tmp0);
99*59d799daSIngo Weinhold
100*59d799daSIngo Weinhold M_restore_stack(5); /* restore the 5 locals we used here */
101*59d799daSIngo Weinhold }
102*59d799daSIngo Weinhold /****************************************************************************/
103