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