1 /*
2 * zfunc - extended precision integral arithmetic non-primitive routines
3 *
4 * Copyright (C) 1999-2007,2021 David I. Bell, Landon Curt Noll
5 * and Ernest Bowen
6 *
7 * Primary author: David I. Bell
8 *
9 * Calc is open software; you can redistribute it and/or modify it under
10 * the terms of the version 2.1 of the GNU Lesser General Public License
11 * as published by the Free Software Foundation.
12 *
13 * Calc is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
16 * Public License for more details.
17 *
18 * A copy of version 2.1 of the GNU Lesser General Public License is
19 * distributed with calc under the filename COPYING-LGPL. You should have
20 * received a copy with calc; if not, write to Free Software Foundation, Inc.
21 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 *
23 * Under source code control: 1990/02/15 01:48:27
24 * File existed as early as: before 1990
25 *
26 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
27 */
28
29
30 #include "zmath.h"
31 #include "alloc.h"
32
33
34 #include "banned.h" /* include after system header <> includes */
35
36
37 ZVALUE _tenpowers_[TEN_MAX+1]; /* table of 10^2^n */
38
39 STATIC long *power10 = NULL;
40 STATIC int max_power10_exp = 0;
41
42 /*
43 * given:
44 *
45 * unsigned long x
46 * or: unsigned long long x
47 * or: long x and x >= 0
48 * or: long long x and x >= 0
49 *
50 * If issq_mod4k[x & 0xfff] == 0, then x cannot be a perfect square
51 * else x might be a perfect square.
52 */
53 STATIC USB8 issq_mod4k[1<<12] = {
54 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
55 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
56 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
57 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
58 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
59 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
60 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
61 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
62 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
63 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
64 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
65 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
66 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
67 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
68 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
69 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
70 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
71 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
72 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
73 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
74 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
75 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
76 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
77 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
78 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
79 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
80 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
81 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
82 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
83 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
84 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
85 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
86 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
87 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
88 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
89 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
90 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
91 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
92 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
93 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
94 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
95 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
96 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
97 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
98 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
99 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
100 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
101 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
102 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
103 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
104 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
105 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
106 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
107 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
108 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
109 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
110 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
111 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
112 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
113 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
114 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
115 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
116 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
117 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
118 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
119 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
120 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
121 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
122 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
123 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
124 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
125 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
126 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
127 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
128 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
129 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
130 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
131 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
132 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
133 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
134 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
135 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
136 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
137 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
138 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
139 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
140 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
141 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
142 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
143 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
144 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
145 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
146 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
147 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
148 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
149 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
150 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
151 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
152 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
153 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
154 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
155 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
156 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
157 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
158 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
159 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
160 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
161 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
162 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
163 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
164 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
165 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
166 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
167 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
168 1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
169 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
170 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
171 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
172 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
173 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
174 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
175 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
176 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
177 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
178 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
179 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
180 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
181 0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
182 };
183
184
185 /*
186 * Compute the factorial of a number.
187 */
188 void
zfact(ZVALUE z,ZVALUE * dest)189 zfact(ZVALUE z, ZVALUE *dest)
190 {
191 long ptwo; /* count of powers of two */
192 long n; /* current multiplication value */
193 long m; /* reduced multiplication value */
194 long mul; /* collected value to multiply by */
195 ZVALUE res, temp;
196
197 if (zisneg(z)) {
198 math_error("Negative argument for factorial");
199 /*NOTREACHED*/
200 }
201 if (zge31b(z)) {
202 math_error("Very large factorial");
203 /*NOTREACHED*/
204 }
205 n = ztolong(z);
206 ptwo = 0;
207 mul = 1;
208 res = _one_;
209 /*
210 * Multiply numbers together, but squeeze out all powers of two.
211 * We will put them back in at the end. Also collect multiple
212 * numbers together until there is a risk of overflow.
213 */
214 for (; n > 1; n--) {
215 for (m = n; ((m & 0x1) == 0); m >>= 1)
216 ptwo++;
217 if (mul <= MAXLONG/m) {
218 mul *= m;
219 continue;
220 }
221 zmuli(res, mul, &temp);
222 zfree(res);
223 res = temp;
224 mul = m;
225 }
226 /*
227 * Multiply by the remaining value, then scale result by
228 * the proper power of two.
229 */
230 if (mul > 1) {
231 zmuli(res, mul, &temp);
232 zfree(res);
233 res = temp;
234 }
235 zshift(res, ptwo, &temp);
236 zfree(res);
237 *dest = temp;
238 }
239
240
241 /*
242 * Compute the permutation function M! / (M - N)!.
243 */
244 void
zperm(ZVALUE z1,ZVALUE z2,ZVALUE * res)245 zperm(ZVALUE z1, ZVALUE z2, ZVALUE *res)
246 {
247 SFULL count;
248 ZVALUE cur, tmp, ans;
249
250 if (zisneg(z1) || zisneg(z2)) {
251 math_error("Negative argument for permutation");
252 /*NOTREACHED*/
253 }
254 if (zrel(z1, z2) < 0) {
255 math_error("Second arg larger than first in permutation");
256 /*NOTREACHED*/
257 }
258 if (zge31b(z2)) {
259 math_error("Very large permutation");
260 /*NOTREACHED*/
261 }
262 count = ztolong(z2);
263 zcopy(z1, &ans);
264 zsub(z1, _one_, &cur);
265 while (--count > 0) {
266 zmul(ans, cur, &tmp);
267 zfree(ans);
268 ans = tmp;
269 zsub(cur, _one_, &tmp);
270 zfree(cur);
271 cur = tmp;
272 }
273 zfree(cur);
274 *res = ans;
275 }
276
277 /*
278 * docomb evaluates binomial coefficient when z1 >= 0, z2 >= 0
279 */
280 S_FUNC int
docomb(ZVALUE z1,ZVALUE z2,ZVALUE * res)281 docomb(ZVALUE z1, ZVALUE z2, ZVALUE *res)
282 {
283 ZVALUE ans;
284 ZVALUE mul, div, temp;
285 FULL count, i;
286 #if BASEB == 16
287 HALF dh[2];
288 #else
289 HALF dh[1];
290 #endif
291
292 if (zrel(z2, z1) > 0)
293 return 0;
294 zsub(z1, z2, &temp);
295
296 if (zge31b(z2) && zge31b(temp)) {
297 zfree(temp);
298 return -2;
299 }
300 if (zrel(temp, z2) < 0)
301 count = ztofull(temp);
302 else
303 count = ztofull(z2);
304 zfree(temp);
305 if (count == 0)
306 return 1;
307 if (count == 1)
308 return 2;
309 div.sign = 0;
310 div.v = dh;
311 div.len = 1;
312 zcopy(z1, &mul);
313 zcopy(z1, &ans);
314 for (i = 2; i <= count; i++) {
315 #if BASEB == 16
316 dh[0] = (HALF)(i & BASE1);
317 dh[1] = (HALF)(i >> BASEB);
318 div.len = 1 + (dh[1] != 0);
319 #else
320 dh[0] = (HALF) i;
321 #endif
322 zsub(mul, _one_, &temp);
323 zfree(mul);
324 mul = temp;
325 zmul(ans, mul, &temp);
326 zfree(ans);
327 zquo(temp, div, &ans, 0);
328 zfree(temp);
329 }
330 zfree(mul);
331 *res = ans;
332 return 3;
333 }
334
335 /*
336 * Compute the combinatorial function M! / ( N! * (M - N)! ).
337 * Returns 0 if result is 0
338 * 1 1
339 * 2 z1
340 * -1 -1
341 * -2 if too complicated
342 * 3 result stored at res
343 */
344 int
zcomb(ZVALUE z1,ZVALUE z2,ZVALUE * res)345 zcomb(ZVALUE z1, ZVALUE z2, ZVALUE *res)
346 {
347 ZVALUE z3, z4;
348 int r;
349
350 if (z2.sign || (!z1.sign && zrel(z2, z1) > 0))
351 return 0;
352 if (zisone(z2))
353 return 2;
354 if (z1.sign) {
355 z1.sign = 0;
356 zsub(z1, _one_, &z3);
357 zadd(z3, z2, &z4);
358 zfree(z3);
359 r = docomb(z4, z2, res);
360 if (r == 2) {
361 *res = z4;
362 r = 3;
363 }
364 else
365 zfree(z4);
366 if (z2.v[0] & 1) {
367 if (r == 1)
368 r = -1;
369 if (r == 3)
370 res->sign = 1;
371 }
372 return r;
373 }
374 return docomb(z1, z2, res);
375 }
376
377
378 /*
379 * Compute the Jacobi function (m / n) for odd n.
380 *
381 * The property of the Jacobi function is: If n>2 is prime then
382 *
383 * the result is 1 if m == x^2 (mod n) for some x.
384 * otherwise the result is -1.
385 *
386 * If n is not prime, then the result does not prove that n is not prime
387 * when the value of the Jacobi is 1.
388 *
389 * Jacobi evaluation of (m / n), where n > 0 is odd AND m > 0 is odd:
390 *
391 * rule 0: (0 / n) == 0
392 * rule 1: (1 / n) == 1
393 * rule 2: (m / n) == (a / n) if m == a % n
394 * rule 3: (m / n) == (2*m / n) if n == 1 % 8 OR n == 7 % 8
395 * rule 4: (m / n) == -(2*m / n) if n != 1 & 8 AND n != 7 % 8
396 * rule 5: (m / n) == (n / m) if m == 3 % 4 AND n == 3 % 4
397 * rule 6: (m / n) == -(n / m) if m != 3 % 4 OR n != 3 % 4
398 *
399 * NOTE: This function returns 0 in invalid Jacobi parameters:
400 * m < 0 OR n is even OR n < 1.
401 */
402 FLAG
zjacobi(ZVALUE z1,ZVALUE z2)403 zjacobi(ZVALUE z1, ZVALUE z2)
404 {
405 ZVALUE p, q, tmp;
406 long lowbit;
407 int val;
408
409 /* firewall */
410 if (ziszero(z1) || zisneg(z1))
411 return 0;
412 if (ziseven(z2) || zisneg(z2))
413 return 0;
414
415 /* assume a value of 1 unless we find otherwise */
416 if (zisone(z1))
417 return 1;
418 val = 1;
419 zcopy(z1, &p);
420 zcopy(z2, &q);
421 for (;;) {
422 zmod(p, q, &tmp, 0);
423 zfree(p);
424 p = tmp;
425 if (ziszero(p)) {
426 zfree(p);
427 zfree(q);
428 return 0;
429 }
430 if (ziseven(p)) {
431 lowbit = zlowbit(p);
432 zshift(p, -lowbit, &tmp);
433 zfree(p);
434 p = tmp;
435 if ((lowbit & 1) && (((*q.v & 0x7) == 3) ||
436 ((*q.v & 0x7) == 5)))
437 val = -val;
438 }
439 if (zisunit(p)) {
440 zfree(p);
441 zfree(q);
442 return val;
443 }
444 if ((*p.v & *q.v & 0x3) == 3)
445 val = -val;
446 tmp = q;
447 q = p;
448 p = tmp;
449 }
450 }
451
452
453 /*
454 * Return the Fibonacci number F(n).
455 * This is evaluated by recursively using the formulas:
456 * F(2N+1) = F(N+1)^2 + F(N)^2
457 * and
458 * F(2N) = F(N+1)^2 - F(N-1)^2
459 */
460 void
zfib(ZVALUE z,ZVALUE * res)461 zfib(ZVALUE z, ZVALUE *res)
462 {
463 long n;
464 int sign;
465 ZVALUE fnm1, fn, fnp1; /* consecutive Fibonacci values */
466 ZVALUE t1, t2, t3;
467 FULL i;
468
469 if (zge31b(z)) {
470 math_error("Very large Fibonacci number");
471 /*NOTREACHED*/
472 }
473 n = ztolong(z);
474 if (n == 0) {
475 *res = _zero_;
476 return;
477 }
478 sign = z.sign && ((n & 0x1) == 0);
479 if (n <= 2) {
480 *res = _one_;
481 res->sign = (BOOL)sign;
482 return;
483 }
484 i = TOPFULL;
485 while ((i & n) == 0)
486 i >>= (FULL)1;
487 i >>= (FULL)1;
488 fnm1 = _zero_;
489 fn = _one_;
490 fnp1 = _one_;
491 while (i) {
492 zsquare(fnm1, &t1);
493 zsquare(fn, &t2);
494 zsquare(fnp1, &t3);
495 zfree(fnm1);
496 zfree(fn);
497 zfree(fnp1);
498 zadd(t2, t3, &fnp1);
499 zsub(t3, t1, &fn);
500 zfree(t1);
501 zfree(t2);
502 zfree(t3);
503 if (i & n) {
504 fnm1 = fn;
505 fn = fnp1;
506 zadd(fnm1, fn, &fnp1);
507 } else {
508 zsub(fnp1, fn, &fnm1);
509 }
510 i >>= (FULL)1;
511 }
512 zfree(fnm1);
513 zfree(fnp1);
514 *res = fn;
515 res->sign = (BOOL)sign;
516 }
517
518
519 /*
520 * Compute the result of raising one number to the power of another
521 * The second number is assumed to be non-negative.
522 * It cannot be too large except for trivial cases.
523 */
524 void
zpowi(ZVALUE z1,ZVALUE z2,ZVALUE * res)525 zpowi(ZVALUE z1, ZVALUE z2, ZVALUE *res)
526 {
527 int sign; /* final sign of number */
528 unsigned long power; /* power to raise to */
529 FULL bit; /* current bit value */
530 long twos; /* count of times 2 is in result */
531 ZVALUE ans, temp;
532
533 sign = (z1.sign && zisodd(z2));
534 z1.sign = 0;
535 z2.sign = 0;
536 if (ziszero(z2) && !ziszero(z1)) { /* number raised to power 0 */
537 *res = _one_;
538 return;
539 }
540 if (zisabsleone(z1)) { /* 0, 1, or -1 raised to a power */
541 ans = _one_;
542 ans.sign = (BOOL)sign;
543 if (*z1.v == 0)
544 ans = _zero_;
545 *res = ans;
546 return;
547 }
548 if (zge31b(z2)) {
549 math_error("Raising to very large power");
550 /*NOTREACHED*/
551 }
552 power = ztoulong(z2);
553 if (zistwo(z1)) { /* two raised to a power */
554 zbitvalue((long) power, res);
555 return;
556 }
557 /*
558 * See if this is a power of ten
559 */
560 if (zistiny(z1) && (*z1.v == 10)) {
561 ztenpow((long) power, res);
562 res->sign = (BOOL)sign;
563 return;
564 }
565 /*
566 * Handle low powers specially
567 */
568 if (power <= 4) {
569 switch ((int) power) {
570 case 1:
571 ans.len = z1.len;
572 ans.v = alloc(ans.len);
573 zcopyval(z1, ans);
574 ans.sign = (BOOL)sign;
575 *res = ans;
576 return;
577 case 2:
578 zsquare(z1, res);
579 return;
580 case 3:
581 zsquare(z1, &temp);
582 zmul(z1, temp, res);
583 zfree(temp);
584 res->sign = (BOOL)sign;
585 return;
586 case 4:
587 zsquare(z1, &temp);
588 zsquare(temp, res);
589 zfree(temp);
590 return;
591 }
592 }
593 /*
594 * Shift out all powers of twos so the multiplies are smaller.
595 * We will shift back the right amount when done.
596 */
597 twos = 0;
598 if (ziseven(z1)) {
599 twos = zlowbit(z1);
600 ans.v = alloc(z1.len);
601 ans.len = z1.len;
602 ans.sign = z1.sign;
603 zcopyval(z1, ans);
604 zshiftr(ans, twos);
605 ztrim(&ans);
606 z1 = ans;
607 twos *= power;
608 }
609 /*
610 * Compute the power by squaring and multiplying.
611 * This uses the left to right method of power raising.
612 */
613 bit = TOPFULL;
614 while ((bit & power) == 0)
615 bit >>= 1;
616 bit >>= 1;
617 zsquare(z1, &ans);
618 if (bit & power) {
619 zmul(ans, z1, &temp);
620 zfree(ans);
621 ans = temp;
622 }
623 bit >>= 1;
624 while (bit) {
625 zsquare(ans, &temp);
626 zfree(ans);
627 ans = temp;
628 if (bit & power) {
629 zmul(ans, z1, &temp);
630 zfree(ans);
631 ans = temp;
632 }
633 bit >>= 1;
634 }
635 /*
636 * Scale back up by proper power of two
637 */
638 if (twos) {
639 zshift(ans, twos, &temp);
640 zfree(ans);
641 ans = temp;
642 zfree(z1);
643 }
644 ans.sign = (BOOL)sign;
645 *res = ans;
646 }
647
648
649 /*
650 * Compute ten to the specified power
651 * This saves some work since the squares of ten are saved.
652 */
653 void
ztenpow(long power,ZVALUE * res)654 ztenpow(long power, ZVALUE *res)
655 {
656 long i;
657 ZVALUE ans;
658 ZVALUE temp;
659
660 if (power <= 0) {
661 *res = _one_;
662 return;
663 }
664 ans = _one_;
665 _tenpowers_[0] = _ten_;
666 for (i = 0; power; i++) {
667 if (_tenpowers_[i].len == 0) {
668 if (i <= TEN_MAX) {
669 zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
670 } else {
671 math_error("cannot compute 10^2^(TEN_MAX+1)");
672 /*NOTREACHED*/
673 }
674 }
675 if (power & 0x1) {
676 zmul(ans, _tenpowers_[i], &temp);
677 zfree(ans);
678 ans = temp;
679 }
680 power /= 2;
681 }
682 *res = ans;
683 }
684
685
686 /*
687 * Calculate modular inverse suppressing unnecessary divisions.
688 * This is based on the Euclidean algorithm for large numbers.
689 * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
690 * Returns TRUE if there is no solution because the numbers
691 * are not relatively prime.
692 */
693 BOOL
zmodinv(ZVALUE u,ZVALUE v,ZVALUE * res)694 zmodinv(ZVALUE u, ZVALUE v, ZVALUE *res)
695 {
696 FULL q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
697 ZVALUE u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
698
699 v.sign = 0;
700 if (zisneg(u) || (zrel(u, v) >= 0))
701 zmod(u, v, &v3, 0);
702 else
703 zcopy(u, &v3);
704 zcopy(v, &u3);
705 u2 = _zero_;
706 v2 = _one_;
707
708 /*
709 * Loop here while the size of the numbers remain above
710 * the size of a HALF. Throughout this loop u3 >= v3.
711 */
712 while ((u3.len > 1) && !ziszero(v3)) {
713 vh = 0;
714 #if LONG_BITS == BASEB
715 uh = u3.v[u3.len - 1];
716 if (v3.len == u3.len)
717 vh = v3.v[v3.len - 1];
718 #else
719 uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
720 if ((v3.len + 1) >= u3.len)
721 vh = v3.v[v3.len - 1];
722 if (v3.len == u3.len)
723 vh = (vh << BASEB) + v3.v[v3.len - 2];
724 #endif
725 A = 1;
726 B = 0;
727 C = 0;
728 D = 1;
729
730 /*
731 * Calculate successive quotients of the continued fraction
732 * expansion using only single precision arithmetic until
733 * greater precision is required.
734 */
735 while ((vh + C) && (vh + D)) {
736 q1 = (uh + A) / (vh + C);
737 q2 = (uh + B) / (vh + D);
738 if (q1 != q2)
739 break;
740 T = A - q1 * C;
741 A = C;
742 C = T;
743 T = B - q1 * D;
744 B = D;
745 D = T;
746 T = uh - q1 * vh;
747 uh = vh;
748 vh = T;
749 }
750
751 /*
752 * If B is zero, then we made no progress because
753 * the calculation requires a very large quotient.
754 * So we must do this step of the calculation in
755 * full precision
756 */
757 if (B == 0) {
758 zquo(u3, v3, &qz, 0);
759 zmul(qz, v2, &tmp1);
760 zsub(u2, tmp1, &tmp2);
761 zfree(tmp1);
762 zfree(u2);
763 u2 = v2;
764 v2 = tmp2;
765 zmul(qz, v3, &tmp1);
766 zsub(u3, tmp1, &tmp2);
767 zfree(tmp1);
768 zfree(u3);
769 u3 = v3;
770 v3 = tmp2;
771 zfree(qz);
772 continue;
773 }
774 /*
775 * Apply the calculated A,B,C,D numbers to the current
776 * values to update them as if the full precision
777 * calculations had been carried out.
778 */
779 zmuli(u2, (long) A, &tmp1);
780 zmuli(v2, (long) B, &tmp2);
781 zadd(tmp1, tmp2, &tmp3);
782 zfree(tmp1);
783 zfree(tmp2);
784 zmuli(u2, (long) C, &tmp1);
785 zmuli(v2, (long) D, &tmp2);
786 zfree(u2);
787 zfree(v2);
788 u2 = tmp3;
789 zadd(tmp1, tmp2, &v2);
790 zfree(tmp1);
791 zfree(tmp2);
792 zmuli(u3, (long) A, &tmp1);
793 zmuli(v3, (long) B, &tmp2);
794 zadd(tmp1, tmp2, &tmp3);
795 zfree(tmp1);
796 zfree(tmp2);
797 zmuli(u3, (long) C, &tmp1);
798 zmuli(v3, (long) D, &tmp2);
799 zfree(u3);
800 zfree(v3);
801 u3 = tmp3;
802 zadd(tmp1, tmp2, &v3);
803 zfree(tmp1);
804 zfree(tmp2);
805 }
806
807 /*
808 * Here when the remaining numbers become single precision in size.
809 * Finish the procedure using single precision calculations.
810 */
811 if (ziszero(v3) && !zisone(u3)) {
812 zfree(u3);
813 zfree(v3);
814 zfree(u2);
815 zfree(v2);
816 return TRUE;
817 }
818 ui3 = ztofull(u3);
819 vi3 = ztofull(v3);
820 zfree(u3);
821 zfree(v3);
822 while (vi3) {
823 q1 = ui3 / vi3;
824 zmuli(v2, (long) q1, &tmp1);
825 zsub(u2, tmp1, &tmp2);
826 zfree(tmp1);
827 zfree(u2);
828 u2 = v2;
829 v2 = tmp2;
830 q2 = ui3 - q1 * vi3;
831 ui3 = vi3;
832 vi3 = q2;
833 }
834 zfree(v2);
835 if (ui3 != 1) {
836 zfree(u2);
837 return TRUE;
838 }
839 if (zisneg(u2)) {
840 zadd(v, u2, res);
841 zfree(u2);
842 return FALSE;
843 }
844 *res = u2;
845 return FALSE;
846 }
847
848
849 /*
850 * Compute the greatest common divisor of a pair of integers.
851 */
852 void
zgcd(ZVALUE z1,ZVALUE z2,ZVALUE * res)853 zgcd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
854 {
855 int h, i, j, k;
856 LEN len, l, m, n, o, p, q;
857 HALF u, v, w, x;
858 HALF *a, *a0, *A, *b, *b0, *B, *c, *d;
859 FULL f, g;
860 ZVALUE gcd;
861 BOOL needw;
862
863 if (zisunit(z1) || zisunit(z2)) {
864 *res = _one_;
865 return;
866 }
867 z1.sign = 0;
868 z2.sign = 0;
869 if (ziszero(z1) || !zcmp(z1, z2)) {
870 zcopy(z2, res);
871 return;
872 }
873 if (ziszero(z2)) {
874 zcopy(z1, res);
875 return;
876 }
877
878 o = 0;
879 while (!(z1.v[o] | z2.v[o])) o++; /* Count common zero digits */
880
881 c = z1.v + o;
882 d = z2.v + o;
883
884 m = z1.len - o;
885 n = z2.len - o;
886 u = *c | *d; /* Count common zero bits */
887 v = 1;
888 p = 0;
889 while (!(u & v)) {
890 v <<= 1;
891 p++;
892 }
893
894 while (!*c) { /* Removing zero digits */
895 c++;
896 m--;
897 }
898
899 while (!*d) {
900 d++;
901 n--;
902 }
903
904
905 u = *d; /* Count zero bits for *d */
906 v = 1;
907 q = 0;
908 while (!(u & v)) {
909 v <<= 1;
910 q++;
911 }
912
913 a0 = A = alloc(m);
914 b0 = B = alloc(n);
915
916 memcpy(A, c, m * sizeof(HALF)); /* Copy c[] to A[] */
917
918 /* Copy d[] to B[], shifting if necessary */
919 if (q) {
920 i = n;
921 b = B + n;
922 d += n;
923 f = 0;
924 while (i--) {
925 f = f << BASEB | *--d;
926 *--b = (HALF) (f >> q);
927 }
928 if (B[n-1] == 0) n--;
929 }
930 else memcpy(B, d, n * sizeof(HALF));
931
932 if (n == 1) { /* One digit case; use Euclid's algorithm */
933 n = m;
934 b0 = A;
935 m = 1;
936 a0 = B;
937 if (m == 1) { /* a has one digit */
938 v = *a0;
939 if (v > 1) { /* Euclid's algorithm */
940 b = b0 + n;
941 i = n;
942 u = 0;
943 while (i--) {
944 f = (FULL) u << BASEB | *--b;
945 u = (HALF) (f % v);
946 }
947 while (u) { w = v % u; v = u; u = w; }
948 }
949 *b0 = v;
950 n = 1;
951 }
952 len = n + o;
953 gcd.v = alloc(len + 1);
954 /* Common zero digits */
955 if (o) memset(gcd.v, 0, o * sizeof(HALF));
956 /* Left shift for common zero bits */
957 if (p) {
958 i = n;
959 f = 0;
960 b = b0;
961 a = gcd.v + o;
962 while (i--) {
963 f = f >> BASEB | (FULL) *b++ << p;
964 *a++ = (HALF) f;
965 }
966 if (f >>= BASEB) {len++; *a = (HALF) f;}
967 } else {
968 memcpy(gcd.v + o, b0, n * sizeof(HALF));
969 }
970 gcd.len = len;
971 gcd.sign = 0;
972 freeh(A);
973 freeh(B);
974 *res = gcd;
975 return;
976 }
977
978 u = B[n-1]; /* Bit count for b */
979 k = (n - 1) * BASEB;
980 while (u >>= 1) k++;
981
982 needw = TRUE;
983
984 w = 0;
985 j = 0;
986 while (m) { /* START OF MAIN LOOP */
987 if (m - n < 2 || needw) {
988 q = 0;
989 u = *a0;
990 v = 1;
991 while (!(u & v)) { /* count zero bits for *a0 */
992 q++;
993 v <<= 1;
994 }
995
996 if (q) { /* right-justify a */
997 a = a0 + m;
998 i = m;
999 f = 0;
1000 while (i--) {
1001 f = f << BASEB | *--a;
1002 *a = (HALF) (f >> q);
1003 }
1004 if (!a0[m-1]) m--; /* top digit vanishes */
1005 }
1006
1007 if (m == 1) break;
1008
1009 u = a0[m-1];
1010 j = (m - 1) * BASEB;
1011 while (u >>= 1) j++; /* counting bits for a */
1012 h = j - k;
1013 if (h < 0) { /* swapping to get h > 0 */
1014 l = m;
1015 m = n;
1016 n = l;
1017 a = a0;
1018 a0 = b0;
1019 b0 = a;
1020 k = j;
1021 h = -h;
1022 needw = TRUE;
1023 }
1024 if (h > 1) {
1025 if (needw) { /* find w = minv(*b0, h0) */
1026 u = 1;
1027 v = *b0;
1028 w = 0;
1029 x = 1;
1030 i = h;
1031 while (i-- && x) {
1032 if (u & x) { u -= v * x; w |= x;}
1033 x <<= 1;
1034 }
1035 needw = FALSE;
1036 }
1037 g = (FULL) (*a0 * w);
1038 if (h < BASEB) {
1039 g &= (1 << h) - 1;
1040 } else {
1041 g &= BASE1;
1042 }
1043 } else {
1044 g = 1;
1045 }
1046 } else {
1047 g = (FULL) (*a0 * w);
1048 }
1049 a = a0;
1050 b = b0;
1051 i = n;
1052 if (g > 1) { /* a - g * b case */
1053 f = 0;
1054 while (i--) {
1055 f = (FULL) *a - g * *b++ - f;
1056 *a++ = (HALF) f;
1057 f >>= BASEB;
1058 f = -f & BASE1;
1059 }
1060 if (f) {
1061 i = m - n;
1062 while (i-- && f) {
1063 f = *a - f;
1064 *a++ = (HALF) f;
1065 f >>= BASEB;
1066 f = -f & BASE1;
1067 }
1068 }
1069 while (m && !*a0) { /* Removing trailing zeros */
1070 m--;
1071 a0++;
1072 }
1073 if (f) { /* a - g * b < 0 */
1074 while (m > 1 && a0[m-1] == BASE1) m--;
1075 *a0 = - *a0;
1076 a = a0;
1077 i = m;
1078 while (--i) {
1079 a++;
1080 *a = ~*a;
1081 }
1082 }
1083 } else { /* abs(a - b) case */
1084 while (i && *a++ == *b++) i--;
1085 q = n - i;
1086 if (m == n) { /* a and b same length */
1087 if (i) { /* a not equal to b */
1088 while (m && a0[m-1] == b0[m-1]) m--;
1089 if (a0[m-1] < b0[m-1]) {
1090 /* Swapping since a < b */
1091 a = a0;
1092 a0 = b0;
1093 b0 = a;
1094 k = j;
1095 }
1096 a = a0 + q;
1097 b = b0 + q;
1098 i = m - q;
1099 f = 0;
1100 while (i--) {
1101 f = (FULL) *a - *b++ - f;
1102 *a++ = (HALF) f;
1103 f >>= BASEB;
1104 f = -f & BASE1;
1105 }
1106 }
1107 } else { /* a has more digits than b */
1108 a = a0 + q;
1109 b = b0 + q;
1110 i = n - q;
1111 f = 0;
1112 while (i--) {
1113 f = (FULL) *a - *b++ - f;
1114 *a++ = (HALF) f;
1115 f >>= BASEB;
1116 f = -f & BASE1;
1117 }
1118 if (f) { while (!*a) *a++ = BASE1;
1119 (*a)--;
1120 }
1121 }
1122 a0 += q;
1123 m -= q;
1124 while (m && !*a0) { /* Removing trailing zeros */
1125 m--;
1126 a0++;
1127 }
1128 }
1129 while (m && !a0[m-1]) m--; /* Removing leading zeros */
1130 }
1131 if (m == 1) { /* a has one digit */
1132 v = *a0;
1133 if (v > 1) { /* Euclid's algorithm */
1134 b = b0 + n;
1135 i = n;
1136 u = 0;
1137 while (i--) {
1138 f = (FULL) u << BASEB | *--b;
1139 u = (HALF) (f % v);
1140 }
1141 while (u) { w = v % u; v = u; u = w; }
1142 }
1143 *b0 = v;
1144 n = 1;
1145 }
1146 len = n + o;
1147 gcd.v = alloc(len + 1);
1148 if (o) memset(gcd.v, 0, o * sizeof(HALF)); /* Common zero digits */
1149 if (p) { /* Left shift for common zero bits */
1150 i = n;
1151 f = 0;
1152 b = b0;
1153 a = gcd.v + o;
1154 while (i--) {
1155 f = (FULL) *b++ << p | f;
1156 *a++ = (HALF) f;
1157 f >>= BASEB;
1158 }
1159 if (f) {
1160 len++; *a = (HALF) f;
1161 }
1162 } else {
1163 memcpy(gcd.v + o, b0, n * sizeof(HALF));
1164 }
1165 gcd.len = len;
1166 gcd.sign = 0;
1167 freeh(A);
1168 freeh(B);
1169 *res = gcd;
1170 return;
1171 }
1172
1173 /*
1174 * Compute the lcm of two integers (least common multiple).
1175 * This is done using the formula: gcd(a,b) * lcm(a,b) = a * b.
1176 */
1177 void
zlcm(ZVALUE z1,ZVALUE z2,ZVALUE * res)1178 zlcm(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1179 {
1180 ZVALUE temp1, temp2;
1181
1182 zgcd(z1, z2, &temp1);
1183 zequo(z1, temp1, &temp2);
1184 zfree(temp1);
1185 zmul(temp2, z2, res);
1186 zfree(temp2);
1187 }
1188
1189
1190 /*
1191 * Return whether or not two numbers are relatively prime to each other.
1192 */
1193 BOOL
zrelprime(ZVALUE z1,ZVALUE z2)1194 zrelprime(ZVALUE z1, ZVALUE z2)
1195 {
1196 FULL rem1, rem2; /* remainders */
1197 ZVALUE rem;
1198 BOOL result;
1199
1200 z1.sign = 0;
1201 z2.sign = 0;
1202 if (ziseven(z1) && ziseven(z2)) /* false if both even */
1203 return FALSE;
1204 if (zisunit(z1) || zisunit(z2)) /* true if either is a unit */
1205 return TRUE;
1206 if (ziszero(z1) || ziszero(z2)) /* false if either is zero */
1207 return FALSE;
1208 if (zistwo(z1) || zistwo(z2)) /* true if either is two */
1209 return TRUE;
1210 /*
1211 * Try reducing each number by the product of the first few odd primes
1212 * to see if any of them are a common factor.
1213 */
1214 rem1 = zmodi(z1, (FULL)3 * 5 * 7 * 11 * 13);
1215 rem2 = zmodi(z2, (FULL)3 * 5 * 7 * 11 * 13);
1216 if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
1217 return FALSE;
1218 if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
1219 return FALSE;
1220 if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
1221 return FALSE;
1222 if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
1223 return FALSE;
1224 if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
1225 return FALSE;
1226 /*
1227 * Try a new batch of primes now
1228 */
1229 rem1 = zmodi(z1, (FULL)17 * 19 * 23);
1230 rem2 = zmodi(z2, (FULL)17 * 19 * 23);
1231 if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
1232 return FALSE;
1233 if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
1234 return FALSE;
1235 if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
1236 return FALSE;
1237 /*
1238 * Yuk, we must actually compute the gcd to know the answer
1239 */
1240 zgcd(z1, z2, &rem);
1241 result = zisunit(rem);
1242 zfree(rem);
1243 return result;
1244 }
1245
1246
1247 /*
1248 * Compute the integer floor of the log of an integer to a specified base.
1249 * The signs of the integers and base are ignored.
1250 * Example: zlog(123456, 10) = 5.
1251 */
1252 long
zlog(ZVALUE z,ZVALUE base)1253 zlog(ZVALUE z, ZVALUE base)
1254 {
1255 ZVALUE *zp; /* current square */
1256 long power; /* current power */
1257 ZVALUE temp; /* temporary */
1258 ZVALUE squares[32]; /* table of squares of base */
1259
1260 /* ignore signs */
1261
1262 z.sign = 0;
1263 base.sign = 0;
1264
1265 /*
1266 * Make sure that the numbers are nonzero and the base is > 1
1267 */
1268 if (ziszero(z) || ziszero(base) || zisone(base)) {
1269 math_error("Zero or too small argument argument for zlog!!!");
1270 /*NOTREACHED*/
1271 }
1272
1273 /*
1274 * Some trivial cases.
1275 */
1276 power = zrel(z, base);
1277 if (power <= 0)
1278 return (power + 1);
1279
1280 /* base - power of two */
1281 if (zisonebit(base))
1282 return (zhighbit(z) / zlowbit(base));
1283
1284 /* base = 10 */
1285 if (base.len == 1 && base.v[0] == 10)
1286 return zlog10(z, NULL);
1287 /*
1288 * Now loop by squaring the base each time, and see whether or
1289 * not each successive square is still smaller than the number.
1290 */
1291 zp = &squares[0];
1292 *zp = base;
1293 while (zp->len * 2 - 1 <= z.len && zrel(z, *zp) > 0) {
1294 /* while square not too large */
1295 zsquare(*zp, zp + 1);
1296 zp++;
1297 }
1298
1299 /*
1300 * Now back down the squares,
1301 */
1302 power = 0;
1303 for (; zp > squares; zp--) {
1304 if (zrel(z, *zp) >= 0) {
1305 zquo(z, *zp, &temp, 0);
1306 if (power)
1307 zfree(z);
1308 z = temp;
1309 power++;
1310 }
1311 zfree(*zp);
1312 power <<= 1;
1313 }
1314 if (zrel(z, *zp) >= 0)
1315 power++;
1316 if (power > 1)
1317 zfree(z);
1318 return power;
1319 }
1320
1321
1322 /*
1323 * Return the integral log base 10 of a number.
1324 *
1325 * If was_10_power != NULL, then this flag is set to TRUE if the
1326 * value was a power of 10, FALSE otherwise.
1327 */
1328 long
zlog10(ZVALUE z,BOOL * was_10_power)1329 zlog10(ZVALUE z, BOOL *was_10_power)
1330 {
1331 ZVALUE *zp; /* current square */
1332 long power; /* current power */
1333 ZVALUE temp; /* temporary */
1334 ZVALUE pow10; /* power of 10 */
1335 FLAG rel; /* relationship */
1336 int i;
1337
1338 if (ziszero(z)) {
1339 math_error("Zero argument argument for zlog10");
1340 /*NOTREACHED*/
1341 }
1342
1343 /* Ignore sign of z */
1344 z.sign = 0;
1345
1346 /* preload power10 table if missing */
1347 if (power10 == NULL) {
1348 long v;
1349
1350 /* determine power10 table size */
1351 for (v=1, max_power10_exp=0;
1352 v <= (long)(MAXLONG/10L);
1353 v *= 10L, ++max_power10_exp) {
1354 }
1355
1356 /* create power10 table */
1357 power10 = malloc(sizeof(long) * (max_power10_exp+1));
1358 if (power10 == NULL) {
1359 math_error("cannot malloc power10 table");
1360 /*NOTREACHED*/
1361 }
1362
1363 /* load power10 table */
1364 for (i=0, v = 1L; i <= max_power10_exp; ++i, v *= 10L) {
1365 power10[i] = v;
1366 }
1367 }
1368
1369 /* assume not a power of ten unless we find out otherwise */
1370 if (was_10_power != NULL) {
1371 *was_10_power = FALSE;
1372 }
1373
1374 /* quick exit for small values */
1375 if (! zgtmaxlong(z)) {
1376 long value = ztolong(z);
1377
1378 for (i=0; i <= max_power10_exp; ++i) {
1379 if (value == power10[i]) {
1380 if (was_10_power != NULL) {
1381 *was_10_power = TRUE;
1382 }
1383 return i;
1384 } else if (value < power10[i]) {
1385 return i-1;
1386 }
1387 }
1388 }
1389
1390 /*
1391 * Loop by squaring the base each time, and see whether or
1392 * not each successive square is still smaller than the number.
1393 */
1394 zp = &_tenpowers_[0];
1395 *zp = _ten_;
1396 while (((zp->len * 2) - 1) <= z.len) { /* while square not too large */
1397 if (zp >= &_tenpowers_[TEN_MAX]) {
1398 math_error("Maximum storable power of 10 reached!");
1399 /*NOTREACHED*/
1400 }
1401 if (zp[1].len == 0)
1402 zsquare(*zp, zp + 1);
1403 zp++;
1404 }
1405
1406 /*
1407 * Now back down the squares, and multiply them together to see
1408 * exactly how many times the base can be raised by.
1409 */
1410 /* find the tenpower table entry < z */
1411 do {
1412 rel = zrel(*zp, z);
1413 if (rel == 0) {
1414 /* quick exit - we match a tenpower entry */
1415 if (was_10_power != NULL) {
1416 *was_10_power = TRUE;
1417 }
1418 return (1L << (zp - _tenpowers_));
1419 }
1420 } while (rel > 0 && --zp >= _tenpowers_);
1421 if (zp < _tenpowers_) {
1422 math_error("fell off bottom of tenpower table!");
1423 /*NOTREACHED*/
1424 }
1425
1426 /* the tenpower value is now our starting comparison value */
1427 zcopy(*zp, &pow10);
1428 power = (1L << (zp - _tenpowers_));
1429
1430 /* try to build up a power of 10 from tenpower table entries */
1431 while (--zp >= _tenpowers_) {
1432
1433 /* try the next lower tenpower value */
1434 zmul(pow10, *zp, &temp);
1435 rel = zrel(temp, z);
1436 if (rel == 0) {
1437 /* exact power of 10 match */
1438 power += (1L << (zp - _tenpowers_));
1439 if (was_10_power != NULL) {
1440 *was_10_power = TRUE;
1441 }
1442 zfree(pow10);
1443 zfree(temp);
1444 return power;
1445
1446 /* ignore this entry if we went too large */
1447 } else if (rel > 0) {
1448 zfree(temp);
1449
1450 /* otherwise increase power and keep going */
1451 } else {
1452 power += (1L << (zp - _tenpowers_));
1453 zfree(pow10);
1454 pow10 = temp;
1455 }
1456 }
1457 zfree(pow10);
1458 return power;
1459 }
1460
1461
1462 /*
1463 * Return the number of times that one number will divide another.
1464 * This works similarly to zlog, except that divisions must be exact.
1465 * For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
1466 */
1467 long
zdivcount(ZVALUE z1,ZVALUE z2)1468 zdivcount(ZVALUE z1, ZVALUE z2)
1469 {
1470 long count; /* number of factors removed */
1471 ZVALUE tmp; /* ignored return value */
1472
1473 if (ziszero(z1) || ziszero(z2) || zisunit(z2))
1474 return 0;
1475 count = zfacrem(z1, z2, &tmp);
1476 zfree(tmp);
1477 return count;
1478 }
1479
1480
1481 /*
1482 * Remove all occurrences of the specified factor from a number.
1483 * Also returns the number of factors removed as a function return value.
1484 * Example: zfacrem(540, 3, &x) returns 3 and sets x to 20.
1485 */
1486 long
zfacrem(ZVALUE z1,ZVALUE z2,ZVALUE * rem)1487 zfacrem(ZVALUE z1, ZVALUE z2, ZVALUE *rem)
1488 {
1489 register ZVALUE *zp; /* current square */
1490 long count; /* total count of divisions */
1491 long worth; /* worth of current square */
1492 long lowbit; /* for zlowbit(z2) */
1493 ZVALUE temp1, temp2, temp3; /* temporaries */
1494 ZVALUE squares[32]; /* table of squares of factor */
1495
1496 z1.sign = 0;
1497 z2.sign = 0;
1498 /*
1499 * Reject trivial cases.
1500 */
1501 if ((z1.len < z2.len) || (zisodd(z1) && ziseven(z2)) ||
1502 ziszero(z2) || zisone(z2) ||
1503 ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
1504 rem->v = alloc(z1.len);
1505 rem->len = z1.len;
1506 rem->sign = 0;
1507 zcopyval(z1, *rem);
1508 return 0;
1509 }
1510 /*
1511 * Handle any power of two special.
1512 */
1513 if (zisonebit(z2)) {
1514 lowbit = zlowbit(z2);
1515 count = zlowbit(z1) / lowbit;
1516 rem->v = alloc(z1.len);
1517 rem->len = z1.len;
1518 rem->sign = 0;
1519 zcopyval(z1, *rem);
1520 zshiftr(*rem, count * lowbit);
1521 ztrim(rem);
1522 return count;
1523 }
1524 /*
1525 * See if the factor goes in even once.
1526 */
1527 zdiv(z1, z2, &temp1, &temp2, 0);
1528 if (!ziszero(temp2)) {
1529 zfree(temp1);
1530 zfree(temp2);
1531 rem->v = alloc(z1.len);
1532 rem->len = z1.len;
1533 rem->sign = 0;
1534 zcopyval(z1, *rem);
1535 return 0;
1536 }
1537 zfree(temp2);
1538 z1 = temp1;
1539 /*
1540 * Now loop by squaring the factor each time, and see whether
1541 * or not each successive square will still divide the number.
1542 */
1543 count = 1;
1544 worth = 1;
1545 zp = &squares[0];
1546 *zp = z2;
1547 while (((zp->len * 2) - 1) <= z1.len) { /* while square not too large */
1548 zsquare(*zp, &temp1);
1549 zdiv(z1, temp1, &temp2, &temp3, 0);
1550 if (!ziszero(temp3)) {
1551 zfree(temp1);
1552 zfree(temp2);
1553 zfree(temp3);
1554 break;
1555 }
1556 zfree(temp3);
1557 zfree(z1);
1558 z1 = temp2;
1559 *++zp = temp1;
1560 worth *= 2;
1561 count += worth;
1562 }
1563
1564 /*
1565 * Now back down the list of squares, and see if the lower powers
1566 * will divide any more times.
1567 */
1568 /*
1569 * We prevent the zp pointer from walking behind squares
1570 * by stopping one short of the end and running the loop one
1571 * more time.
1572 *
1573 * We could stop the loop with just zp >= squares, but stopping
1574 * short and running the loop one last time manually helps make
1575 * code checkers such as insure happy.
1576 */
1577 for (; zp > squares; zp--, worth /= 2) {
1578 if (zp->len <= z1.len) {
1579 zdiv(z1, *zp, &temp1, &temp2, 0);
1580 if (ziszero(temp2)) {
1581 temp3 = z1;
1582 z1 = temp1;
1583 temp1 = temp3;
1584 count += worth;
1585 }
1586 zfree(temp1);
1587 zfree(temp2);
1588 }
1589 if (zp != squares)
1590 zfree(*zp);
1591 }
1592 /* run the loop manually one last time */
1593 if (zp == squares) {
1594 if (zp->len <= z1.len) {
1595 zdiv(z1, *zp, &temp1, &temp2, 0);
1596 if (ziszero(temp2)) {
1597 temp3 = z1;
1598 z1 = temp1;
1599 temp1 = temp3;
1600 count += worth;
1601 }
1602 zfree(temp1);
1603 zfree(temp2);
1604 }
1605 if (zp != squares)
1606 zfree(*zp);
1607 }
1608
1609 *rem = z1;
1610 return count;
1611 }
1612
1613
1614 /*
1615 * Keep dividing a number by the gcd of it with another number until the
1616 * result is relatively prime to the second number. Returns the number
1617 * of divisions made, and if this is positive, stores result at res.
1618 */
1619 long
zgcdrem(ZVALUE z1,ZVALUE z2,ZVALUE * res)1620 zgcdrem(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1621 {
1622 ZVALUE tmp1, tmp2, tmp3, tmp4;
1623 long count, onecount;
1624 long sh;
1625
1626 if (ziszero(z1) || ziszero(z2)) {
1627 math_error("Zero argument in call to zgcdrem!!!");
1628 /*NOTREACHED*/
1629 }
1630 /*
1631 * Begin by taking the gcd for the first time.
1632 * If the number is already relatively prime, then we are done.
1633 */
1634 z1.sign = 0;
1635 z2.sign = 0;
1636 if (zisone(z2))
1637 return 0;
1638 if (zisonebit(z2)) {
1639 sh = zlowbit(z1);
1640 if (sh == 0)
1641 return 0;
1642 zshift(z1, -sh, res);
1643 return 1 + (sh - 1)/zlowbit(z2);
1644 }
1645 if (zisonebit(z1)) {
1646 if (zisodd(z2))
1647 return 0;
1648 *res = _one_;
1649 return zlowbit(z1);
1650 }
1651
1652 zgcd(z1, z2, &tmp1);
1653 if (zisunit(tmp1) || ziszero(tmp1)) {
1654 zfree(tmp1);
1655 return 0;
1656 }
1657 zequo(z1, tmp1, &tmp2);
1658 count = 1;
1659 z1 = tmp2;
1660 z2 = tmp1;
1661 /*
1662 * Now keep alternately taking the gcd and removing factors until
1663 * the gcd becomes one.
1664 */
1665 while (!zisunit(z2)) {
1666 onecount = zfacrem(z1, z2, &tmp3);
1667 if (onecount) {
1668 count += onecount;
1669 zfree(z1);
1670 z1 = tmp3;
1671 } else {
1672 zfree(tmp3);
1673 }
1674 zgcd(z1, z2, &tmp4);
1675 zfree(z2);
1676 z2 = tmp4;
1677 }
1678 zfree(z2);
1679 *res = z1;
1680 return count;
1681 }
1682
1683
1684 /*
1685 * Return the number of digits (base 10) in a number, ignoring the sign.
1686 */
1687 long
zdigits(ZVALUE z1)1688 zdigits(ZVALUE z1)
1689 {
1690 long count, val;
1691
1692 z1.sign = 0;
1693 if (!zge16b(z1)) { /* do small numbers ourself */
1694 count = 1;
1695 val = 10;
1696 while (*z1.v >= (HALF)val) {
1697 count++;
1698 val *= 10;
1699 }
1700 return count;
1701 }
1702 return (zlog10(z1, NULL) + 1);
1703 }
1704
1705
1706 /*
1707 * Return the single digit at the specified decimal place of a number,
1708 * where 0 means the rightmost digit. Example: zdigit(1234, 1) = 3.
1709 */
1710 long
zdigit(ZVALUE z1,long n)1711 zdigit(ZVALUE z1, long n)
1712 {
1713 ZVALUE tmp1, tmp2;
1714 long res;
1715
1716 z1.sign = 0;
1717 if (ziszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
1718 return 0;
1719 if (n == 0)
1720 return zmodi(z1, 10L);
1721 if (n == 1)
1722 return zmodi(z1, 100L) / 10;
1723 if (n == 2)
1724 return zmodi(z1, 1000L) / 100;
1725 if (n == 3)
1726 return zmodi(z1, 10000L) / 1000;
1727 ztenpow(n, &tmp1);
1728 zquo(z1, tmp1, &tmp2, 0);
1729 res = zmodi(tmp2, 10L);
1730 zfree(tmp1);
1731 zfree(tmp2);
1732 return res;
1733 }
1734
1735
1736 /*
1737 * z is to be a nonnegative integer
1738 * If z is the square of a integer stores at dest the square root of z;
1739 * otherwise stores at z an integer differing from the square root
1740 * by less than 1. Returns the sign of the true square root minus
1741 * the calculated integer. Type of rounding is determined by
1742 * rnd as follows: rnd = 0 gives round down, rnd = 1
1743 * rounds up, rnd = 8 rounds to even integer, rnd = 9 rounds to odd
1744 * integer, rnd = 16 rounds to nearest integer.
1745 */
1746 FLAG
zsqrt(ZVALUE z,ZVALUE * dest,long rnd)1747 zsqrt(ZVALUE z, ZVALUE *dest, long rnd)
1748 {
1749 HALF *a, *A, *b, *a0, u;
1750 int i, j, j1, j2, k, k1, m, m0, m1, n, n0, o;
1751 FULL d, e, f, g, h, s, t, x, topbit;
1752 int remsign;
1753 BOOL up, onebit;
1754 ZVALUE sqrt;
1755
1756 if (z.sign) {
1757 math_error("Square root of negative number");
1758 /*NOTREACHED*/
1759 }
1760 if (ziszero(z)) {
1761 *dest = _zero_;
1762 return 0;
1763 }
1764 m0 = z.len;
1765 o = m0 & 1;
1766 m = m0 + o; /* m is smallest even number >= z.len */
1767 n0 = n = m / 2;
1768 f = z.v[z.len - 1];
1769 k = 1;
1770 while (f >>= 2)
1771 k++;
1772 if (!o)
1773 k += BASEB/2;
1774 j = BASEB - k;
1775 m1 = m;
1776 if (k == BASEB) {
1777 m1 += 2;
1778 n0++;
1779 }
1780 A = alloc(m1);
1781 A[m1] = 0;
1782 a0 = A + n0;
1783 memcpy(A, z.v, m0 * sizeof(HALF));
1784 if (o)
1785 A[m - 1] = 0;
1786 if (n == 1) {
1787 if (j)
1788 f = (FULL) A[1] << j | A[0] >> k;
1789 else
1790 f = A[1];
1791 g = (FULL) A[0] << (j + BASEB);
1792 d = e = topbit = (FULL)1 << (k - 1);
1793 } else {
1794 if (j)
1795 f = (FULL) A[m-1] << (j + BASEB) | (FULL) A[m-2] << j |
1796 A[m-3] >> k;
1797 else
1798 f = (FULL) A[m-1] << BASEB | A[m-2];
1799 g = (FULL) A[m-3] << (j + BASEB) | (FULL) A[m-4] << j;
1800 d = e = topbit = (FULL)1 << (BASEB + k - 1);
1801 }
1802
1803 s = (f & topbit);
1804 f <<= 1;
1805 if (g & TOPFULL)
1806 f++;
1807 g <<= 1;
1808 if (s) {
1809 f -= 4 * d;
1810 e = 2 * d - 1;
1811 }
1812 else
1813 f -= d;
1814 while (d >>= 1) {
1815 if (!(s | f | g))
1816 break;
1817 while (d && (f & topbit) == s) {
1818 d >>= 1;
1819 f <<= 1;
1820 if (g & TOPFULL)
1821 f++;
1822 g <<= 1;
1823 }
1824 if (d == 0)
1825 break;
1826 if (s)
1827 f += e + 1;
1828 else
1829 f -= e;
1830 t = f & topbit;
1831 f <<= 1;
1832 if (g & TOPFULL)
1833 f++;
1834 g <<= 1;
1835 if (t == 0 && f < d)
1836 t = topbit;
1837 f -= d;
1838 if (s)
1839 e -= d - !t;
1840 else
1841 e += d - (t > 0);
1842 s = t;
1843 }
1844 if (n0 == 1) {
1845 A[1] = (HALF)e;
1846 A[0] = (HALF)f;
1847 m = 1;
1848 goto done;
1849 }
1850 if (n0 == 2) {
1851 A[3] = (HALF)(e >> BASEB);
1852 A[2] = (HALF)e;
1853 A[1] = (HALF)(f >> BASEB);
1854 A[0] = (HALF)f;
1855 m = 2;
1856 goto done;
1857 }
1858 u = (HALF)(s ? BASE1 : 0);
1859 if (k < BASEB) {
1860 A[m1 - 1] = (HALF)(e >> (BASEB - 1));
1861 A[m1 - 2] = ((HALF)(e << 1) | (HALF)(s > 0));
1862 A[m1 - 3] = (HALF)(f >> BASEB);
1863 A[m1 - 4] = (HALF)f;
1864 m = m1 - 2;
1865 k1 = k + 1;
1866 } else {
1867 A[m1 - 1] = 1;
1868 A[m1 - 2] = (HALF)(e >> (BASEB - 1));
1869 A[m1 - 3] = ((HALF)(e << 1) | (HALF)(s > 0));
1870 A[m1 - 4] = u;
1871 A[m1 - 5] = (HALF)(f >> BASEB);
1872 A[m1 - 6] = (HALF)f;
1873 m = m1 - 3;
1874 k1 = 1;
1875 }
1876 h = e >> k;
1877 onebit = ((e & ((FULL)1 << (k - 1))) ? TRUE : FALSE);
1878 j2 = BASEB - k1;
1879 j1 = BASEB + j2;
1880 while (m > n0) {
1881 a = A + m - 1;
1882 if (j2)
1883 f = (FULL) *a << j1 | (FULL) a[-1] << j2 | a[-2] >> k1;
1884 else
1885 f = (FULL) *a << BASEB | a[-1];
1886 if (u)
1887 f = ~f;
1888 x = f / h;
1889 if (x) {
1890 if (onebit && x > 2 * (f % h) + 2)
1891 x--;
1892 b = a + 1;
1893 i = m1 - m;
1894 a -= i + 1;
1895 if (u) {
1896 f = *a + x * (BASE - x);
1897 *a++ = (HALF)f;
1898 u = (HALF)(f >> BASEB);
1899 while (i--) {
1900 f = *a + x * *b++ + u;
1901 *a++ = (HALF)f;
1902 u = (HALF)(f >> BASEB);
1903 }
1904 u += *a;
1905 x = ~x + !u;
1906 if (!(x & TOPHALF))
1907 a[1] -= 1;
1908 } else {
1909 f = *a - x * x;
1910 *a++ = (HALF)f;
1911 u = -(HALF)(f >> BASEB);
1912 while (i--) {
1913 f = *a - x * *b++ - u;
1914 *a++ = (HALF)f;
1915 u = -(HALF)(f >> BASEB);
1916 }
1917 u = *a - u;
1918 x = x + u;
1919 if (x & TOPHALF)
1920 a[1] |= 1;
1921 }
1922 *a = ((HALF)(x << 1) | (HALF)(u > 0));
1923 } else {
1924 *a = u;
1925 }
1926 m--;
1927 if (*--a == u) {
1928 while (m > 1 && *--a == u)
1929 m--;
1930 }
1931 }
1932 i = n;
1933 a = a0;
1934 while (i--) {
1935 *a >>= 1;
1936 if (a[1] & 1) *a |= TOPHALF;
1937 a++;
1938 }
1939 s = u;
1940 done: if (s == 0) {
1941 while (m > 0 && A[m - 1] == 0)
1942 m--;
1943 if (m == 0) {
1944 remsign = 0;
1945 sqrt.v = alloc(n);
1946 sqrt.len = n;
1947 sqrt.sign = 0;
1948 memcpy(sqrt.v, a0, n * sizeof(HALF));
1949 freeh(A);
1950 *dest = sqrt;
1951 return remsign;
1952 }
1953 }
1954 if (rnd & 16) {
1955 if (s == 0) {
1956 if (m != n) {
1957 up = (m > n);
1958 } else {
1959 i = n;
1960 b = a0 + n;
1961 a = A + n;
1962 while (i > 0 && *--a == *--b)
1963 i--;
1964 up = (i > 0 && *a > *b);
1965 }
1966 } else {
1967 while (m > 1 && A[m - 1] == BASE1)
1968 m--;
1969 if (m != n) {
1970 up = (m < n);
1971 } else {
1972 i = n;
1973 b = a0 + n;
1974 a = A + n;
1975 while (i > 0 && *--a + *--b == BASE1)
1976 i--;
1977 up = ((FULL) *a + *b >= BASE);
1978 }
1979 }
1980 }
1981 else
1982 if (rnd & 8)
1983 up = (((rnd ^ *a0) & 1) ? TRUE : FALSE);
1984 else
1985 up = ((rnd & 1) ? TRUE : FALSE);
1986 if (up) {
1987 remsign = -1;
1988 i = n;
1989 a = a0;
1990 while (i-- && *a == BASE1)
1991 *a++ = 0;
1992 if (i >= 0) {
1993 (*a)++;
1994 } else {
1995 n++;
1996 *a = 1;
1997 }
1998 } else {
1999 remsign = 1;
2000 }
2001 sqrt.v = alloc(n);
2002 sqrt.len = n;
2003 sqrt.sign = 0;
2004 memcpy(sqrt.v, a0, n * sizeof(HALF));
2005 freeh(A);
2006 *dest = sqrt;
2007 return remsign;
2008
2009 }
2010
2011 /*
2012 * Take an arbitrary root of a number (to the greatest integer).
2013 * This uses the following iteration to get the K-th root of N:
2014 * x = ((K-1) * x + N / x^(K-1)) / K
2015 */
2016 void
zroot(ZVALUE z1,ZVALUE z2,ZVALUE * dest)2017 zroot(ZVALUE z1, ZVALUE z2, ZVALUE *dest)
2018 {
2019 ZVALUE ztry, quo, old, temp, temp2;
2020 ZVALUE k1; /* holds k - 1 */
2021 int sign;
2022 long i;
2023 LEN highbit, k;
2024 SIUNION sival;
2025
2026 sign = z1.sign;
2027 if (sign && ziseven(z2)) {
2028 math_error("Even root of negative number");
2029 /*NOTREACHED*/
2030 }
2031 if (ziszero(z2) || zisneg(z2)) {
2032 math_error("Non-positive root");
2033 /*NOTREACHED*/
2034 }
2035 if (ziszero(z1)) { /* root of zero */
2036 *dest = _zero_;
2037 return;
2038 }
2039 if (zisunit(z2)) { /* first root */
2040 zcopy(z1, dest);
2041 return;
2042 }
2043 if (zge31b(z2)) { /* humongous root */
2044 *dest = _one_;
2045 dest->sign = (BOOL)((HALF)sign);
2046 return;
2047 }
2048 k = (LEN)ztolong(z2);
2049 highbit = zhighbit(z1);
2050 if (highbit < k) { /* too high a root */
2051 *dest = _one_;
2052 dest->sign = (BOOL)((HALF)sign);
2053 return;
2054 }
2055 sival.ivalue = k - 1;
2056 k1.v = &sival.silow;
2057 /* ignore Saber-C warning #112 - get ushort from uint */
2058 /* OK to ignore on name zroot`sival */
2059 k1.len = 1 + (sival.sihigh != 0);
2060 k1.sign = 0;
2061 z1.sign = 0;
2062 /*
2063 * Allocate the numbers to use for the main loop.
2064 * The size and high bits of the final result are correctly set here.
2065 * Notice that the remainder of the test value is rubbish, but this
2066 * is unimportant.
2067 */
2068 highbit = (highbit + k - 1) / k;
2069 ztry.len = (highbit / BASEB) + 1;
2070 ztry.v = alloc(ztry.len);
2071 zclearval(ztry);
2072 ztry.v[ztry.len-1] = ((HALF)1 << (highbit % BASEB));
2073 ztry.sign = 0;
2074 old.v = alloc(ztry.len);
2075 old.len = 1;
2076 zclearval(old);
2077 old.sign = 0;
2078 /*
2079 * Main divide and average loop
2080 */
2081 for (;;) {
2082 zpowi(ztry, k1, &temp);
2083 zquo(z1, temp, &quo, 0);
2084 zfree(temp);
2085 i = zrel(ztry, quo);
2086 if (i <= 0) {
2087 /*
2088 * Current try is less than or equal to the root since
2089 * it is less than the quotient. If the quotient is
2090 * equal to the try, we are all done. Also, if the
2091 * try is equal to the old value, we are done since
2092 * no improvement occurred. If not, save the improved
2093 * value and loop some more.
2094 */
2095 if ((i == 0) || (zcmp(old, ztry) == 0)) {
2096 zfree(quo);
2097 zfree(old);
2098 ztry.sign = (BOOL)((HALF)sign);
2099 zquicktrim(ztry);
2100 *dest = ztry;
2101 return;
2102 }
2103 old.len = ztry.len;
2104 zcopyval(ztry, old);
2105 }
2106 /* average current try and quotient for the new try */
2107 zmul(ztry, k1, &temp);
2108 zfree(ztry);
2109 zadd(quo, temp, &temp2);
2110 zfree(temp);
2111 zfree(quo);
2112 zquo(temp2, z2, &ztry, 0);
2113 zfree(temp2);
2114 }
2115 }
2116
2117
2118 /*
2119 * Test to see if a number is an exact square or not.
2120 */
2121 BOOL
zissquare(ZVALUE z)2122 zissquare(ZVALUE z)
2123 {
2124 long n;
2125 ZVALUE tmp;
2126
2127 /* negative values are never perfect squares */
2128 if (zisneg(z)) {
2129 return FALSE;
2130 }
2131
2132 /* ignore trailing zero words */
2133 while ((z.len > 1) && (*z.v == 0)) {
2134 z.len--;
2135 z.v++;
2136 }
2137
2138 /* zero or one is a perfect square */
2139 if (zisabsleone(z)) {
2140 return TRUE;
2141 }
2142
2143 /* check mod 4096 values */
2144 if (issq_mod4k[(int)(*z.v & 0xfff)] == 0) {
2145 return FALSE;
2146 }
2147
2148 /* must do full square root test now */
2149 n = !zsqrt(z, &tmp, 0);
2150 zfree(tmp);
2151 return (n ? TRUE : FALSE);
2152 }
2153