1/*
2 * surd - calculate using quadratic surds of the form: a + b * sqrt(D).
3 *
4 * Copyright (C) 1999  David I. Bell
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	1990/02/15 01:50:38
21 * File existed as early as:	before 1990
22 *
23 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
24 */
25
26
27obj surd {a, b};		/* definition of the surd object */
28
29global surd_type = -1;		/* type of surd (value of D) */
30static obj surd surd__;		/* example surd for testing against */
31
32
33define surd(a,b)
34{
35	local x;
36
37	obj surd x;
38	x.a = a;
39	x.b = b;
40	return x;
41}
42
43
44define surd_print(a)
45{
46	print "surd(" : a.a : ", " : a.b : ")" :;
47}
48
49
50define surd_conj(a)
51{
52	local	x;
53
54	obj surd x;
55	x.a = a.a;
56	x.b = -a.b;
57	return x;
58}
59
60
61define surd_norm(a)
62{
63	return a.a^2 + abs(surd_type) * a.b^2;
64}
65
66
67define surd_value(a, xepsilon)
68{
69	local	epsilon;
70
71	epsilon = xepsilon;
72	if (isnull(epsilon))
73		epsilon = epsilon();
74	return a.a + a.b * sqrt(surd_type, epsilon);
75}
76
77define surd_add(a, b)
78{
79	local obj surd	x;
80
81	if (!istype(b, x)) {
82		x.a = a.a + b;
83		x.b = a.b;
84		return x;
85	}
86	if (!istype(a, x)) {
87		x.a = a + b.a;
88		x.b = b.b;
89		return x;
90	}
91	x.a = a.a + b.a;
92	x.b = a.b + b.b;
93	if (x.b)
94		return x;
95	return x.a;
96}
97
98
99define surd_sub(a, b)
100{
101	local obj surd	x;
102
103	if (!istype(b, x)) {
104		x.a = a.a - b;
105		x.b = a.b;
106		return x;
107	}
108	if (!istype(a, x)) {
109		x.a = a - b.a;
110		x.b = -b.b;
111		return x;
112	}
113	x.a = a.a - b.a;
114	x.b = a.b - b.b;
115	if (x.b)
116		return x;
117	return x.a;
118}
119
120
121define surd_inc(a)
122{
123	local	x;
124
125	x = a;
126	x.a++;
127	return x;
128}
129
130
131define surd_dec(a)
132{
133	local	x;
134
135	x = a;
136	x.a--;
137	return x;
138}
139
140
141define surd_neg(a)
142{
143	local obj surd	x;
144
145	x.a = -a.a;
146	x.b = -a.b;
147	return x;
148}
149
150
151define surd_mul(a, b)
152{
153	local obj surd	x;
154
155	if (!istype(b, x)) {
156		x.a = a.a * b;
157		x.b = a.b * b;
158	} else if (!istype(a, x)) {
159		x.a = b.a * a;
160		x.b = b.b * a;
161	} else {
162		x.a = a.a * b.a + surd_type * a.b * b.b;
163		x.b = a.a * b.b + a.b * b.a;
164	}
165	if (x.b)
166		return x;
167	return x.a;
168}
169
170
171define surd_square(a)
172{
173	local obj surd	x;
174
175	x.a = a.a^2 + a.b^2 * surd_type;
176	x.b = a.a * a.b * 2;
177	if (x.b)
178		return x;
179	return x.a;
180}
181
182
183define surd_scale(a, b)
184{
185	local obj surd	x;
186
187	x.a = scale(a.a, b);
188	x.b = scale(a.b, b);
189	return x;
190}
191
192
193define surd_shift(a, b)
194{
195	local obj surd	x;
196
197	x.a = a.a << b;
198	x.b = a.b << b;
199	if (x.b)
200		return x;
201	return x.a;
202}
203
204
205define surd_div(a, b)
206{
207	local x, y;
208
209	if ((a == 0) && b)
210		return 0;
211	obj surd x;
212	if (!istype(b, x)) {
213		x.a = a.a / b;
214		x.b = a.b / b;
215		return x;
216	}
217	y = b;
218	y.b = -b.b;
219	return (a * y) / (b.a^2 - surd_type * b.b^2);
220}
221
222
223define surd_inv(a)
224{
225	return 1 / a;
226}
227
228
229define surd_sgn(a)
230{
231	if (surd_type < 0)
232		quit "Taking sign of complex surd";
233	if (a.a == 0)
234		return sgn(a.b);
235	if (a.b == 0)
236		return sgn(a.a);
237	if ((a.a > 0) && (a.b > 0))
238		return 1;
239	if ((a.a < 0) && (a.b < 0))
240		return -1;
241	return sgn(a.a^2 - a.b^2 * surd_type) * sgn(a.a);
242}
243
244
245define surd_cmp(a, b)
246{
247	if (!istype(a, surd__))
248		return ((b.b != 0) || (a != b.a));
249	if (!istype(b, surd__))
250		return ((a.b != 0) || (b != a.a));
251	return ((a.a != b.a) || (a.b != b.b));
252}
253
254
255define surd_rel(a, b)
256{
257	local x, y;
258
259	if (surd_type < 0)
260		quit "Relative comparison of complex surds";
261	if (!istype(a, surd__)) {
262		x = a - b.a;
263		y = -b.b;
264	} else if (!istype(b, surd__)) {
265		x = a.a - b;
266		y = a.b;
267	} else {
268		x = a.a - b.a;
269		y = a.b - b.b;
270	}
271	if (y == 0)
272		return sgn(x);
273	if (x == 0)
274		return sgn(y);
275	if ((x < 0) && (y < 0))
276		return -1;
277	if ((x > 0) && (y > 0))
278		return 1;
279	return sgn(x^2 - y^2 * surd_type) * sgn(x);
280}
281
282if (config("resource_debug") & 3) {
283    print "obj surd {a, b} defined";
284    print "surd_type defined";
285    print "set surd_type as needed";
286}
287