1 /****************************************************************/
2 /* file analysis.c
3
4 ARIBAS interpreter for Arithmetic
5 Copyright (C) 1996-2010 O.Forster
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20
21 Address of the author
22
23 Otto Forster
24 Math. Institut der LMU
25 Theresienstr. 39
26 D-80333 Muenchen, Germany
27
28 Email forster@mathematik.uni-muenchen.de
29 */
30 /****************************************************************/
31 /*
32 ** analysis.c
33 ** transcendental functions
34 **
35 ** date of last change
36 ** 1995-02-21
37 ** 1997-04-14: changed inipi
38 ** 2010-01-23: bugfix in sin0 and cos0 with floatprec > 2048
39 */
40
41 #include "common.h"
42
43 PUBLIC void inianalys _((void));
44 PUBLIC int lognum _((int prec, numdata *nptr, word2 *hilf));
45 PUBLIC int expnum _((int prec, numdata *nptr, word2 *hilf));
46
47 /*--------------------------------------------------------*/
48 PRIVATE truc expsym, logsym, sqrtsym, sinsym, cossym;
49 PRIVATE truc tansym, atansym, atan2sym, asinsym, acossym;
50 PRIVATE truc pisym;
51
52 PRIVATE truc inipi _((int prec));
53 PRIVATE int Gget1flt _((truc symb, numdata *nptr));
54 PRIVATE truc Fsqrt _((void));
55 PRIVATE truc Fexp _((void));
56 PRIVATE truc Flog _((void));
57 PRIVATE truc Fsin _((void));
58 PRIVATE truc Fcos _((void));
59 PRIVATE truc Ftan _((void));
60 PRIVATE truc Gtrig _((truc symb));
61 PRIVATE truc Fatan _((void));
62 PRIVATE truc Fatan2 _((void));
63 PRIVATE truc Fasin _((void));
64 PRIVATE truc Facos _((void));
65 PRIVATE truc Garcus _((truc symb));
66
67 PRIVATE int atannum _((int prec, numdata *nptr1, numdata *nptr2,
68 word2 *hilf));
69 PRIVATE int atanprep _((int prec, numdata *nptr1, numdata *nptr2,
70 word2 *x, int *segptr));
71 PRIVATE int trignum _((int prec, numdata *nptr, word2 *hilf, truc symb));
72 PRIVATE int expovfl _((numdata *nptr, word2 *hilf));
73 PRIVATE long redmod _((int prec, numdata *nptr, word2 *modul, int modlen,
74 word2 *hilf));
75 PRIVATE int exp0 _((int prec, word2 *x, int n, word2 *z, word2 *hilf));
76 PRIVATE int exp0aux _((word2 *x, int n, unsigned a, int k, word2 *temp));
77 PRIVATE int exp1aux _((word2 *x, int n, unsigned a, int k, word2 *temp));
78 PRIVATE int sin0 _((int prec, word2 *x, int n, word2 *z, word2 *hilf));
79 PRIVATE int cos0 _((int prec, word2 *x, int n, word2 *z, word2 *hilf));
80 PRIVATE int log0 _((int prec, word2 *x, int n, word2 *z, word2 *hilf));
81 PRIVATE unsigned log1_16 _((unsigned x));
82 PRIVATE int atan0 _((int prec, word2 *x, int n, word2 *z, word2 *hilf));
83
84 PRIVATE int curfltprec;
85
86 PRIVATE word2 LOG2DAT[] = /* log(2) */
87 #ifdef FPREC_HIGH
88 {322, 0x154C, 0xB783, 0x64F1, 0xEC3F, 0x53DA, 0x501E, 0xF281, 0x9316,
89 0xDB4A, 0xF949, 0x56C9, 0xC921, 0xBE2E, 0x341C, 0x94F0, 0xF58D, 0xCA8, 0x4A00,
90 0xD287, 0x3D7, 0x554B, 0xE00C, 0x5497, 0x75DF, 0xFB0C, 0x2D06, 0xECA4, 0x850,
91 0xEE6E, 0xEC2F, 0xEF22, 0x5B8A, 0x364F, 0x3C9F, 0x78B6, 0x39CE, 0x897A, 0x8438,
92 0x1E23, 0x3316, 0x52AB, 0xC60C, 0xA6C4, 0x1A63, 0x62B, 0xEDD, 0xE8F7, 0x449F,
93 0x3EA8, 0xC51C, 0x26FA, 0xA415, 0x6425, 0x84E0, 0xF958, 0x767D, 0xC5E5, 0x23FA,
94 0x8A0E, 0xB31D, 0xC0B1, 0xBD0D, 0x3A49, 0x6AB0, 0x85DB, 0xADD8, 0xC8DA, 0xB4AF,
95 0x175E, 0x374E, 0xA892, 0xFFF3, 0xF07A, 0x891E, 0xDEA, 0x2625, 0x8F68, 0x339D,
96 0x9C38, 0x72F1, 0xCECB, 0x45AE, 0xAC9F, 0x7CEB, 0x5F6F, 0x15C0, 0xE761, 0x2096,
97 0x6C47, 0x9D42, 0xFBBD, 0xD18B, 0x972C, 0xC724, 0xBD67, 0x11BB, 0xAB1, 0x38B9,
98 0xA0C2, 0x26FD, 0x4738, 0xAEBD, 0xD24A, 0x696D, 0x61C1, 0xD5E3, 0x2413, 0xC29,
99 0x156E, 0x7487, 0xDC4E, 0x4460, 0x9518, 0x646A, 0x901E, 0x2658, 0xD762, 0x3958,
100 0xD737, 0xCE2, 0xEF2F, 0x207C, 0xC4E9, 0xB61C, 0x2AC5, 0x7D05, 0xBEBA, 0x9BA2,
101 0x5733, 0x1A0C, 0x839, 0xE499, 0x60, 0x302, 0x6AF5, 0x6319, 0x6213, 0xD2F9,
102 0x3D0B, 0x28D5, 0x5C1, 0x86B9, 0xCEE8, 0x2B20, 0x36E0, 0x49F2, 0xF3D9, 0x16FA,
103 0xBBB, 0x2109, 0xC994, 0x83ED, 0x4221, 0xD3C5, 0x8C66, 0x22B8, 0x5E92, 0xA3CF,
104 0x6B1C, 0xFD44, 0x61AF, 0xB982, 0x9538, 0x5C1F, 0x268A, 0x755, 0xFBCF, 0x5177,
105 0x8D6F, 0x4EF9, 0x228A, 0x93D1, 0xA172, 0xDC8E, 0x731C, 0x2554, 0x44A0, 0x889B,
106 0x30AF, 0xE6D3, 0x96D4, 0x9834, 0x8F96, 0xB6C6, 0x5570, 0x73EE, 0x1AE2, 0xA195,
107 0x7598, 0x853D, 0xB365, 0x2DB3, 0x4D16, 0xC18B, 0x5064, 0xB518, 0x5F50, 0xB31B,
108 0x1B2D, 0x735D, 0x78F, 0x6CB1, 0x6C60, 0x3CDB, 0xAE31, 0x7B9D, 0xB1E1, 0x5179,
109 0x955D, 0xD2C, 0x1735, 0xA54, 0xC48, 0x7AA3, 0x5CFE, 0xB601, 0x74D, 0x8E82,
110 0x5E14, 0x7F8A, 0x6A9C, 0xA337, 0x3564, 0x9B33, 0x2566, 0x95D, 0xD1D6, 0x1E0B,
111 0x4C1A, 0x514C, 0x9393, 0x4E65, 0xCCCC, 0xCD33, 0xB479, 0xE732, 0xC943, 0x90E5,
112 0xDB89, 0x775, 0x1746, 0xB396, 0x1400, 0x23DE, 0x7D2E, 0xFA15, 0xFC1E, 0x9D6D,
113 0xEE56, 0x51A2, 0x8FE5, 0x30F8, 0x610D, 0xFB90, 0xFB5B, 0xCA11,
114 #else
115 {66,
116 #endif
117 0x7F4, 0xD5C6,
118 0xF3F, 0x97C5, 0xDA2D, 0xE3A2, 0x2F20, 0xA187, 0x655F, 0x3248, 0x3830, 0xA6BD,
119 0xF5DF, 0x48CA, 0x9D65, 0x87B1, 0x72CE, 0xF74B, 0x7657, 0xA0EC, 0x256F, 0x603B,
120 0xB136, 0x9BC3, 0xB9EA, 0x387E, 0x317C, 0xDA11, 0x1ACB, 0xE8C5, 0x224A, 0xCA16,
121 0x3E96, 0xB825, 0x1169, 0x3B29, 0x2757, 0x2144, 0xC138, 0xAE35, 0xED2E, 0x1B10,
122 0x4AFA, 0x52FB, 0x5595, 0xAC98, 0x6DEB, 0x7620, 0xE7B8, 0xFA2B, 0x8BAA, 0x175B,
123 0x8A0D, 0xB62D, 0x7298, 0x4326, 0x40F3, 0xF6AF, 0x3F2, 0xB398, 0xC9E3, 0x79AB,
124 0xD1CF, 0x17F7, 0xB172};
125
126 PRIVATE word2 PI4THDAT[] = /* pi/4 */
127 #ifdef FPREC_HIGH
128 {322, 0x23A9, 0xE8F3, 0xBEC7, 0xC97F, 0x59E7, 0x1C9E, 0x900B, 0x4031,
129 0xB5A8, 0xC82, 0x4698, 0x702F, 0xD55E, 0xFEF6, 0x6E74, 0xD7CE, 0xF482, 0x1D03,
130 0xD172, 0xEA15, 0xF032, 0x92EC, 0xC64B, 0xCA01, 0x5983, 0xD2BF, 0x378C, 0xF401,
131 0x6FB8, 0xAF42, 0x2BD7, 0x5151, 0x3320, 0x254B, 0xE6CC, 0x1447, 0xDB7F, 0xBB1B,
132 0xCED4, 0x6CBA, 0x44CE, 0x14ED, 0xCF9B, 0xDBEB, 0xDA3E, 0x8918, 0x865A, 0x27B0,
133 0x1797, 0xD831, 0x9027, 0x53ED, 0xB06A, 0x1AE, 0x4130, 0x382F, 0xE5DB, 0x530E,
134 0xAD9E, 0x9406, 0xF8FF, 0x37BD, 0x3DBA, 0x1E76, 0xC975, 0x46DE, 0x6026, 0xDCB2,
135 0xC1D4, 0x7026, 0xD27C, 0xFAB4, 0x36C3, 0x8492, 0x3402, 0x35C9, 0x4DF4, 0xC08F,
136 0x90A6, 0xB7DC, 0x86FF, 0xDDC1, 0x8D8F, 0xEA98, 0x93B4, 0x5AA9, 0xD5B0, 0x9127,
137 0xD006, 0x481C, 0x2170, 0xDD76, 0xB81B, 0xD7AF, 0xCEE2, 0x2970, 0x1F61, 0xE7ED,
138 0x515B, 0xA186, 0x233B, 0xC3A2, 0xA090, 0x964F, 0x99B2, 0xC05D, 0x4E6B, 0x5947,
139 0x287C, 0xCAA6, 0x1FBE, 0xFC14, 0x2E8E, 0x8EF9, 0x4DE, 0xC2DB, 0xDBBB, 0x4CE8,
140 0x2AD4, 0xE9CA, 0x2583, 0xBDA, 0xB615, 0x6834, 0x1A94, 0xE23C, 0x6AF4, 0x2718,
141 0x99C3, 0x5B26, 0xBDBA, 0x9A10, 0x8871, 0xE6D7, 0xA787, 0x3C12, 0x1A72, 0x801,
142 0xA921, 0xD120, 0x4B82, 0x108E, 0xE0FD, 0x5BFC, 0x43DB, 0xAB31, 0x74E5, 0x4FA0,
143 0x8E2, 0x46E2, 0xBAD9, 0x88C0, 0x7709, 0x5D6C, 0x7A61, 0x1757, 0xBBE1, 0x200C,
144 0x177B, 0x2B18, 0x521F, 0x6A64, 0x3EC8, 0x273, 0xD876, 0x864, 0xD98A, 0xFA06,
145 0xF12F, 0xEE6B, 0x1AD2, 0xD226, 0xCEE3, 0x619D, 0x4A25, 0x94E0, 0x1E8C, 0x33D7,
146 0xDB09, 0xAE8C, 0xABF5, 0xE4C7, 0xA6E1, 0xF85, 0xB397, 0xC7D, 0x5D06, 0x7157,
147 0x8AEA, 0xEF0A, 0x58DB, 0x8504, 0xECFB, 0xBA64, 0xDF1C, 0x21AB, 0xA855, 0x7A33,
148 0x450, 0x170D, 0xAD33, 0xC42D, 0x8AAA, 0x8E5A, 0x1572, 0x510, 0x98FA, 0x2618,
149 0x15D2, 0x6AE5, 0xEA95, 0x497C, 0x3995, 0x1718, 0x9558, 0xCBF6, 0xDE2B, 0x52C9,
150 0x6F4C, 0x5DF0, 0xB5C5, 0xA28F, 0xEC07, 0x83A2, 0x9B27, 0x8603, 0x180E, 0x772C,
151 0xE39E, 0xCE3B, 0x2E36, 0x5E46, 0x3290, 0x217C, 0xCA18, 0x6C08, 0xF174, 0x9804,
152 0x4ABC, 0x354E, 0x670C, 0x966D, 0x7096, 0x2907, 0x9ED5, 0x52BB, 0x2085, 0xF356,
153 0x1C62, 0xAD96, 0xDCA3, 0x5D23, 0x8365, 0xCF5F, 0xFD24, 0x3FA8,
154 #else
155 {66,
156 #endif
157 0x6916, 0xD39A,
158 0x1C55, 0x4836, 0x98DA, 0xBF05, 0xA163, 0x7CB8, 0xC200, 0x5B3D, 0xECE4, 0x6651,
159 0x4928, 0x1FE6, 0x7C4B, 0x2411, 0xAE9F, 0x9FA5, 0x5A89, 0x6BFB, 0xEE38, 0xB7ED,
160 0xF406, 0x5CB6, 0xBFF, 0xED6B, 0xA637, 0x42E9, 0xF44C, 0x7EC6, 0x625E, 0xB576,
161 0xE485, 0xC245, 0x6D51, 0x356D, 0x4FE1, 0x1437, 0xF25F, 0xA6D, 0x302B, 0x431B,
162 0xCD3A, 0x19B3, 0xEF95, 0x4DD, 0x8E34, 0x879, 0x514A, 0x9B22, 0x3B13, 0xBEA6,
163 0x20B, 0xCC74, 0x8A67, 0x4E08, 0x2902, 0x1CD1, 0x80DC, 0x628B, 0xC4C6, 0xC234,
164 0x2168, 0xDAA2, 0xC90F};
165
166 PRIVATE word2 ATAN4DAT[] = /* atan(1/4) */
167 #ifdef FPREC_HIGH
168 {322, 0x6B3F, 0x5870, 0xADCF, 0xE549, 0x1C7F, 0x82B2, 0x9F6F, 0x5450,
169 0xE8F, 0xB2E1, 0x95F5, 0x9CE5, 0x998F, 0x64AE, 0x9F2F, 0xDE06, 0xF67, 0xF39,
170 0x111B, 0x3858, 0x25DF, 0x7580, 0x7910, 0xC3A6, 0x5FC3, 0xD1A8, 0xD87A, 0xE0C5,
171 0x559D, 0xDDFE, 0x68B1, 0x81C0, 0x1970, 0x7B17, 0xE38E, 0x6D8A, 0xCB0F, 0x5873,
172 0x6156, 0x89FD, 0xD3FF, 0x5798, 0xD222, 0x62A3, 0x2B89, 0x5A2C, 0x260F, 0xCCE9,
173 0x293D, 0xE516, 0x4EFB, 0xBD90, 0x4872, 0x20F, 0x4EE8, 0xE6DA, 0x95AF, 0x4E41,
174 0x3EB, 0xACB3, 0x4BB, 0xFA89, 0x7C7B, 0xE015, 0xB5E, 0xD3B4, 0xE52, 0x8AAD,
175 0x23A1, 0xC362, 0x18B1, 0x55A6, 0x43E0, 0xD662, 0x5797, 0x3973, 0x13D8, 0x973E,
176 0xF08F, 0x3D2F, 0x90F, 0x24A3, 0x1BEF, 0xE4D5, 0xBD6F, 0x59FF, 0x62B7, 0xDBB9,
177 0x1D86, 0xBB9E, 0xAFC2, 0x44FD, 0xA3B, 0x2F5C, 0x6A0E, 0x683, 0x9A47, 0xFCE5,
178 0xA135, 0x4C26, 0x53DD, 0x2C34, 0x24E0, 0xD837, 0xB3AE, 0xBC2B, 0xD2D4, 0x734B,
179 0x1532, 0x7380, 0x5C33, 0x575, 0xC233, 0x2399, 0xBAE6, 0x7055, 0xF5A3, 0x4C75,
180 0x25A3, 0x2990, 0x818C, 0xAD7F, 0x8904, 0x4750, 0xFF60, 0x5C39, 0xFAA1, 0x9C67,
181 0xF168, 0xE82D, 0x6B10, 0xD0C5, 0xD5B4, 0x9B3A, 0x9454, 0x29C8, 0x4BE0, 0x81CB,
182 0x4EAB, 0xEAF8, 0x579A, 0x338B, 0x4F0D, 0x3013, 0xDBAF, 0xD9A3, 0xF5E0, 0x6F27,
183 0x5EDE, 0xB593, 0xCBCE, 0x3DD7, 0x5BCE, 0x2D5D, 0xA47C, 0x9F21, 0x2301, 0xCBE5,
184 0x11E3, 0x9164, 0xCBDC, 0x79E8, 0x4BA8, 0xCEF1, 0x17C3, 0x56EC, 0x494, 0x7880,
185 0x80E9, 0x31E6, 0x23C7, 0xEA04, 0xE55F, 0xFA48, 0x7482, 0xE6EF, 0x8AEF, 0x170F,
186 0xAB2D, 0xC9E8, 0x9A8F, 0x6030, 0x5040, 0xC199, 0xF2AD, 0x8FA1, 0x8187, 0xF447,
187 0x7EED, 0x6403, 0xC38D, 0xC50D, 0x32FA, 0xEF27, 0x2E51, 0xCE00, 0x3E5A, 0x8682,
188 0x2519, 0xEA4, 0x2861, 0x3E51, 0xB4BB, 0x6267, 0xF41, 0x3D03, 0xFA02, 0x9927,
189 0x4748, 0x5F53, 0x2703, 0x4487, 0xBCA2, 0x63CE, 0x7F56, 0x950C, 0x492E, 0x1875,
190 0xEEBA, 0xE60, 0xDA37, 0x9A95, 0x1906, 0x77BF, 0xE0DF, 0xB528, 0xFCFC, 0xA1D,
191 0xDBFB, 0x2714, 0xFB05, 0x8AA1, 0x6928, 0xB682, 0x65A9, 0x50C7, 0x38B1, 0x4DD3,
192 0x1BED, 0x4652, 0xD4F7, 0xE21E, 0xA8AB, 0xE37, 0x8EC9, 0xAEA6, 0x1541, 0x16BE,
193 0x3E03, 0xA59D, 0x5E8, 0x5574, 0x77B4, 0xA43E, 0x3638, 0x2903,
194 #else
195 {66,
196 #endif
197 0x6F3D, 0xC10E,
198 0xF52C, 0x8AB6, 0x7616, 0xFD3E, 0x172D, 0x2301, 0x3520, 0x463F, 0x9382, 0xA65,
199 0xC31B, 0x1204, 0xD9EA, 0x1DE7, 0x6EE5, 0xF2CB, 0xCF96, 0x738F, 0xF848, 0x9C81,
200 0xEB38, 0x48D5, 0x234, 0xF724, 0x5756, 0x7231, 0xFDA2, 0xA695, 0xF8FC, 0xADD0,
201 0x67A2, 0x60DE, 0xC0E6, 0x70EA, 0x83D0, 0x93E6, 0xC1C7, 0xF089, 0x8F0A, 0xC68F,
202 0xE960, 0xA316, 0xC16F, 0x5D94, 0x9A39, 0x4945, 0x64AE, 0x69D9, 0x2512, 0x9D9F,
203 0xDE8E, 0xE0DA, 0xE22C, 0xEA40, 0x6A9F, 0x85F9, 0x7DE8, 0xE7BD, 0x5B71, 0xBAC5,
204 0x5901, 0xEBF2, 0x3EB6};
205
206 #define LOG2(prec) LOG2DAT + (*LOG2DAT - (prec))
207 #define PI4TH(prec) PI4THDAT + (*PI4THDAT - (prec))
208 #define ATAN4(prec) ATAN4DAT + (*ATAN4DAT - (prec))
209
210 /*----------------------------------------------------------------*/
inianalys()211 PUBLIC void inianalys()
212 {
213 expsym = newsymsig("exp", sFBINARY, (wtruc)Fexp, s_rr);
214 sqrtsym = newsymsig("sqrt", sFBINARY, (wtruc)Fsqrt,s_rr);
215 sinsym = newsymsig("sin", sFBINARY, (wtruc)Fsin, s_rr);
216 cossym = newsymsig("cos", sFBINARY, (wtruc)Fcos, s_rr);
217 tansym = newsymsig("tan", sFBINARY, (wtruc)Ftan, s_rr);
218 logsym = newsymsig("log", sFBINARY, (wtruc)Flog, s_rr);
219 atansym = newsymsig("arctan",sFBINARY, (wtruc)Fatan,s_rr);
220 atan2sym = newsymsig("arctan2",sFBINARY, (wtruc)Fatan2,s_rrr);
221 asinsym = newsymsig("arcsin",sFBINARY, (wtruc)Fasin,s_rr);
222 acossym = newsymsig("arccos",sFBINARY, (wtruc)Facos,s_rr);
223
224 pisym = newsym("pi", sSCONSTANT, inipi(FltPrec[MaxFltLevel]));
225 }
226 /*------------------------------------------------------------------*/
inipi(prec)227 PRIVATE truc inipi(prec)
228 int prec;
229 {
230 numdata acc;
231
232 acc.sign = 0;
233 acc.digits = PI4TH(prec);
234 acc.len = prec;
235 acc.expo = -(prec<<4) + 2; /* (pi/4)*4 */
236 return(mk0float(&acc));
237 }
238 /*----------------------------------------------------------------*/
239 /*
240 ** Setzt nptr->digits = AriBuf
241 ** und holt float aus argStkPtr nach nptr;
242 ** setzt curfltprec und gibt prec = curfltprec + 1 zurueck.
243 ** Im Fehlerfall Rueckgabewert = ERROR
244 */
Gget1flt(symb,nptr)245 PRIVATE int Gget1flt(symb,nptr)
246 truc symb;
247 numdata *nptr;
248 {
249 int prec, type;
250
251 type = chknum(symb,argStkPtr);
252 if(type == aERROR)
253 return(aERROR);
254
255 curfltprec = deffltprec();
256 prec = curfltprec + 1;
257 nptr->digits = AriBuf;
258 getnumtrunc(prec,argStkPtr,nptr);
259 return(prec);
260 }
261 /*----------------------------------------------------------------*/
Fsqrt()262 PRIVATE truc Fsqrt()
263 {
264 numdata acc;
265 word2 *hilf, *x;
266 long m;
267 int prec;
268 int sh, len, rlen;
269
270 prec = Gget1flt(sqrtsym,&acc);
271 if(prec == aERROR)
272 return(brkerr());
273 if(acc.sign) {
274 error(sqrtsym,err_p0num,*argStkPtr);
275 return(brkerr());
276 }
277 if((len = acc.len)) {
278 sh = (curfltprec << 5) + 8;
279 sh -= bitlen(*(AriBuf + len - 1)) + ((len - 1) << 4) - 1;
280 len = shiftarr(AriBuf,len,sh);
281 m = acc.expo - sh;
282 if(m & 1) {
283 len = shlarr(AriBuf,len,1);
284 m -= 1;
285 }
286 x = AriScratch;
287 hilf = AriScratch + aribufSize;
288 cpyarr(AriBuf,len,x);
289 acc.len = bigsqrt(x,len,AriBuf,&rlen,hilf);
290 acc.expo = (m >> 1);
291 }
292 return(mkfloat(curfltprec,&acc));
293 }
294 /*----------------------------------------------------------------*/
Fexp()295 PRIVATE truc Fexp()
296 {
297 numdata acc;
298 int prec;
299 int ret;
300
301 prec = Gget1flt(expsym,&acc);
302 if(prec == aERROR)
303 return(brkerr());
304 ret = expnum(prec,&acc,AriScratch);
305 if(ret == aERROR) {
306 error(expsym,err_ovfl,voidsym);
307 return(brkerr());
308 }
309 return(mkfloat(curfltprec,&acc));
310 }
311 /*----------------------------------------------------------------*/
Flog()312 PRIVATE truc Flog()
313 {
314 numdata acc;
315 int prec;
316 int ret;
317
318 prec = Gget1flt(logsym,&acc);
319 if(prec == aERROR)
320 return(brkerr());
321 ret = lognum(prec,&acc,AriScratch);
322 if(ret == aERROR) {
323 error(logsym,err_pnum,*argStkPtr);
324 return(brkerr());
325 }
326 return(mkfloat(curfltprec,&acc));
327 }
328 /*----------------------------------------------------------------*/
Fsin()329 PRIVATE truc Fsin()
330 {
331 return(Gtrig(sinsym));
332 }
333 /*----------------------------------------------------------------*/
Fcos()334 PRIVATE truc Fcos()
335 {
336 return(Gtrig(cossym));
337 }
338 /*----------------------------------------------------------------*/
Ftan()339 PRIVATE truc Ftan()
340 {
341 return(Gtrig(tansym));
342 }
343 /*----------------------------------------------------------------*/
Gtrig(symb)344 PRIVATE truc Gtrig(symb)
345 truc symb;
346 {
347 numdata acc, acc2;
348 int prec;
349 int ret;
350
351 prec = Gget1flt(symb,&acc);
352 if(prec == aERROR)
353 return(brkerr());
354
355 if(symb == sinsym || symb == cossym)
356 ret = trignum(prec,&acc,AriScratch,symb);
357 else { /* symb == tansym */
358 acc2.digits = AriScratch + aribufSize;
359 cpynumdat(&acc,&acc2);
360 ret = trignum(prec,&acc2,AriScratch,cossym);
361 ret = trignum(prec,&acc,AriScratch,sinsym);
362 ret = divtrunc(prec,&acc,&acc2,AriScratch);
363 }
364 if(ret == aERROR) {
365 error(symb,err_ovfl,voidsym);
366 return(brkerr());
367 }
368 return(mkfloat(curfltprec,&acc));
369 }
370 /*----------------------------------------------------------------*/
Fatan()371 PRIVATE truc Fatan()
372 {
373 truc res;
374
375 ARGpush(constone);
376 res = Fatan2();
377 ARGpop();
378 return(res);
379 }
380 /*----------------------------------------------------------------*/
Fatan2()381 PRIVATE truc Fatan2()
382 {
383 numdata acc1, acc2;
384 word2 *hilf;
385 int type, prec;
386 int ret;
387
388 type = chknums(atan2sym,argStkPtr-1,2);
389 if(type == aERROR)
390 return(brkerr());
391 acc1.digits = AriBuf;
392 acc2.digits = AriScratch;
393 hilf = AriScratch + aribufSize;
394 curfltprec = fltprec(type);
395 prec = curfltprec + 1;
396
397 getnumtrunc(prec,argStkPtr-1,&acc1);
398 getnumtrunc(prec,argStkPtr,&acc2);
399 ret = atannum(prec,&acc1,&acc2,hilf);
400 if(ret == aERROR) {
401 error(atansym,err_ovfl,voidsym);
402 return(brkerr());
403 }
404 return(mkfloat(curfltprec,&acc1));
405 }
406 /*----------------------------------------------------------------*/
Fasin()407 PRIVATE truc Fasin()
408 {
409 return(Garcus(asinsym));
410 }
411 /*----------------------------------------------------------------*/
Facos()412 PRIVATE truc Facos()
413 {
414 return(Garcus(acossym));
415 }
416 /*----------------------------------------------------------------*/
Garcus(symb)417 PRIVATE truc Garcus(symb)
418 truc symb;
419 {
420 numdata acc1, acc2;
421 word2 *x, *y, *z, *hilf;
422 int prec, prec2, ret, cmp, rlen;
423 int n, m;
424
425 prec = Gget1flt(symb,&acc1);
426 if(prec == aERROR)
427 return(brkerr());
428 prec2 = prec + prec;
429
430 x = acc1.digits;
431 acc2.digits = y = AriScratch;
432 z = AriScratch + aribufSize;
433 hilf = z + prec2 + 2;
434 setarr(z,prec2,0);
435 z[prec2] = 1;
436
437 n = alignfix(prec,&acc1);
438 if(n == aERROR || (cmp = cmparr(x,n,z+prec,prec+1)) > 0) {
439 error(symb,err_range,*argStkPtr);
440 return(brkerr());
441 }
442 if(cmp == 0) /* abs(x) = 1 */
443 int2numdat(0,&acc2);
444 else if(n == 0) /* x = 0 */
445 int2numdat(1,&acc2);
446 else {
447 m = multbig(x,n,x,n,y,hilf);
448 m = sub1arr(y,m,z,prec2+1); /* z = 1 - x*x */
449 m = bigsqrt(y,m,z,&rlen,hilf);
450 cpyarr(z,m,y); /* y = sqrt(1 - x*x) */
451 acc2.len = m;
452 acc2.sign = 0;
453 acc2.expo = -(prec<<4);
454 }
455 if(symb == asinsym)
456 ret = atannum(prec,&acc1,&acc2,hilf);
457 else {
458 ret = atannum(prec,&acc2,&acc1,hilf);
459 cpynumdat(&acc2,&acc1);
460 }
461 if(ret == aERROR) {
462 error(symb,err_ovfl,voidsym);
463 return(brkerr());
464 }
465 return(mkfloat(curfltprec,&acc1));
466 }
467 /*----------------------------------------------------------------*/
468 /*
469 ** Ersetzt nptr1 durch den Arcus-Tangens des Quotienten
470 ** aus nptr1 und nptr2
471 */
atannum(prec,nptr1,nptr2,hilf)472 PRIVATE int atannum(prec,nptr1,nptr2,hilf)
473 int prec;
474 numdata *nptr1, *nptr2;
475 word2 *hilf;
476 {
477 word2 *x, *z;
478 int seg, s, sign, m, n;
479
480 x = hilf;
481 hilf += prec << 1;
482 z = nptr1->digits;
483
484 n = atanprep(prec,nptr1,nptr2,x,&seg);
485 if(n < 0)
486 return(n); /* aERROR */
487 m = atan0(prec,x,n,z,hilf);
488 nptr1->sign = sign = (seg < 0 ? MINUSBYTE : 0);
489 if(sign)
490 seg = -seg - 2;
491 s = ((seg + 2) >> 1) & 0xFFFE;
492 if(s) {
493 n = multarr(PI4TH(prec),prec,s,hilf);
494 if(seg & 2) {
495 m = sub1arr(z,m,hilf,n);
496 }
497 else
498 m = addarr(z,m,hilf,n);
499 }
500 nptr1->len = m;
501 nptr1->expo = -(prec << 4);
502 return(m);
503 }
504 /*----------------------------------------------------------------*/
505 /*
506 ** Falls Zahl nptr2 groesser als Zahl in nptr1,
507 ** wird nptr1 durch nptr2 dividiert,
508 ** andernfalls wird nptr2 durch nptr1 dividiert.
509 ** Ist len der Rueckgabewert so erhaelt man
510 ** Ergebnis = (x,len) * (2**16)**(-prec)
511 ** Platz x muss genuegend lang sein
512 ** Arbeitet destruktiv auf nptr1 und nptr2 !!!
513 */
atanprep(prec,nptr1,nptr2,x,segp)514 PRIVATE int atanprep(prec,nptr1,nptr2,x,segp)
515 int prec;
516 numdata *nptr1, *nptr2;
517 word2 *x;
518 int *segp;
519 {
520 numdata *temp;
521 long diff, sh;
522 int n, m, sign1, sign2;
523 int cmp;
524
525 sign1 = nptr1->sign;
526 sign2 = nptr2->sign;
527 n = alignfloat(prec+1,nptr1);
528 m = alignfloat(prec+1,nptr2);
529 if(!n) {
530 if(!m)
531 return(aERROR);
532 else {
533 *segp = (sign2 ? 8 : 0);
534 return(0);
535 }
536 }
537 else if(!m) {
538 *segp = (sign1 ? -4 : 4);
539 return(0);
540 }
541 if(!sign1)
542 *segp = (sign2 ? 4 : 0);
543 else
544 *segp = (sign2 ? -8 : -4);
545 if(sign1 != sign2) {
546 temp = nptr1;
547 nptr1 = nptr2;
548 nptr2 = temp;
549 }
550
551 diff = nptr1->expo - nptr2->expo;
552 if(diff == 0)
553 cmp = (nptr1->digits[prec] > nptr2->digits[prec]);
554 else
555 cmp = (diff > 0);
556
557 if(cmp) { /* nptr1 groesser */
558 *segp += 2;
559 temp = nptr1;
560 nptr1 = nptr2;
561 nptr2 = temp;
562 }
563 n = divtrunc(prec,nptr1,nptr2,x);
564 cpyarr(nptr1->digits,n,x);
565 sh = (prec << 4) + nptr1->expo;
566 n = lshiftarr(x,n,sh);
567 return(n);
568 }
569 /*----------------------------------------------------------------*/
trignum(prec,nptr,hilf,symb)570 PRIVATE int trignum(prec,nptr,hilf,symb)
571 int prec;
572 numdata *nptr;
573 word2 *hilf;
574 truc symb;
575 {
576 word2 *x;
577 int m, n;
578
579 m = redmod(prec,nptr,PI4TH(prec),prec,hilf);
580 if(m & 1) {
581 nptr->len =
582 sub1arr(nptr->digits,nptr->len,PI4TH(prec),prec);
583 }
584 x = hilf;
585 hilf += prec;
586 n = nptr->len;
587 cpyarr(nptr->digits,n,x);
588 if(symb == cossym)
589 m += 2;
590 if((m+1) & 2)
591 nptr->len = cos0(prec,x,n,nptr->digits,hilf);
592 else
593 nptr->len = sin0(prec,x,n,nptr->digits,hilf);
594 nptr->sign = (m & 4 ? MINUSBYTE : 0);
595 nptr->expo = -(prec << 4);
596 return(nptr->len);
597 }
598 /*----------------------------------------------------------------*/
599 /*
600 ** Die durch nptr dargestellte Zahl z wird ersetzt durch exp(z)
601 ** Bei overflow wird aERROR zurueckgegeben
602 */
expnum(prec,nptr,hilf)603 PUBLIC int expnum(prec,nptr,hilf)
604 int prec;
605 numdata *nptr;
606 word2 *hilf;
607 {
608 word2 *x;
609 long m;
610 int ovfl, n;
611
612 ovfl = expovfl(nptr,hilf);
613 if(ovfl > 0)
614 return(aERROR);
615 else if(ovfl < 0) {
616 int2numdat(0,nptr);
617 return(0);
618 }
619 m = redmod(prec+1,nptr,LOG2(prec+1),prec+1,hilf);
620 x = hilf;
621 hilf += prec;
622 n = nptr->len - 1;
623 cpyarr(nptr->digits+1,n,x);
624 nptr->len = exp0(prec,x,n,nptr->digits,hilf);
625 nptr->sign = 0;
626 nptr->expo = m - (prec << 4);
627 return(nptr->len);
628 }
629 /*----------------------------------------------------------------*/
expovfl(nptr,hilf)630 PRIVATE int expovfl(nptr,hilf)
631 numdata *nptr;
632 word2 *hilf;
633 {
634 int n;
635
636 if(nptr->expo <= 1) {
637 cpyarr(nptr->digits,nptr->len,hilf);
638 n = lshiftarr(hilf,nptr->len,nptr->expo);
639 if(n == 0 || (n <= 2 && big2long(hilf,n) < exprange))
640 return(0);
641 }
642 return(nptr->sign ? -1 : 1);
643 }
644 /*----------------------------------------------------------------*/
645 /*
646 ** Die durch nptr gegebene Zahl wird destruktiv dargestellt als
647 ** ((nptr->digits,nptr->len) + ret*(modul,modlen)) * (2**16)**(-prec),
648 ** wobei ret der Rueckgabewert ist.
649 ** Die Zahl (nptr->digits,nptr->len) ist nicht negativ und < (2**16)**prec
650 ** ret kann auch negativ sein
651 */
redmod(prec,nptr,modul,modlen,hilf)652 PRIVATE long redmod(prec,nptr,modul,modlen,hilf)
653 int prec, modlen;
654 numdata *nptr;
655 word2 *modul, *hilf;
656 {
657 word2 *x, *quot;
658 word4 u;
659 long ret;
660 int n, len, rlen;
661
662 x = nptr->digits;
663 len = lshiftarr(x,nptr->len,(prec << 4) + nptr->expo);
664 quot = hilf + prec + 1;
665 n = divbig(x,len,modul,modlen,quot,&rlen,hilf);
666 if(n <= 2)
667 u = big2long(quot,n);
668 if(n > 2 || u >= 0x80000000) {
669 error(scratch("redmod"),err_ovfl,voidsym);
670 return(LONGERROR);
671 }
672 else
673 ret = u;
674 if(nptr->sign) {
675 ret = -ret;
676 if(rlen) {
677 rlen = sub1arr(x,rlen,modul,modlen);
678 ret--;
679 }
680 }
681 nptr->len = rlen;
682 return(ret);
683 }
684 /*----------------------------------------------------------------*/
lognum(prec,nptr,hilf)685 PUBLIC int lognum(prec,nptr,hilf)
686 int prec;
687 numdata *nptr;
688 word2 *hilf;
689 {
690 word2 *x, *z;
691 word2 aa[2];
692 word4 u;
693 long expo;
694 int m, n, len;
695
696 if(nptr->sign || nptr->len == 0)
697 return(aERROR);
698 x = nptr->digits;
699 z = hilf;
700 hilf += prec + 2;
701
702 normfloat(prec,nptr);
703 n = shlarr(x,prec,1);
704 expo = nptr->expo + (prec << 4) - 1;
705 len = log0(prec,x,n,z,hilf);
706 cpyarr(z,len,x);
707 if(expo) {
708 u = (expo > 0 ? expo : -expo);
709 m = long2big(u,aa);
710 n = multbig(LOG2(prec),prec,aa,m,z,hilf);
711 if(expo > 0) {
712 len = addarr(x,len,z,n);
713 nptr->sign = 0;
714 }
715 else if(cmparr(x,len,z,n) >= 0) {
716 len = subarr(x,len,z,n);
717 nptr->sign = 0;
718 }
719 else {
720 len = sub1arr(x,len,z,n);
721 nptr->sign = MINUSBYTE;
722 }
723 }
724 if(len == 0)
725 int2numdat(0,nptr);
726 else {
727 nptr->len = len;
728 nptr->expo = -(prec << 4);
729 }
730 return(len);
731 }
732 /*----------------------------------------------------------------*/
733 /*
734 ** Berechnet die Exponentialfunktion von (x,n) * (2**16)**(-prec)
735 ** Ist len der Rueckgabewert, so erhaelt man
736 ** Resultat = (z,len) * (2**16)**(-prec)
737 ** Es wird vorausgesetzt, dass n <= prec ist
738 ** Platz hilf muss mindestens prec + 2 lang sein
739 ** Platz z muss mindestens prec + 1 lang sein
740 */
741
exp0(prec,x,n,z,hilf)742 PRIVATE int exp0(prec,x,n,z,hilf)
743 int prec, n;
744 word2 *x, *z, *hilf;
745 {
746 int len;
747
748 setarr(z,prec,0);
749 z[prec] = 1;
750 len = prec + 1;
751 while(--n >= 0)
752 len = exp0aux(z,len,x[n],prec-n,hilf);
753 return(len);
754 }
755 /*----------------------------------------------------------------*/
756 /*
757 ** Multipliziert (x,n) mit exp(a * (2**16)**(-k))
758 ** Platz temp muss mindestens n + 2 lang sein
759 ** Arbeitet destruktiv auf x !!!
760 */
exp0aux(x,n,a,k,temp)761 PRIVATE int exp0aux(x,n,a,k,temp)
762 word2 *x, *temp;
763 unsigned a;
764 int n, k;
765 {
766 int i, m;
767 word2 rest;
768
769 if(a == 0)
770 return(n);
771 temp++;
772 cpyarr(x,n,temp);
773 m = n;
774 for(i=1; m>k; i++) {
775 m = multarr(temp+k-1,m-k+1,a,temp-1) - 1;
776 m = divarr(temp,m,i,&rest);
777 n = addarr(x,n,temp,m);
778 }
779 return(n);
780 }
781 /*----------------------------------------------------------------*/
782 /*
783 ** Dividiert (x,n) durch exp(a * (2**16)**(-k))
784 ** Platz temp muss mindestens n + 2 lang sein
785 ** Arbeitet destruktiv auf x !!!
786 */
exp1aux(x,n,a,k,temp)787 PRIVATE int exp1aux(x,n,a,k,temp)
788 word2 *x, *temp;
789 unsigned a;
790 int n, k;
791 {
792 int i, m;
793 word2 rest;
794
795 if(a == 0)
796 return(n);
797 temp++;
798 cpyarr(x,n,temp);
799 m = n;
800 for(i=1; m>k; i++) {
801 m = multarr(temp+k-1,m-k+1,a,temp-1) - 1;
802 m = divarr(temp,m,i,&rest);
803 n = (i&1 ? subarr(x,n,temp,m) : addarr(x,n,temp,m));
804 }
805 return(n);
806 }
807 /*----------------------------------------------------------------*/
808 /*
809 ** Berechnet die Funktion sin(x) von (x,n) * (2**16)**(-prec)
810 ** Ist len der Rueckgabewert so erhaelt man
811 ** Resultat = (z,len) * (2**16)**(-prec)
812 */
sin0(prec,x,n,z,hilf)813 PRIVATE int sin0(prec,x,n,z,hilf)
814 int prec, n;
815 word2 *x, *z, *hilf;
816 {
817 word2 *temp, *temp1, *x2;
818 unsigned i;
819 int len, m, m2;
820
821 m = (prec + 1) << 1;
822 temp = hilf + m;
823 temp1 = temp + m;
824 x2 = temp1 + m;
825 cpyarr(x,n,temp);
826 m = n;
827 cpyarr(temp,m,z);
828 len = m;
829 m2 = multfix(prec,x,n,x,n,x2,hilf);
830 for(i=2; m>0; i+=2) {
831 m = multfix(prec,x2,m2,temp,m,temp1,hilf);
832 cpyarr(temp1,m,temp);
833 if(i<=255)
834 m = divarr(temp,m,i*(i+1),hilf);
835 else {
836 #ifdef M64_32
837 m = div4arr(temp,m,i*(i+1),hilf);
838 #else
839 m = divarr(temp,m,i,hilf);
840 m = divarr(temp,m,i+1,hilf);
841 #endif
842 }
843 if(i & 2)
844 len = subarr(z,len,temp,m);
845 else
846 len = addarr(z,len,temp,m);
847 }
848 return(len);
849 }
850 /*----------------------------------------------------------------*/
851 /*
852 ** Berechnet die Funktion cos(x) von (x,n) * (2**16)**(-prec)
853 ** Ist len der Rueckgabewert so erhaelt man
854 ** Resultat = (z,len) * (2**16)**(-prec)
855 */
cos0(prec,x,n,z,hilf)856 PRIVATE int cos0(prec,x,n,z,hilf)
857 int prec, n;
858 word2 *x, *z, *hilf;
859 {
860 word2 *temp, *temp1, *x2;
861 int len, m, m2;
862 unsigned i;
863
864 m = (prec + 1) << 1;
865 temp = hilf + m;
866 temp1 = temp + m;
867 x2 = temp1 + m;
868 for(i=0; i<prec; i++)
869 z[i] = 0;
870 z[prec] = 1;
871 len = prec + 1;
872 m2 = multfix(prec,x,n,x,n,x2,hilf);
873 cpyarr(x2,m2,temp);
874 m = m2;
875 for(i=1; ; i+=2) {
876 if(i<=255)
877 m = divarr(temp,m,i*(i+1),hilf);
878 else {
879 #ifdef M64_32
880 m = div4arr(temp,m,i*(i+1),hilf);
881 #else
882 m = divarr(temp,m,i,hilf);
883 m = divarr(temp,m,i+1,hilf);
884 #endif
885 }
886 if(i & 2)
887 len = addarr(z,len,temp,m);
888 else
889 len = subarr(z,len,temp,m);
890 if(m <= 1)
891 break;
892 m = multfix(prec,x2,m2,temp,m,temp1,hilf);
893 cpyarr(temp1,m,temp);
894 }
895 return(len);
896 }
897 /*----------------------------------------------------------------*/
898 /*
899 ** Berechnet die Funktion log(x) von (x,n) * (2**16)**(-prec)
900 ** Ist len der Rueckgabewert so erhaelt man
901 ** Resultat = (z,len) * (2**16)**(-prec)
902 ** Es wird vorausgesetzt, dass (2**16)**prec <= (x,n) < 2*(2**16)**prec
903 ** Platz hilf muss prec + 2 lang sein
904 ** Arbeitet destruktiv auf x !!!
905 */
log0(prec,x,n,z,hilf)906 PRIVATE int log0(prec,x,n,z,hilf)
907 int prec, n;
908 word2 *x, *z, *hilf;
909 {
910 unsigned u;
911 int k, len;
912
913 setarr(z,prec,0);
914
915 k = prec - 1;
916 u = log1_16(x[k]);
917 len = 0;
918
919 while(1) {
920 n = exp1aux(x,n,u,prec-k,hilf);
921 while(n <= prec) {
922 u--;
923 n = exp0aux(x,n,1,prec-k,hilf);
924 }
925 len = incarr(z+k,len,u);
926 u = x[k];
927 if(u == 0)
928 u = x[k-1];
929 else if(u == 1)
930 u = 0xFFFF;
931 else
932 continue;
933 if(--k < 0)
934 break;
935 if(len)
936 len++;
937 }
938 while(len > 0 && z[len-1] == 0)
939 len--;
940 return(len);
941 }
942 /*----------------------------------------------------------------*/
943 /*
944 ** berechnet 2**16 * log(1 + x*(2**-16))
945 ** gerundet auf ganze Zahl
946 */
log1_16(x)947 PRIVATE unsigned log1_16(x)
948 unsigned x;
949 {
950 static word4 logtab[] = {0x49A58844,0x108598B5,0x04081596,0x01008055,
951 0x00400801,0x00100080,0x00040008,0x00010000};
952 /* log(2**m/(2**m - 1)) * 2**32 fuer m = 2,4,...,16 */
953
954 word4 *logptr;
955 word4 xx,d,z;
956 int m;
957
958 logptr = logtab;
959 z = 0;
960 xx = x;
961 xx <<= 16;
962 m = 1;
963 while(m < 16) {
964 d = xx;
965 d >>= 1;
966 d += 0x80000000;
967 d >>= m;
968 if(xx >= d) {
969 xx -= d;
970 z += *logptr;
971 }
972 else {
973 m += 2;
974 logptr++;
975 }
976 }
977 return(z >> 16);
978 }
979 /*----------------------------------------------------------------*/
980 /*
981 ** Berechnet die Funktion atan(x) von (x,n) * (2**16)**(-prec)
982 ** Ist len der Rueckgabewert so erhaelt man
983 ** Resultat = (z,len) * (2**16)**(-prec)
984 ** Es wird vorausgesetzt, dass (x,n) < (1/2)*(2**16)**prec
985 ** Platz hilf muss 8*(prec+1) lang sein
986 */
atan0(prec,x,n,z,hilf)987 PRIVATE int atan0(prec,x,n,z,hilf)
988 int prec, n;
989 word2 *x, *z, *hilf;
990 {
991 word2 *u, *v, *y, *temp, *temp1, *x2, *arctan;
992 int i, k, m, m1, m2, len;
993
994 i = (prec+1) << 1;
995 u = temp = hilf;
996 v = temp1 = temp + i;
997 y = x2 = temp1 + i;
998 hilf = x2 + i;
999
1000 len = 0; /* z = 0 */
1001 setarr(u,prec-1,0);
1002 u[prec-1] = 0x4000; /* u = 1/4 */
1003 while(cmparr(x,n,u,prec) >= 0) { /* x = (x-u)/(1+u*x) */
1004 arctan = ATAN4(prec);
1005 cpyarr(x,n,v);
1006 k = shrarr(v,n,2);
1007 setarr(v+k,prec-k,0);
1008 v[prec] = 1;
1009 n = subarr(x,n,u,prec);
1010 n = divfix(prec,x,n,v,prec+1,y,hilf);
1011 cpyarr(y,n,x);
1012 len = addarr(z,len,arctan,prec);
1013 }
1014 len = addarr(z,len,x,n);
1015 cpyarr(x,n,temp);
1016 m = n;
1017 m2 = multfix(prec,x,n,x,n,x2,hilf);
1018 for(i=3; m>1; i+=2) {
1019 m = multfix(prec,x2,m2,temp,m,temp1,hilf);
1020 cpyarr(temp1,m,temp);
1021 m1 = divarr(temp1,m,i,hilf);
1022 if(i & 2)
1023 len = subarr(z,len,temp1,m1);
1024 else
1025 len = addarr(z,len,temp1,m1);
1026 }
1027 return(len);
1028 }
1029 /******************************************************************/
1030