1 A new version of mpfr is expected which isn't in Solaris yet, so we |
|
2 are removing that functionality till we can get mpfr updated and |
|
3 then will revert back these changes. |
|
4 |
|
5 Not suitable for upstream |
|
6 |
|
7 --- a/data/org.gnome.calculator.gschema.xml Fri Jan 29 15:39:20 2016 |
|
8 +++ b/data/org.gnome.calculator.gschema.xml Fri Jan 29 15:40:20 2016 |
|
9 @@ -22,6 +22,7 @@ |
|
10 <schema path="/org/gnome/calculator/" id="org.gnome.calculator" gettext-domain="calculator"> |
|
11 <key type="i" name="accuracy"> |
|
12 <default>9</default> |
|
13 + <range min="0" max="9"/> |
|
14 <summary>Accuracy value</summary> |
|
15 <description>The number of digits displayed after the numeric point</description> |
|
16 </key> |
|
17 @@ -82,10 +83,5 @@ |
|
18 <summary>Target units</summary> |
|
19 <description>Units to convert the current calculation into</description> |
|
20 </key> |
|
21 - <key type="i" name="precision"> |
|
22 - <default>2000</default> |
|
23 - <summary>Internal precision</summary> |
|
24 - <description>The internal precision used with the MPFR library</description> |
|
25 - </key> |
|
26 </schema> |
|
27 </schemalist> |
|
28 --- a/src/Makefile.am Fri Jan 29 15:40:46 2016 |
|
29 +++ b/src/Makefile.am Fri Jan 29 15:41:34 2016 |
|
30 @@ -33,8 +33,7 @@ |
|
31 --pkg gtk+-3.0 \ |
|
32 --pkg gtksourceview-3.0 \ |
|
33 --pkg libxml-2.0 \ |
|
34 - $(top_builddir)/lib/libcalculator.vapi \ |
|
35 - $(top_builddir)/lib/mpfr.vapi |
|
36 + $(top_builddir)/lib/libcalculator.vapi |
|
37 |
|
38 gnome_calculator_CPPFLAGS = \ |
|
39 $(AM_CPPFLAGS) \ |
|
40 @@ -54,8 +53,7 @@ |
|
41 --pkg gio-2.0 \ |
|
42 --pkg gtksourceview-3.0 \ |
|
43 --pkg libxml-2.0 \ |
|
44 - $(top_builddir)/lib/libcalculator.vapi \ |
|
45 - $(top_builddir)/lib/mpfr.vapi |
|
46 + $(top_builddir)/lib/libcalculator.vapi |
|
47 |
|
48 gcalccmd_LDADD = \ |
|
49 $(top_builddir)/lib/libcalculator.la |
|
50 --- a/lib/Makefile.am Fri Jan 29 15:42:13 2016 |
|
51 +++ b/lib/Makefile.am Fri Jan 29 15:42:24 2016 |
|
52 @@ -10,7 +10,6 @@ |
|
53 |
|
54 libcalculator_la_SOURCES = \ |
|
55 config.vapi \ |
|
56 - mpfr.vapi \ |
|
57 currency.vala \ |
|
58 equation.vala \ |
|
59 equation-lexer.vala \ |
|
60 @@ -37,8 +36,7 @@ |
|
61 libcalculator_la_LIBADD = \ |
|
62 $(LIBCALCULATOR_LIBS) \ |
|
63 -lgmp \ |
|
64 - -lm \ |
|
65 - -lmpfr |
|
66 + -lm |
|
67 |
|
68 EXTRA_DIST = \ |
|
69 libcalculator.h \ |
|
70 --- a/lib/equation-parser.vala Fri Jan 29 15:47:35 2016 |
|
71 +++ b/lib/equation-parser.vala Fri Jan 29 16:00:51 2016 |
|
72 @@ -78,20 +78,8 @@ |
|
73 var r = right.solve (); |
|
74 if (r == null) |
|
75 return null; |
|
76 - var z = solve_r (r); |
|
77 + return solve_r (r); |
|
78 |
|
79 - /* check for errors */ |
|
80 - Number.check_flags (); |
|
81 - if (Number.error != null) |
|
82 - { |
|
83 - var tmpleft = right; |
|
84 - var tmpright = right; |
|
85 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
86 - while (tmpright.right != null) tmpright = tmpright.right; |
|
87 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
88 - Number.error = null; |
|
89 - } |
|
90 - return z; |
|
91 } |
|
92 |
|
93 public abstract Number solve_r (Number r); |
|
94 @@ -110,20 +98,7 @@ |
|
95 var r = right.solve (); |
|
96 if (l == null || r == null) |
|
97 return null; |
|
98 - var z = solve_lr (l, r); |
|
99 - |
|
100 - /* check for errors */ |
|
101 - Number.check_flags (); |
|
102 - if (Number.error != null) |
|
103 - { |
|
104 - var tmpleft = left; |
|
105 - var tmpright = right; |
|
106 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
107 - while (tmpright.right != null) tmpright = tmpright.right; |
|
108 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
109 - Number.error = null; |
|
110 - } |
|
111 - return z; |
|
112 + return solve_lr (l, r); |
|
113 } |
|
114 |
|
115 public abstract Number solve_lr (Number left, Number r); |
|
116 @@ -261,18 +236,6 @@ |
|
117 value = value.multiply (t); |
|
118 } |
|
119 |
|
120 - /* check for errors */ |
|
121 - Number.check_flags (); |
|
122 - if (Number.error != null) |
|
123 - { |
|
124 - var tmpleft = left; |
|
125 - var tmpright = right; |
|
126 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
127 - while (tmpright.right != null) tmpright = tmpright.right; |
|
128 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
129 - Number.error = null; |
|
130 - } |
|
131 - |
|
132 return value; |
|
133 } |
|
134 } |
|
135 @@ -392,14 +355,6 @@ |
|
136 if (tmp != null) |
|
137 tmp = tmp.xpowy_integer (pow); |
|
138 |
|
139 - /* check for errors */ |
|
140 - Number.check_flags (); |
|
141 - if (Number.error != null) |
|
142 - { |
|
143 - parser.set_error (ErrorCode.MP, Number.error); |
|
144 - Number.error = null; |
|
145 - } |
|
146 - |
|
147 return tmp; |
|
148 } |
|
149 } |
|
150 @@ -574,17 +529,7 @@ |
|
151 |
|
152 public override Number solve_lr (Number l, Number r) |
|
153 { |
|
154 - var z = l.divide (r); |
|
155 - if (Number.error != null) |
|
156 - { |
|
157 - var tmpleft = right; |
|
158 - var tmpright = right; |
|
159 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
160 - while (tmpright.right != null) tmpright = tmpright.right; |
|
161 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
162 - Number.error = null; |
|
163 - } |
|
164 - return z; |
|
165 + return l.divide (r); |
|
166 } |
|
167 } |
|
168 |
|
169 @@ -604,21 +549,8 @@ |
|
170 var mod = right.solve (); |
|
171 if (base_value == null || exponent == null || mod == null) |
|
172 return null; |
|
173 - var z = base_value.modular_exponentiation (exponent, mod); |
|
174 + return base_value.modular_exponentiation (exponent, mod); |
|
175 |
|
176 - /* check for errors */ |
|
177 - Number.check_flags (); |
|
178 - if (Number.error != null) |
|
179 - { |
|
180 - var tmpleft = left; |
|
181 - var tmpright = right; |
|
182 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
183 - while (tmpright.right != null) tmpright = tmpright.right; |
|
184 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
185 - Number.error = null; |
|
186 - } |
|
187 - |
|
188 - return z; |
|
189 } |
|
190 else |
|
191 { |
|
192 @@ -626,21 +558,8 @@ |
|
193 var r = right.solve (); |
|
194 if (l == null || r == null) |
|
195 return null; |
|
196 - var z = solve_lr (l, r); |
|
197 + return solve_lr (l, r); |
|
198 |
|
199 - /* check for errors */ |
|
200 - Number.check_flags (); |
|
201 - if (Number.error != null) |
|
202 - { |
|
203 - var tmpleft = left; |
|
204 - var tmpright = right; |
|
205 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
206 - while (tmpright.right != null) tmpright = tmpright.right; |
|
207 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
208 - Number.error = null; |
|
209 - } |
|
210 - |
|
211 - return z; |
|
212 } |
|
213 } |
|
214 |
|
215 @@ -709,21 +628,8 @@ |
|
216 if (val == null) |
|
217 return null; |
|
218 |
|
219 - var z = val.xpowy_integer (pow); |
|
220 + return val.xpowy_integer (pow); |
|
221 |
|
222 - /* check for errors */ |
|
223 - Number.check_flags (); |
|
224 - if (Number.error != null) |
|
225 - { |
|
226 - var tmpleft = left; |
|
227 - var tmpright = right; |
|
228 - while (tmpleft.left != null) tmpleft = tmpleft.left; |
|
229 - while (tmpright.right != null) tmpright = tmpright.right; |
|
230 - parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index, tmpright.token.end_index); |
|
231 - Number.error = null; |
|
232 - } |
|
233 - |
|
234 - return z; |
|
235 } |
|
236 } |
|
237 |
|
238 @@ -1036,10 +942,10 @@ |
|
239 } |
|
240 |
|
241 representation_base = this.representation_base; |
|
242 - error_code = this.error; |
|
243 - error_token = this.error_token; |
|
244 - error_start = this.error_token_start; |
|
245 - error_end = this.error_token_end; |
|
246 + error_code = ErrorCode.NONE; |
|
247 + error_token = null; |
|
248 + error_start = 0; |
|
249 + error_end = 0; |
|
250 return ans; |
|
251 } |
|
252 |
|
253 --- a/lib/equation.vala Fri Jan 29 16:01:11 2016 |
|
254 +++ b/lib/equation.vala Fri Jan 29 16:02:09 2016 |
|
255 @@ -117,17 +117,15 @@ |
|
256 public new Number? parse (out uint representation_base = null, out ErrorCode error_code = null, out string? error_token = null, out uint? error_start = null, out uint? error_end = null) |
|
257 { |
|
258 var parser = new EquationParser (this, expression); |
|
259 - Number.error = null; |
|
260 + mp_clear_error(); |
|
261 |
|
262 var z = parser.parse (out representation_base, out error_code, out error_token, out error_start, out error_end); |
|
263 |
|
264 /* Error during parsing */ |
|
265 if (error_code != ErrorCode.NONE) |
|
266 - { |
|
267 return null; |
|
268 - } |
|
269 |
|
270 - if (Number.error != null) |
|
271 + if (mp_get_error() != null) |
|
272 { |
|
273 error_code = ErrorCode.MP; |
|
274 return null; |
|
275 --- a/src/gnome-calculator.vala Fri Jan 29 16:05:17 2016 |
|
276 +++ b/src/gnome-calculator.vala Fri Jan 29 16:06:27 2016 |
|
277 @@ -56,7 +56,6 @@ |
|
278 var target_currency = settings.get_string ("target-currency"); |
|
279 var source_units = settings.get_string ("source-units"); |
|
280 var target_units = settings.get_string ("target-units"); |
|
281 - var precision = settings.get_int ("precision"); |
|
282 |
|
283 var equation = new MathEquation (); |
|
284 equation.accuracy = accuracy; |
|
285 @@ -69,7 +68,6 @@ |
|
286 equation.target_currency = target_currency; |
|
287 equation.source_units = source_units; |
|
288 equation.target_units = target_units; |
|
289 - Number.precision = precision; |
|
290 |
|
291 add_action_entries (app_entries, this); |
|
292 |
|
293 @@ -173,7 +171,7 @@ |
|
294 } |
|
295 else if (error == ErrorCode.MP) |
|
296 { |
|
297 - stderr.printf ("Error: %s\n", Number.error); |
|
298 + stderr.printf ("Error: %s\n", mp_get_error()); |
|
299 return Posix.EXIT_FAILURE; |
|
300 } |
|
301 else |
|
302 --- a/src/gcalccmd.vala Fri Jan 29 16:03:42 2016 |
|
303 +++ b/src/gcalccmd.vala Fri Jan 29 16:04:56 2016 |
|
304 @@ -34,18 +34,9 @@ |
|
305 |
|
306 result_serializer.set_representation_base (representation_base); |
|
307 if (z != null) |
|
308 - { |
|
309 - var str = result_serializer.to_string (z); |
|
310 - if (result_serializer.error != null) |
|
311 - { |
|
312 - stderr.printf ("%s\n", result_serializer.error); |
|
313 - result_serializer.error = null; |
|
314 - } |
|
315 - else |
|
316 - stdout.printf ("%s\n", str); |
|
317 - } |
|
318 + stdout.printf ("%s\n", result_serializer.to_string (z)); |
|
319 else if (ret == ErrorCode.MP) |
|
320 - stderr.printf ("Error %s\n", Number.error); |
|
321 + stderr.printf ("Error %s\n", mp_get_error()); |
|
322 else |
|
323 stderr.printf ("Error %d\n", ret); |
|
324 } |
|
325 --- a/lib/math-equation.vala Fri Jan 29 16:33:36 2016 |
|
326 +++ b/lib/math-equation.vala Fri Jan 29 16:39:56 2016 |
|
327 @@ -786,11 +786,6 @@ |
|
328 ans_end_mark = create_mark (null, end, true); |
|
329 apply_tag (ans_tag, start, end); |
|
330 |
|
331 - if (serializer.error != null) |
|
332 - { |
|
333 - status = serializer.error; |
|
334 - serializer.error = null; |
|
335 - } |
|
336 } |
|
337 |
|
338 public new void insert (string text) |
|
339 @@ -964,13 +959,11 @@ |
|
340 break; |
|
341 |
|
342 case ErrorCode.MP: |
|
343 - if (Number.error != null) // LEGACY, should never be run |
|
344 + if (mp_get_error () != null) |
|
345 + solvedata.error = mp_get_error (); |
|
346 + else if (error_token != null) |
|
347 { |
|
348 - solvedata.error = Number.error; |
|
349 - } |
|
350 - else if (error_token != null) // should always be run |
|
351 - { |
|
352 - solvedata.error = _("%s").printf (error_token); |
|
353 + solvedata.error = _("Malformed expression at token '%s'").printf(error_token); |
|
354 solvedata.error_start = error_start; |
|
355 solvedata.error_end = error_end; |
|
356 } |
|
357 @@ -1015,8 +1008,6 @@ |
|
358 |
|
359 /* Fix thousand separator offsets in the start and end offsets of error token. */ |
|
360 error_token_fix_thousands_separator (); |
|
361 - /* Fix missing Parenthesis before the start and after the end offsets of error token */ |
|
362 - error_token_fix_parenthesis (); |
|
363 |
|
364 /* Notify the GUI about the change in error token locations. */ |
|
365 notify_property ("error-token-end"); |
|
366 @@ -1087,70 +1078,6 @@ |
|
367 end.forward_chars (length); |
|
368 } |
|
369 } |
|
370 - |
|
371 - /* Fix the offsets to consider starting and ending parenthesis */ |
|
372 - private void error_token_fix_parenthesis () |
|
373 - { |
|
374 - unichar c; |
|
375 - int count = 0; |
|
376 - int real_end = display.index_of_nth_char (error_token_end); |
|
377 - int real_start = display.index_of_nth_char (error_token_start); |
|
378 - |
|
379 - /* checks if there are more opening/closing parenthesis than closing/opening parenthesis */ |
|
380 - for (int i = real_start; display.get_next_char (ref i, out c) && i <= real_end;) |
|
381 - { |
|
382 - if (c.to_string () == "(") count++; |
|
383 - if (c.to_string () == ")") count--; |
|
384 - } |
|
385 - |
|
386 - /* if there are more opening than closing parenthesis and there are closing parenthesis |
|
387 - after the end offset, include those in the offsets */ |
|
388 - for (int i = real_end; display.get_next_char (ref i, out c) && count > 0;) |
|
389 - { |
|
390 - if (c.to_string () == ")") |
|
391 - { |
|
392 - state.error_token_end++; |
|
393 - count--; |
|
394 - } |
|
395 - else |
|
396 - { |
|
397 - break; |
|
398 - } |
|
399 - } |
|
400 - |
|
401 - /* the same for closing parenthesis */ |
|
402 - for (int i = real_start; display.get_prev_char (ref i, out c) && count < 0;) |
|
403 - { |
|
404 - if (c.to_string () == "(") |
|
405 - { |
|
406 - state.error_token_start--; |
|
407 - count++; |
|
408 - } |
|
409 - else |
|
410 - { |
|
411 - break; |
|
412 - } |
|
413 - } |
|
414 - |
|
415 - real_end = display.index_of_nth_char (error_token_end); |
|
416 - real_start = display.index_of_nth_char (error_token_start); |
|
417 - |
|
418 - unichar d; |
|
419 - |
|
420 - /* if there are opening parenthesis directly before aswell as closing parenthesis directly after the offsets, include those aswell */ |
|
421 - while (display.get_next_char (ref real_end, out d) && display.get_prev_char (ref real_start, out c)) |
|
422 - { |
|
423 - if (c.to_string () == "(" && d.to_string () == ")") |
|
424 - { |
|
425 - state.error_token_start--; |
|
426 - state.error_token_end++; |
|
427 - } |
|
428 - else |
|
429 - { |
|
430 - break; |
|
431 - } |
|
432 - } |
|
433 - } |
|
434 |
|
435 private void* factorize_real () |
|
436 { |
|
437 --- a/src/math-preferences.vala Fri Jan 29 16:40:37 2016 |
|
438 +++ b/src/math-preferences.vala Fri Jan 29 16:43:30 2016 |
|
439 @@ -90,7 +90,7 @@ |
|
440 var string = _("Show %d decimal _places"); |
|
441 var tokens = string.split ("%d", 2); |
|
442 |
|
443 - var decimal_places_adjustment = new Gtk.Adjustment (0.0, 0.0, 100.0, 1.0, 1.0, 0.0); |
|
444 + var decimal_places_adjustment = new Gtk.Adjustment (0.0, 0.0, 9.0, 1.0, 1.0, 0.0); |
|
445 decimal_places_spin = new Gtk.SpinButton (decimal_places_adjustment, 0.0, 0); |
|
446 |
|
447 if (tokens.length > 0) |
|
448 --- a/lib/serializer.vala Mon Feb 1 10:45:23 2016 |
|
449 +++ b/lib/serializer.vala Mon Feb 1 10:46:37 2016 |
|
450 @@ -32,9 +32,6 @@ |
|
451 private unichar tsep; /* Locale specific thousands separator. */ |
|
452 private int tsep_count; /* Number of digits between separator. */ |
|
453 |
|
454 - /* is set when an error (for example precision error while converting) occurs */ |
|
455 - public string? error { get; set; default = null; } |
|
456 - |
|
457 public Serializer (DisplayFormat format, int number_base, int trailing_digits) |
|
458 { |
|
459 var radix_string = Posix.nl_langinfo (Posix.NLItem.RADIXCHAR); |
|
460 @@ -322,17 +319,8 @@ |
|
461 |
|
462 var d = t3.to_integer (); |
|
463 |
|
464 - if (d < 16 && d >= 0) |
|
465 - { |
|
466 - string.prepend_c (digits[d]); |
|
467 - } |
|
468 - else |
|
469 - { |
|
470 - string.prepend_c ('?'); |
|
471 - error = _("Precision error"); |
|
472 - string.assign ("0"); |
|
473 - break; |
|
474 - } |
|
475 + string.prepend_c (d < 16 ? digits[d] : '?'); |
|
476 + |
|
477 n_digits++; |
|
478 |
|
479 temp = t; |
|
480 --- a/tests/test-number.vala Mon Feb 1 10:48:08 2016 |
|
481 +++ b/tests/test-number.vala Mon Feb 1 10:49:35 2016 |
|
482 @@ -72,6 +72,22 @@ |
|
483 pass (); |
|
484 } |
|
485 |
|
486 +private void test_float () |
|
487 +{ |
|
488 + for (var a = -10.0f; a <= 10.0f; a += 0.5f) |
|
489 + { |
|
490 + var z = new Number.float (a); |
|
491 + if (z.to_float () != a) |
|
492 + { |
|
493 + fail ("Number.float (%f).to_float () -> %f, expected %f".printf (a, z.to_float (), a)); |
|
494 + return; |
|
495 + } |
|
496 + } |
|
497 + |
|
498 + pass (); |
|
499 +} |
|
500 + |
|
501 + |
|
502 private void test_double () |
|
503 { |
|
504 for (var a = -10.0; a <= 10.0; a += 0.5) |
|
505 @@ -1074,6 +1090,7 @@ |
|
506 test_integer (); |
|
507 test_unsigned_integer (); |
|
508 test_fraction (); |
|
509 + test_float (); |
|
510 test_double (); |
|
511 test_complex (); |
|
512 test_polar (); |
|
513 --- a/tests/test-equation.c Mon Feb 1 12:22:33 2016 |
|
514 +++ b/tests/test-equation.c Mon Feb 1 12:22:47 2016 |
|
515 @@ -95,7 +95,7 @@ |
|
516 const gchar* _tmp1_ = NULL; |
|
517 const gchar* _tmp2_ = NULL; |
|
518 gchar* _tmp3_ = NULL; |
|
519 - _tmp1_ = number_get_error (); |
|
520 + _tmp1_ = mp_get_error (); |
|
521 _tmp2_ = _tmp1_; |
|
522 _tmp3_ = g_strdup_printf ("ErrorCode.MP(\"%s\")", _tmp2_); |
|
523 result = _tmp3_; |
|
524 |
|
525 --- a/lib/number.vala Fri Jan 29 16:44:36 2016 |
|
526 +++ b/lib/number.vala Mon Feb 1 12:02:40 2016 |
|
527 @@ -27,8 +27,15 @@ |
|
528 * THE MP USERS GUIDE. |
|
529 */ |
|
530 |
|
531 -using MPFR; |
|
532 |
|
533 +/* Size of the multiple precision values */ |
|
534 +private const int SIZE = 1000; |
|
535 + |
|
536 +/* Base for numbers */ |
|
537 +private const int BASE = 10000; |
|
538 + |
|
539 +private const int T = 100; |
|
540 + |
|
541 private delegate int BitwiseFunc (int v1, int v2); |
|
542 |
|
543 public enum AngleUnit |
|
544 @@ -41,37 +48,69 @@ |
|
545 /* Object for a high precision floating point number representation */ |
|
546 public class Number : Object |
|
547 { |
|
548 - /* real and imaginary part of a Number */ |
|
549 - private MPFloat re_num { get; set; } |
|
550 - private MPFloat im_num { get; set; } |
|
551 + /* Sign (+1, -1) or 0 for the value zero */ |
|
552 + public int re_sign; |
|
553 + public int im_sign; |
|
554 |
|
555 - public static ulong precision { get; set; default = 1000; } |
|
556 + /* Exponent (to base BASE) */ |
|
557 + public int re_exponent; |
|
558 + public int im_exponent; |
|
559 |
|
560 - /* Stores the error msg if an error occurs during calculation. Otherwise should be null */ |
|
561 - public static string? error { get; set; default = null; } |
|
562 + /* Normalized fraction */ |
|
563 + public int re_fraction[1000]; // SIZE |
|
564 + public int im_fraction[1000]; // SIZE |
|
565 |
|
566 public Number.integer (int64 value) |
|
567 { |
|
568 - var tmp = MPFloat.init2 ((Precision) precision); |
|
569 - tmp.set_signed_integer ((long)value, Round.NEAREST); |
|
570 - re_num = tmp; |
|
571 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
572 - tmp2.set_unsigned_integer (0, Round.NEAREST); |
|
573 - im_num = tmp2; |
|
574 + if (value < 0) |
|
575 + { |
|
576 + value = -value; |
|
577 + re_sign = -1; |
|
578 + } |
|
579 + else if (value > 0) |
|
580 + re_sign = 1; |
|
581 + else |
|
582 + re_sign = 0; |
|
583 + |
|
584 + while (value != 0) |
|
585 + { |
|
586 + re_fraction[re_exponent] = (int) (value % BASE); |
|
587 + re_exponent++; |
|
588 + value /= BASE; |
|
589 + } |
|
590 + for (var i = 0; i < re_exponent / 2; i++) |
|
591 + { |
|
592 + int t = re_fraction[i]; |
|
593 + re_fraction[i] = re_fraction[re_exponent - i - 1]; |
|
594 + re_fraction[re_exponent - i - 1] = t; |
|
595 + } |
|
596 } |
|
597 |
|
598 public Number.unsigned_integer (uint64 x) |
|
599 { |
|
600 - var tmp = MPFloat.init2 ((Precision) precision); |
|
601 - tmp.set_unsigned_integer ((ulong)x, Round.NEAREST); |
|
602 - re_num = tmp; |
|
603 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
604 - tmp2.set_unsigned_integer (0, Round.NEAREST); |
|
605 - im_num = tmp2; |
|
606 + if (x == 0) |
|
607 + re_sign = 0; |
|
608 + else |
|
609 + re_sign = 1; |
|
610 + |
|
611 + while (x != 0) |
|
612 + { |
|
613 + re_fraction[re_exponent] = (int) (x % BASE); |
|
614 + x = x / BASE; |
|
615 + re_exponent++; |
|
616 + } |
|
617 + for (var i = 0; i < re_exponent / 2; i++) |
|
618 + { |
|
619 + int t = re_fraction[i]; |
|
620 + re_fraction[i] = re_fraction[re_exponent - i - 1]; |
|
621 + re_fraction[re_exponent - i - 1] = t; |
|
622 + } |
|
623 } |
|
624 |
|
625 public Number.fraction (int64 numerator, int64 denominator) |
|
626 { |
|
627 + mp_gcd (ref numerator, ref denominator); |
|
628 + |
|
629 if (denominator < 0) |
|
630 { |
|
631 numerator = -numerator; |
|
632 @@ -81,39 +120,203 @@ |
|
633 Number.integer (numerator); |
|
634 if (denominator != 1) |
|
635 { |
|
636 - var tmp = re_num; |
|
637 - tmp.divide_signed_integer (re_num, (long) denominator, Round.NEAREST); |
|
638 - re_num = tmp; |
|
639 + var z = divide_integer (denominator); |
|
640 + re_sign = z.re_sign; |
|
641 + im_sign = z.im_sign; |
|
642 + re_exponent = z.re_exponent; |
|
643 + im_exponent = z.im_exponent; |
|
644 + for (var i = 0; i < z.re_fraction.length; i++) |
|
645 + { |
|
646 + re_fraction[i] = z.re_fraction[i]; |
|
647 + im_fraction[i] = z.im_fraction[i]; |
|
648 + } |
|
649 } |
|
650 } |
|
651 |
|
652 - /* Helper constructor. Creates new Number from already existing MPFloat. */ |
|
653 - public Number.mpfloat (MPFloat value) |
|
654 + public Number.float (float value) |
|
655 { |
|
656 - re_num = value; |
|
657 - var tmp = MPFloat.init2 ((Precision) precision); |
|
658 - tmp.set_unsigned_integer (0, Round.NEAREST); |
|
659 - im_num = tmp; |
|
660 + var z = new Number.integer (0); |
|
661 + |
|
662 + if (value != 0) |
|
663 + { |
|
664 + /* CHECK SIGN */ |
|
665 + var rj = 0f; |
|
666 + if (value < 0.0f) |
|
667 + { |
|
668 + z.re_sign = -1; |
|
669 + rj = -value; |
|
670 + } |
|
671 + else if (value > 0.0f) |
|
672 + { |
|
673 + z.re_sign = 1; |
|
674 + rj = value; |
|
675 + } |
|
676 + |
|
677 + /* INCREASE IE AND DIVIDE RJ BY 16. */ |
|
678 + var ie = 0; |
|
679 + while (rj >= 1.0f) |
|
680 + { |
|
681 + ie++; |
|
682 + rj *= 0.0625f; |
|
683 + } |
|
684 + while (rj < 0.0625f) |
|
685 + { |
|
686 + ie--; |
|
687 + rj *= 16.0f; |
|
688 + } |
|
689 + |
|
690 + /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16. |
|
691 + * SET re_exponent TO 0 |
|
692 + */ |
|
693 + z.re_exponent = 0; |
|
694 + |
|
695 + /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */ |
|
696 + for (var i = 0; i < T + 4; i++) |
|
697 + { |
|
698 + rj *= BASE; |
|
699 + z.re_fraction[i] = (int) rj; |
|
700 + rj -= z.re_fraction[i]; |
|
701 + } |
|
702 + |
|
703 + /* NORMALIZE RESULT */ |
|
704 + mp_normalize (ref z); |
|
705 + |
|
706 + /* Computing MAX */ |
|
707 + var ib = int.max (BASE * 7 * BASE, 32767) / 16; |
|
708 + var tp = 1; |
|
709 + |
|
710 + /* NOW MULTIPLY BY 16**IE */ |
|
711 + if (ie < 0) |
|
712 + { |
|
713 + var k = -ie; |
|
714 + for (var i = 1; i <= k; i++) |
|
715 + { |
|
716 + tp <<= 4; |
|
717 + if (tp <= ib && tp != BASE && i < k) |
|
718 + continue; |
|
719 + z = z.divide_integer (tp); |
|
720 + tp = 1; |
|
721 + } |
|
722 + } |
|
723 + else if (ie > 0) |
|
724 + { |
|
725 + for (var i = 1; i <= ie; i++) |
|
726 + { |
|
727 + tp <<= 4; |
|
728 + if (tp <= ib && tp != BASE && i < ie) |
|
729 + continue; |
|
730 + z = z.multiply_integer (tp); |
|
731 + tp = 1; |
|
732 + } |
|
733 + } |
|
734 + } |
|
735 + |
|
736 + re_sign = z.re_sign; |
|
737 + im_sign = z.im_sign; |
|
738 + re_exponent = z.re_exponent; |
|
739 + im_exponent = z.im_exponent; |
|
740 + for (var i = 0; i < z.re_fraction.length; i++) |
|
741 + { |
|
742 + re_fraction[i] = z.re_fraction[i]; |
|
743 + im_fraction[i] = z.im_fraction[i]; |
|
744 + } |
|
745 } |
|
746 |
|
747 public Number.double (double value) |
|
748 { |
|
749 - var tmp = MPFloat.init2 ((Precision) precision); |
|
750 - tmp.set_double (value, Round.NEAREST); |
|
751 - re_num = tmp; |
|
752 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
753 - tmp2.set_unsigned_integer (0, Round.NEAREST); |
|
754 - im_num = tmp2; |
|
755 + var z = new Number.integer (0); |
|
756 + |
|
757 + if (value != 0) |
|
758 + { |
|
759 + /* CHECK SIGN */ |
|
760 + var dj = 0.0; |
|
761 + if (value < 0.0) |
|
762 + { |
|
763 + z.re_sign = -1; |
|
764 + dj = -value; |
|
765 + } |
|
766 + else if (value > 0.0) |
|
767 + { |
|
768 + z.re_sign = 1; |
|
769 + dj = value; |
|
770 + } |
|
771 + |
|
772 + /* INCREASE IE AND DIVIDE DJ BY 16. */ |
|
773 + var ie = 0; |
|
774 + for (ie = 0; dj >= 1.0; ie++) |
|
775 + dj *= 1.0/16.0; |
|
776 + |
|
777 + for ( ; dj < 1.0/16.0; ie--) |
|
778 + dj *= 16.0; |
|
779 + |
|
780 + /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16 |
|
781 + * SET re_exponent TO 0 |
|
782 + */ |
|
783 + z.re_exponent = 0; |
|
784 + |
|
785 + /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */ |
|
786 + for (var i = 0; i < T + 4; i++) |
|
787 + { |
|
788 + dj *= (double) BASE; |
|
789 + z.re_fraction[i] = (int) dj; |
|
790 + dj -= (double) z.re_fraction[i]; |
|
791 + } |
|
792 + |
|
793 + /* NORMALIZE RESULT */ |
|
794 + mp_normalize (ref z); |
|
795 + |
|
796 + /* Computing MAX */ |
|
797 + var ib = int.max (BASE * 7 * BASE, 32767) / 16; |
|
798 + var tp = 1; |
|
799 + |
|
800 + /* NOW MULTIPLY BY 16**IE */ |
|
801 + if (ie < 0) |
|
802 + { |
|
803 + var k = -ie; |
|
804 + for (var i = 1; i <= k; ++i) |
|
805 + { |
|
806 + tp <<= 4; |
|
807 + if (tp <= ib && tp != BASE && i < k) |
|
808 + continue; |
|
809 + z = z.divide_integer (tp); |
|
810 + tp = 1; |
|
811 + } |
|
812 + } |
|
813 + else if (ie > 0) |
|
814 + { |
|
815 + for (var i = 1; i <= ie; ++i) |
|
816 + { |
|
817 + tp <<= 4; |
|
818 + if (tp <= ib && tp != BASE && i < ie) |
|
819 + continue; |
|
820 + z = z.multiply_integer (tp); |
|
821 + tp = 1; |
|
822 + } |
|
823 + } |
|
824 + } |
|
825 + |
|
826 + re_sign = z.re_sign; |
|
827 + im_sign = z.im_sign; |
|
828 + re_exponent = z.re_exponent; |
|
829 + im_exponent = z.im_exponent; |
|
830 + for (var i = 0; i < z.re_fraction.length; i++) |
|
831 + { |
|
832 + re_fraction[i] = z.re_fraction[i]; |
|
833 + im_fraction[i] = z.im_fraction[i]; |
|
834 + } |
|
835 } |
|
836 |
|
837 public Number.complex (Number x, Number y) |
|
838 { |
|
839 - var tmp = MPFloat.init2 ((Precision) precision); |
|
840 - tmp.@set (x.re_num, Round.NEAREST); |
|
841 - re_num = tmp; |
|
842 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
843 - tmp2.@set (y.re_num, Round.NEAREST); |
|
844 - im_num = tmp2; |
|
845 + re_sign = x.re_sign; |
|
846 + re_exponent = x.re_exponent; |
|
847 + for (var i = 0; i < im_fraction.length; i++) |
|
848 + re_fraction[i] = x.re_fraction[i]; |
|
849 + |
|
850 + im_sign = y.re_sign; |
|
851 + im_exponent = y.re_exponent; |
|
852 + for (var i = 0; i < im_fraction.length; i++) |
|
853 + im_fraction[i] = y.re_fraction[i]; |
|
854 } |
|
855 |
|
856 public Number.polar (Number r, Number theta, AngleUnit unit = AngleUnit.RADIANS) |
|
857 @@ -125,35 +328,28 @@ |
|
858 |
|
859 public Number.eulers () |
|
860 { |
|
861 - var tmp = MPFloat.init2 ((Precision) precision); |
|
862 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
863 - tmp2.set_unsigned_integer (1, Round.NEAREST); |
|
864 - /* e^1, since mpfr doesn't have a function to return e */ |
|
865 - tmp.exp (tmp2, Round.NEAREST); |
|
866 - re_num = tmp; |
|
867 - var tmp3 = MPFloat.init2 ((Precision) precision); |
|
868 - tmp3.set_unsigned_integer (0, Round.NEAREST); |
|
869 - im_num = tmp3; |
|
870 + var z = new Number.integer (1).epowy (); |
|
871 + re_sign = z.re_sign; |
|
872 + im_sign = z.im_sign; |
|
873 + re_exponent = z.re_exponent; |
|
874 + im_exponent = z.im_exponent; |
|
875 + for (var i = 0; i < z.re_fraction.length; i++) |
|
876 + { |
|
877 + re_fraction[i] = z.re_fraction[i]; |
|
878 + im_fraction[i] = z.im_fraction[i]; |
|
879 + } |
|
880 } |
|
881 |
|
882 public Number.i () |
|
883 { |
|
884 - var tmp = MPFloat.init2 ((Precision) precision); |
|
885 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
886 - tmp.set_unsigned_integer (0, Round.NEAREST); |
|
887 - tmp2.set_unsigned_integer (1, Round.NEAREST); |
|
888 - re_num = tmp; |
|
889 - im_num = tmp2; |
|
890 + im_sign = 1; |
|
891 + im_exponent = 1; |
|
892 + im_fraction[0] = 1; |
|
893 } |
|
894 |
|
895 public Number.pi () |
|
896 { |
|
897 - var tmp = MPFloat.init2 ((Precision) precision); |
|
898 - tmp.const_pi (Round.NEAREST); |
|
899 - re_num = tmp; |
|
900 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
901 - tmp2.set_unsigned_integer (0, Round.NEAREST); |
|
902 - im_num = tmp2; |
|
903 + Number.double (Math.PI); |
|
904 } |
|
905 |
|
906 /* Sets z to be a uniform random number in the range [0, 1] */ |
|
907 @@ -162,42 +358,123 @@ |
|
908 Number.double (Random.next_double ()); |
|
909 } |
|
910 |
|
911 - ~Number () |
|
912 - { |
|
913 - re_num.clear (); |
|
914 - im_num.clear (); |
|
915 - } |
|
916 - |
|
917 public int64 to_integer () |
|
918 { |
|
919 - return re_num.get_signed_integer (Round.NEAREST); |
|
920 + int64 z = 0; |
|
921 + |
|
922 + /* |x| <= 1 */ |
|
923 + if (re_sign == 0 || re_exponent <= 0) |
|
924 + return 0; |
|
925 + |
|
926 + /* Multiply digits together */ |
|
927 + for (var i = 0; i < re_exponent; i++) |
|
928 + { |
|
929 + var t = z; |
|
930 + z = z * BASE + re_fraction[i]; |
|
931 + |
|
932 + /* Check for overflow */ |
|
933 + if (z <= t) |
|
934 + return 0; |
|
935 + } |
|
936 + |
|
937 + /* Validate result */ |
|
938 + var v = z; |
|
939 + for (var i = re_exponent - 1; i >= 0; i--) |
|
940 + { |
|
941 + /* Get last digit */ |
|
942 + var digit = v - (v / BASE) * BASE; |
|
943 + if (re_fraction[i] != digit) |
|
944 + return 0; |
|
945 + |
|
946 + v /= BASE; |
|
947 + } |
|
948 + if (v != 0) |
|
949 + return 0; |
|
950 + |
|
951 + return re_sign * z; |
|
952 } |
|
953 |
|
954 public uint64 to_unsigned_integer () |
|
955 { |
|
956 - return re_num.get_unsigned_integer (Round.NEAREST); |
|
957 + /* x <= 1 */ |
|
958 + if (re_sign <= 0 || re_exponent <= 0) |
|
959 + return 0; |
|
960 + |
|
961 + /* Multiply digits together */ |
|
962 + uint64 z = 0; |
|
963 + for (var i = 0; i < re_exponent; i++) |
|
964 + { |
|
965 + var t = z; |
|
966 + z = z * BASE + re_fraction[i]; |
|
967 + |
|
968 + /* Check for overflow */ |
|
969 + if (z <= t) |
|
970 + return 0; |
|
971 + } |
|
972 + |
|
973 + /* Validate result */ |
|
974 + var v = z; |
|
975 + for (var i = re_exponent - 1; i >= 0; i--) |
|
976 + { |
|
977 + /* Get last digit */ |
|
978 + var digit = (uint64) v - (v / BASE) * BASE; |
|
979 + if (re_fraction[i] != digit) |
|
980 + return 0; |
|
981 + |
|
982 + v /= BASE; |
|
983 + } |
|
984 + if (v != 0) |
|
985 + return 0; |
|
986 + |
|
987 + return z; |
|
988 } |
|
989 |
|
990 public float to_float () |
|
991 { |
|
992 - return re_num.get_float (Round.NEAREST); |
|
993 + if (is_zero ()) |
|
994 + return 0f; |
|
995 + |
|
996 + var z = 0f; |
|
997 + for (var i = 0; i < T; i++) |
|
998 + { |
|
999 + if (re_fraction[i] != 0) |
|
1000 + z += re_fraction[i] * Math.powf (BASE, re_exponent - i - 1); |
|
1001 + } |
|
1002 + |
|
1003 + if (re_sign < 0) |
|
1004 + return -z; |
|
1005 + else |
|
1006 + return z; |
|
1007 } |
|
1008 |
|
1009 public double to_double () |
|
1010 { |
|
1011 - return re_num.get_double (Round.NEAREST); |
|
1012 + if (is_zero ()) |
|
1013 + return 0d; |
|
1014 + |
|
1015 + var z = 0d; |
|
1016 + for (var i = 0; i < T; i++) |
|
1017 + { |
|
1018 + if (re_fraction[i] != 0) |
|
1019 + z += re_fraction[i] * Math.pow (BASE, re_exponent - i - 1); |
|
1020 + } |
|
1021 + |
|
1022 + if (re_sign < 0) |
|
1023 + return -z; |
|
1024 + else |
|
1025 + return z; |
|
1026 } |
|
1027 |
|
1028 /* Return true if the value is x == 0 */ |
|
1029 public bool is_zero () |
|
1030 { |
|
1031 - return re_num.is_zero () && im_num.is_zero (); |
|
1032 + return re_sign == 0 && im_sign == 0; |
|
1033 } |
|
1034 |
|
1035 /* Return true if x < 0 */ |
|
1036 public bool is_negative () |
|
1037 { |
|
1038 - return re_num.sgn () < 0; |
|
1039 + return re_sign < 0; |
|
1040 } |
|
1041 |
|
1042 /* Return true if x is integer */ |
|
1043 @@ -205,8 +482,35 @@ |
|
1044 { |
|
1045 if (is_complex ()) |
|
1046 return false; |
|
1047 + /* Multiplication and division by 10000 is used to get around a |
|
1048 + * limitation to the "fix" for Sun bugtraq bug #4006391 in the |
|
1049 + * floor () routine in mp.c, when the re_exponent is less than 1. |
|
1050 + */ |
|
1051 + var t3 = new Number.integer (10000); |
|
1052 + var t1 = multiply (t3); |
|
1053 + t1 = t1.divide (t3); |
|
1054 + var t2 = t1.floor (); |
|
1055 + return t1.equals (t2); |
|
1056 |
|
1057 - return re_num.is_integer () != 0; |
|
1058 + /* Correct way to check for integer */ |
|
1059 + /* |
|
1060 + |
|
1061 + // Zero is an integer |
|
1062 + if (is_zero ()) |
|
1063 + return true; |
|
1064 + |
|
1065 + // fractional |
|
1066 + if (re_exponent <= 0) |
|
1067 + return false; |
|
1068 + |
|
1069 + // Look for fractional components |
|
1070 + for (var i = re_exponent; i < SIZE; i++) |
|
1071 + { |
|
1072 + if (re_fraction[i] != 0) |
|
1073 + return false; |
|
1074 + } |
|
1075 + |
|
1076 + return true;*/ |
|
1077 } |
|
1078 |
|
1079 /* Return true if x is a positive integer */ |
|
1080 @@ -215,7 +519,7 @@ |
|
1081 if (is_complex ()) |
|
1082 return false; |
|
1083 else |
|
1084 - return re_num.sgn () >= 0 && is_integer (); |
|
1085 + return re_sign >= 0 && is_integer (); |
|
1086 } |
|
1087 |
|
1088 /* Return true if x is a natural number (an integer ≥ 0) */ |
|
1089 @@ -224,34 +528,19 @@ |
|
1090 if (is_complex ()) |
|
1091 return false; |
|
1092 else |
|
1093 - return re_num.sgn () > 0 && is_integer (); |
|
1094 + return re_sign > 0 && is_integer (); |
|
1095 } |
|
1096 |
|
1097 /* Return true if x has an imaginary component */ |
|
1098 public bool is_complex () |
|
1099 { |
|
1100 - return !im_num.is_zero (); |
|
1101 + return im_sign != 0; |
|
1102 } |
|
1103 |
|
1104 - /* Return error if overflow or underflow */ |
|
1105 - public static void check_flags () |
|
1106 - { |
|
1107 - if (mpfr_is_underflow () != 0) |
|
1108 - { |
|
1109 - /* Translators: Error displayed when underflow error occured */ |
|
1110 - error = _("Underflow error"); |
|
1111 - } |
|
1112 - else if (mpfr_is_overflow () != 0) |
|
1113 - { |
|
1114 - /* Translators: Error displayed when overflow error occured */ |
|
1115 - error = _("Overflow error"); |
|
1116 - } |
|
1117 - } |
|
1118 - |
|
1119 /* Return true if x == y */ |
|
1120 public bool equals (Number y) |
|
1121 { |
|
1122 - return re_num.is_equal (y.re_num); |
|
1123 + return compare (y) == 0; |
|
1124 } |
|
1125 |
|
1126 /* Returns: |
|
1127 @@ -261,25 +550,62 @@ |
|
1128 */ |
|
1129 public int compare (Number y) |
|
1130 { |
|
1131 - return re_num.cmp (y.re_num); |
|
1132 + if (re_sign != y.re_sign) |
|
1133 + { |
|
1134 + if (re_sign > y.re_sign) |
|
1135 + return 1; |
|
1136 + else |
|
1137 + return -1; |
|
1138 + } |
|
1139 + |
|
1140 + /* x = y = 0 */ |
|
1141 + if (is_zero ()) |
|
1142 + return 0; |
|
1143 + |
|
1144 + /* See if numbers are of different magnitude */ |
|
1145 + if (re_exponent != y.re_exponent) |
|
1146 + { |
|
1147 + if (re_exponent > y.re_exponent) |
|
1148 + return re_sign; |
|
1149 + else |
|
1150 + return -re_sign; |
|
1151 + } |
|
1152 + |
|
1153 + /* Compare fractions */ |
|
1154 + for (var i = 0; i < SIZE; i++) |
|
1155 + { |
|
1156 + if (re_fraction[i] == y.re_fraction[i]) |
|
1157 + continue; |
|
1158 + |
|
1159 + if (re_fraction[i] > y.re_fraction[i]) |
|
1160 + return re_sign; |
|
1161 + else |
|
1162 + return -re_sign; |
|
1163 + } |
|
1164 + |
|
1165 + /* x = y */ |
|
1166 + return 0; |
|
1167 } |
|
1168 |
|
1169 /* Sets z = sgn (x) */ |
|
1170 public Number sgn () |
|
1171 { |
|
1172 - var z = new Number.integer (re_num.sgn ()); |
|
1173 - return z; |
|
1174 + if (is_zero ()) |
|
1175 + return new Number.integer (0); |
|
1176 + else if (is_negative ()) |
|
1177 + return new Number.integer (-1); |
|
1178 + else |
|
1179 + return new Number.integer (1); |
|
1180 } |
|
1181 |
|
1182 /* Sets z = −x */ |
|
1183 public Number invert_sign () |
|
1184 { |
|
1185 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1186 - tmp.neg (re_num, Round.NEAREST); |
|
1187 - var z = new Number.mpfloat (tmp); |
|
1188 - var tmp_im = z.im_num; |
|
1189 - tmp_im.neg (im_num, Round.NEAREST); |
|
1190 - z.im_num = tmp_im; |
|
1191 + var z = copy (); |
|
1192 + |
|
1193 + z.re_sign = -z.re_sign; |
|
1194 + z.im_sign = -z.im_sign; |
|
1195 + |
|
1196 return z; |
|
1197 } |
|
1198 |
|
1199 @@ -294,13 +620,13 @@ |
|
1200 x_real = x_real.multiply (x_real); |
|
1201 x_im = x_im.multiply (x_im); |
|
1202 var z = x_real.add (x_im); |
|
1203 - return z.sqrt (); |
|
1204 + return z.root (2); |
|
1205 } |
|
1206 else |
|
1207 { |
|
1208 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1209 - tmp.abs (re_num, Round.NEAREST); |
|
1210 - var z = new Number.mpfloat (tmp); |
|
1211 + var z = copy (); |
|
1212 + if (z.re_sign < 0) |
|
1213 + z.re_sign = -z.re_sign; |
|
1214 return z; |
|
1215 } |
|
1216 } |
|
1217 @@ -311,7 +637,7 @@ |
|
1218 if (is_zero ()) |
|
1219 { |
|
1220 /* Translators: Error display when attempting to take argument of zero */ |
|
1221 - error = _("Argument not defined for zero"); |
|
1222 + mperr (_("Argument not defined for zero")); |
|
1223 return new Number.integer (0); |
|
1224 } |
|
1225 |
|
1226 @@ -355,12 +681,8 @@ |
|
1227 /* Sets z = ‾̅x */ |
|
1228 public Number conjugate () |
|
1229 { |
|
1230 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1231 - tmp.neg (im_num, Round.NEAREST); |
|
1232 var z = copy (); |
|
1233 - var tmp2 = z.im_num; |
|
1234 - tmp2.clear (); |
|
1235 - z.im_num = tmp; |
|
1236 + z.im_sign = -z.im_sign; |
|
1237 return z; |
|
1238 } |
|
1239 |
|
1240 @@ -368,11 +690,13 @@ |
|
1241 public Number real_component () |
|
1242 { |
|
1243 var z = copy (); |
|
1244 - var tmp = z.im_num; |
|
1245 - tmp.clear (); |
|
1246 - tmp = MPFloat.init2 ((Precision) precision); |
|
1247 - tmp.set_unsigned_integer (0, Round.NEAREST); |
|
1248 - z.im_num = tmp; |
|
1249 + |
|
1250 + /* Clear imaginary component */ |
|
1251 + z.im_sign = 0; |
|
1252 + z.im_exponent = 0; |
|
1253 + for (var i = 0; i < z.im_fraction.length; i++) |
|
1254 + z.im_fraction[i] = 0; |
|
1255 + |
|
1256 return z; |
|
1257 } |
|
1258 |
|
1259 @@ -380,30 +704,65 @@ |
|
1260 public Number imaginary_component () |
|
1261 { |
|
1262 /* Copy imaginary component to real component */ |
|
1263 - var z = copy (); |
|
1264 - var tmp = z.re_num; |
|
1265 - tmp.clear (); |
|
1266 - z.re_num = z.im_num; |
|
1267 - tmp = MPFloat.init2 ((Precision) precision); |
|
1268 - tmp.set_unsigned_integer (0, Round.NEAREST); |
|
1269 - z.im_num = tmp; |
|
1270 + var z = new Number (); |
|
1271 + z.re_sign = im_sign; |
|
1272 + z.re_exponent = im_exponent; |
|
1273 + |
|
1274 + for (var i = 0; i < z.im_fraction.length; i++) |
|
1275 + z.re_fraction[i] = im_fraction[i]; |
|
1276 + |
|
1277 + /* Clear (old) imaginary component */ |
|
1278 + z.im_sign = 0; |
|
1279 + z.im_exponent = 0; |
|
1280 + for (var i = 0; i < z.im_fraction.length; i++) |
|
1281 + z.im_fraction[i] = 0; |
|
1282 + |
|
1283 return z; |
|
1284 } |
|
1285 |
|
1286 public Number integer_component () |
|
1287 { |
|
1288 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1289 - tmp.trunc (re_num); |
|
1290 - var z = new Number.mpfloat (tmp); |
|
1291 + /* Clear re_fraction */ |
|
1292 + var z = copy (); |
|
1293 + for (var i = z.re_exponent; i < SIZE; i++) |
|
1294 + z.re_fraction[i] = 0; |
|
1295 + z.im_sign = 0; |
|
1296 + z.im_exponent = 0; |
|
1297 + for (var i = 0; i < z.im_fraction.length; i++) |
|
1298 + z.im_fraction[i] = 0; |
|
1299 + |
|
1300 return z; |
|
1301 + |
|
1302 } |
|
1303 |
|
1304 /* Sets z = x mod 1 */ |
|
1305 public Number fractional_component () |
|
1306 { |
|
1307 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1308 - tmp.frac (re_num, Round.NEAREST); |
|
1309 - var z = new Number.mpfloat (tmp); |
|
1310 + /* fractional component of zero is 0 */ |
|
1311 + if (is_zero ()) |
|
1312 + return new Number.integer (0); |
|
1313 + |
|
1314 + /* All fractional */ |
|
1315 + if (re_exponent <= 0) |
|
1316 + return this; |
|
1317 + |
|
1318 + /* Shift fractional component */ |
|
1319 + var shift = re_exponent; |
|
1320 + for (var i = shift; i < SIZE && re_fraction[i] == 0; i++) |
|
1321 + shift++; |
|
1322 + var z = new Number.integer (0); |
|
1323 + z.re_sign = re_sign; |
|
1324 + z.re_exponent = re_exponent - shift; |
|
1325 + for (var i = 0; i < SIZE; i++) |
|
1326 + { |
|
1327 + if (i + shift >= SIZE) |
|
1328 + z.re_fraction[i] = 0; |
|
1329 + else |
|
1330 + z.re_fraction[i] = re_fraction[i + shift]; |
|
1331 + } |
|
1332 + if (z.re_fraction[0] == 0) |
|
1333 + z.re_sign = 0; |
|
1334 + |
|
1335 return z; |
|
1336 } |
|
1337 |
|
1338 @@ -416,9 +775,36 @@ |
|
1339 /* Sets z = ⌊x⌋ */ |
|
1340 public Number floor () |
|
1341 { |
|
1342 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1343 - tmp.floor (re_num); |
|
1344 - var z = new Number.mpfloat (tmp); |
|
1345 + /* Integer component of zero = 0 */ |
|
1346 + if (is_zero ()) |
|
1347 + return this; |
|
1348 + |
|
1349 + /* If all fractional then no integer component */ |
|
1350 + if (re_exponent <= 0) |
|
1351 + { |
|
1352 + if (is_negative ()) |
|
1353 + return new Number.integer (-1); |
|
1354 + else |
|
1355 + return new Number.integer (0); |
|
1356 + } |
|
1357 + |
|
1358 + /* Clear fractional component */ |
|
1359 + var z = copy (); |
|
1360 + var have_fraction = false; |
|
1361 + for (var i = z.re_exponent; i < SIZE; i++) |
|
1362 + { |
|
1363 + if (z.re_fraction[i] != 0) |
|
1364 + have_fraction = true; |
|
1365 + z.re_fraction[i] = 0; |
|
1366 + } |
|
1367 + z.im_sign = 0; |
|
1368 + z.im_exponent = 0; |
|
1369 + for (var i = 0; i < z.im_fraction.length; i++) |
|
1370 + z.im_fraction[i] = 0; |
|
1371 + |
|
1372 + if (have_fraction && is_negative ()) |
|
1373 + z = z.add (new Number.integer (-1)); |
|
1374 + |
|
1375 return z; |
|
1376 } |
|
1377 |
|
1378 @@ -425,19 +811,29 @@ |
|
1379 /* Sets z = ⌈x⌉ */ |
|
1380 public Number ceiling () |
|
1381 { |
|
1382 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1383 - tmp.ceil (re_num); |
|
1384 - var z = new Number.mpfloat (tmp); |
|
1385 - return z; |
|
1386 + var z = floor (); |
|
1387 + var f = fractional_component (); |
|
1388 + if (f.is_zero ()) |
|
1389 + return z; |
|
1390 + return z.add (new Number.integer (1)); |
|
1391 } |
|
1392 |
|
1393 /* Sets z = [x] */ |
|
1394 public Number round () |
|
1395 { |
|
1396 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1397 - tmp.round (re_num); |
|
1398 - var z = new Number.mpfloat (tmp); |
|
1399 - return z; |
|
1400 + var do_floor = !is_negative (); |
|
1401 + |
|
1402 + var f = fractional_component (); |
|
1403 + f = f.multiply_integer (2); |
|
1404 + f = f.abs (); |
|
1405 + if (f.compare (new Number.integer (1)) >= 0) |
|
1406 + do_floor = !do_floor; |
|
1407 + |
|
1408 + if (do_floor) |
|
1409 + return floor (); |
|
1410 + else |
|
1411 + return ceiling (); |
|
1412 + |
|
1413 } |
|
1414 |
|
1415 /* Sets z = 1 ÷ x */ |
|
1416 @@ -482,24 +878,10 @@ |
|
1417 /* Sets z = x^y */ |
|
1418 public Number xpowy (Number y) |
|
1419 { |
|
1420 - /* 0^-n invalid */ |
|
1421 - if (is_zero () && y.is_negative ()) |
|
1422 + if (y.is_integer ()) |
|
1423 + return xpowy_integer (y.to_integer ()); |
|
1424 + else |
|
1425 { |
|
1426 - /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */ |
|
1427 - error = _("The power of zero is undefined for a negative exponent"); |
|
1428 - return new Number.integer (0); |
|
1429 - } |
|
1430 - |
|
1431 - /* 0^0 is indeterminate */ |
|
1432 - if (is_zero () && y.is_zero ()) |
|
1433 - { |
|
1434 - /* Translators: Error displayed when attempted to raise 0 to power of zero */ |
|
1435 - error = _("Zero raised to zero is undefined"); |
|
1436 - return new Number.integer (0); |
|
1437 - } |
|
1438 - |
|
1439 - if (!y.is_integer ()) |
|
1440 - { |
|
1441 var reciprocal = y.reciprocal (); |
|
1442 if (reciprocal.is_integer ()) |
|
1443 return root (reciprocal.to_integer ()); |
|
1444 @@ -506,29 +888,6 @@ |
|
1445 else |
|
1446 return pwr (y); |
|
1447 } |
|
1448 - |
|
1449 - Number t; |
|
1450 - Number t2; |
|
1451 - if (y.is_negative ()) |
|
1452 - { |
|
1453 - t = reciprocal (); |
|
1454 - t2 = y.invert_sign (); |
|
1455 - } |
|
1456 - else |
|
1457 - { |
|
1458 - t = this; |
|
1459 - t2 = y; |
|
1460 - } |
|
1461 - |
|
1462 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1463 - tmp.power (t.re_num, t2.re_num, Round.NEAREST); |
|
1464 - var z = new Number.mpfloat (tmp); |
|
1465 - var tmp2 = z.im_num; |
|
1466 - tmp2.clear (); |
|
1467 - tmp = MPFloat.init2 ((Precision) precision); |
|
1468 - tmp.@set (t.im_num, Round.NEAREST); |
|
1469 - z.im_num = tmp; |
|
1470 - return z; |
|
1471 } |
|
1472 |
|
1473 /* Sets z = x^y */ |
|
1474 @@ -538,7 +897,7 @@ |
|
1475 if (is_zero () && n < 0) |
|
1476 { |
|
1477 /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */ |
|
1478 - error = _("The power of zero is undefined for a negative exponent"); |
|
1479 + mperr (_("The power of zero is undefined for a negative exponent")); |
|
1480 return new Number.integer (0); |
|
1481 } |
|
1482 |
|
1483 @@ -546,10 +905,22 @@ |
|
1484 if (is_zero () && n == 0) |
|
1485 { |
|
1486 /* Translators: Error displayed when attempted to raise 0 to power of zero */ |
|
1487 - error = _("Zero raised to zero is undefined"); |
|
1488 + mperr (_("Zero raised to zero is undefined")); |
|
1489 return new Number.integer (0); |
|
1490 } |
|
1491 |
|
1492 + /* x^0 = 1 */ |
|
1493 + if (n == 0) |
|
1494 + return new Number.integer (1); |
|
1495 + |
|
1496 + /* 0^n = 0 */ |
|
1497 + if (is_zero ()) |
|
1498 + return new Number.integer (0); |
|
1499 + |
|
1500 + /* x^1 = x */ |
|
1501 + if (n == 1) |
|
1502 + return this; |
|
1503 + |
|
1504 Number t; |
|
1505 if (n < 0) |
|
1506 { |
|
1507 @@ -557,44 +928,17 @@ |
|
1508 n = -n; |
|
1509 } |
|
1510 else |
|
1511 - { |
|
1512 t = this; |
|
1513 - } |
|
1514 |
|
1515 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1516 - tmp.power_signed_integer (t.re_num, (long) n, Round.NEAREST); |
|
1517 - var z = new Number.mpfloat (tmp); |
|
1518 - var tmp2 = z.im_num; |
|
1519 - tmp2.clear (); |
|
1520 - tmp = MPFloat.init2 ((Precision) precision); |
|
1521 - tmp.@set (t.im_num, Round.NEAREST); |
|
1522 - z.im_num = tmp; |
|
1523 - return z; |
|
1524 - } |
|
1525 - |
|
1526 - private Number pwr (Number y) |
|
1527 - { |
|
1528 - /* (-x)^y imaginary */ |
|
1529 - /* FIXME: Make complex numbers optional */ |
|
1530 - /*if (re_sign < 0) |
|
1531 + var z = new Number.integer (1); |
|
1532 + while (n != 0) |
|
1533 { |
|
1534 - mperr (_("The power of negative numbers is only defined for integer exponents")); |
|
1535 - return new Number.integer (0); |
|
1536 - }*/ |
|
1537 - |
|
1538 - /* 0^y = 0, 0^-y undefined */ |
|
1539 - if (is_zero ()) |
|
1540 - { |
|
1541 - if (y.is_negative ()) |
|
1542 - error = _("The power of zero is undefined for a negative exponent"); |
|
1543 - return new Number.integer (0); |
|
1544 + if (n % 2 == 1) |
|
1545 + z = z.multiply (t); |
|
1546 + t = t.multiply (t); |
|
1547 + n = n / 2; |
|
1548 } |
|
1549 - |
|
1550 - /* x^0 = 1 */ |
|
1551 - if (y.is_zero ()) |
|
1552 - return new Number.integer (1); |
|
1553 - |
|
1554 - return y.multiply (ln ()).epowy (); |
|
1555 + return z; |
|
1556 } |
|
1557 |
|
1558 /* Sets z = n√x */ |
|
1559 @@ -623,10 +967,22 @@ |
|
1560 /* Sets z = √x */ |
|
1561 public Number sqrt () |
|
1562 { |
|
1563 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1564 - tmp.sqrt (re_num, Round.NEAREST); |
|
1565 - var z = new Number.mpfloat (tmp); |
|
1566 - return z; |
|
1567 + if (is_zero ()) |
|
1568 + return this; |
|
1569 + |
|
1570 + /* FIXME: Make complex numbers optional */ |
|
1571 + /*if (re_sign < 0) |
|
1572 + { |
|
1573 + mperr (_("Square root is undefined for negative values")); |
|
1574 + return new Number.integer (0); |
|
1575 + }*/ |
|
1576 + else |
|
1577 + { |
|
1578 + var t = root (-2); |
|
1579 + var z = multiply (t); |
|
1580 + return z.ext (t.re_fraction[0], z.re_fraction[0]); |
|
1581 + } |
|
1582 + |
|
1583 } |
|
1584 |
|
1585 /* Sets z = ln x */ |
|
1586 @@ -636,7 +992,7 @@ |
|
1587 if (is_zero ()) |
|
1588 { |
|
1589 /* Translators: Error displayed when attempting to take logarithm of zero */ |
|
1590 - error = _("Logarithm of zero is undefined"); |
|
1591 + mperr (_("Logarithm of zero is undefined")); |
|
1592 return new Number.integer (0); |
|
1593 } |
|
1594 |
|
1595 @@ -669,7 +1025,7 @@ |
|
1596 if (is_zero ()) |
|
1597 { |
|
1598 /* Translators: Error displayed when attempting to take logarithm of zero */ |
|
1599 - error = _("Logarithm of zero is undefined"); |
|
1600 + mperr (_("Logarithm of zero is undefined")); |
|
1601 return new Number.integer (0); |
|
1602 } |
|
1603 |
|
1604 @@ -691,17 +1047,17 @@ |
|
1605 if (is_negative () || is_complex ()) |
|
1606 { |
|
1607 /* Translators: Error displayed when attempted take the factorial of a negative or complex number */ |
|
1608 - error = _("Factorial is only defined for non-negative real numbers"); |
|
1609 + mperr (_("Factorial is only defined for non-negative real numbers")); |
|
1610 return new Number.integer (0); |
|
1611 } |
|
1612 |
|
1613 - var tmp = add (new Number.integer (1)); |
|
1614 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
1615 + var val = to_double (); |
|
1616 |
|
1617 /* Factorial(x) = Gamma(x+1) - This is the formula used to calculate Factorial.*/ |
|
1618 - tmp2.gamma (tmp.re_num, Round.NEAREST); |
|
1619 + var fact = Math.tgamma(val+1); |
|
1620 |
|
1621 - return new Number.mpfloat (tmp2); |
|
1622 + return new Number.double(fact); |
|
1623 + |
|
1624 } |
|
1625 |
|
1626 /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */ |
|
1627 @@ -716,41 +1072,23 @@ |
|
1628 /* Sets z = x + y */ |
|
1629 public Number add (Number y) |
|
1630 { |
|
1631 - if (is_complex () || y.is_complex ()) |
|
1632 - { |
|
1633 - Number real_z, im_z; |
|
1634 - |
|
1635 - var real_x = real_component (); |
|
1636 - var im_x = imaginary_component (); |
|
1637 - var real_y = y.real_component (); |
|
1638 - var im_y = y.imaginary_component (); |
|
1639 - |
|
1640 - real_z = real_x.add_real (real_y); |
|
1641 - im_z = im_x.add_real (im_y); |
|
1642 - |
|
1643 - return new Number.complex (real_z, im_z); |
|
1644 - } |
|
1645 - else |
|
1646 - return add_real (y); |
|
1647 + return add_with_sign (1, y); |
|
1648 } |
|
1649 |
|
1650 - public Number add_real (Number y) |
|
1651 - { |
|
1652 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1653 - tmp.add (re_num, y.re_num, Round.NEAREST); |
|
1654 - var z = new Number.mpfloat (tmp); |
|
1655 - return z; |
|
1656 - } |
|
1657 - |
|
1658 /* Sets z = x − y */ |
|
1659 public Number subtract (Number y) |
|
1660 { |
|
1661 - return add (y.invert_sign ()); |
|
1662 + return add_with_sign (-1, y); |
|
1663 } |
|
1664 |
|
1665 /* Sets z = x × y */ |
|
1666 public Number multiply (Number y) |
|
1667 { |
|
1668 + /* x*0 = 0*y = 0 */ |
|
1669 + if (is_zero () || y.is_zero ()) |
|
1670 + return new Number.integer (0); |
|
1671 + |
|
1672 + /* (a+bi)(c+di) = (ac-bd)+(ad+bc)i */ |
|
1673 if (is_complex () || y.is_complex ()) |
|
1674 { |
|
1675 Number t1, t2, real_z, im_z; |
|
1676 @@ -774,23 +1112,17 @@ |
|
1677 return multiply_real (y); |
|
1678 } |
|
1679 |
|
1680 - public Number multiply_real (Number y) |
|
1681 - { |
|
1682 - var z = new Number.integer (0); |
|
1683 - var tmp = z.re_num; |
|
1684 - tmp.multiply (re_num, y.re_num, Round.NEAREST); |
|
1685 - z.re_num = tmp; |
|
1686 - return z; |
|
1687 - } |
|
1688 - |
|
1689 /* Sets z = x × y */ |
|
1690 public Number multiply_integer (int64 y) |
|
1691 { |
|
1692 - var z = new Number.integer (0); |
|
1693 - var tmp = z.re_num; |
|
1694 - tmp.multiply_signed_integer (re_num, (long) y, Round.NEAREST); |
|
1695 - z.re_num = tmp; |
|
1696 - return z; |
|
1697 + if (is_complex ()) |
|
1698 + { |
|
1699 + var re_z = real_component ().multiply_integer_real (y);; |
|
1700 + var im_z = imaginary_component ().multiply_integer_real (y); |
|
1701 + return new Number.complex (re_z, im_z); |
|
1702 + } |
|
1703 + else |
|
1704 + return multiply_integer_real (y); |
|
1705 } |
|
1706 |
|
1707 /* Sets z = x ÷ y */ |
|
1708 @@ -799,31 +1131,23 @@ |
|
1709 if (y.is_zero ()) |
|
1710 { |
|
1711 /* Translators: Error displayed attempted to divide by zero */ |
|
1712 - error = _("Division by zero is undefined"); |
|
1713 + mperr (_("Division by zero is undefined")); |
|
1714 return new Number.integer (0); |
|
1715 } |
|
1716 + /* 0/y = 0 */ |
|
1717 + if (is_zero ()) |
|
1718 + return this; |
|
1719 |
|
1720 - if (is_complex () || y.is_complex ()) |
|
1721 - { |
|
1722 - var a = real_component (); |
|
1723 - var b = imaginary_component (); |
|
1724 - var c = y.real_component (); |
|
1725 - var d = y.imaginary_component (); |
|
1726 + /* z = x × y⁻¹ */ |
|
1727 + /* FIXME: Set re_exponent to zero to avoid overflow in multiply??? */ |
|
1728 + var t = y.reciprocal (); |
|
1729 + var ie = t.re_exponent; |
|
1730 + t.re_exponent = 0; |
|
1731 + var i = t.re_fraction[0]; |
|
1732 + var z = multiply (t); |
|
1733 + z = z.ext (i, z.re_fraction[0]); |
|
1734 + z.re_exponent += ie; |
|
1735 |
|
1736 - var tmp = a.multiply (c).add (b.multiply (d)); |
|
1737 - var tmp_2 = c.xpowy_integer (2).add (d.xpowy_integer (2)); |
|
1738 - var z_1 = tmp.divide (tmp_2); |
|
1739 - |
|
1740 - tmp = b.multiply (c).subtract (a.multiply (d)); |
|
1741 - var z_2 = tmp.divide (tmp_2); |
|
1742 - |
|
1743 - var z = new Number.complex (z_1, z_2); |
|
1744 - return z; |
|
1745 - } |
|
1746 - |
|
1747 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1748 - tmp.divide (re_num, y.re_num, Round.NEAREST); |
|
1749 - var z = new Number.mpfloat (tmp); |
|
1750 return z; |
|
1751 } |
|
1752 |
|
1753 @@ -830,7 +1154,14 @@ |
|
1754 /* Sets z = x ÷ y */ |
|
1755 public Number divide_integer (int64 y) |
|
1756 { |
|
1757 - return divide (new Number.integer (y)); |
|
1758 + if (is_complex ()) |
|
1759 + { |
|
1760 + var re_z = real_component ().divide_integer_real (y); |
|
1761 + var im_z = imaginary_component ().divide_integer_real (y); |
|
1762 + return new Number.complex (re_z, im_z); |
|
1763 + } |
|
1764 + else |
|
1765 + return divide_integer_real (y); |
|
1766 } |
|
1767 |
|
1768 /* Sets z = x mod y */ |
|
1769 @@ -839,7 +1170,7 @@ |
|
1770 if (!is_integer () || !y.is_integer ()) |
|
1771 { |
|
1772 /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */ |
|
1773 - error = _("Modulus division is only defined for integers"); |
|
1774 + mperr (_("Modulus division is only defined for integers")); |
|
1775 return new Number.integer (0); |
|
1776 } |
|
1777 |
|
1778 @@ -857,7 +1188,7 @@ |
|
1779 /* Sets z = x ^ y mod p */ |
|
1780 public Number modular_exponentiation (Number exp, Number mod) |
|
1781 { |
|
1782 - var base_value = copy (); |
|
1783 + var base_value = this.copy (); |
|
1784 if (exp.is_negative ()) |
|
1785 base_value = base_value.reciprocal (); |
|
1786 var exp_value = exp.abs (); |
|
1787 @@ -927,51 +1258,87 @@ |
|
1788 public Number tan (AngleUnit unit = AngleUnit.RADIANS) |
|
1789 { |
|
1790 /* Check for undefined values */ |
|
1791 - var x_radians = to_radians (unit); |
|
1792 - var check = x_radians.subtract (new Number.pi ().divide_integer (2)).divide (new Number.pi ()); |
|
1793 - |
|
1794 - if (check.is_integer ()) |
|
1795 + var cos_x = cos (unit); |
|
1796 + if (cos_x.is_zero ()) |
|
1797 { |
|
1798 /* Translators: Error displayed when tangent value is undefined */ |
|
1799 - error = _("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)"); |
|
1800 + mperr (_("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)")); |
|
1801 return new Number.integer (0); |
|
1802 } |
|
1803 |
|
1804 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1805 - tmp.tan (x_radians.re_num, Round.NEAREST); |
|
1806 - var z = new Number.mpfloat (tmp); |
|
1807 - return z; |
|
1808 + /* tan (x) = sin (x) / cos (x) */ |
|
1809 + return sin (unit).divide (cos_x); |
|
1810 } |
|
1811 |
|
1812 /* Sets z = sin⁻¹ x */ |
|
1813 public Number asin (AngleUnit unit = AngleUnit.RADIANS) |
|
1814 { |
|
1815 - if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0) |
|
1816 - { |
|
1817 - /* Translators: Error displayed when inverse sine value is undefined */ |
|
1818 - error = _("Inverse sine is undefined for values outside [-1, 1]"); |
|
1819 + /* asin⁻¹(0) = 0 */ |
|
1820 + if (is_zero ()) |
|
1821 return new Number.integer (0); |
|
1822 + |
|
1823 + /* sin⁻¹(x) = tan⁻¹(x / √(1 - x²)), |x| < 1 */ |
|
1824 + if (re_exponent <= 0) |
|
1825 + { |
|
1826 + var t1 = new Number.integer (1); |
|
1827 + var t2 = t1; |
|
1828 + t1 = t1.subtract (this); |
|
1829 + t2 = t2.add (this); |
|
1830 + t2 = t1.multiply (t2); |
|
1831 + t2 = t2.root (-2); |
|
1832 + var z = multiply (t2); |
|
1833 + z = z.atan (unit); |
|
1834 + return z; |
|
1835 } |
|
1836 |
|
1837 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1838 - tmp.asin (re_num, Round.NEAREST); |
|
1839 - var z = new Number.mpfloat (tmp); |
|
1840 - return z.from_radians (unit); |
|
1841 + /* sin⁻¹(1) = π/2, sin⁻¹(-1) = -π/2 */ |
|
1842 + var t2 = new Number.integer (re_sign); |
|
1843 + if (equals (t2)) |
|
1844 + { |
|
1845 + var z = new Number.pi ().divide_integer (2 * t2.re_sign); |
|
1846 + return z.from_radians (unit); |
|
1847 + } |
|
1848 + |
|
1849 + /* Translators: Error displayed when inverse sine value is undefined */ |
|
1850 + mperr (_("Inverse sine is undefined for values outside [-1, 1]")); |
|
1851 + return new Number.integer (0); |
|
1852 } |
|
1853 |
|
1854 /* Sets z = cos⁻¹ x */ |
|
1855 public Number acos (AngleUnit unit = AngleUnit.RADIANS) |
|
1856 { |
|
1857 - if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0) |
|
1858 + var pi = new Number.pi (); |
|
1859 + var t1 = new Number.integer (1); |
|
1860 + var n1 = new Number.integer (-1); |
|
1861 + |
|
1862 + Number z; |
|
1863 + if (compare (t1) > 0 || compare (n1) < 0) |
|
1864 { |
|
1865 /* Translators: Error displayed when inverse cosine value is undefined */ |
|
1866 - error = _("Inverse cosine is undefined for values outside [-1, 1]"); |
|
1867 - return new Number.integer (0); |
|
1868 + mperr (_("Inverse cosine is undefined for values outside [-1, 1]")); |
|
1869 + z = new Number.integer (0); |
|
1870 } |
|
1871 + else if (is_zero ()) |
|
1872 + z = pi.divide_integer (2); |
|
1873 + else if (equals (t1)) |
|
1874 + z = new Number.integer (0); |
|
1875 + else if (equals (n1)) |
|
1876 + z = pi; |
|
1877 + else |
|
1878 + { |
|
1879 + /* cos⁻¹(x) = tan⁻¹(√(1 - x²) / x) */ |
|
1880 + Number y; |
|
1881 + var t2 = multiply (this); |
|
1882 + t2 = t1.subtract (t2); |
|
1883 + t2 = t2.sqrt (); |
|
1884 + t2 = t2.divide (this); |
|
1885 + y = t2.atan (AngleUnit.RADIANS); |
|
1886 + if (re_sign > 0) |
|
1887 + z = y; |
|
1888 + else |
|
1889 + z = y.add (pi); |
|
1890 + } |
|
1891 |
|
1892 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1893 - tmp.acos (re_num, Round.NEAREST); |
|
1894 - var z = new Number.mpfloat (tmp); |
|
1895 return z.from_radians (unit); |
|
1896 } |
|
1897 |
|
1898 @@ -978,9 +1345,63 @@ |
|
1899 /* Sets z = tan⁻¹ x */ |
|
1900 public Number atan (AngleUnit unit = AngleUnit.RADIANS) |
|
1901 { |
|
1902 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1903 - tmp.atan (re_num, Round.NEAREST); |
|
1904 - var z = new Number.mpfloat (tmp); |
|
1905 + if (is_zero ()) |
|
1906 + return new Number.integer (0); |
|
1907 + |
|
1908 + var t2 = this; |
|
1909 + var rx = 0f; |
|
1910 + if (re_exponent.abs () <= 2) |
|
1911 + rx = to_float (); |
|
1912 + |
|
1913 + /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */ |
|
1914 + var q = 1; |
|
1915 + Number z; |
|
1916 + while (t2.re_exponent >= 0) |
|
1917 + { |
|
1918 + if (t2.re_exponent == 0 && 2 * (t2.re_fraction[0] + 1) <= BASE) |
|
1919 + break; |
|
1920 + |
|
1921 + q *= 2; |
|
1922 + |
|
1923 + /* t = t / (√(t² + 1) + 1) */ |
|
1924 + z = t2.multiply (t2); |
|
1925 + z = z.add (new Number.integer (1)); |
|
1926 + z = z.sqrt (); |
|
1927 + z = z.add (new Number.integer (1)); |
|
1928 + t2 = t2.divide (z); |
|
1929 + } |
|
1930 + |
|
1931 + /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */ |
|
1932 + z = t2; |
|
1933 + var t1 = t2.multiply (t2); |
|
1934 + |
|
1935 + /* SERIES LOOP. REDUCE T IF POSSIBLE. */ |
|
1936 + for (var i = 1; ; i += 2) |
|
1937 + { |
|
1938 + if (T + 2 + t2.re_exponent <= 1) |
|
1939 + break; |
|
1940 + |
|
1941 + t2 = t2.multiply (t1).multiply_integer (-i).divide_integer (i + 2); |
|
1942 + |
|
1943 + z = z.add (t2); |
|
1944 + if (t2.is_zero ()) |
|
1945 + break; |
|
1946 + } |
|
1947 + |
|
1948 + /* CORRECT FOR ARGUMENT REDUCTION */ |
|
1949 + z = z.multiply_integer (q); |
|
1950 + |
|
1951 + /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS re_exponent |
|
1952 + * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK) |
|
1953 + */ |
|
1954 + if (re_exponent.abs () <= 2) |
|
1955 + { |
|
1956 + float ry = z.to_float (); |
|
1957 + /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */ |
|
1958 + if (Math.fabs (ry - Math.atan (rx)) >= Math.fabs (ry) * 0.01) |
|
1959 + mperr ("*** ERROR OCCURRED IN ATAN, RESULT INCORRECT ***"); |
|
1960 + } |
|
1961 + |
|
1962 return z.from_radians (unit); |
|
1963 } |
|
1964 |
|
1965 @@ -987,27 +1408,88 @@ |
|
1966 /* Sets z = sinh x */ |
|
1967 public Number sinh () |
|
1968 { |
|
1969 - var tmp = MPFloat.init2 ((Precision) precision); |
|
1970 - tmp.sinh (re_num, Round.NEAREST); |
|
1971 - var z = new Number.mpfloat (tmp); |
|
1972 - return z; |
|
1973 + /* sinh (0) = 0 */ |
|
1974 + if (is_zero ()) |
|
1975 + return new Number.integer (0); |
|
1976 + |
|
1977 + /* WORK WITH ABS (X) */ |
|
1978 + var abs_x = abs (); |
|
1979 + |
|
1980 + /* If |x| < 1 USE EXP TO AVOID CANCELLATION, otherwise IF TOO LARGE EPOWY GIVES ERROR MESSAGE */ |
|
1981 + Number z; |
|
1982 + if (abs_x.re_exponent <= 0) |
|
1983 + { |
|
1984 + /* ((e^|x| + 1) * (e^|x| - 1)) / e^|x| */ |
|
1985 + // FIXME: Solves to e^|x| - e^-|x|, why not lower branch always? */ |
|
1986 + var exp_x = abs_x.epowy (); |
|
1987 + var a = exp_x.add (new Number.integer (1)); |
|
1988 + var b = exp_x.add (new Number.integer (-1)); |
|
1989 + z = a.multiply (b); |
|
1990 + z = z.divide (exp_x); |
|
1991 + } |
|
1992 + else |
|
1993 + { |
|
1994 + /* e^|x| - e^-|x| */ |
|
1995 + var exp_x = abs_x.epowy (); |
|
1996 + z = exp_x.reciprocal (); |
|
1997 + z = exp_x.subtract (z); |
|
1998 + } |
|
1999 + |
|
2000 + /* DIVIDE BY TWO AND RESTORE re_sign */ |
|
2001 + z = z.divide_integer (2); |
|
2002 + return z.multiply_integer (re_sign); |
|
2003 } |
|
2004 |
|
2005 /* Sets z = cosh x */ |
|
2006 public Number cosh () |
|
2007 { |
|
2008 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2009 - tmp.cosh (re_num, Round.NEAREST); |
|
2010 - var z = new Number.mpfloat (tmp); |
|
2011 - return z; |
|
2012 + /* cosh (0) = 1 */ |
|
2013 + if (is_zero ()) |
|
2014 + return new Number.integer (1); |
|
2015 + |
|
2016 + /* cosh (x) = (e^x + e^-x) / 2 */ |
|
2017 + var t = abs (); |
|
2018 + t = t.epowy (); |
|
2019 + var z = t.reciprocal (); |
|
2020 + z = t.add (z); |
|
2021 + return z.divide_integer (2); |
|
2022 } |
|
2023 |
|
2024 /* Sets z = tanh x */ |
|
2025 public Number tanh () |
|
2026 { |
|
2027 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2028 - tmp.tanh (re_num, Round.NEAREST); |
|
2029 - var z = new Number.mpfloat (tmp); |
|
2030 + /* tanh (0) = 0 */ |
|
2031 + if (is_zero ()) |
|
2032 + return new Number.integer (0); |
|
2033 + |
|
2034 + var t = abs (); |
|
2035 + |
|
2036 + /* SEE IF ABS (X) SO LARGE THAT RESULT IS +-1 */ |
|
2037 + var r__1 = (float) T * 0.5 * Math.log ((float) BASE); |
|
2038 + var z = new Number.double (r__1); |
|
2039 + if (t.compare (z) > 0) |
|
2040 + return new Number.integer (re_sign); |
|
2041 + |
|
2042 + /* If |x| >= 1/2 use ?, otherwise use ? to avoid cancellation */ |
|
2043 + /* |tanh (x)| = (e^|2x| - 1) / (e^|2x| + 1) */ |
|
2044 + t = t.multiply_integer (2); |
|
2045 + if (t.re_exponent > 0) |
|
2046 + { |
|
2047 + t = t.epowy (); |
|
2048 + z = t.add (new Number.integer (-1)); |
|
2049 + t = t.add (new Number.integer (1)); |
|
2050 + z = z.divide (t); |
|
2051 + } |
|
2052 + else |
|
2053 + { |
|
2054 + t = t.epowy (); |
|
2055 + z = t.add (new Number.integer (1)); |
|
2056 + t = t.add (new Number.integer (-1)); |
|
2057 + z = t.divide (z); |
|
2058 + } |
|
2059 + |
|
2060 + /* Restore re_sign */ |
|
2061 + z.re_sign = re_sign * z.re_sign; |
|
2062 return z; |
|
2063 } |
|
2064 |
|
2065 @@ -1014,10 +1496,12 @@ |
|
2066 /* Sets z = sinh⁻¹ x */ |
|
2067 public Number asinh () |
|
2068 { |
|
2069 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2070 - tmp.asinh (re_num, Round.NEAREST); |
|
2071 - var z = new Number.mpfloat (tmp); |
|
2072 - return z; |
|
2073 + /* sinh⁻¹(x) = ln (x + √(x² + 1)) */ |
|
2074 + var t = multiply (this); |
|
2075 + t = t.add (new Number.integer (1)); |
|
2076 + t = t.sqrt (); |
|
2077 + t = add (t); |
|
2078 + return t.ln (); |
|
2079 } |
|
2080 |
|
2081 /* Sets z = cosh⁻¹ x */ |
|
2082 @@ -1028,14 +1512,16 @@ |
|
2083 if (compare (t) < 0) |
|
2084 { |
|
2085 /* Translators: Error displayed when inverse hyperbolic cosine value is undefined */ |
|
2086 - error = _("Inverse hyperbolic cosine is undefined for values less than one"); |
|
2087 + mperr (_("Inverse hyperbolic cosine is undefined for values less than one")); |
|
2088 return new Number.integer (0); |
|
2089 } |
|
2090 |
|
2091 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2092 - tmp.acosh (re_num, Round.NEAREST); |
|
2093 - var z = new Number.mpfloat (tmp); |
|
2094 - return z; |
|
2095 + /* cosh⁻¹(x) = ln (x + √(x² - 1)) */ |
|
2096 + t = multiply (this); |
|
2097 + t = t.add (new Number.integer (-1)); |
|
2098 + t = t.sqrt (); |
|
2099 + t = add (t); |
|
2100 + return t.ln (); |
|
2101 } |
|
2102 |
|
2103 /* Sets z = tanh⁻¹ x */ |
|
2104 @@ -1045,14 +1531,17 @@ |
|
2105 if (compare (new Number.integer (1)) >= 0 || compare (new Number.integer (-1)) <= 0) |
|
2106 { |
|
2107 /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */ |
|
2108 - error = _("Inverse hyperbolic tangent is undefined for values outside [-1, 1]"); |
|
2109 + mperr (_("Inverse hyperbolic tangent is undefined for values outside [-1, 1]")); |
|
2110 return new Number.integer (0); |
|
2111 } |
|
2112 |
|
2113 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2114 - tmp.atanh (re_num, Round.NEAREST); |
|
2115 - var z = new Number.mpfloat (tmp); |
|
2116 - return z; |
|
2117 + /* atanh (x) = 0.5 * ln ((1 + x) / (1 - x)) */ |
|
2118 + var n = add (new Number.integer (1)); |
|
2119 + var d = invert_sign (); |
|
2120 + d = d.add (new Number.integer (1)); |
|
2121 + var z = n.divide (d); |
|
2122 + z = z.ln (); |
|
2123 + return z.divide_integer (2); |
|
2124 } |
|
2125 |
|
2126 /* Sets z = boolean AND for each bit in x and z */ |
|
2127 @@ -1062,7 +1551,7 @@ |
|
2128 is_positive_integer () || !y.is_positive_integer ()) |
|
2129 { |
|
2130 /* Translators: Error displayed when boolean AND attempted on non-integer values */ |
|
2131 - error = _("Boolean AND is only defined for positive integers"); |
|
2132 + mperr (_("Boolean AND is only defined for positive integers")); |
|
2133 } |
|
2134 |
|
2135 return bitwise (y, (v1, v2) => { return v1 & v2; }, 0); |
|
2136 @@ -1074,7 +1563,7 @@ |
|
2137 if (!is_positive_integer () || !y.is_positive_integer ()) |
|
2138 { |
|
2139 /* Translators: Error displayed when boolean OR attempted on non-integer values */ |
|
2140 - error = _("Boolean OR is only defined for positive integers"); |
|
2141 + mperr (_("Boolean OR is only defined for positive integers")); |
|
2142 } |
|
2143 |
|
2144 return bitwise (y, (v1, v2) => { return v1 | v2; }, 0); |
|
2145 @@ -1086,7 +1575,7 @@ |
|
2146 if (!is_positive_integer () || !y.is_positive_integer ()) |
|
2147 { |
|
2148 /* Translators: Error displayed when boolean XOR attempted on non-integer values */ |
|
2149 - error = _("Boolean XOR is only defined for positive integers"); |
|
2150 + mperr (_("Boolean XOR is only defined for positive integers")); |
|
2151 } |
|
2152 |
|
2153 return bitwise (y, (v1, v2) => { return v1 ^ v2; }, 0); |
|
2154 @@ -1098,7 +1587,7 @@ |
|
2155 if (!is_positive_integer ()) |
|
2156 { |
|
2157 /* Translators: Error displayed when boolean XOR attempted on non-integer values */ |
|
2158 - error = _("Boolean NOT is only defined for positive integers"); |
|
2159 + mperr (_("Boolean NOT is only defined for positive integers")); |
|
2160 } |
|
2161 |
|
2162 return bitwise (new Number.integer (0), (v1, v2) => { return v1 ^ 0xF; }, wordlen); |
|
2163 @@ -1121,7 +1610,7 @@ |
|
2164 if (!is_integer ()) |
|
2165 { |
|
2166 /* Translators: Error displayed when bit shift attempted on non-integer values */ |
|
2167 - error = _("Shift is only possible on integer values"); |
|
2168 + mperr (_("Shift is only possible on integer values")); |
|
2169 return new Number.integer (0); |
|
2170 } |
|
2171 |
|
2172 @@ -1253,60 +1742,1056 @@ |
|
2173 private Number copy () |
|
2174 { |
|
2175 var z = new Number (); |
|
2176 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2177 - var tmp2 = MPFloat.init2 ((Precision) precision); |
|
2178 - tmp.@set(re_num, Round.NEAREST); |
|
2179 - tmp2.@set(im_num, Round.NEAREST); |
|
2180 - z.re_num = tmp; |
|
2181 - z.im_num = tmp2; |
|
2182 + z.re_sign = re_sign; |
|
2183 + z.im_sign = im_sign; |
|
2184 + z.re_exponent = re_exponent; |
|
2185 + z.im_exponent = im_exponent; |
|
2186 + for (var i = 0; i < re_fraction.length; i++) |
|
2187 + { |
|
2188 + z.re_fraction[i] = re_fraction[i]; |
|
2189 + z.im_fraction[i] = im_fraction[i]; |
|
2190 + } |
|
2191 + |
|
2192 return z; |
|
2193 } |
|
2194 |
|
2195 + private Number add_with_sign (int y_sign, Number y) |
|
2196 + { |
|
2197 + if (is_complex () || y.is_complex ()) |
|
2198 + { |
|
2199 + Number real_z, im_z; |
|
2200 + |
|
2201 + var real_x = real_component (); |
|
2202 + var im_x = imaginary_component (); |
|
2203 + var real_y = y.real_component (); |
|
2204 + var im_y = y.imaginary_component (); |
|
2205 + |
|
2206 + real_z = real_x.add_real (y_sign * y.re_sign, real_y); |
|
2207 + im_z = im_x.add_real (y_sign * y.im_sign, im_y); |
|
2208 + |
|
2209 + return new Number.complex (real_z, im_z); |
|
2210 + } |
|
2211 + else |
|
2212 + return add_real (y_sign * y.re_sign, y); |
|
2213 + } |
|
2214 + |
|
2215 + |
|
2216 private Number epowy_real () |
|
2217 { |
|
2218 - var z = copy (); |
|
2219 - var tmp = z.re_num; |
|
2220 - tmp.exp (re_num, Round.NEAREST); |
|
2221 - z.re_num = tmp; |
|
2222 + /* e^0 = 1 */ |
|
2223 + if (is_zero ()) |
|
2224 + return new Number.integer (1); |
|
2225 + |
|
2226 + /* If |x| < 1 use exp */ |
|
2227 + if (re_exponent <= 0) |
|
2228 + return exp (); |
|
2229 + |
|
2230 + /* NOW SAFE TO CONVERT X TO REAL */ |
|
2231 + var rx = to_double (); |
|
2232 + |
|
2233 + /* SAVE re_sign AND WORK WITH ABS (X) */ |
|
2234 + var xs = re_sign; |
|
2235 + var t2 = abs (); |
|
2236 + |
|
2237 + /* GET fractional AND INTEGER PARTS OF ABS (X) */ |
|
2238 + var ix = t2.to_integer (); |
|
2239 + t2 = t2.fractional_component (); |
|
2240 + |
|
2241 + /* ATTACH re_sign TO fractional PART AND COMPUTE EXP OF IT */ |
|
2242 + t2.re_sign *= xs; |
|
2243 + var z = t2.exp (); |
|
2244 + |
|
2245 + /* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS (X) LARGE |
|
2246 + * (BUT ONLY ONE EXTRA DIGIT IF T < 4) |
|
2247 + */ |
|
2248 + var tss = 0; |
|
2249 + if (T < 4) |
|
2250 + tss = T + 1; |
|
2251 + else |
|
2252 + tss = T + 2; |
|
2253 + |
|
2254 + /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */ |
|
2255 + /* Computing MIN */ |
|
2256 + var t1 = new Number.integer (xs); |
|
2257 + |
|
2258 + t2.re_sign = 0; |
|
2259 + for (var i = 2 ; ; i++) |
|
2260 + { |
|
2261 + if (int.min (tss, tss + 2 + t1.re_exponent) <= 2) |
|
2262 + break; |
|
2263 + |
|
2264 + t1 = t1.divide_integer (i * xs); |
|
2265 + t2 = t2.add (t1); |
|
2266 + if (t1.is_zero ()) |
|
2267 + break; |
|
2268 + } |
|
2269 + |
|
2270 + /* RAISE E OR 1/E TO POWER IX */ |
|
2271 + if (xs > 0) |
|
2272 + t2 = t2.add (new Number.integer (2)); |
|
2273 + t2 = t2.xpowy_integer (ix); |
|
2274 + |
|
2275 + /* MULTIPLY EXPS OF INTEGER AND fractional PARTS */ |
|
2276 + z = z.multiply (t2); |
|
2277 + |
|
2278 + /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS (X) LARGE |
|
2279 + * (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW) |
|
2280 + */ |
|
2281 + if (Math.fabs (rx) > 10.0f) |
|
2282 + return z; |
|
2283 + var rz = z.to_double (); |
|
2284 + var r__1 = rz - Math.exp (rx); |
|
2285 + if (Math.fabs (r__1) < rz * 0.01f) |
|
2286 + return z; |
|
2287 + |
|
2288 + /* THE FOLLOWING MESSAGE MAY INDICATE THAT |
|
2289 + * B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE |
|
2290 + * RESULT UNDERFLOWED. |
|
2291 + */ |
|
2292 + mperr ("*** ERROR OCCURRED IN EPOWY, RESULT INCORRECT ***"); |
|
2293 return z; |
|
2294 + |
|
2295 } |
|
2296 |
|
2297 + /* Return e^x for |x| < 1 USING AN o (SQRt (T).m (T)) ALGORITHM |
|
2298 + * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE- |
|
2299 + * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM |
|
2300 + * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165). |
|
2301 + * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL |
|
2302 + * UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL. |
|
2303 + */ |
|
2304 + private Number exp () |
|
2305 + { |
|
2306 + /* e^0 = 1 */ |
|
2307 + if (is_zero ()) |
|
2308 + return new Number.integer (1); |
|
2309 + |
|
2310 + /* Only defined for |x| < 1 */ |
|
2311 + if (re_exponent > 0) |
|
2312 + { |
|
2313 + mperr ("*** ABS (X) NOT LESS THAN 1 IN CALL TO MP_EXP ***"); |
|
2314 + return new Number.integer (0); |
|
2315 + } |
|
2316 + |
|
2317 + var t1 = this; |
|
2318 + var rlb = Math.log (BASE); |
|
2319 + |
|
2320 + /* Compute approximately optimal q (and divide x by 2^q) */ |
|
2321 + var q = (int) (Math.sqrt (T * 0.48 * rlb) + re_exponent * 1.44 * rlb); |
|
2322 + |
|
2323 + /* HALVE Q TIMES */ |
|
2324 + if (q > 0) |
|
2325 + { |
|
2326 + var ib = BASE << 2; |
|
2327 + var ic = 1; |
|
2328 + for (var i = 1; i <= q; i++) |
|
2329 + { |
|
2330 + ic *= 2; |
|
2331 + if (ic < ib && ic != BASE && i < q) |
|
2332 + continue; |
|
2333 + t1 = t1.divide_integer (ic); |
|
2334 + ic = 1; |
|
2335 + } |
|
2336 + } |
|
2337 + |
|
2338 + if (t1.is_zero ()) |
|
2339 + return new Number.integer (0); |
|
2340 + |
|
2341 + /* Sum series, reducing t where possible */ |
|
2342 + var z = t1.copy (); |
|
2343 + var t2 = t1; |
|
2344 + for (var i = 2; T + t2.re_exponent - z.re_exponent > 0; i++) |
|
2345 + { |
|
2346 + t2 = t1.multiply (t2); |
|
2347 + t2 = t2.divide_integer (i); |
|
2348 + z = t2.add (z); |
|
2349 + if (t2.is_zero ()) |
|
2350 + break; |
|
2351 + } |
|
2352 + |
|
2353 + /* Apply (x+1)^2 - 1 = x (2 + x) for q iterations */ |
|
2354 + for (var i = 1; i <= q; i++) |
|
2355 + { |
|
2356 + t1 = z.add (new Number.integer (2)); |
|
2357 + z = t1.multiply (z); |
|
2358 + } |
|
2359 + |
|
2360 + return z.add (new Number.integer (1)); |
|
2361 + } |
|
2362 + |
|
2363 + private Number pwr (Number y) |
|
2364 + { |
|
2365 + /* (-x)^y imaginary */ |
|
2366 + /* FIXME: Make complex numbers optional */ |
|
2367 + /*if (re_sign < 0) |
|
2368 + { |
|
2369 + mperr (_("The power of negative numbers is only defined for integer exponents")); |
|
2370 + return new Number.integer (0); |
|
2371 + }*/ |
|
2372 + |
|
2373 + /* 0^y = 0, 0^-y undefined */ |
|
2374 + if (is_zero ()) |
|
2375 + { |
|
2376 + if (y.re_sign < 0) |
|
2377 + mperr (_("The power of zero is undefined for a negative exponent")); |
|
2378 + return new Number.integer (0); |
|
2379 + } |
|
2380 + |
|
2381 + /* x^0 = 1 */ |
|
2382 + if (y.is_zero ()) |
|
2383 + return new Number.integer (1); |
|
2384 + |
|
2385 + return y.multiply (ln ()).epowy (); |
|
2386 + } |
|
2387 + |
|
2388 private Number root_real (int64 n) |
|
2389 { |
|
2390 + /* x^(1/1) = x */ |
|
2391 + if (n == 1) |
|
2392 + return this; |
|
2393 + |
|
2394 + /* x^(1/0) invalid */ |
|
2395 + if (n == 0) |
|
2396 + { |
|
2397 + mperr (_("Root must be non-zero")); |
|
2398 + return new Number.integer (0); |
|
2399 + } |
|
2400 + |
|
2401 + var np = n.abs (); |
|
2402 + |
|
2403 + /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */ |
|
2404 + if (np > int.max (BASE, 64)) |
|
2405 + { |
|
2406 + mperr ("*** ABS (N) TOO LARGE IN CALL TO ROOT ***"); |
|
2407 + return new Number.integer (0); |
|
2408 + } |
|
2409 + |
|
2410 + /* 0^(1/n) = 0 for positive n */ |
|
2411 + if (is_zero ()) |
|
2412 + { |
|
2413 + if (n <= 0) |
|
2414 + mperr (_("Negative root of zero is undefined")); |
|
2415 + return new Number.integer (0); |
|
2416 + } |
|
2417 + |
|
2418 + // FIXME: Imaginary root |
|
2419 + if (re_sign < 0 && np % 2 == 0) |
|
2420 + { |
|
2421 + mperr (_("nth root of negative number is undefined for even n")); |
|
2422 + return new Number.integer (0); |
|
2423 + } |
|
2424 + |
|
2425 + /* DIVIDE re_exponent BY NP */ |
|
2426 + var ex = re_exponent / np; |
|
2427 + |
|
2428 + /* Get initial approximation */ |
|
2429 + var t1 = copy (); |
|
2430 + t1.re_exponent = 0; |
|
2431 + var approximation = Math.exp (((np * ex - re_exponent) * Math.log (BASE) - Math.log (Math.fabs (t1.to_float ()))) / np); |
|
2432 + t1 = new Number.double (approximation); |
|
2433 + t1.re_sign = re_sign; |
|
2434 + t1.re_exponent -= (int) ex; |
|
2435 + |
|
2436 + /* MAIN ITERATION LOOP */ |
|
2437 + var it0 = 3; |
|
2438 + var t = it0; |
|
2439 + Number t2; |
|
2440 + while (true) |
|
2441 + { |
|
2442 + /* t1 = t1 - ((t1 * ((x * t1^np) - 1)) / np) */ |
|
2443 + t2 = t1.xpowy_integer (np); |
|
2444 + t2 = multiply (t2); |
|
2445 + t2 = t2.add (new Number.integer (-1)); |
|
2446 + t2 = t1.multiply (t2); |
|
2447 + t2 = t2.divide_integer (np); |
|
2448 + t1 = t1.subtract (t2); |
|
2449 + |
|
2450 + /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE |
|
2451 + * NEWTONS METHOD HAS 2ND ORDER CONVERGENCE). |
|
2452 + */ |
|
2453 + if (t >= T) |
|
2454 + break; |
|
2455 + |
|
2456 + var ts3 = t; |
|
2457 + var ts2 = t; |
|
2458 + t = T; |
|
2459 + do |
|
2460 + { |
|
2461 + ts2 = t; |
|
2462 + t = (t + it0) / 2; |
|
2463 + } while (t > ts3); |
|
2464 + t = int.min (ts2, T); |
|
2465 + } |
|
2466 + |
|
2467 + /* NOW r (I2) IS X**(-1/NP) |
|
2468 + * CHECK THAT NEWTON ITERATION WAS CONVERGING |
|
2469 + */ |
|
2470 + if (t2.re_sign != 0 && (t1.re_exponent - t2.re_exponent) << 1 < T - it0) |
|
2471 + { |
|
2472 + /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL, |
|
2473 + * OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP |
|
2474 + * IS NOT ACCURATE ENOUGH. |
|
2475 + */ |
|
2476 + mperr ("*** ERROR OCCURRED IN ROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***"); |
|
2477 + } |
|
2478 + |
|
2479 + if (n >= 0) |
|
2480 + { |
|
2481 + t1 = t1.xpowy_integer (n - 1); |
|
2482 + return multiply (t1); |
|
2483 + } |
|
2484 + |
|
2485 + return t1; |
|
2486 + |
|
2487 + } |
|
2488 + |
|
2489 + /* ROUTINE CALLED BY DIVIDE AND SQRT TO ENSURE THAT |
|
2490 + * RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY |
|
2491 + * CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS. |
|
2492 + */ |
|
2493 + private Number ext (int i, int j) |
|
2494 + { |
|
2495 + if (is_zero () || T <= 2 || i == 0) |
|
2496 + return this; |
|
2497 + |
|
2498 + /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */ |
|
2499 + var q = (j + 1) / i + 1; |
|
2500 + var s = BASE * re_fraction[T - 2] + re_fraction[T - 1]; |
|
2501 + |
|
2502 + /* SET LAST TWO DIGITS TO ZERO */ |
|
2503 var z = copy (); |
|
2504 - var tmp = z.re_num; |
|
2505 - tmp.root (re_num, (ulong) n, Round.NEAREST); |
|
2506 - z.re_num = tmp; |
|
2507 - return z; |
|
2508 + if (s <= q) |
|
2509 + { |
|
2510 + z.re_fraction[T - 2] = 0; |
|
2511 + z.re_fraction[T - 1] = 0; |
|
2512 + return z; |
|
2513 + } |
|
2514 + |
|
2515 + if (s + q < BASE * BASE) |
|
2516 + return z; |
|
2517 + |
|
2518 + /* ROUND UP HERE */ |
|
2519 + z.re_fraction[T - 2] = BASE - 1; |
|
2520 + z.re_fraction[T - 1] = BASE; |
|
2521 + |
|
2522 + /* NORMALIZE X (LAST DIGIT B IS OK IN MULTIPLY_INTEGER) */ |
|
2523 + return z.multiply_integer (1); |
|
2524 } |
|
2525 |
|
2526 + |
|
2527 private Number ln_real () |
|
2528 { |
|
2529 - var tmp = MPFloat.init2 ((Precision) precision); |
|
2530 - tmp.log (re_num, Round.NEAREST); |
|
2531 - var z = new Number.mpfloat (tmp); |
|
2532 + // ln(e^1) = 1, fixes precision loss |
|
2533 + if (equals (new Number.eulers ())) |
|
2534 + return new Number.integer (1); |
|
2535 + |
|
2536 + /* LOOP TO GET APPROXIMATE Ln (X) USING SINGLE-PRECISION */ |
|
2537 + var t1 = copy (); |
|
2538 + var z = new Number.integer (0); |
|
2539 + for (var k = 0; k < 10; k++) |
|
2540 + { |
|
2541 + /* COMPUTE FINAL CORRECTION ACCURATELY USING LNS */ |
|
2542 + var t2 = t1.add (new Number.integer (-1)); |
|
2543 + if (t2.is_zero () || t2.re_exponent + 1 <= 0) |
|
2544 + return z.add (t2.lns ()); |
|
2545 + |
|
2546 + /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */ |
|
2547 + var e = t1.re_exponent; |
|
2548 + t1.re_exponent = 0; |
|
2549 + var rx = t1.to_double (); |
|
2550 + t1.re_exponent = e; |
|
2551 + var rlx = Math.log (rx) + e * Math.log (BASE); |
|
2552 + t2 = new Number.double (-rlx); |
|
2553 + |
|
2554 + /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */ |
|
2555 + z = z.subtract (t2); |
|
2556 + t2 = t2.epowy (); |
|
2557 + |
|
2558 + /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */ |
|
2559 + t1 = t1.multiply (t2); |
|
2560 + } |
|
2561 + |
|
2562 + mperr ("*** ERROR IN LN, ITERATION NOT CONVERGING ***"); |
|
2563 return z; |
|
2564 + |
|
2565 } |
|
2566 |
|
2567 + /* RETURNS MP Y = Ln (1+X) IF X IS AN MP NUMBER SATISFYING THE |
|
2568 + * CONDITION ABS (X) < 1/B, ERROR OTHERWISE. |
|
2569 + * USES NEWTONS METHOD TO SOLVE THE EQUATION |
|
2570 + * EXP1(-Y) = X, THEN REVERSES re_sign OF Y. |
|
2571 + */ |
|
2572 + private Number lns () |
|
2573 + { |
|
2574 + /* ln (1+0) = 0 */ |
|
2575 + if (is_zero ()) |
|
2576 + return this; |
|
2577 + |
|
2578 + /* Get starting approximation -ln (1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */ |
|
2579 + var t2 = copy (); |
|
2580 + var t1 = divide_integer (4); |
|
2581 + t1 = t1.add (new Number.fraction (-1, 3)); |
|
2582 + t1 = multiply (t1); |
|
2583 + t1 = t1.add (new Number.fraction (1, 2)); |
|
2584 + t1 = multiply (t1); |
|
2585 + t1 = t1.add (new Number.integer (-1)); |
|
2586 + var z = multiply (t1); |
|
2587 + |
|
2588 + /* Solve using Newtons method */ |
|
2589 + var it0 = 5; |
|
2590 + var t = it0; |
|
2591 + Number t3; |
|
2592 + while (true) |
|
2593 + { |
|
2594 + /* t3 = (e^t3 - 1) */ |
|
2595 + /* z = z - (t2 + t3 + (t2 * t3)) */ |
|
2596 + t3 = z.epowy (); |
|
2597 + t3 = t3.add (new Number.integer (-1)); |
|
2598 + t1 = t2.multiply (t3); |
|
2599 + t3 = t3.add (t1); |
|
2600 + t3 = t2.add (t3); |
|
2601 + z = z.subtract (t3); |
|
2602 + if (t >= T) |
|
2603 + break; |
|
2604 + |
|
2605 + /* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE. |
|
2606 + * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE, |
|
2607 + * WE CAN ALMOST DOUBLE T EACH TIME. |
|
2608 + */ |
|
2609 + var ts3 = t; |
|
2610 + var ts2 = t; |
|
2611 + t = T; |
|
2612 + do |
|
2613 + { |
|
2614 + ts2 = t; |
|
2615 + t = (t + it0) / 2; |
|
2616 + } while (t > ts3); |
|
2617 + t = ts2; |
|
2618 + } |
|
2619 + |
|
2620 + /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */ |
|
2621 + if (t3.re_sign != 0 && t3.re_exponent << 1 > it0 - T) |
|
2622 + mperr ("*** ERROR OCCURRED IN LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***"); |
|
2623 + |
|
2624 + z.re_sign = -z.re_sign; |
|
2625 + |
|
2626 + return z; |
|
2627 + } |
|
2628 + |
|
2629 + private Number add_real (int y_sign, Number y) |
|
2630 + { |
|
2631 + var re_sign_prod = y_sign * re_sign; |
|
2632 + |
|
2633 + /* 0 + y = y */ |
|
2634 + if (is_zero ()) |
|
2635 + { |
|
2636 + if (y_sign != y.re_sign) |
|
2637 + return y.invert_sign (); |
|
2638 + else |
|
2639 + return y; |
|
2640 + } |
|
2641 + /* x + 0 = x */ |
|
2642 + else if (y.is_zero ()) |
|
2643 + return this; |
|
2644 + |
|
2645 + var exp_diff = re_exponent - y.re_exponent; |
|
2646 + var med = exp_diff.abs (); |
|
2647 + var x_largest = false; |
|
2648 + if (exp_diff < 0) |
|
2649 + x_largest = false; |
|
2650 + else if (exp_diff > 0) |
|
2651 + x_largest = true; |
|
2652 + else |
|
2653 + { |
|
2654 + /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */ |
|
2655 + if (re_sign_prod < 0) |
|
2656 + { |
|
2657 + /* signs are not equal. find out which mantissa is larger. */ |
|
2658 + int j; |
|
2659 + for (j = 0; j < T; j++) |
|
2660 + { |
|
2661 + int i = re_fraction[j] - y.re_fraction[j]; |
|
2662 + if (i != 0) |
|
2663 + { |
|
2664 + if (i < 0) |
|
2665 + x_largest = false; |
|
2666 + else if (i > 0) |
|
2667 + x_largest = true; |
|
2668 + break; |
|
2669 + } |
|
2670 + } |
|
2671 + |
|
2672 + /* Both mantissas equal, so result is zero. */ |
|
2673 + if (j >= T) |
|
2674 + return new Number.integer (0); |
|
2675 + } |
|
2676 + } |
|
2677 + |
|
2678 + var z = new Number.integer (0); |
|
2679 + int[] big_fraction, small_fraction; |
|
2680 + if (x_largest) |
|
2681 + { |
|
2682 + z.re_sign = re_sign; |
|
2683 + z.re_exponent = re_exponent; |
|
2684 + big_fraction = re_fraction; |
|
2685 + small_fraction = y.re_fraction; |
|
2686 + } |
|
2687 + else |
|
2688 + { |
|
2689 + z.re_sign = y_sign; |
|
2690 + z.re_exponent = y.re_exponent; |
|
2691 + big_fraction = y.re_fraction; |
|
2692 + small_fraction = re_fraction; |
|
2693 + } |
|
2694 + |
|
2695 + /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */ |
|
2696 + for (var i = 3; i >= med; i--) |
|
2697 + z.re_fraction[T + i] = 0; |
|
2698 + |
|
2699 + if (re_sign_prod >= 0) |
|
2700 + { |
|
2701 + /* This is probably insufficient overflow detection, but it makes us not crash at least. */ |
|
2702 + if (T + 3 < med) |
|
2703 + { |
|
2704 + mperr (_("Overflow: the result couldn't be calculated")); |
|
2705 + return new Number.integer (0); |
|
2706 + } |
|
2707 + |
|
2708 + /* HERE DO ADDITION, re_exponent (Y) >= re_exponent (X) */ |
|
2709 + var i = 0; |
|
2710 + for (i = T + 3; i >= T; i--) |
|
2711 + z.re_fraction[i] = small_fraction[i - med]; |
|
2712 + |
|
2713 + var c = 0; |
|
2714 + for (; i >= med; i--) |
|
2715 + { |
|
2716 + c = big_fraction[i] + small_fraction[i - med] + c; |
|
2717 + |
|
2718 + if (c < BASE) |
|
2719 + { |
|
2720 + /* NO CARRY GENERATED HERE */ |
|
2721 + z.re_fraction[i] = c; |
|
2722 + c = 0; |
|
2723 + } |
|
2724 + else |
|
2725 + { |
|
2726 + /* CARRY GENERATED HERE */ |
|
2727 + z.re_fraction[i] = c - BASE; |
|
2728 + c = 1; |
|
2729 + } |
|
2730 + } |
|
2731 + |
|
2732 + for (; i >= 0; i--) |
|
2733 + { |
|
2734 + c = big_fraction[i] + c; |
|
2735 + if (c < BASE) |
|
2736 + { |
|
2737 + z.re_fraction[i] = c; |
|
2738 + i--; |
|
2739 + |
|
2740 + /* NO CARRY POSSIBLE HERE */ |
|
2741 + for (; i >= 0; i--) |
|
2742 + z.re_fraction[i] = big_fraction[i]; |
|
2743 + |
|
2744 + c = 0; |
|
2745 + break; |
|
2746 + } |
|
2747 + |
|
2748 + z.re_fraction[i] = 0; |
|
2749 + c = 1; |
|
2750 + } |
|
2751 + |
|
2752 + /* MUST SHIFT RIGHT HERE AS CARRY OFF END */ |
|
2753 + if (c != 0) |
|
2754 + { |
|
2755 + for (var j = T + 3; j > 0; j--) |
|
2756 + z.re_fraction[j] = z.re_fraction[j - 1]; |
|
2757 + z.re_fraction[0] = 1; |
|
2758 + z.re_exponent++; |
|
2759 + } |
|
2760 + } |
|
2761 + else |
|
2762 + { |
|
2763 + var c = 0; |
|
2764 + var i = 0; |
|
2765 + for (i = T + med - 1; i >= T; i--) |
|
2766 + { |
|
2767 + /* HERE DO SUBTRACTION, ABS (Y) > ABS (X) */ |
|
2768 + z.re_fraction[i] = c - small_fraction[i - med]; |
|
2769 + c = 0; |
|
2770 + |
|
2771 + /* BORROW GENERATED HERE */ |
|
2772 + if (z.re_fraction[i] < 0) |
|
2773 + { |
|
2774 + c = -1; |
|
2775 + z.re_fraction[i] += BASE; |
|
2776 + } |
|
2777 + } |
|
2778 + |
|
2779 + for (; i >= med; i--) |
|
2780 + { |
|
2781 + c = big_fraction[i] + c - small_fraction[i - med]; |
|
2782 + if (c >= 0) |
|
2783 + { |
|
2784 + /* NO BORROW GENERATED HERE */ |
|
2785 + z.re_fraction[i] = c; |
|
2786 + c = 0; |
|
2787 + } |
|
2788 + else |
|
2789 + { |
|
2790 + /* BORROW GENERATED HERE */ |
|
2791 + z.re_fraction[i] = c + BASE; |
|
2792 + c = -1; |
|
2793 + } |
|
2794 + } |
|
2795 + |
|
2796 + for (; i >= 0; i--) |
|
2797 + { |
|
2798 + c = big_fraction[i] + c; |
|
2799 + |
|
2800 + if (c >= 0) |
|
2801 + { |
|
2802 + z.re_fraction[i] = c; |
|
2803 + i--; |
|
2804 + |
|
2805 + /* NO CARRY POSSIBLE HERE */ |
|
2806 + for (; i >= 0; i--) |
|
2807 + z.re_fraction[i] = big_fraction[i]; |
|
2808 + |
|
2809 + break; |
|
2810 + } |
|
2811 + |
|
2812 + z.re_fraction[i] = c + BASE; |
|
2813 + c = -1; |
|
2814 + } |
|
2815 + } |
|
2816 + |
|
2817 + mp_normalize (ref z); |
|
2818 + |
|
2819 + return z; |
|
2820 + } |
|
2821 + |
|
2822 + private Number multiply_real (Number y) |
|
2823 + { |
|
2824 + /* x*0 = 0*y = 0 */ |
|
2825 + if (re_sign == 0 || y.re_sign == 0) |
|
2826 + return new Number.integer (0); |
|
2827 + |
|
2828 + var z = new Number.integer (0); |
|
2829 + z.re_sign = re_sign * y.re_sign; |
|
2830 + z.re_exponent = re_exponent + y.re_exponent; |
|
2831 + |
|
2832 + var r = new Number.integer (0); |
|
2833 + |
|
2834 + /* PERFORM MULTIPLICATION */ |
|
2835 + var c = 8; |
|
2836 + for (var i = 0; i < T; i++) |
|
2837 + { |
|
2838 + var xi = re_fraction[i]; |
|
2839 + |
|
2840 + /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */ |
|
2841 + if (xi == 0) |
|
2842 + continue; |
|
2843 + |
|
2844 + /* Computing MIN */ |
|
2845 + for (var j = 0; j < int.min (T, T + 3 - i); j++) |
|
2846 + r.re_fraction[i+j+1] += xi * y.re_fraction[j]; |
|
2847 + c--; |
|
2848 + if (c > 0) |
|
2849 + continue; |
|
2850 + |
|
2851 + /* CHECK FOR LEGAL BASE B DIGIT */ |
|
2852 + if (xi < 0 || xi >= BASE) |
|
2853 + { |
|
2854 + mperr ("*** ILLEGAL BASE B DIGIT IN CALL TO MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***"); |
|
2855 + return new Number.integer (0); |
|
2856 + } |
|
2857 + |
|
2858 + /* PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME, |
|
2859 + * FASTER THAN DOING IT EVERY TIME. |
|
2860 + */ |
|
2861 + for (var j = T + 3; j >= 0; j--) |
|
2862 + { |
|
2863 + int ri = r.re_fraction[j] + c; |
|
2864 + if (ri < 0) |
|
2865 + { |
|
2866 + mperr ("*** INTEGER OVERFLOW IN MULTIPLY, B TOO LARGE ***"); |
|
2867 + return new Number.integer (0); |
|
2868 + } |
|
2869 + c = ri / BASE; |
|
2870 + r.re_fraction[j] = ri - BASE * c; |
|
2871 + } |
|
2872 + if (c != 0) |
|
2873 + { |
|
2874 + mperr ("*** ILLEGAL BASE B DIGIT IN CALL TO MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***"); |
|
2875 + return new Number.integer (0); |
|
2876 + } |
|
2877 + c = 8; |
|
2878 + } |
|
2879 + |
|
2880 + if (c != 8) |
|
2881 + { |
|
2882 + c = 0; |
|
2883 + for (var i = T + 3; i >= 0; i--) |
|
2884 + { |
|
2885 + int ri = r.re_fraction[i] + c; |
|
2886 + if (ri < 0) |
|
2887 + { |
|
2888 + mperr ("*** INTEGER OVERFLOW IN MULTIPLY, B TOO LARGE ***"); |
|
2889 + return new Number.integer (0); |
|
2890 + } |
|
2891 + c = ri / BASE; |
|
2892 + r.re_fraction[i] = ri - BASE * c; |
|
2893 + } |
|
2894 + |
|
2895 + if (c != 0) |
|
2896 + { |
|
2897 + mperr ("*** ILLEGAL BASE B DIGIT IN CALL TO MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***"); |
|
2898 + return new Number.integer (0); |
|
2899 + } |
|
2900 + } |
|
2901 + |
|
2902 + /* Clear complex part */ |
|
2903 + z.im_sign = 0; |
|
2904 + z.im_exponent = 0; |
|
2905 + for (var i = 0; i < z.im_fraction.length; i++) |
|
2906 + z.im_fraction[i] = 0; |
|
2907 + |
|
2908 + /* NORMALIZE AND ROUND RESULT */ |
|
2909 + // FIXME: Use stack variable because of mp_normalize brokeness |
|
2910 + for (var i = 0; i < SIZE; i++) |
|
2911 + z.re_fraction[i] = r.re_fraction[i]; |
|
2912 + mp_normalize (ref z); |
|
2913 + |
|
2914 + return z; |
|
2915 + } |
|
2916 + |
|
2917 + private Number multiply_integer_real (int64 y) |
|
2918 + { |
|
2919 + /* x*0 = 0*y = 0 */ |
|
2920 + if (is_zero () || y == 0) |
|
2921 + return new Number.integer (0); |
|
2922 + |
|
2923 + /* x*1 = x, x*-1 = -x */ |
|
2924 + // FIXME: Why is this not working? mp_ext is using this function to do a normalization |
|
2925 + /*if (y == 1 || y == -1) |
|
2926 + { |
|
2927 + if (y < 0) |
|
2928 + z = invert_sign (); |
|
2929 + else |
|
2930 + z = this; |
|
2931 + return z; |
|
2932 + }*/ |
|
2933 + |
|
2934 + /* Copy x as z may also refer to x */ |
|
2935 + var z = new Number.integer (0); |
|
2936 + |
|
2937 + if (y < 0) |
|
2938 + { |
|
2939 + y = -y; |
|
2940 + z.re_sign = -re_sign; |
|
2941 + } |
|
2942 + else |
|
2943 + z.re_sign = re_sign; |
|
2944 + z.re_exponent = re_exponent + 4; |
|
2945 + |
|
2946 + /* FORM PRODUCT IN ACCUMULATOR */ |
|
2947 + int64 c = 0; |
|
2948 + |
|
2949 + /* IF y*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE |
|
2950 + * DOUBLE-PRECISION MULTIPLICATION. |
|
2951 + */ |
|
2952 + |
|
2953 + /* Computing MAX */ |
|
2954 + if (y >= int.max (BASE << 3, 32767 / BASE)) |
|
2955 + { |
|
2956 + /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */ |
|
2957 + var j1 = y / BASE; |
|
2958 + var j2 = y - j1 * BASE; |
|
2959 + |
|
2960 + /* FORM PRODUCT */ |
|
2961 + for (var i = T + 3; i >= 0; i--) |
|
2962 + { |
|
2963 + var c1 = c / BASE; |
|
2964 + var c2 = c - BASE * c1; |
|
2965 + var ix = 0; |
|
2966 + if (i > 3) |
|
2967 + ix = re_fraction[i - 4]; |
|
2968 + |
|
2969 + var t = j2 * ix + c2; |
|
2970 + var is = t / BASE; |
|
2971 + c = j1 * ix + c1 + is; |
|
2972 + z.re_fraction[i] = (int) (t - BASE * is); |
|
2973 + } |
|
2974 + } |
|
2975 + else |
|
2976 + { |
|
2977 + int64 ri = 0; |
|
2978 + for (var i = T + 3; i >= 4; i--) |
|
2979 + { |
|
2980 + ri = y * re_fraction[i - 4] + c; |
|
2981 + c = ri / BASE; |
|
2982 + z.re_fraction[i] = (int) (ri - BASE * c); |
|
2983 + } |
|
2984 + |
|
2985 + /* CHECK FOR INTEGER OVERFLOW */ |
|
2986 + if (ri < 0) |
|
2987 + { |
|
2988 + mperr ("*** INTEGER OVERFLOW IN multiply_integer, B TOO LARGE ***"); |
|
2989 + return new Number.integer (0); |
|
2990 + } |
|
2991 + |
|
2992 + /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */ |
|
2993 + for (var i = 3; i >= 0; i--) |
|
2994 + { |
|
2995 + var t = c; |
|
2996 + c = t / BASE; |
|
2997 + z.re_fraction[i] = (int) (t - BASE * c); |
|
2998 + } |
|
2999 + } |
|
3000 + |
|
3001 + /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */ |
|
3002 + while (c != 0) |
|
3003 + { |
|
3004 + for (var i = T + 3; i >= 1; i--) |
|
3005 + z.re_fraction[i] = z.re_fraction[i - 1]; |
|
3006 + var t = c; |
|
3007 + c = t / BASE; |
|
3008 + z.re_fraction[0] = (int) (t - BASE * c); |
|
3009 + z.re_exponent++; |
|
3010 + } |
|
3011 + |
|
3012 + if (c < 0) |
|
3013 + { |
|
3014 + mperr ("*** INTEGER OVERFLOW IN multiply_integer, B TOO LARGE ***"); |
|
3015 + return new Number.integer (0); |
|
3016 + } |
|
3017 + |
|
3018 + z.im_sign = 0; |
|
3019 + z.im_exponent = 0; |
|
3020 + for (var i = 0; i < z.im_fraction.length; i++) |
|
3021 + z.im_fraction[i] = 0; |
|
3022 + mp_normalize (ref z); |
|
3023 + return z; |
|
3024 + } |
|
3025 + |
|
3026 private Number reciprocal_real () |
|
3027 { |
|
3028 /* 1/0 invalid */ |
|
3029 if (is_zero ()) |
|
3030 { |
|
3031 - error = _("Reciprocal of zero is undefined"); |
|
3032 + mperr (_("Reciprocal of zero is undefined")); |
|
3033 return new Number.integer (0); |
|
3034 } |
|
3035 |
|
3036 - var tmp = MPFloat.init2 ((Precision) precision); |
|
3037 - tmp.set_unsigned_integer (1, Round.NEAREST); |
|
3038 - tmp.divide (tmp, re_num, Round.NEAREST); |
|
3039 - var z = copy (); |
|
3040 - z.re_num.clear (); |
|
3041 - z.re_num = tmp; |
|
3042 - return z; |
|
3043 + /* Start by approximating value using floating point */ |
|
3044 + var t1 = copy (); |
|
3045 + t1.re_exponent = 0; |
|
3046 + t1 = new Number.double (1.0 / t1.to_double ()); |
|
3047 + t1.re_exponent -= re_exponent; |
|
3048 + |
|
3049 + var t = 3; |
|
3050 + var it0 = t; |
|
3051 + Number t2; |
|
3052 + while (true) |
|
3053 + { |
|
3054 + /* t1 = t1 - (t1 * ((x * t1) - 1)) (2*t1 - t1^2*x) */ |
|
3055 + t2 = multiply (t1); |
|
3056 + t2 = t2.add (new Number.integer (-1)); |
|
3057 + t2 = t1.multiply (t2); |
|
3058 + t1 = t1.subtract (t2); |
|
3059 + if (t >= T) |
|
3060 + break; |
|
3061 + |
|
3062 + /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE |
|
3063 + * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE). |
|
3064 + */ |
|
3065 + var ts3 = t; |
|
3066 + var ts2 = 0; |
|
3067 + t = T; |
|
3068 + do |
|
3069 + { |
|
3070 + ts2 = t; |
|
3071 + t = (t + it0) / 2; |
|
3072 + } while (t > ts3); |
|
3073 + t = int.min (ts2, T); |
|
3074 + } |
|
3075 + |
|
3076 + /* RETURN IF NEWTON ITERATION WAS CONVERGING */ |
|
3077 + if (t2.re_sign != 0 && (t1.re_exponent - t2.re_exponent) << 1 < T - it0) |
|
3078 + { |
|
3079 + /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL, |
|
3080 + * OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH. |
|
3081 + */ |
|
3082 + mperr ("*** ERROR OCCURRED IN RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***"); |
|
3083 + } |
|
3084 + |
|
3085 + return t1; |
|
3086 } |
|
3087 |
|
3088 + private Number divide_integer_real (int64 y) |
|
3089 + { |
|
3090 + /* x/0 */ |
|
3091 + if (y == 0) |
|
3092 + { |
|
3093 + /* Translators: Error displayed attempted to divide by zero */ |
|
3094 + mperr (_("Division by zero is undefined")); |
|
3095 + return new Number.integer (0); |
|
3096 + } |
|
3097 |
|
3098 + /* 0/y = 0 */ |
|
3099 + if (is_zero ()) |
|
3100 + return new Number.integer (0); |
|
3101 + |
|
3102 + /* Division by -1 or 1 just changes re_sign */ |
|
3103 + if (y == 1 || y == -1) |
|
3104 + { |
|
3105 + if (y < 0) |
|
3106 + return invert_sign (); |
|
3107 + else |
|
3108 + return this; |
|
3109 + } |
|
3110 + |
|
3111 + var z = new Number.integer (0); |
|
3112 + if (y < 0) |
|
3113 + { |
|
3114 + y = -y; |
|
3115 + z.re_sign = -re_sign; |
|
3116 + } |
|
3117 + else |
|
3118 + z.re_sign = re_sign; |
|
3119 + z.re_exponent = re_exponent; |
|
3120 + |
|
3121 + int64 c = 0; |
|
3122 + int64 i = 0; |
|
3123 + |
|
3124 + /* IF y*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE |
|
3125 + * LONG DIVISION. ASSUME AT LEAST 16-BIT WORD. |
|
3126 + */ |
|
3127 + |
|
3128 + /* Computing MAX */ |
|
3129 + var b2 = int.max (BASE << 3, 32767 / BASE); |
|
3130 + if (y < b2) |
|
3131 + { |
|
3132 + /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */ |
|
3133 + int64 r1 = 0; |
|
3134 + do |
|
3135 + { |
|
3136 + c = BASE * c; |
|
3137 + if (i < T) |
|
3138 + c += re_fraction[i]; |
|
3139 + i++; |
|
3140 + r1 = c / y; |
|
3141 + if (r1 < 0) |
|
3142 + { |
|
3143 + mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***"); |
|
3144 + return new Number.integer (0); |
|
3145 + } |
|
3146 + } while (r1 == 0); |
|
3147 + |
|
3148 + /* ADJUST re_exponent AND GET T+4 DIGITS IN QUOTIENT */ |
|
3149 + z.re_exponent += (int) (1 - i); |
|
3150 + z.re_fraction[0] = (int) r1; |
|
3151 + c = BASE * (c - y * r1); |
|
3152 + int64 kh = 1; |
|
3153 + if (i < T) |
|
3154 + { |
|
3155 + kh = T + 1 - i; |
|
3156 + for (var k = 1; k < kh; k++) |
|
3157 + { |
|
3158 + c += re_fraction[i]; |
|
3159 + z.re_fraction[k] = (int) (c / y); |
|
3160 + c = BASE * (c - y * z.re_fraction[k]); |
|
3161 + i++; |
|
3162 + } |
|
3163 + if (c < 0) |
|
3164 + { |
|
3165 + mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***"); |
|
3166 + return new Number.integer (0); |
|
3167 + } |
|
3168 + } |
|
3169 + |
|
3170 + for (var k = kh; k < T + 4; k++) |
|
3171 + { |
|
3172 + z.re_fraction[k] = (int) (c / y); |
|
3173 + c = BASE * (c - y * z.re_fraction[k]); |
|
3174 + } |
|
3175 + if (c < 0) |
|
3176 + { |
|
3177 + mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***"); |
|
3178 + return new Number.integer (0); |
|
3179 + } |
|
3180 + |
|
3181 + mp_normalize (ref z); |
|
3182 + return z; |
|
3183 + } |
|
3184 + |
|
3185 + /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */ |
|
3186 + var j1 = y / BASE; |
|
3187 + var j2 = y - j1 * BASE; |
|
3188 + |
|
3189 + /* LOOK FOR FIRST NONZERO DIGIT */ |
|
3190 + var c2 = 0; |
|
3191 + do |
|
3192 + { |
|
3193 + c = BASE * c + c2; |
|
3194 + c2 = i < T ? re_fraction[i] : 0; |
|
3195 + i++; |
|
3196 + } while (c < j1 || (c == j1 && c2 < j2)); |
|
3197 + |
|
3198 + /* COMPUTE T+4 QUOTIENT DIGITS */ |
|
3199 + z.re_exponent += (int) (1 - i); |
|
3200 + i--; |
|
3201 + |
|
3202 + /* MAIN LOOP FOR LARGE ABS (y) CASE */ |
|
3203 + for (var k = 1; k <= T + 4; k++) |
|
3204 + { |
|
3205 + /* GET APPROXIMATE QUOTIENT FIRST */ |
|
3206 + var ir = c / (j1 + 1); |
|
3207 + |
|
3208 + /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */ |
|
3209 + var iq = c - ir * j1; |
|
3210 + if (iq >= b2) |
|
3211 + { |
|
3212 + /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */ |
|
3213 + ir++; |
|
3214 + iq -= j1; |
|
3215 + } |
|
3216 + |
|
3217 + iq = iq * BASE - ir * j2; |
|
3218 + if (iq < 0) |
|
3219 + { |
|
3220 + /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */ |
|
3221 + ir--; |
|
3222 + iq += y; |
|
3223 + } |
|
3224 + |
|
3225 + if (i < T) |
|
3226 + iq += re_fraction[i]; |
|
3227 + i++; |
|
3228 + var iqj = iq / y; |
|
3229 + |
|
3230 + /* r (K) = QUOTIENT, C = REMAINDER */ |
|
3231 + z.re_fraction[k - 1] = (int) (iqj + ir); |
|
3232 + c = iq - y * iqj; |
|
3233 + |
|
3234 + if (c < 0) |
|
3235 + { |
|
3236 + /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */ |
|
3237 + mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***"); |
|
3238 + return new Number.integer (0); |
|
3239 + } |
|
3240 + } |
|
3241 + |
|
3242 + mp_normalize (ref z); |
|
3243 + |
|
3244 + /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */ |
|
3245 + mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***"); |
|
3246 + return new Number.integer (0); |
|
3247 + } |
|
3248 + |
|
3249 + |
|
3250 + |
|
3251 private Number from_radians (AngleUnit unit) |
|
3252 { |
|
3253 switch (unit) |
|
3254 @@ -1340,23 +2825,146 @@ |
|
3255 } |
|
3256 } |
|
3257 |
|
3258 + /* z = sin (x) -1 >= x >= 1, do_sin = 1 |
|
3259 + * z = cos (x) -1 >= x >= 1, do_sin = 0 |
|
3260 + */ |
|
3261 + private Number sin1 (bool do_sin) |
|
3262 + { |
|
3263 + /* sin (0) = 0, cos (0) = 1 */ |
|
3264 + if (is_zero ()) |
|
3265 + { |
|
3266 + if (do_sin) |
|
3267 + return new Number.integer (0); |
|
3268 + else |
|
3269 + return new Number.integer (1); |
|
3270 + } |
|
3271 + |
|
3272 + var t2 = multiply (this); |
|
3273 + if (t2.compare (new Number.integer (1)) > 0) |
|
3274 + mperr ("*** ABS (X) > 1 IN CALL TO SIN1 ***"); |
|
3275 + |
|
3276 + Number t1; |
|
3277 + int i; |
|
3278 + Number z; |
|
3279 + if (do_sin) |
|
3280 + { |
|
3281 + t1 = this; |
|
3282 + z = t1; |
|
3283 + i = 2; |
|
3284 + } |
|
3285 + else |
|
3286 + { |
|
3287 + t1 = new Number.integer (1); |
|
3288 + z = new Number.integer (0); |
|
3289 + i = 1; |
|
3290 + } |
|
3291 + |
|
3292 + /* Taylor series */ |
|
3293 + /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */ |
|
3294 + var b2 = 2 * int.max (BASE, 64); |
|
3295 + do |
|
3296 + { |
|
3297 + if (T + t1.re_exponent <= 0) |
|
3298 + break; |
|
3299 + |
|
3300 + /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING |
|
3301 + * DIVISION BY I*(I+1) HAS TO BE SPLIT UP. |
|
3302 + */ |
|
3303 + t1 = t2.multiply (t1); |
|
3304 + if (i > b2) |
|
3305 + { |
|
3306 + t1 = t1.divide_integer (-i); |
|
3307 + t1 = t1.divide_integer (i + 1); |
|
3308 + } |
|
3309 + else |
|
3310 + t1 = t1.divide_integer (-i * (i + 1)); |
|
3311 + z = t1.add (z); |
|
3312 + |
|
3313 + i += 2; |
|
3314 + } while (t1.re_sign != 0); |
|
3315 + |
|
3316 + if (!do_sin) |
|
3317 + z = z.add (new Number.integer (1)); |
|
3318 + |
|
3319 + return z; |
|
3320 + } |
|
3321 + |
|
3322 + |
|
3323 private Number sin_real (AngleUnit unit) |
|
3324 { |
|
3325 + /* sin (0) = 0 */ |
|
3326 + if (is_zero ()) |
|
3327 + return new Number.integer (0); |
|
3328 + |
|
3329 var x_radians = to_radians (unit); |
|
3330 - var z = new Number.integer (0); |
|
3331 - var tmp = z.re_num; |
|
3332 - tmp.sin (x_radians.re_num, Round.NEAREST); |
|
3333 - z.re_num = tmp; |
|
3334 + |
|
3335 + var xs = x_radians.re_sign; |
|
3336 + x_radians = x_radians.abs (); |
|
3337 + |
|
3338 + /* USE SIN1 IF ABS (X) <= 1 */ |
|
3339 + Number z; |
|
3340 + if (x_radians.compare (new Number.integer (1)) <= 0) |
|
3341 + z = x_radians.sin1 (true); |
|
3342 + /* FIND ABS (X) MODULO 2PI */ |
|
3343 + else |
|
3344 + { |
|
3345 + z = new Number.pi ().divide_integer (4); |
|
3346 + x_radians = x_radians.divide (z); |
|
3347 + x_radians = x_radians.divide_integer (8); |
|
3348 + x_radians = x_radians.fractional_component (); |
|
3349 + |
|
3350 + /* SUBTRACT 1/2, SAVE re_sign AND TAKE ABS */ |
|
3351 + x_radians = x_radians.add (new Number.fraction (-1, 2)); |
|
3352 + xs = -xs * x_radians.re_sign; |
|
3353 + if (xs == 0) |
|
3354 + return new Number.integer (0); |
|
3355 + |
|
3356 + x_radians.re_sign = 1; |
|
3357 + x_radians = x_radians.multiply_integer (4); |
|
3358 + |
|
3359 + /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */ |
|
3360 + if (x_radians.re_exponent > 0) |
|
3361 + x_radians = x_radians.add (new Number.integer (-2)); |
|
3362 + |
|
3363 + if (x_radians.is_zero ()) |
|
3364 + return new Number.integer (0); |
|
3365 + |
|
3366 + x_radians.re_sign = 1; |
|
3367 + x_radians = x_radians.multiply_integer (2); |
|
3368 + |
|
3369 + /* NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE |
|
3370 + * POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT |
|
3371 + */ |
|
3372 + if (x_radians.re_exponent > 0) |
|
3373 + { |
|
3374 + x_radians = x_radians.add (new Number.integer (-2)); |
|
3375 + x_radians = x_radians.multiply (z); |
|
3376 + z = x_radians.sin1 (false); |
|
3377 + } |
|
3378 + else |
|
3379 + { |
|
3380 + x_radians = x_radians.multiply (z); |
|
3381 + z = x_radians.sin1 (true); |
|
3382 + } |
|
3383 + } |
|
3384 + |
|
3385 + z.re_sign = xs; |
|
3386 return z; |
|
3387 } |
|
3388 |
|
3389 private Number cos_real (AngleUnit unit) |
|
3390 { |
|
3391 - var x_radians = to_radians (unit); |
|
3392 - var tmp = MPFloat.init2 ((Precision) precision); |
|
3393 - tmp.cos (x_radians.re_num, Round.NEAREST); |
|
3394 - var z = new Number.mpfloat (tmp); |
|
3395 - return z; |
|
3396 + /* cos (0) = 1 */ |
|
3397 + if (is_zero ()) |
|
3398 + return new Number.integer (1); |
|
3399 + |
|
3400 + /* Use power series if |x| <= 1 */ |
|
3401 + var z = to_radians (unit).abs (); |
|
3402 + if (z.compare (new Number.integer (1)) <= 0) |
|
3403 + return z.sin1 (false); |
|
3404 + else |
|
3405 + /* cos (x) = sin (π/2 - |x|) */ |
|
3406 + return new Number.pi ().divide_integer (2).subtract (z).sin (AngleUnit.RADIANS); |
|
3407 } |
|
3408 |
|
3409 private Number bitwise (Number y, BitwiseFunc bitwise_operator, int wordlen) |
|
3410 @@ -1370,7 +2978,7 @@ |
|
3411 offset_out = offset1 > offset2 ? offset1 : offset2; |
|
3412 if (offset_out > 0 && (offset_out < offset1 || offset_out < offset2)) |
|
3413 { |
|
3414 - error = ("Overflow. Try a bigger word size"); |
|
3415 + mperr ("Overflow. Try a bigger word size"); |
|
3416 return new Number.integer (0); |
|
3417 } |
|
3418 |
|
3419 @@ -1420,6 +3028,29 @@ |
|
3420 |
|
3421 // FIXME: Re-add overflow and underflow detection |
|
3422 |
|
3423 +static string? mp_error = null; |
|
3424 + |
|
3425 +/* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND |
|
3426 + * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR. |
|
3427 + */ |
|
3428 +public void mperr (string text) |
|
3429 +{ |
|
3430 + mp_error = text; |
|
3431 +} |
|
3432 + |
|
3433 +/* Returns error string or null if no error */ |
|
3434 +// FIXME: Global variable |
|
3435 +public string mp_get_error () |
|
3436 +{ |
|
3437 + return mp_error; |
|
3438 +} |
|
3439 + |
|
3440 +/* Clear any current error */ |
|
3441 +public void mp_clear_error () |
|
3442 +{ |
|
3443 + mp_error = null; |
|
3444 +} |
|
3445 + |
|
3446 /* Sets z from a string representation in 'text'. */ |
|
3447 public Number? mp_set_from_string (string str, int default_base = 10) |
|
3448 { |
|
3449 @@ -1468,7 +3099,6 @@ |
|
3450 |
|
3451 /* Convert integer part */ |
|
3452 var z = new Number.integer (0); |
|
3453 - |
|
3454 while (str.get_next_char (ref index, out c)) |
|
3455 { |
|
3456 var i = char_val (c, number_base); |
|
3457 @@ -1600,6 +3230,73 @@ |
|
3458 return null; |
|
3459 } |
|
3460 |
|
3461 +/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE |
|
3462 + * GREATEST COMMON DIVISOR OF K AND L. |
|
3463 + * SAVE INPUT PARAMETERS IN LOCAL VARIABLES |
|
3464 + */ |
|
3465 +public void mp_gcd (ref int64 k, ref int64 l) |
|
3466 +{ |
|
3467 + var i = k.abs (); |
|
3468 + var j = l.abs (); |
|
3469 + if (j == 0) |
|
3470 + { |
|
3471 + /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */ |
|
3472 + k = 1; |
|
3473 + l = 0; |
|
3474 + if (i == 0) |
|
3475 + k = 0; |
|
3476 + return; |
|
3477 + } |
|
3478 + |
|
3479 + /* EUCLIDEAN ALGORITHM LOOP */ |
|
3480 + do |
|
3481 + { |
|
3482 + i %= j; |
|
3483 + if (i == 0) |
|
3484 + { |
|
3485 + k = k / j; |
|
3486 + l = l / j; |
|
3487 + return; |
|
3488 + } |
|
3489 + j %= i; |
|
3490 + } while (j != 0); |
|
3491 + |
|
3492 + /* HERE J IS THE GCD OF K AND L */ |
|
3493 + k = k / i; |
|
3494 + l = l / i; |
|
3495 +} |
|
3496 + |
|
3497 +// FIXME: Is r.re_fraction large enough? It seems to be in practise but it may be T+4 instead of T |
|
3498 +// FIXME: There is some sort of stack corruption/use of unitialised variables here. Some functions are |
|
3499 +// using stack variables as x otherwise there are corruption errors. e.g. "Cos (45) - 1/Sqrt (2) = -0" |
|
3500 +// (try in scientific mode) |
|
3501 +public void mp_normalize (ref Number x) |
|
3502 +{ |
|
3503 + int start_index; |
|
3504 + |
|
3505 + /* Find first non-zero digit */ |
|
3506 + for (start_index = 0; start_index < SIZE && x.re_fraction[start_index] == 0; start_index++); |
|
3507 + |
|
3508 + /* Mark as zero */ |
|
3509 + if (start_index >= SIZE) |
|
3510 + { |
|
3511 + x.re_sign = 0; |
|
3512 + x.re_exponent = 0; |
|
3513 + return; |
|
3514 + } |
|
3515 + |
|
3516 + /* Shift left so first digit is non-zero */ |
|
3517 + if (start_index > 0) |
|
3518 + { |
|
3519 + x.re_exponent -= start_index; |
|
3520 + var i = 0; |
|
3521 + for (; (i + start_index) < SIZE; i++) |
|
3522 + x.re_fraction[i] = x.re_fraction[i + start_index]; |
|
3523 + for (; i < SIZE; i++) |
|
3524 + x.re_fraction[i] = 0; |
|
3525 + } |
|
3526 +} |
|
3527 + |
|
3528 /* Returns true if x is cannot be represented in a binary word of length 'wordlen' */ |
|
3529 public bool mp_is_overflow (Number x, int wordlen) |
|
3530 { |
|
3531 |
|
3532 --- gnome-calculator-3.18.2/configure.ac 2016-08-09 14:02:40.743464258 +0000 |
|
3533 +++ gnome-calculator-3.18.2/configure.ac 2016-08-09 14:02:52.623794159 +0000 |
|
3534 @@ -25,9 +25,6 @@ |
|
3535 |
|
3536 AC_SUBST([GLIB_REQUIRED]) |
|
3537 |
|
3538 -# FIXME We link to it too, so this check needs to be better.... |
|
3539 -AC_CHECK_HEADER([mpfr.h], [], [AC_MSG_ERROR([The mpfr header is missing. Please, install mpfr])]) |
|
3540 - |
|
3541 PKG_CHECK_MODULES(LIBCALCULATOR, [ |
|
3542 glib-2.0 >= $GLIB_REQUIRED |
|
3543 gio-2.0 >= $GLIB_REQUIRED |
|