[gnome-calculator] ported to MPFR
- From: Arth Patel <arthpatel src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnome-calculator] ported to MPFR
- Date: Mon, 22 Dec 2014 18:42:36 +0000 (UTC)
commit 0e863eeb121d670d8f5266a82b6152e98b2428ea
Author: Daniel Renninghoff <daniel renninghoff gmail com>
Date: Wed Dec 17 21:21:57 2014 +0100
ported to MPFR
Signed-off-by: PioneerAxon <arth svnit gmail com>
data/org.gnome.calculator.gschema.xml | 6 +-
src/Makefile.am | 21 +-
src/currency.vala | 2 +-
src/equation-lexer.vala | 2 +-
src/equation-parser.vala | 125 ++-
src/equation.vala | 6 +-
src/gcalccmd.vala | 10 +-
src/gnome-calculator.vala | 6 +-
src/math-converter.vala | 2 +-
src/math-equation.vala | 85 ++-
src/math-preferences.vala | 2 +-
src/math-window.vala | 2 +-
src/mpfr.vapi | 133 ++
src/number.vala | 2478 ++++++---------------------------
src/serializer.vala | 16 +-
src/test-equation.vala | 32 +-
src/test-number.vala | 16 -
src/unit.vala | 4 +-
18 files changed, 788 insertions(+), 2160 deletions(-)
---
diff --git a/data/org.gnome.calculator.gschema.xml b/data/org.gnome.calculator.gschema.xml
index ee6684f..2a21b4e 100644
--- a/data/org.gnome.calculator.gschema.xml
+++ b/data/org.gnome.calculator.gschema.xml
@@ -21,7 +21,6 @@
<schema path="/org/gnome/calculator/" id="org.gnome.calculator" gettext-domain="calculator">
<key type="i" name="accuracy">
<default>9</default>
- <range min="0" max="9"/>
<summary>Accuracy value</summary>
<description>The number of digits displayed after the numeric point</description>
</key>
@@ -82,5 +81,10 @@
<summary>Target units</summary>
<description>Units to convert the current calculation into</description>
</key>
+ <key type="i" name="precision">
+ <default>2000</default>
+ <summary>Internal precision</summary>
+ <description>The internal precision used with the MPFR library</description>
+ </key>
</schema>
</schemalist>
diff --git a/src/Makefile.am b/src/Makefile.am
index d8a408e..d3a52d7 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -18,6 +18,7 @@ resources.c: $(top_srcdir)/data/gnome-calculator.gresource.xml $(shell $(GLIB_CO
gnome_calculator_SOURCES = \
config.vapi \
+ mpfr.vapi \
gnome-calculator.vala \
currency.vala \
equation.vala \
@@ -50,10 +51,12 @@ gnome_calculator_VALAFLAGS = \
gnome_calculator_LDADD = \
$(GNOME_CALCULATOR_LIBS) \
- -lm
+ -lmpfr \
+ -lgmp
gcalccmd_SOURCES = \
config.vapi \
+ mpfr.vapi \
gcalccmd.vala \
currency.vala \
equation.vala \
@@ -73,10 +76,12 @@ gcalccmd_VALAFLAGS = \
gcalccmd_LDADD = \
$(GCALCCMD_LIBS) \
- -lm
+ -lmpfr \
+ -lgmp
test_number_SOURCES = \
config.vapi \
+ mpfr.vapi \
test-number.vala \
number.vala \
serializer.vala
@@ -89,10 +94,13 @@ test_number_VALAFLAGS = \
test_number_LDADD = \
$(GCALCCMD_LIBS) \
- -lm
+ -lm \
+ -lmpfr \
+ -lgmp
test_equation_SOURCES = \
config.vapi \
+ mpfr.vapi \
test-equation.vala \
currency.vala \
equation.vala \
@@ -112,10 +120,12 @@ test_equation_VALAFLAGS = \
test_equation_LDADD = \
$(GCALCCMD_LIBS) \
- -lm
+ -lmpfr \
+ -lgmp
test_serializer_SOURCES = \
config.vapi \
+ mpfr.vapi \
test-serializer.vala \
number.vala \
serializer.vala
@@ -128,6 +138,7 @@ test_serializer_VALAFLAGS = \
test_serializer_LDADD = \
$(GCALCCMD_LIBS) \
- -lm
+ -lmpfr \
+ -lgmp
-include $(top_srcdir)/git.mk
diff --git a/src/currency.vala b/src/currency.vala
index 6435599..db66ee6 100644
--- a/src/currency.vala
+++ b/src/currency.vala
@@ -444,7 +444,7 @@ public class CurrencyManager : Object
{
warning ("Couldn't download IMF currency rate file: %s", e.message);
}
-
+
downloading_imf_rates = false;
load_rates ();
}
diff --git a/src/equation-lexer.vala b/src/equation-lexer.vala
index 652769e..20731af 100644
--- a/src/equation-lexer.vala
+++ b/src/equation-lexer.vala
@@ -120,7 +120,7 @@ public class PreLexer : Object
if (c == ',' || c == '.')
return LexerTokenType.PL_DECIMAL;
-
+
if (c.isdigit ())
return LexerTokenType.PL_DIGIT;
diff --git a/src/equation-parser.vala b/src/equation-parser.vala
index d06485b..29f66bd 100644
--- a/src/equation-parser.vala
+++ b/src/equation-parser.vala
@@ -78,7 +78,20 @@ public abstract class RNode : ParseNode
var r = right.solve ();
if (r == null)
return null;
- return solve_r (r);
+ var z = solve_r (r);
+
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ var tmpleft = right;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+ return z;
}
public abstract Number solve_r (Number r);
@@ -97,7 +110,20 @@ public abstract class LRNode : ParseNode
var r = right.solve ();
if (l == null || r == null)
return null;
- return solve_lr (l, r);
+ var z = solve_lr (l, r);
+
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ var tmpleft = left;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+ return z;
}
public abstract Number solve_lr (Number left, Number r);
@@ -109,7 +135,7 @@ public class ConstantNode : ParseNode
{
base (parser, token, precedence, associativity);
}
-
+
public override Number? solve ()
{
return mp_set_from_string (token.text, parser.number_base);
@@ -118,7 +144,7 @@ public class ConstantNode : ParseNode
public class AssignNode : RNode
{
- public AssignNode (Parser parser, LexerToken? token, uint precedence, Associativity associativity)
+ public AssignNode (Parser parser, LexerToken? token, uint precedence, Associativity associativity)
{
base (parser, token, precedence, associativity);
}
@@ -170,7 +196,7 @@ public class VariableNode : ParseNode
}
public override Number? solve ()
- {
+ {
/* If defined, then get the variable */
var ans = parser.get_variable (token.text);
if (ans != null)
@@ -235,6 +261,18 @@ public class VariableWithPowerNode : ParseNode
value = value.multiply (t);
}
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ var tmpleft = left;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+
return value;
}
}
@@ -353,6 +391,15 @@ public class FunctionNode : ParseNode
if (tmp != null)
tmp = tmp.xpowy_integer (pow);
+
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ parser.set_error (ErrorCode.MP, Number.error);
+ Number.error = null;
+ }
+
return tmp;
}
}
@@ -527,7 +574,17 @@ public class DivideNode : LRNode
public override Number solve_lr (Number l, Number r)
{
- return l.divide (r);
+ var z = l.divide (r);
+ if (Number.error != null)
+ {
+ var tmpleft = right;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+ return z;
}
}
@@ -547,7 +604,21 @@ public class ModulusDivideNode : LRNode
var mod = right.solve ();
if (base_value == null || exponent == null || mod == null)
return null;
- return base_value.modular_exponentiation (exponent, mod);
+ var z = base_value.modular_exponentiation (exponent, mod);
+
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ var tmpleft = left;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+
+ return z;
}
else
{
@@ -555,7 +626,21 @@ public class ModulusDivideNode : LRNode
var r = right.solve ();
if (l == null || r == null)
return null;
- return solve_lr (l, r);
+ var z = solve_lr (l, r);
+
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ var tmpleft = left;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+
+ return z;
}
}
@@ -624,7 +709,21 @@ public class XPowYIntegerNode : ParseNode
if (val == null)
return null;
- return val.xpowy_integer (pow);
+ var z = val.xpowy_integer (pow);
+
+ /* check for errors */
+ Number.check_flags ();
+ if (Number.error != null)
+ {
+ var tmpleft = left;
+ var tmpright = right;
+ while (tmpleft.left != null) tmpleft = tmpleft.left;
+ while (tmpright.right != null) tmpright = tmpright.right;
+ parser.set_error (ErrorCode.MP, Number.error, tmpleft.token.start_index,
tmpright.token.end_index);
+ Number.error = null;
+ }
+
+ return z;
}
}
@@ -937,10 +1036,10 @@ public class Parser
}
representation_base = this.representation_base;
- error_code = ErrorCode.NONE;
- error_token = null;
- error_start = 0;
- error_end = 0;
+ error_code = this.error;
+ error_token = this.error_token;
+ error_start = this.error_token_start;
+ error_end = this.error_token_end;
return ans;
}
diff --git a/src/equation.vala b/src/equation.vala
index f1a0afb..3c49610 100644
--- a/src/equation.vala
+++ b/src/equation.vala
@@ -117,15 +117,17 @@ public class Equation : Object
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)
{
var parser = new EquationParser (this, expression);
- mp_clear_error ();
+ Number.error = null;
var z = parser.parse (out representation_base, out error_code, out error_token, out error_start, out
error_end);
/* Error during parsing */
if (error_code != ErrorCode.NONE)
+ {
return null;
+ }
- if (mp_get_error () != null)
+ if (Number.error != null)
{
error_code = ErrorCode.MP;
return null;
diff --git a/src/gcalccmd.vala b/src/gcalccmd.vala
index 88bf3f6..dc8783e 100644
--- a/src/gcalccmd.vala
+++ b/src/gcalccmd.vala
@@ -34,9 +34,15 @@ static void solve (string equation)
result_serializer.set_representation_base (representation_base);
if (z != null)
- stdout.printf ("%s\n", result_serializer.to_string (z));
+ {
+ var str = result_serializer.to_string (z);
+ if (result_serializer.error != null)
+ stderr.printf ("%s\n", result_serializer.error);
+ else
+ stdout.printf ("%s\n", str);
+ }
else if (ret == ErrorCode.MP)
- stderr.printf ("Error %s\n", mp_get_error ());
+ stderr.printf ("Error %s\n", Number.error);
else
stderr.printf ("Error %d\n", ret);
}
diff --git a/src/gnome-calculator.vala b/src/gnome-calculator.vala
index 694d384..dc8df49 100644
--- a/src/gnome-calculator.vala
+++ b/src/gnome-calculator.vala
@@ -31,7 +31,7 @@ public class Calculator : Gtk.Application
{ "about", about_cb, null, null, null },
{ "quit", quit_cb, null, null, null },
};
-
+
public Calculator ()
{
Object (flags : ApplicationFlags.NON_UNIQUE);
@@ -56,6 +56,7 @@ public class Calculator : Gtk.Application
var target_currency = settings.get_string ("target-currency");
var source_units = settings.get_string ("source-units");
var target_units = settings.get_string ("target-units");
+ var precision = settings.get_int ("precision");
var equation = new MathEquation ();
equation.accuracy = accuracy;
@@ -68,6 +69,7 @@ public class Calculator : Gtk.Application
equation.target_currency = target_currency;
equation.source_units = source_units;
equation.target_units = target_units;
+ Number.precision = precision;
add_action_entries (app_entries, this);
@@ -171,7 +173,7 @@ public class Calculator : Gtk.Application
}
else if (error == ErrorCode.MP)
{
- stderr.printf ("Error: %s\n", mp_get_error ());
+ stderr.printf ("Error: %s\n", Number.error);
return Posix.EXIT_FAILURE;
}
else
diff --git a/src/math-converter.vala b/src/math-converter.vala
index 15eba6f..1a69dc8 100644
--- a/src/math-converter.vala
+++ b/src/math-converter.vala
@@ -18,7 +18,7 @@ public class MathConverter : Gtk.Box
private Gtk.ComboBox to_combo;
private Gtk.Label result_label;
-
+
public signal void changed ();
public MathConverter (MathEquation equation)
diff --git a/src/math-equation.vala b/src/math-equation.vala
index cb61057..fa60685 100644
--- a/src/math-equation.vala
+++ b/src/math-equation.vala
@@ -47,7 +47,7 @@ public class MathEquation : Gtk.SourceBuffer
{
private Gtk.TextTag ans_tag;
- /* Word size in bits */
+ /* Word size in bits */
private int _word_size;
public int word_size
{
@@ -781,6 +781,11 @@ public class MathEquation : Gtk.SourceBuffer
ans_start_mark = create_mark (null, start, false);
ans_end_mark = create_mark (null, end, true);
apply_tag (ans_tag, start, end);
+
+ if (serializer.error != null)
+ {
+ status = serializer.error;
+ }
}
public new void insert (string text)
@@ -954,11 +959,13 @@ public class MathEquation : Gtk.SourceBuffer
break;
case ErrorCode.MP:
- if (mp_get_error () != null)
- solvedata.error = mp_get_error ();
- else if (error_token != null) /* Uncategorized error. Show error token to user */
+ if (Number.error != null) // LEGACY, should never be run
{
- solvedata.error = _("Malformed expression at token '%s'").printf (error_token);
+ solvedata.error = Number.error;
+ }
+ else if (error_token != null) // should always be run
+ {
+ solvedata.error = _("%s").printf (error_token);
solvedata.error_start = error_start;
solvedata.error_end = error_end;
}
@@ -1003,6 +1010,8 @@ public class MathEquation : Gtk.SourceBuffer
/* Fix thousand separator offsets in the start and end offsets of error token. */
error_token_fix_thousands_separator ();
+ /* Fix missing Parenthesis before the start and after the end offsets of error token */
+ error_token_fix_parenthesis ();
/* Notify the GUI about the change in error token locations. */
notify_property ("error-token-end");
@@ -1074,6 +1083,70 @@ public class MathEquation : Gtk.SourceBuffer
}
}
+ /* Fix the offsets to consider starting and ending parenthesis */
+ private void error_token_fix_parenthesis ()
+ {
+ unichar c;
+ int count = 0;
+ int real_end = display.index_of_nth_char (error_token_end);
+ int real_start = display.index_of_nth_char (error_token_start);
+
+ /* checks if there are more opening/closing parenthesis than closing/opening parenthesis */
+ for (int i = real_start; display.get_next_char (ref i, out c) && i <= real_end;)
+ {
+ if (c.to_string () == "(") count++;
+ if (c.to_string () == ")") count--;
+ }
+
+ /* if there are more opening than closing parenthesis and there are closing parenthesis
+ after the end offset, include those in the offsets */
+ for (int i = real_end; display.get_next_char (ref i, out c) && count > 0;)
+ {
+ if (c.to_string () == ")")
+ {
+ state.error_token_end++;
+ count--;
+ }
+ else
+ {
+ break;
+ }
+ }
+
+ /* the same for closing parenthesis */
+ for (int i = real_start; display.get_prev_char (ref i, out c) && count < 0;)
+ {
+ if (c.to_string () == "(")
+ {
+ state.error_token_start--;
+ count++;
+ }
+ else
+ {
+ break;
+ }
+ }
+
+ real_end = display.index_of_nth_char (error_token_end);
+ real_start = display.index_of_nth_char (error_token_start);
+
+ unichar d;
+
+ /* if there are opening parenthesis directly before aswell as closing parenthesis directly after the
offsets, include those aswell */
+ while (display.get_next_char (ref real_end, out d) && display.get_prev_char (ref real_start, out c))
+ {
+ if (c.to_string () == "(" && d.to_string () == ")")
+ {
+ state.error_token_start--;
+ state.error_token_end++;
+ }
+ else
+ {
+ break;
+ }
+ }
+ }
+
private void* factorize_real ()
{
var x = number;
@@ -1159,7 +1232,7 @@ public class MathEquation : Gtk.SourceBuffer
var z = number;
if (z == null)
{
- /* This message is displayed in the status bar when a bit shift operation is performed and the
display does not contain a number */
+ /* This message is displayed in the status bar when a bit shift operation is performed and the
display does not contain a number */
status = _("No sane value to bitwise shift");
return;
}
diff --git a/src/math-preferences.vala b/src/math-preferences.vala
index 525082d..d1256a1 100644
--- a/src/math-preferences.vala
+++ b/src/math-preferences.vala
@@ -90,7 +90,7 @@ public class MathPreferencesDialog : Gtk.Dialog
var string = _("Show %d decimal _places");
var tokens = string.split ("%d", 2);
- var decimal_places_adjustment = new Gtk.Adjustment (0.0, 0.0, 9.0, 1.0, 1.0, 0.0);
+ var decimal_places_adjustment = new Gtk.Adjustment (0.0, 0.0, 100.0, 1.0, 1.0, 0.0);
decimal_places_spin = new Gtk.SpinButton (decimal_places_adjustment, 0.0, 0);
if (tokens.length > 0)
diff --git a/src/math-window.vala b/src/math-window.vala
index 9690807..4ba7b07 100644
--- a/src/math-window.vala
+++ b/src/math-window.vala
@@ -90,7 +90,7 @@ public class MathWindow : Gtk.ApplicationWindow
var vbox = new Gtk.Box (Gtk.Orientation.VERTICAL, 6);
vbox.border_width = 6;
- main_vbox.pack_start (vbox, true, true, 0);
+ main_vbox.pack_start (vbox, true, true, 0);
vbox.show ();
var scrolled_window = new Gtk.ScrolledWindow (null, null);
diff --git a/src/mpfr.vapi b/src/mpfr.vapi
new file mode 100644
index 0000000..8c08cab
--- /dev/null
+++ b/src/mpfr.vapi
@@ -0,0 +1,133 @@
+/*
+ * Copyright (C) 2014 Daniel Renninghoff
+ *
+ * This program is free software: you can redistribute it and/or modify it under
+ * the terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 2 of the License, or (at your option) any later
+ * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
+ * license.
+ */
+
+/* A small guide on using MPFR with gnome-calculator:
+ * Using it is pretty much self-explanatory when you look at the code in
+ * number.vala.
+ * C syntax: mpfr_add (result, op1, op2, MPFR_RNDN);
+ * Vala syntax: result.add (op1, op2, Round.NEAREST);
+ *
+ * The result has to be initialized with init2() before using it (same in C).
+ * Since MPFR is a C library you have to do manual memory management. This means
+ * that after init2()ing something, you have to clear() it at the end. Clearing
+ * re_num and im_num is taken are of by the destructor of Number. Just make sure
+ * to manually clear any temporary helper variables you use.
+ */
+
+[CCode (cheader_filename="mpfr.h")]
+namespace MPFR {
+ [SimpleType]
+ [IntegerType]
+ [CCode (cname = "mpfr_prec_t", has_type_id = false)]
+ public struct Precision {}
+
+ [CCode (cname = "mpfr_rnd_t", has_type_id = false)]
+ public enum Round {
+ [CCode (cname = "MPFR_RNDN")]
+ NEAREST,
+ [CCode (cname = "MPFR_RNDZ")]
+ ZERO,
+ [CCode (cname = "MPFR_RNDU")]
+ UP,
+ [CCode (cname = "MPFR_RNDD")]
+ DOWN,
+ [CCode (cname = "MPFR_RNDA")]
+ AWAY,
+ [CCode (cname = "MPFR_RNDF")]
+ FAITHFUL,
+ [CCode (cname = "MPFR_RNDNA")]
+ NEARESTAWAY
+ }
+
+ [CCode (cname = "__mpfr_struct", cprefix = "mpfr_", has_destroy_function = false, copy_function = "",
has_type_id = false)]
+ public struct MPFloat {
+ [CCode (cname="mpfr_init2")]
+ public MPFloat.init2 (Precision prec);
+
+ [CCode (cname="mpfr_set_ui")]
+ public int set_unsigned_integer (ulong op, Round rnd);
+ [CCode (cname="mpfr_set_si")]
+ public int set_signed_integer (long op, Round rnd);
+ [CCode (cname="mpfr_set_flt")]
+ public int set_float (float op, Round rnd);
+ [CCode (cname="mpfr_set_d")]
+ public int set_double (double op, Round rnd);
+ [CCode (cname="mpfr_set")]
+ public int set (MPFloat op, Round rnd);
+
+ public void clear ();
+
+ public int add (MPFloat op1, MPFloat op2, Round rnd);
+ [CCode (cname="mpfr_sub")]
+ public int subtract (MPFloat op1, MPFloat op2, Round rnd);
+ [CCode (cname="mpfr_mul")]
+ public int multiply (MPFloat op1, MPFloat op2, Round rnd);
+ [CCode (cname="mpfr_div")]
+ public int divide (MPFloat op1, MPFloat op2, Round rnd);
+
+ [CCode (cname="mpfr_get_si")]
+ public long get_signed_integer (Round rnd);
+ [CCode (cname="mpfr_get_ui")]
+ public ulong get_unsigned_integer (Round rnd);
+ [CCode (cname="mpfr_get_flt")]
+ public float get_float (Round rnd);
+ [CCode (cname="mpfr_get_d")]
+ public double get_double (Round rnd);
+
+ public int const_pi (Round rnd);
+ [CCode (cname="mpfr_zero_p")]
+ public bool is_zero ();
+ public int sgn ();
+ [CCode (cname="mpfr_equal_p")]
+ public bool is_equal (MPFloat op2);
+ public int cmp (MPFloat op2);
+ public int sqrt (MPFloat op, Round rnd);
+ public int neg (MPFloat op, Round rnd);
+ [CCode (cname="mpfr_pow_si")]
+ public int power_signed_integer (MPFloat op1, long op2, Round rnd);
+ [CCode (cname="mpfr_mul_si")]
+ public int multiply_signed_integer (MPFloat op1, long op2, Round rnd);
+ [CCode (cname="mpfr_div_si")]
+ public int divide_signed_integer (MPFloat op1, long op2, Round rnd);
+ public int floor (MPFloat op);
+ [CCode (cname="mpfr_pow")]
+ public int power (MPFloat op1, MPFloat op2, Round rnd);
+ public int exp (MPFloat op, Round rnd);
+ public int log (MPFloat op, Round rnd);
+ public int sin (MPFloat op, Round rnd);
+ public int cos (MPFloat op, Round rnd);
+ public int tan (MPFloat op, Round rnd);
+ public int asin (MPFloat op, Round rnd);
+ public int acos (MPFloat op, Round rnd);
+ public int atan (MPFloat op, Round rnd);
+ public int sinh (MPFloat op, Round rnd);
+ public int cosh (MPFloat op, Round rnd);
+ public int tanh (MPFloat op, Round rnd);
+ public int asinh (MPFloat op, Round rnd);
+ public int acosh (MPFloat op, Round rnd);
+ public int atanh (MPFloat op, Round rnd);
+ public int abs (MPFloat op, Round rnd);
+ public int root (MPFloat op, ulong k, Round rnd);
+ public int rint (MPFloat op, Round rnd);
+ public int frac (MPFloat op, Round rnd);
+ public int ceil (MPFloat op);
+ public int trunc (MPFloat op);
+ public int round (MPFloat op);
+ [CCode (cname="mpfr_integer_p")]
+ public int is_integer ();
+ public int gamma (MPFloat op, Round rnd);
+ }
+
+ [CCode (cname = "mpfr_underflow_p")]
+ public int mpfr_is_underflow ();
+
+ [CCode (cname = "mpfr_overflow_p")]
+ public int mpfr_is_overflow ();
+}
diff --git a/src/number.vala b/src/number.vala
index 3d9c01f..d6b419c 100644
--- a/src/number.vala
+++ b/src/number.vala
@@ -27,20 +27,7 @@
* THE MP USERS GUIDE.
*/
-/* Size of the multiple precision values */
-private const int SIZE = 1000;
-
-/* Base for numbers */
-private const int BASE = 10000;
-
-//2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT
-//MP.t = (int) ((float) (accuracy) * Math.log ((float)10.) / Math.log ((float) BASE) + (float) 2.0);
-//if (MP.t > SIZE)
-//{
-// mperr ("SIZE TOO SMALL IN CALL TO MPSET, INCREASE SIZE AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d
***", MP.t);
-// MP.t = SIZE;
-//}
-private const int T = 100;
+using MPFR;
private delegate int BitwiseFunc (int v1, int v2);
@@ -51,75 +38,40 @@ public enum AngleUnit
GRADIANS
}
-/* Object for a high precision floating point number representation
- *
- * x = sign * (BASE^(exponent-1) + BASE^(exponent-2) + ...)
- */
+/* Object for a high precision floating point number representation */
public class Number : Object
{
- /* Sign (+1, -1) or 0 for the value zero */
- public int re_sign;
- public int im_sign;
+ /* real and imaginary part of a Number */
+ private MPFloat re_num { get; set; }
+ private MPFloat im_num { get; set; }
- /* Exponent (to base BASE) */
- public int re_exponent;
- public int im_exponent;
+ public static ulong precision { get; set; default = 1000; }
- /* Normalized fraction */
- public int re_fraction[1000]; // SIZE
- public int im_fraction[1000]; // SIZE
+ /* Stores the error msg if an error occurs during calculation. Otherwise should be null */
+ public static string? error { get; set; default = null; }
public Number.integer (int64 value)
{
- if (value < 0)
- {
- value = -value;
- re_sign = -1;
- }
- else if (value > 0)
- re_sign = 1;
- else
- re_sign = 0;
-
- while (value != 0)
- {
- re_fraction[re_exponent] = (int) (value % BASE);
- re_exponent++;
- value /= BASE;
- }
- for (var i = 0; i < re_exponent / 2; i++)
- {
- int t = re_fraction[i];
- re_fraction[i] = re_fraction[re_exponent - i - 1];
- re_fraction[re_exponent - i - 1] = t;
- }
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_signed_integer ((long)value, Round.NEAREST);
+ re_num = tmp;
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp2.set_unsigned_integer (0, Round.NEAREST);
+ im_num = tmp2;
}
public Number.unsigned_integer (uint64 x)
{
- if (x == 0)
- re_sign = 0;
- else
- re_sign = 1;
-
- while (x != 0)
- {
- re_fraction[re_exponent] = (int) (x % BASE);
- x = x / BASE;
- re_exponent++;
- }
- for (var i = 0; i < re_exponent / 2; i++)
- {
- int t = re_fraction[i];
- re_fraction[i] = re_fraction[re_exponent - i - 1];
- re_fraction[re_exponent - i - 1] = t;
- }
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_unsigned_integer ((ulong)x, Round.NEAREST);
+ re_num = tmp;
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp2.set_unsigned_integer (0, Round.NEAREST);
+ im_num = tmp2;
}
public Number.fraction (int64 numerator, int64 denominator)
{
- mp_gcd (ref numerator, ref denominator);
-
if (denominator < 0)
{
numerator = -numerator;
@@ -129,203 +81,39 @@ public class Number : Object
Number.integer (numerator);
if (denominator != 1)
{
- var z = divide_integer (denominator);
- re_sign = z.re_sign;
- im_sign = z.im_sign;
- re_exponent = z.re_exponent;
- im_exponent = z.im_exponent;
- for (var i = 0; i < z.re_fraction.length; i++)
- {
- re_fraction[i] = z.re_fraction[i];
- im_fraction[i] = z.im_fraction[i];
- }
+ var tmp = re_num;
+ tmp.divide_signed_integer (re_num, (long) denominator, Round.NEAREST);
+ re_num = tmp;
}
}
- public Number.float (float value)
+ /* Helper constructor. Creates new Number from already existing MPFloat. */
+ public Number.mpfloat (MPFloat value)
{
- var z = new Number.integer (0);
-
- if (value != 0)
- {
- /* CHECK SIGN */
- var rj = 0f;
- if (value < 0.0f)
- {
- z.re_sign = -1;
- rj = -value;
- }
- else if (value > 0.0f)
- {
- z.re_sign = 1;
- rj = value;
- }
-
- /* INCREASE IE AND DIVIDE RJ BY 16. */
- var ie = 0;
- while (rj >= 1.0f)
- {
- ie++;
- rj *= 0.0625f;
- }
- while (rj < 0.0625f)
- {
- ie--;
- rj *= 16.0f;
- }
-
- /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
- * SET re_exponent TO 0
- */
- z.re_exponent = 0;
-
- /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
- for (var i = 0; i < T + 4; i++)
- {
- rj *= BASE;
- z.re_fraction[i] = (int) rj;
- rj -= z.re_fraction[i];
- }
-
- /* NORMALIZE RESULT */
- mp_normalize (ref z);
-
- /* Computing MAX */
- var ib = int.max (BASE * 7 * BASE, 32767) / 16;
- var tp = 1;
-
- /* NOW MULTIPLY BY 16**IE */
- if (ie < 0)
- {
- var k = -ie;
- for (var i = 1; i <= k; i++)
- {
- tp <<= 4;
- if (tp <= ib && tp != BASE && i < k)
- continue;
- z = z.divide_integer (tp);
- tp = 1;
- }
- }
- else if (ie > 0)
- {
- for (var i = 1; i <= ie; i++)
- {
- tp <<= 4;
- if (tp <= ib && tp != BASE && i < ie)
- continue;
- z = z.multiply_integer (tp);
- tp = 1;
- }
- }
- }
-
- re_sign = z.re_sign;
- im_sign = z.im_sign;
- re_exponent = z.re_exponent;
- im_exponent = z.im_exponent;
- for (var i = 0; i < z.re_fraction.length; i++)
- {
- re_fraction[i] = z.re_fraction[i];
- im_fraction[i] = z.im_fraction[i];
- }
+ re_num = value;
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_unsigned_integer (0, Round.NEAREST);
+ im_num = tmp;
}
-
+
public Number.double (double value)
{
- var z = new Number.integer (0);
-
- if (value != 0)
- {
- /* CHECK SIGN */
- var dj = 0.0;
- if (value < 0.0)
- {
- z.re_sign = -1;
- dj = -value;
- }
- else if (value > 0.0)
- {
- z.re_sign = 1;
- dj = value;
- }
-
- /* INCREASE IE AND DIVIDE DJ BY 16. */
- var ie = 0;
- for (ie = 0; dj >= 1.0; ie++)
- dj *= 1.0/16.0;
-
- for ( ; dj < 1.0/16.0; ie--)
- dj *= 16.0;
-
- /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
- * SET re_exponent TO 0
- */
- z.re_exponent = 0;
-
- /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
- for (var i = 0; i < T + 4; i++)
- {
- dj *= (double) BASE;
- z.re_fraction[i] = (int) dj;
- dj -= (double) z.re_fraction[i];
- }
-
- /* NORMALIZE RESULT */
- mp_normalize (ref z);
-
- /* Computing MAX */
- var ib = int.max (BASE * 7 * BASE, 32767) / 16;
- var tp = 1;
-
- /* NOW MULTIPLY BY 16**IE */
- if (ie < 0)
- {
- var k = -ie;
- for (var i = 1; i <= k; ++i)
- {
- tp <<= 4;
- if (tp <= ib && tp != BASE && i < k)
- continue;
- z = z.divide_integer (tp);
- tp = 1;
- }
- }
- else if (ie > 0)
- {
- for (var i = 1; i <= ie; ++i)
- {
- tp <<= 4;
- if (tp <= ib && tp != BASE && i < ie)
- continue;
- z = z.multiply_integer (tp);
- tp = 1;
- }
- }
- }
-
- re_sign = z.re_sign;
- im_sign = z.im_sign;
- re_exponent = z.re_exponent;
- im_exponent = z.im_exponent;
- for (var i = 0; i < z.re_fraction.length; i++)
- {
- re_fraction[i] = z.re_fraction[i];
- im_fraction[i] = z.im_fraction[i];
- }
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_double (value, Round.NEAREST);
+ re_num = tmp;
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp2.set_unsigned_integer (0, Round.NEAREST);
+ im_num = tmp2;
}
public Number.complex (Number x, Number y)
{
- re_sign = x.re_sign;
- re_exponent = x.re_exponent;
- for (var i = 0; i < im_fraction.length; i++)
- re_fraction[i] = x.re_fraction[i];
-
- im_sign = y.re_sign;
- im_exponent = y.re_exponent;
- for (var i = 0; i < im_fraction.length; i++)
- im_fraction[i] = y.re_fraction[i];
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp set (x.re_num, Round.NEAREST);
+ re_num = tmp;
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp2 set (y.re_num, Round.NEAREST);
+ im_num = tmp2;
}
public Number.polar (Number r, Number theta, AngleUnit unit = AngleUnit.RADIANS)
@@ -336,30 +124,36 @@ public class Number : Object
}
public Number.eulers ()
- {
- var z = new Number.integer (1).epowy ();
- re_sign = z.re_sign;
- im_sign = z.im_sign;
- re_exponent = z.re_exponent;
- im_exponent = z.im_exponent;
- for (var i = 0; i < z.re_fraction.length; i++)
- {
- re_fraction[i] = z.re_fraction[i];
- im_fraction[i] = z.im_fraction[i];
- }
+ {
+ var tmp = MPFloat.init2 ((Precision) precision);
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp2.set_unsigned_integer (1, Round.NEAREST);
+ /* e^1, since mpfr doesn't have a function to return e */
+ tmp.exp (tmp2, Round.NEAREST);
+ re_num = tmp;
+ var tmp3 = MPFloat.init2 ((Precision) precision);
+ tmp3.set_unsigned_integer (0, Round.NEAREST);
+ im_num = tmp3;
}
public Number.i ()
{
- im_sign = 1;
- im_exponent = 1;
- im_fraction[0] = 1;
+ var tmp = MPFloat.init2 ((Precision) precision);
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp.set_unsigned_integer (0, Round.NEAREST);
+ tmp2.set_unsigned_integer (1, Round.NEAREST);
+ re_num = tmp;
+ im_num = tmp2;
}
public Number.pi ()
{
- // FIXME: Should generate PI to required accuracy
- Number.double (Math.PI);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.const_pi (Round.NEAREST);
+ re_num = tmp;
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp2.set_unsigned_integer (0, Round.NEAREST);
+ im_num = tmp2;
}
/* Sets z to be a uniform random number in the range [0, 1] */
@@ -368,123 +162,42 @@ public class Number : Object
Number.double (Random.next_double ());
}
- public int64 to_integer ()
+ ~Number ()
{
- int64 z = 0;
-
- /* |x| <= 1 */
- if (re_sign == 0 || re_exponent <= 0)
- return 0;
-
- /* Multiply digits together */
- for (var i = 0; i < re_exponent; i++)
- {
- var t = z;
- z = z * BASE + re_fraction[i];
-
- /* Check for overflow */
- if (z <= t)
- return 0;
- }
-
- /* Validate result */
- var v = z;
- for (var i = re_exponent - 1; i >= 0; i--)
- {
- /* Get last digit */
- var digit = v - (v / BASE) * BASE;
- if (re_fraction[i] != digit)
- return 0;
-
- v /= BASE;
- }
- if (v != 0)
- return 0;
+ re_num.clear ();
+ im_num.clear ();
+ }
- return re_sign * z;
+ public int64 to_integer ()
+ {
+ return re_num.get_signed_integer (Round.NEAREST);
}
public uint64 to_unsigned_integer ()
{
- /* x <= 1 */
- if (re_sign <= 0 || re_exponent <= 0)
- return 0;
-
- /* Multiply digits together */
- uint64 z = 0;
- for (var i = 0; i < re_exponent; i++)
- {
- var t = z;
- z = z * BASE + re_fraction[i];
-
- /* Check for overflow */
- if (z <= t)
- return 0;
- }
-
- /* Validate result */
- var v = z;
- for (var i = re_exponent - 1; i >= 0; i--)
- {
- /* Get last digit */
- var digit = (uint64) v - (v / BASE) * BASE;
- if (re_fraction[i] != digit)
- return 0;
-
- v /= BASE;
- }
- if (v != 0)
- return 0;
-
- return z;
+ return re_num.get_unsigned_integer (Round.NEAREST);
}
public float to_float ()
{
- if (is_zero ())
- return 0f;
-
- var z = 0f;
- for (var i = 0; i < T; i++)
- {
- if (re_fraction[i] != 0)
- z += re_fraction[i] * Math.powf (BASE, re_exponent - i - 1);
- }
-
- if (re_sign < 0)
- return -z;
- else
- return z;
+ return re_num.get_float (Round.NEAREST);
}
public double to_double ()
{
- if (is_zero ())
- return 0d;
-
- var z = 0d;
- for (var i = 0; i < T; i++)
- {
- if (re_fraction[i] != 0)
- z += re_fraction[i] * Math.pow (BASE, re_exponent - i - 1);
- }
-
- if (re_sign < 0)
- return -z;
- else
- return z;
+ return re_num.get_double (Round.NEAREST);
}
/* Return true if the value is x == 0 */
public bool is_zero ()
{
- return re_sign == 0 && im_sign == 0;
+ return re_num.is_zero () && im_num.is_zero ();
}
/* Return true if x < 0 */
public bool is_negative ()
{
- return re_sign < 0;
+ return re_num.sgn () < 0;
}
/* Return true if x is integer */
@@ -493,36 +206,7 @@ public class Number : Object
if (is_complex ())
return false;
- /* This fix is required for 1/3 repiprocal not being detected as an integer */
- /* Multiplication and division by 10000 is used to get around a
- * limitation to the "fix" for Sun bugtraq bug #4006391 in the
- * floor () routine in mp.c, when the re_exponent is less than 1.
- */
- var t3 = new Number.integer (10000);
- var t1 = multiply (t3);
- t1 = t1.divide (t3);
- var t2 = t1.floor ();
- return t1.equals (t2);
-
- /* Correct way to check for integer */
- /*
-
- // Zero is an integer
- if (is_zero ())
- return true;
-
- // fractional
- if (re_exponent <= 0)
- return false;
-
- // Look for fractional components
- for (var i = re_exponent; i < SIZE; i++)
- {
- if (re_fraction[i] != 0)
- return false;
- }
-
- return true;*/
+ return re_num.is_integer () != 0;
}
/* Return true if x is a positive integer */
@@ -531,7 +215,7 @@ public class Number : Object
if (is_complex ())
return false;
else
- return re_sign >= 0 && is_integer ();
+ return re_num.sgn () >= 0 && is_integer ();
}
/* Return true if x is a natural number (an integer ≥ 0) */
@@ -540,19 +224,34 @@ public class Number : Object
if (is_complex ())
return false;
else
- return re_sign > 0 && is_integer ();
+ return re_num.sgn () > 0 && is_integer ();
}
/* Return true if x has an imaginary component */
public bool is_complex ()
{
- return im_sign != 0;
+ return !im_num.is_zero ();
+ }
+
+ /* Return error if overflow or underflow */
+ public static void check_flags ()
+ {
+ if (mpfr_is_underflow () != 0)
+ {
+ /* Translators: Error displayed when underflow error occured */
+ error = _("Underflow error");
+ }
+ else if (mpfr_is_overflow () != 0)
+ {
+ /* Translators: Error displayed when overflow error occured */
+ error = _("Overflow error");
+ }
}
/* Return true if x == y */
public bool equals (Number y)
{
- return compare (y) == 0;
+ return re_num.is_equal (y.re_num);
}
/* Returns:
@@ -562,62 +261,25 @@ public class Number : Object
*/
public int compare (Number y)
{
- if (re_sign != y.re_sign)
- {
- if (re_sign > y.re_sign)
- return 1;
- else
- return -1;
- }
-
- /* x = y = 0 */
- if (is_zero ())
- return 0;
-
- /* See if numbers are of different magnitude */
- if (re_exponent != y.re_exponent)
- {
- if (re_exponent > y.re_exponent)
- return re_sign;
- else
- return -re_sign;
- }
-
- /* Compare fractions */
- for (var i = 0; i < SIZE; i++)
- {
- if (re_fraction[i] == y.re_fraction[i])
- continue;
-
- if (re_fraction[i] > y.re_fraction[i])
- return re_sign;
- else
- return -re_sign;
- }
-
- /* x = y */
- return 0;
+ return re_num.cmp (y.re_num);
}
/* Sets z = sgn (x) */
public Number sgn ()
{
- if (is_zero ())
- return new Number.integer (0);
- else if (is_negative ())
- return new Number.integer (-1);
- else
- return new Number.integer (1);
+ var z = new Number.integer (re_num.sgn ());
+ return z;
}
/* Sets z = −x */
public Number invert_sign ()
{
- var z = copy ();
-
- z.re_sign = -z.re_sign;
- z.im_sign = -z.im_sign;
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.neg (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ var tmp_im = z.im_num;
+ tmp_im.neg (im_num, Round.NEAREST);
+ z.im_num = tmp_im;
return z;
}
@@ -632,13 +294,13 @@ public class Number : Object
x_real = x_real.multiply (x_real);
x_im = x_im.multiply (x_im);
var z = x_real.add (x_im);
- return z.root (2);
+ return z.sqrt ();
}
else
{
- var z = copy ();
- if (z.re_sign < 0)
- z.re_sign = -z.re_sign;
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.abs (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z;
}
}
@@ -649,7 +311,7 @@ public class Number : Object
if (is_zero ())
{
/* Translators: Error display when attempting to take argument of zero */
- mperr (_("Argument not defined for zero"));
+ error = _("Argument not defined for zero");
return new Number.integer (0);
}
@@ -693,8 +355,12 @@ public class Number : Object
/* Sets z = ‾̅x */
public Number conjugate ()
{
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.neg (im_num, Round.NEAREST);
var z = copy ();
- z.im_sign = -z.im_sign;
+ var tmp2 = z.im_num;
+ tmp2.clear ();
+ z.im_num = tmp;
return z;
}
@@ -702,13 +368,11 @@ public class Number : Object
public Number real_component ()
{
var z = copy ();
-
- /* Clear imaginary component */
- z.im_sign = 0;
- z.im_exponent = 0;
- for (var i = 0; i < z.im_fraction.length; i++)
- z.im_fraction[i] = 0;
-
+ var tmp = z.im_num;
+ tmp.clear ();
+ tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_unsigned_integer (0, Round.NEAREST);
+ z.im_num = tmp;
return z;
}
@@ -716,64 +380,30 @@ public class Number : Object
public Number imaginary_component ()
{
/* Copy imaginary component to real component */
- var z = new Number ();
- z.re_sign = im_sign;
- z.re_exponent = im_exponent;
-
- for (var i = 0; i < z.im_fraction.length; i++)
- z.re_fraction[i] = im_fraction[i];
-
- /* Clear (old) imaginary component */
- z.im_sign = 0;
- z.im_exponent = 0;
- for (var i = 0; i < z.im_fraction.length; i++)
- z.im_fraction[i] = 0;
-
+ var z = copy ();
+ var tmp = z.re_num;
+ tmp.clear ();
+ z.re_num = z.im_num;
+ tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_unsigned_integer (0, Round.NEAREST);
+ z.im_num = tmp;
return z;
}
public Number integer_component ()
{
- /* Clear re_fraction */
- var z = copy ();
- for (var i = z.re_exponent; i < SIZE; i++)
- z.re_fraction[i] = 0;
- z.im_sign = 0;
- z.im_exponent = 0;
- for (var i = 0; i < z.im_fraction.length; i++)
- z.im_fraction[i] = 0;
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.trunc (re_num);
+ var z = new Number.mpfloat (tmp);
return z;
}
/* Sets z = x mod 1 */
public Number fractional_component ()
{
- /* fractional component of zero is 0 */
- if (is_zero ())
- return new Number.integer (0);
-
- /* All fractional */
- if (re_exponent <= 0)
- return this;
-
- /* Shift fractional component */
- var shift = re_exponent;
- for (var i = shift; i < SIZE && re_fraction[i] == 0; i++)
- shift++;
- var z = new Number.integer (0);
- z.re_sign = re_sign;
- z.re_exponent = re_exponent - shift;
- for (var i = 0; i < SIZE; i++)
- {
- if (i + shift >= SIZE)
- z.re_fraction[i] = 0;
- else
- z.re_fraction[i] = re_fraction[i + shift];
- }
- if (z.re_fraction[0] == 0)
- z.re_sign = 0;
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.frac (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z;
}
@@ -786,64 +416,28 @@ public class Number : Object
/* Sets z = ⌊x⌋ */
public Number floor ()
{
- /* Integer component of zero = 0 */
- if (is_zero ())
- return this;
-
- /* If all fractional then no integer component */
- if (re_exponent <= 0)
- {
- if (is_negative ())
- return new Number.integer (-1);
- else
- return new Number.integer (0);
- }
-
- /* Clear fractional component */
- var z = copy ();
- var have_fraction = false;
- for (var i = z.re_exponent; i < SIZE; i++)
- {
- if (z.re_fraction[i] != 0)
- have_fraction = true;
- z.re_fraction[i] = 0;
- }
- z.im_sign = 0;
- z.im_exponent = 0;
- for (var i = 0; i < z.im_fraction.length; i++)
- z.im_fraction[i] = 0;
-
- if (have_fraction && is_negative ())
- z = z.add (new Number.integer (-1));
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.floor (re_num);
+ var z = new Number.mpfloat (tmp);
return z;
}
/* Sets z = ⌈x⌉ */
public Number ceiling ()
{
- var z = floor ();
- var f = fractional_component ();
- if (f.is_zero ())
- return z;
- return z.add (new Number.integer (1));
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.ceil (re_num);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = [x] */
public Number round ()
{
- var do_floor = !is_negative ();
-
- var f = fractional_component ();
- f = f.multiply_integer (2);
- f = f.abs ();
- if (f.compare (new Number.integer (1)) >= 0)
- do_floor = !do_floor;
-
- if (do_floor)
- return floor ();
- else
- return ceiling ();
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.round (re_num);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = 1 ÷ x */
@@ -868,6 +462,7 @@ public class Number : Object
/* Sets z = e^x */
public Number epowy ()
{
+
/* e^0 = 1 */
if (is_zero ())
return new Number.integer (1);
@@ -887,9 +482,23 @@ public class Number : Object
/* Sets z = x^y */
public Number xpowy (Number y)
{
- if (y.is_integer ())
- return xpowy_integer (y.to_integer ());
- else
+ /* 0^-n invalid */
+ if (is_zero () && y.is_negative ())
+ {
+ /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */
+ error = _("The power of zero is undefined for a negative exponent");
+ return new Number.integer (0);
+ }
+
+ /* 0^0 is indeterminate */
+ if (is_zero () && y.is_zero ())
+ {
+ /* Translators: Error displayed when attempted to raise 0 to power of zero */
+ error = _("Zero raised to zero is undefined");
+ return new Number.integer (0);
+ }
+
+ if (!y.is_integer ())
{
var reciprocal = y.reciprocal ();
if (reciprocal.is_integer ())
@@ -897,6 +506,29 @@ public class Number : Object
else
return pwr (y);
}
+
+ Number t;
+ Number t2;
+ if (y.is_negative ())
+ {
+ t = reciprocal ();
+ t2 = y.invert_sign ();
+ }
+ else
+ {
+ t = this;
+ t2 = y;
+ }
+
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.power (t.re_num, t2.re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ var tmp2 = z.im_num;
+ tmp2.clear ();
+ tmp = MPFloat.init2 ((Precision) precision);
+ tmp set (t.im_num, Round.NEAREST);
+ z.im_num = tmp;
+ return z;
}
/* Sets z = x^y */
@@ -906,7 +538,7 @@ public class Number : Object
if (is_zero () && n < 0)
{
/* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */
- mperr (_("The power of zero is undefined for a negative exponent"));
+ error = _("The power of zero is undefined for a negative exponent");
return new Number.integer (0);
}
@@ -914,22 +546,10 @@ public class Number : Object
if (is_zero () && n == 0)
{
/* Translators: Error displayed when attempted to raise 0 to power of zero */
- mperr (_("Zero raised to zero is undefined"));
+ error = _("Zero raised to zero is undefined");
return new Number.integer (0);
}
- /* x^0 = 1 */
- if (n == 0)
- return new Number.integer (1);
-
- /* 0^n = 0 */
- if (is_zero ())
- return new Number.integer (0);
-
- /* x^1 = x */
- if (n == 1)
- return this;
-
Number t;
if (n < 0)
{
@@ -937,17 +557,44 @@ public class Number : Object
n = -n;
}
else
+ {
t = this;
+ }
- var z = new Number.integer (1);
- while (n != 0)
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.power_signed_integer (t.re_num, (long) n, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ var tmp2 = z.im_num;
+ tmp2.clear ();
+ tmp = MPFloat.init2 ((Precision) precision);
+ tmp set (t.im_num, Round.NEAREST);
+ z.im_num = tmp;
+ return z;
+ }
+
+ private Number pwr (Number y)
+ {
+ /* (-x)^y imaginary */
+ /* FIXME: Make complex numbers optional */
+ /*if (re_sign < 0)
{
- if (n % 2 == 1)
- z = z.multiply (t);
- t = t.multiply (t);
- n = n / 2;
+ mperr (_("The power of negative numbers is only defined for integer exponents"));
+ return new Number.integer (0);
+ }*/
+
+ /* 0^y = 0, 0^-y undefined */
+ if (is_zero ())
+ {
+ if (y.is_negative ())
+ error = _("The power of zero is undefined for a negative exponent");
+ return new Number.integer (0);
}
- return z;
+
+ /* x^0 = 1 */
+ if (y.is_zero ())
+ return new Number.integer (1);
+
+ return y.multiply (ln ()).epowy ();
}
/* Sets z = n√x */
@@ -976,21 +623,10 @@ public class Number : Object
/* Sets z = √x */
public Number sqrt ()
{
- if (is_zero ())
- return this;
-
- /* FIXME: Make complex numbers optional */
- /*if (re_sign < 0)
- {
- mperr (_("Square root is undefined for negative values"));
- return new Number.integer (0);
- }*/
- else
- {
- var t = root (-2);
- var z = multiply (t);
- return z.ext (t.re_fraction[0], z.re_fraction[0]);
- }
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.sqrt (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = ln x */
@@ -1000,7 +636,7 @@ public class Number : Object
if (is_zero ())
{
/* Translators: Error displayed when attempting to take logarithm of zero */
- mperr (_("Logarithm of zero is undefined"));
+ error = _("Logarithm of zero is undefined");
return new Number.integer (0);
}
@@ -1033,7 +669,7 @@ public class Number : Object
if (is_zero ())
{
/* Translators: Error displayed when attempting to take logarithm of zero */
- mperr (_("Logarithm of zero is undefined"));
+ error = _("Logarithm of zero is undefined");
return new Number.integer (0);
}
@@ -1052,19 +688,20 @@ public class Number : Object
{
/* Factorial Not defined for Complex or for negative numbers */
- if(is_negative () || is_complex ())
+ if (is_negative () || is_complex ())
{
/* Translators: Error displayed when attempted take the factorial of a negative or complex
number */
- mperr (_("Factorial is only defined for non-negative real numbers"));
+ error = _("Factorial is only defined for non-negative real numbers");
return new Number.integer (0);
}
- var val = to_double ();
+ var tmp = add (new Number.integer (1));
+ var tmp2 = MPFloat.init2 ((Precision) precision);
/* Factorial(x) = Gamma(x+1) - This is the formula used to calculate Factorial.*/
- var fact = Math.tgamma(val+1);
+ tmp2.gamma (tmp.re_num, Round.NEAREST);
- return new Number.double(fact);
+ return new Number.mpfloat (tmp2);
}
/* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
@@ -1079,23 +716,41 @@ public class Number : Object
/* Sets z = x + y */
public Number add (Number y)
{
- return add_with_sign (1, y);
+ if (is_complex () || y.is_complex ())
+ {
+ Number real_z, im_z;
+
+ var real_x = real_component ();
+ var im_x = imaginary_component ();
+ var real_y = y.real_component ();
+ var im_y = y.imaginary_component ();
+
+ real_z = real_x.add_real (real_y);
+ im_z = im_x.add_real (im_y);
+
+ return new Number.complex (real_z, im_z);
+ }
+ else
+ return add_real (y);
+ }
+
+ public Number add_real (Number y)
+ {
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.add (re_num, y.re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = x − y */
public Number subtract (Number y)
{
- return add_with_sign (-1, y);
+ return add (y.invert_sign ());
}
/* Sets z = x × y */
public Number multiply (Number y)
{
- /* x*0 = 0*y = 0 */
- if (is_zero () || y.is_zero ())
- return new Number.integer (0);
-
- /* (a+bi)(c+di) = (ac-bd)+(ad+bc)i */
if (is_complex () || y.is_complex ())
{
Number t1, t2, real_z, im_z;
@@ -1119,58 +774,63 @@ public class Number : Object
return multiply_real (y);
}
+ public Number multiply_real (Number y)
+ {
+ var z = new Number.integer (0);
+ var tmp = z.re_num;
+ tmp.multiply (re_num, y.re_num, Round.NEAREST);
+ z.re_num = tmp;
+ return z;
+ }
+
/* Sets z = x × y */
public Number multiply_integer (int64 y)
{
- if (is_complex ())
- {
- var re_z = real_component ().multiply_integer_real (y);;
- var im_z = imaginary_component ().multiply_integer_real (y);
- return new Number.complex (re_z, im_z);
- }
- else
- return multiply_integer_real (y);
+ var z = new Number.integer (0);
+ var tmp = z.re_num;
+ tmp.multiply_signed_integer (re_num, (long) y, Round.NEAREST);
+ z.re_num = tmp;
+ return z;
}
/* Sets z = x ÷ y */
public Number divide (Number y)
{
- /* x/0 */
if (y.is_zero ())
{
/* Translators: Error displayed attempted to divide by zero */
- mperr (_("Division by zero is undefined"));
+ error = _("Division by zero is undefined");
return new Number.integer (0);
}
- /* 0/y = 0 */
- if (is_zero ())
- return this;
+ if (is_complex () || y.is_complex ())
+ {
+ var a = real_component ();
+ var b = imaginary_component ();
+ var c = y.real_component ();
+ var d = y.imaginary_component ();
+
+ var tmp = a.multiply (c).add (b.multiply (d));
+ var tmp_2 = c.xpowy_integer (2).add (d.xpowy_integer (2));
+ var z_1 = tmp.divide (tmp_2);
- /* z = x × y⁻¹ */
- /* FIXME: Set re_exponent to zero to avoid overflow in multiply??? */
- var t = y.reciprocal ();
- var ie = t.re_exponent;
- t.re_exponent = 0;
- var i = t.re_fraction[0];
- var z = multiply (t);
- z = z.ext (i, z.re_fraction[0]);
- z.re_exponent += ie;
+ tmp = b.multiply (c).subtract (a.multiply (d));
+ var z_2 = tmp.divide (tmp_2);
+ var z = new Number.complex (z_1, z_2);
+ return z;
+ }
+
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.divide (re_num, y.re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z;
}
/* Sets z = x ÷ y */
public Number divide_integer (int64 y)
{
- if (is_complex ())
- {
- var re_z = real_component ().divide_integer_real (y);
- var im_z = imaginary_component ().divide_integer_real (y);
- return new Number.complex (re_z, im_z);
- }
- else
- return divide_integer_real (y);
+ return divide (new Number.integer (y));
}
/* Sets z = x mod y */
@@ -1179,7 +839,7 @@ public class Number : Object
if (!is_integer () || !y.is_integer ())
{
/* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
- mperr (_("Modulus division is only defined for integers"));
+ error = _("Modulus division is only defined for integers");
return new Number.integer (0);
}
@@ -1197,7 +857,7 @@ public class Number : Object
/* Sets z = x ^ y mod p */
public Number modular_exponentiation (Number exp, Number mod)
{
- var base_value = this.copy ();
+ var base_value = copy ();
if (exp.is_negative ())
base_value = base_value.reciprocal ();
var exp_value = exp.abs ();
@@ -1267,250 +927,97 @@ public class Number : Object
public Number tan (AngleUnit unit = AngleUnit.RADIANS)
{
/* Check for undefined values */
- var cos_x = cos (unit);
- if (cos_x.is_zero ())
+ var x_radians = to_radians (unit);
+ var check = x_radians.subtract (new Number.pi ().divide_integer (2)).divide (new Number.pi ());
+
+ if (check.is_integer ())
{
/* Translators: Error displayed when tangent value is undefined */
- mperr (_("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)"));
+ error = _("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)");
return new Number.integer (0);
}
- /* tan (x) = sin (x) / cos (x) */
- return sin (unit).divide (cos_x);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.tan (x_radians.re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = sin⁻¹ x */
public Number asin (AngleUnit unit = AngleUnit.RADIANS)
{
- /* asin⁻¹(0) = 0 */
- if (is_zero ())
- return new Number.integer (0);
-
- /* sin⁻¹(x) = tan⁻¹(x / √(1 - x²)), |x| < 1 */
- if (re_exponent <= 0)
- {
- var t1 = new Number.integer (1);
- var t2 = t1;
- t1 = t1.subtract (this);
- t2 = t2.add (this);
- t2 = t1.multiply (t2);
- t2 = t2.root (-2);
- var z = multiply (t2);
- z = z.atan (unit);
- return z;
- }
-
- /* sin⁻¹(1) = π/2, sin⁻¹(-1) = -π/2 */
- var t2 = new Number.integer (re_sign);
- if (equals (t2))
+ if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0)
{
- var z = new Number.pi ().divide_integer (2 * t2.re_sign);
- return z.from_radians (unit);
+ /* Translators: Error displayed when inverse sine value is undefined */
+ error = _("Inverse sine is undefined for values outside [-1, 1]");
+ return new Number.integer (0);
}
- /* Translators: Error displayed when inverse sine value is undefined */
- mperr (_("Inverse sine is undefined for values outside [-1, 1]"));
- return new Number.integer (0);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.asin (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z.from_radians (unit);
}
/* Sets z = cos⁻¹ x */
public Number acos (AngleUnit unit = AngleUnit.RADIANS)
{
- var pi = new Number.pi ();
- var t1 = new Number.integer (1);
- var n1 = new Number.integer (-1);
-
- Number z;
- if (compare (t1) > 0 || compare (n1) < 0)
+ if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0)
{
/* Translators: Error displayed when inverse cosine value is undefined */
- mperr (_("Inverse cosine is undefined for values outside [-1, 1]"));
- z = new Number.integer (0);
- }
- else if (is_zero ())
- z = pi.divide_integer (2);
- else if (equals (t1))
- z = new Number.integer (0);
- else if (equals (n1))
- z = pi;
- else
- {
- /* cos⁻¹(x) = tan⁻¹(√(1 - x²) / x) */
- Number y;
- var t2 = multiply (this);
- t2 = t1.subtract (t2);
- t2 = t2.sqrt ();
- t2 = t2.divide (this);
- y = t2.atan (AngleUnit.RADIANS);
- if (re_sign > 0)
- z = y;
- else
- z = y.add (pi);
+ error = _("Inverse cosine is undefined for values outside [-1, 1]");
+ return new Number.integer (0);
}
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.acos (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z.from_radians (unit);
}
/* Sets z = tan⁻¹ x */
public Number atan (AngleUnit unit = AngleUnit.RADIANS)
{
- if (is_zero ())
- return new Number.integer (0);
-
- var t2 = this;
- var rx = 0f;
- if (re_exponent.abs () <= 2)
- rx = to_float ();
-
- /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
- var q = 1;
- Number z;
- while (t2.re_exponent >= 0)
- {
- if (t2.re_exponent == 0 && 2 * (t2.re_fraction[0] + 1) <= BASE)
- break;
-
- q *= 2;
-
- /* t = t / (√(t² + 1) + 1) */
- z = t2.multiply (t2);
- z = z.add (new Number.integer (1));
- z = z.sqrt ();
- z = z.add (new Number.integer (1));
- t2 = t2.divide (z);
- }
-
- /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
- z = t2;
- var t1 = t2.multiply (t2);
-
- /* SERIES LOOP. REDUCE T IF POSSIBLE. */
- for (var i = 1; ; i += 2)
- {
- if (T + 2 + t2.re_exponent <= 1)
- break;
-
- t2 = t2.multiply (t1).multiply_integer (-i).divide_integer (i + 2);
-
- z = z.add (t2);
- if (t2.is_zero ())
- break;
- }
-
- /* CORRECT FOR ARGUMENT REDUCTION */
- z = z.multiply_integer (q);
-
- /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS re_exponent
- * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
- */
- if (re_exponent.abs () <= 2)
- {
- float ry = z.to_float ();
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
- if (Math.fabs (ry - Math.atan (rx)) >= Math.fabs (ry) * 0.01)
- mperr ("*** ERROR OCCURRED IN ATAN, RESULT INCORRECT ***");
- }
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.atan (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z.from_radians (unit);
}
/* Sets z = sinh x */
public Number sinh ()
{
- /* sinh (0) = 0 */
- if (is_zero ())
- return new Number.integer (0);
-
- /* WORK WITH ABS (X) */
- var abs_x = abs ();
-
- /* If |x| < 1 USE EXP TO AVOID CANCELLATION, otherwise IF TOO LARGE EPOWY GIVES ERROR MESSAGE */
- Number z;
- if (abs_x.re_exponent <= 0)
- {
- /* ((e^|x| + 1) * (e^|x| - 1)) / e^|x| */
- // FIXME: Solves to e^|x| - e^-|x|, why not lower branch always? */
- var exp_x = abs_x.epowy ();
- var a = exp_x.add (new Number.integer (1));
- var b = exp_x.add (new Number.integer (-1));
- z = a.multiply (b);
- z = z.divide (exp_x);
- }
- else
- {
- /* e^|x| - e^-|x| */
- var exp_x = abs_x.epowy ();
- z = exp_x.reciprocal ();
- z = exp_x.subtract (z);
- }
-
- /* DIVIDE BY TWO AND RESTORE re_sign */
- z = z.divide_integer (2);
- return z.multiply_integer (re_sign);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.sinh (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = cosh x */
public Number cosh ()
{
- /* cosh (0) = 1 */
- if (is_zero ())
- return new Number.integer (1);
-
- /* cosh (x) = (e^x + e^-x) / 2 */
- var t = abs ();
- t = t.epowy ();
- var z = t.reciprocal ();
- z = t.add (z);
- return z.divide_integer (2);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.cosh (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = tanh x */
public Number tanh ()
{
- /* tanh (0) = 0 */
- if (is_zero ())
- return new Number.integer (0);
-
- var t = abs ();
-
- /* SEE IF ABS (X) SO LARGE THAT RESULT IS +-1 */
- var r__1 = (float) T * 0.5 * Math.log ((float) BASE);
- var z = new Number.double (r__1);
- if (t.compare (z) > 0)
- return new Number.integer (re_sign);
-
- /* If |x| >= 1/2 use ?, otherwise use ? to avoid cancellation */
- /* |tanh (x)| = (e^|2x| - 1) / (e^|2x| + 1) */
- t = t.multiply_integer (2);
- if (t.re_exponent > 0)
- {
- t = t.epowy ();
- z = t.add (new Number.integer (-1));
- t = t.add (new Number.integer (1));
- z = z.divide (t);
- }
- else
- {
- t = t.epowy ();
- z = t.add (new Number.integer (1));
- t = t.add (new Number.integer (-1));
- z = t.divide (z);
- }
-
- /* Restore re_sign */
- z.re_sign = re_sign * z.re_sign;
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.tanh (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z;
}
/* Sets z = sinh⁻¹ x */
public Number asinh ()
{
- /* sinh⁻¹(x) = ln (x + √(x² + 1)) */
- var t = multiply (this);
- t = t.add (new Number.integer (1));
- t = t.sqrt ();
- t = add (t);
- return t.ln ();
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.asinh (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = cosh⁻¹ x */
@@ -1521,16 +1028,14 @@ public class Number : Object
if (compare (t) < 0)
{
/* Translators: Error displayed when inverse hyperbolic cosine value is undefined */
- mperr (_("Inverse hyperbolic cosine is undefined for values less than one"));
+ error = _("Inverse hyperbolic cosine is undefined for values less than one");
return new Number.integer (0);
}
- /* cosh⁻¹(x) = ln (x + √(x² - 1)) */
- t = multiply (this);
- t = t.add (new Number.integer (-1));
- t = t.sqrt ();
- t = add (t);
- return t.ln ();
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.acosh (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = tanh⁻¹ x */
@@ -1540,26 +1045,24 @@ public class Number : Object
if (compare (new Number.integer (1)) >= 0 || compare (new Number.integer (-1)) <= 0)
{
/* Translators: Error displayed when inverse hyperbolic tangent value is undefined */
- mperr (_("Inverse hyperbolic tangent is undefined for values outside [-1, 1]"));
+ error = _("Inverse hyperbolic tangent is undefined for values outside [-1, 1]");
return new Number.integer (0);
}
- /* atanh (x) = 0.5 * ln ((1 + x) / (1 - x)) */
- var n = add (new Number.integer (1));
- var d = invert_sign ();
- d = d.add (new Number.integer (1));
- var z = n.divide (d);
- z = z.ln ();
- return z.divide_integer (2);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.atanh (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
/* Sets z = boolean AND for each bit in x and z */
public Number and (Number y)
{
- if (!is_positive_integer () || !y.is_positive_integer ())
+ if (!
+ is_positive_integer () || !y.is_positive_integer ())
{
/* Translators: Error displayed when boolean AND attempted on non-integer values */
- mperr (_("Boolean AND is only defined for positive integers"));
+ error = _("Boolean AND is only defined for positive integers");
}
return bitwise (y, (v1, v2) => { return v1 & v2; }, 0);
@@ -1571,7 +1074,7 @@ public class Number : Object
if (!is_positive_integer () || !y.is_positive_integer ())
{
/* Translators: Error displayed when boolean OR attempted on non-integer values */
- mperr (_("Boolean OR is only defined for positive integers"));
+ error = _("Boolean OR is only defined for positive integers");
}
return bitwise (y, (v1, v2) => { return v1 | v2; }, 0);
@@ -1583,7 +1086,7 @@ public class Number : Object
if (!is_positive_integer () || !y.is_positive_integer ())
{
/* Translators: Error displayed when boolean XOR attempted on non-integer values */
- mperr (_("Boolean XOR is only defined for positive integers"));
+ error = _("Boolean XOR is only defined for positive integers");
}
return bitwise (y, (v1, v2) => { return v1 ^ v2; }, 0);
@@ -1595,7 +1098,7 @@ public class Number : Object
if (!is_positive_integer ())
{
/* Translators: Error displayed when boolean XOR attempted on non-integer values */
- mperr (_("Boolean NOT is only defined for positive integers"));
+ error = _("Boolean NOT is only defined for positive integers");
}
return bitwise (new Number.integer (0), (v1, v2) => { return v1 ^ 0xF; }, wordlen);
@@ -1618,7 +1121,7 @@ public class Number : Object
if (!is_integer ())
{
/* Translators: Error displayed when bit shift attempted on non-integer values */
- mperr (_("Shift is only possible on integer values"));
+ error = _("Shift is only possible on integer values");
return new Number.integer (0);
}
@@ -1750,834 +1253,38 @@ public class Number : Object
private Number copy ()
{
var z = new Number ();
- z.re_sign = re_sign;
- z.im_sign = im_sign;
- z.re_exponent = re_exponent;
- z.im_exponent = im_exponent;
- for (var i = 0; i < re_fraction.length; i++)
- {
- z.re_fraction[i] = re_fraction[i];
- z.im_fraction[i] = im_fraction[i];
- }
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ var tmp2 = MPFloat.init2 ((Precision) precision);
+ tmp set(re_num, Round.NEAREST);
+ tmp2 set(im_num, Round.NEAREST);
+ z.re_num = tmp;
+ z.im_num = tmp2;
return z;
}
- private Number add_with_sign (int y_sign, Number y)
- {
- if (is_complex () || y.is_complex ())
- {
- Number real_z, im_z;
-
- var real_x = real_component ();
- var im_x = imaginary_component ();
- var real_y = y.real_component ();
- var im_y = y.imaginary_component ();
-
- real_z = real_x.add_real (y_sign * y.re_sign, real_y);
- im_z = im_x.add_real (y_sign * y.im_sign, im_y);
-
- return new Number.complex (real_z, im_z);
- }
- else
- return add_real (y_sign * y.re_sign, y);
- }
-
private Number epowy_real ()
{
- /* e^0 = 1 */
- if (is_zero ())
- return new Number.integer (1);
-
- /* If |x| < 1 use exp */
- if (re_exponent <= 0)
- return exp ();
-
- /* NOW SAFE TO CONVERT X TO REAL */
- var rx = to_double ();
-
- /* SAVE re_sign AND WORK WITH ABS (X) */
- var xs = re_sign;
- var t2 = abs ();
-
- /* GET fractional AND INTEGER PARTS OF ABS (X) */
- var ix = t2.to_integer ();
- t2 = t2.fractional_component ();
-
- /* ATTACH re_sign TO fractional PART AND COMPUTE EXP OF IT */
- t2.re_sign *= xs;
- var z = t2.exp ();
-
- /* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS (X) LARGE
- * (BUT ONLY ONE EXTRA DIGIT IF T < 4)
- */
- var tss = 0;
- if (T < 4)
- tss = T + 1;
- else
- tss = T + 2;
-
- /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
- /* Computing MIN */
- var t1 = new Number.integer (xs);
-
- t2.re_sign = 0;
- for (var i = 2 ; ; i++)
- {
- if (int.min (tss, tss + 2 + t1.re_exponent) <= 2)
- break;
-
- t1 = t1.divide_integer (i * xs);
- t2 = t2.add (t1);
- if (t1.is_zero ())
- break;
- }
-
- /* RAISE E OR 1/E TO POWER IX */
- if (xs > 0)
- t2 = t2.add (new Number.integer (2));
- t2 = t2.xpowy_integer (ix);
-
- /* MULTIPLY EXPS OF INTEGER AND fractional PARTS */
- z = z.multiply (t2);
-
- /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS (X) LARGE
- * (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
- */
- if (Math.fabs (rx) > 10.0f)
- return z;
-
- var rz = z.to_double ();
- var r__1 = rz - Math.exp (rx);
- if (Math.fabs (r__1) < rz * 0.01f)
- return z;
-
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT
- * B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE
- * RESULT UNDERFLOWED.
- */
- mperr ("*** ERROR OCCURRED IN EPOWY, RESULT INCORRECT ***");
+ var z = copy ();
+ var tmp = z.re_num;
+ tmp.exp (re_num, Round.NEAREST);
+ z.re_num = tmp;
return z;
}
- /* Return e^x for |x| < 1 USING AN o (SQRt (T).m (T)) ALGORITHM
- * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
- * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
- * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
- * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
- * UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL.
- */
- private Number exp ()
- {
- /* e^0 = 1 */
- if (is_zero ())
- return new Number.integer (1);
-
- /* Only defined for |x| < 1 */
- if (re_exponent > 0)
- {
- mperr ("*** ABS (X) NOT LESS THAN 1 IN CALL TO MP_EXP ***");
- return new Number.integer (0);
- }
-
- var t1 = this;
- var rlb = Math.log (BASE);
-
- /* Compute approximately optimal q (and divide x by 2^q) */
- var q = (int) (Math.sqrt (T * 0.48 * rlb) + re_exponent * 1.44 * rlb);
-
- /* HALVE Q TIMES */
- if (q > 0)
- {
- var ib = BASE << 2;
- var ic = 1;
- for (var i = 1; i <= q; i++)
- {
- ic *= 2;
- if (ic < ib && ic != BASE && i < q)
- continue;
- t1 = t1.divide_integer (ic);
- ic = 1;
- }
- }
-
- if (t1.is_zero ())
- return new Number.integer (0);
-
- /* Sum series, reducing t where possible */
- var z = t1.copy ();
- var t2 = t1;
- for (var i = 2; T + t2.re_exponent - z.re_exponent > 0; i++)
- {
- t2 = t1.multiply (t2);
- t2 = t2.divide_integer (i);
- z = t2.add (z);
- if (t2.is_zero ())
- break;
- }
-
- /* Apply (x+1)^2 - 1 = x (2 + x) for q iterations */
- for (var i = 1; i <= q; i++)
- {
- t1 = z.add (new Number.integer (2));
- z = t1.multiply (z);
- }
-
- return z.add (new Number.integer (1));
- }
-
- private Number pwr (Number y)
- {
- /* (-x)^y imaginary */
- /* FIXME: Make complex numbers optional */
- /*if (re_sign < 0)
- {
- mperr (_("The power of negative numbers is only defined for integer exponents"));
- return new Number.integer (0);
- }*/
-
- /* 0^y = 0, 0^-y undefined */
- if (is_zero ())
- {
- if (y.re_sign < 0)
- mperr (_("The power of zero is undefined for a negative exponent"));
- return new Number.integer (0);
- }
-
- /* x^0 = 1 */
- if (y.is_zero ())
- return new Number.integer (1);
-
- return y.multiply (ln ()).epowy ();
- }
-
private Number root_real (int64 n)
{
- /* x^(1/1) = x */
- if (n == 1)
- return this;
-
- /* x^(1/0) invalid */
- if (n == 0)
- {
- mperr (_("Root must be non-zero"));
- return new Number.integer (0);
- }
-
- var np = n.abs ();
-
- /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */
- if (np > int.max (BASE, 64))
- {
- mperr ("*** ABS (N) TOO LARGE IN CALL TO ROOT ***");
- return new Number.integer (0);
- }
-
- /* 0^(1/n) = 0 for positive n */
- if (is_zero ())
- {
- if (n <= 0)
- mperr (_("Negative root of zero is undefined"));
- return new Number.integer (0);
- }
-
- // FIXME: Imaginary root
- if (re_sign < 0 && np % 2 == 0)
- {
- mperr (_("nth root of negative number is undefined for even n"));
- return new Number.integer (0);
- }
-
- /* DIVIDE re_exponent BY NP */
- var ex = re_exponent / np;
-
- /* Get initial approximation */
- var t1 = copy ();
- t1.re_exponent = 0;
- var approximation = Math.exp (((np * ex - re_exponent) * Math.log (BASE) - Math.log (Math.fabs
(t1.to_float ()))) / np);
- t1 = new Number.double (approximation);
- t1.re_sign = re_sign;
- t1.re_exponent -= (int) ex;
-
- /* MAIN ITERATION LOOP */
- var it0 = 3;
- var t = it0;
- Number t2;
- while (true)
- {
- /* t1 = t1 - ((t1 * ((x * t1^np) - 1)) / np) */
- t2 = t1.xpowy_integer (np);
- t2 = multiply (t2);
- t2 = t2.add (new Number.integer (-1));
- t2 = t1.multiply (t2);
- t2 = t2.divide_integer (np);
- t1 = t1.subtract (t2);
-
- /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
- * NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
- */
- if (t >= T)
- break;
-
- var ts3 = t;
- var ts2 = t;
- t = T;
- do
- {
- ts2 = t;
- t = (t + it0) / 2;
- } while (t > ts3);
- t = int.min (ts2, T);
- }
-
- /* NOW r (I2) IS X**(-1/NP)
- * CHECK THAT NEWTON ITERATION WAS CONVERGING
- */
- if (t2.re_sign != 0 && (t1.re_exponent - t2.re_exponent) << 1 < T - it0)
- {
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
- * OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
- * IS NOT ACCURATE ENOUGH.
- */
- mperr ("*** ERROR OCCURRED IN ROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
- }
-
- if (n >= 0)
- {
- t1 = t1.xpowy_integer (n - 1);
- return multiply (t1);
- }
-
- return t1;
- }
-
- /* ROUTINE CALLED BY DIVIDE AND SQRT TO ENSURE THAT
- * RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
- * CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.
- */
- private Number ext (int i, int j)
- {
- if (is_zero () || T <= 2 || i == 0)
- return this;
-
- /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
- var q = (j + 1) / i + 1;
- var s = BASE * re_fraction[T - 2] + re_fraction[T - 1];
-
- /* SET LAST TWO DIGITS TO ZERO */
var z = copy ();
- if (s <= q)
- {
- z.re_fraction[T - 2] = 0;
- z.re_fraction[T - 1] = 0;
- return z;
- }
-
- if (s + q < BASE * BASE)
- return z;
-
- /* ROUND UP HERE */
- z.re_fraction[T - 2] = BASE - 1;
- z.re_fraction[T - 1] = BASE;
-
- /* NORMALIZE X (LAST DIGIT B IS OK IN MULTIPLY_INTEGER) */
- return z.multiply_integer (1);
- }
-
- private Number ln_real ()
- {
- // ln(e^1) = 1, fixes precision loss
- if (equals (new Number.eulers ()))
- return new Number.integer (1);
-
- /* LOOP TO GET APPROXIMATE Ln (X) USING SINGLE-PRECISION */
- var t1 = copy ();
- var z = new Number.integer (0);
- for (var k = 0; k < 10; k++)
- {
- /* COMPUTE FINAL CORRECTION ACCURATELY USING LNS */
- var t2 = t1.add (new Number.integer (-1));
- if (t2.is_zero () || t2.re_exponent + 1 <= 0)
- return z.add (t2.lns ());
-
- /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
- var e = t1.re_exponent;
- t1.re_exponent = 0;
- var rx = t1.to_double ();
- t1.re_exponent = e;
- var rlx = Math.log (rx) + e * Math.log (BASE);
- t2 = new Number.double (-rlx);
-
- /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
- z = z.subtract (t2);
- t2 = t2.epowy ();
-
- /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
- t1 = t1.multiply (t2);
- }
-
- mperr ("*** ERROR IN LN, ITERATION NOT CONVERGING ***");
- return z;
- }
-
- /* RETURNS MP Y = Ln (1+X) IF X IS AN MP NUMBER SATISFYING THE
- * CONDITION ABS (X) < 1/B, ERROR OTHERWISE.
- * USES NEWTONS METHOD TO SOLVE THE EQUATION
- * EXP1(-Y) = X, THEN REVERSES re_sign OF Y.
- */
- private Number lns ()
- {
- /* ln (1+0) = 0 */
- if (is_zero ())
- return this;
-
- /* Get starting approximation -ln (1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */
- var t2 = copy ();
- var t1 = divide_integer (4);
- t1 = t1.add (new Number.fraction (-1, 3));
- t1 = multiply (t1);
- t1 = t1.add (new Number.fraction (1, 2));
- t1 = multiply (t1);
- t1 = t1.add (new Number.integer (-1));
- var z = multiply (t1);
-
- /* Solve using Newtons method */
- var it0 = 5;
- var t = it0;
- Number t3;
- while (true)
- {
- /* t3 = (e^t3 - 1) */
- /* z = z - (t2 + t3 + (t2 * t3)) */
- t3 = z.epowy ();
- t3 = t3.add (new Number.integer (-1));
- t1 = t2.multiply (t3);
- t3 = t3.add (t1);
- t3 = t2.add (t3);
- z = z.subtract (t3);
- if (t >= T)
- break;
-
- /* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
- * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
- * WE CAN ALMOST DOUBLE T EACH TIME.
- */
- var ts3 = t;
- var ts2 = t;
- t = T;
- do
- {
- ts2 = t;
- t = (t + it0) / 2;
- } while (t > ts3);
- t = ts2;
- }
-
- /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
- if (t3.re_sign != 0 && t3.re_exponent << 1 > it0 - T)
- mperr ("*** ERROR OCCURRED IN LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
-
- z.re_sign = -z.re_sign;
-
- return z;
- }
-
- private Number add_real (int y_sign, Number y)
- {
- var re_sign_prod = y_sign * re_sign;
-
- /* 0 + y = y */
- if (is_zero ())
- {
- if (y_sign != y.re_sign)
- return y.invert_sign ();
- else
- return y;
- }
- /* x + 0 = x */
- else if (y.is_zero ())
- return this;
-
- var exp_diff = re_exponent - y.re_exponent;
- var med = exp_diff.abs ();
- var x_largest = false;
- if (exp_diff < 0)
- x_largest = false;
- else if (exp_diff > 0)
- x_largest = true;
- else
- {
- /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
- if (re_sign_prod < 0)
- {
- /* signs are not equal. find out which mantissa is larger. */
- int j;
- for (j = 0; j < T; j++)
- {
- int i = re_fraction[j] - y.re_fraction[j];
- if (i != 0)
- {
- if (i < 0)
- x_largest = false;
- else if (i > 0)
- x_largest = true;
- break;
- }
- }
-
- /* Both mantissas equal, so result is zero. */
- if (j >= T)
- return new Number.integer (0);
- }
- }
-
- var z = new Number.integer (0);
-
- int[] big_fraction, small_fraction;
- if (x_largest)
- {
- z.re_sign = re_sign;
- z.re_exponent = re_exponent;
- big_fraction = re_fraction;
- small_fraction = y.re_fraction;
- }
- else
- {
- z.re_sign = y_sign;
- z.re_exponent = y.re_exponent;
- big_fraction = y.re_fraction;
- small_fraction = re_fraction;
- }
-
- /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
- for (var i = 3; i >= med; i--)
- z.re_fraction[T + i] = 0;
-
- if (re_sign_prod >= 0)
- {
- /* This is probably insufficient overflow detection, but it makes us not crash at least. */
- if (T + 3 < med)
- {
- mperr (_("Overflow: the result couldn't be calculated"));
- return new Number.integer (0);
- }
-
- /* HERE DO ADDITION, re_exponent (Y) >= re_exponent (X) */
- var i = 0;
- for (i = T + 3; i >= T; i--)
- z.re_fraction[i] = small_fraction[i - med];
-
- var c = 0;
- for (; i >= med; i--)
- {
- c = big_fraction[i] + small_fraction[i - med] + c;
-
- if (c < BASE)
- {
- /* NO CARRY GENERATED HERE */
- z.re_fraction[i] = c;
- c = 0;
- }
- else
- {
- /* CARRY GENERATED HERE */
- z.re_fraction[i] = c - BASE;
- c = 1;
- }
- }
-
- for (; i >= 0; i--)
- {
- c = big_fraction[i] + c;
- if (c < BASE)
- {
- z.re_fraction[i] = c;
- i--;
-
- /* NO CARRY POSSIBLE HERE */
- for (; i >= 0; i--)
- z.re_fraction[i] = big_fraction[i];
-
- c = 0;
- break;
- }
-
- z.re_fraction[i] = 0;
- c = 1;
- }
-
- /* MUST SHIFT RIGHT HERE AS CARRY OFF END */
- if (c != 0)
- {
- for (var j = T + 3; j > 0; j--)
- z.re_fraction[j] = z.re_fraction[j - 1];
- z.re_fraction[0] = 1;
- z.re_exponent++;
- }
- }
- else
- {
- var c = 0;
- var i = 0;
- for (i = T + med - 1; i >= T; i--)
- {
- /* HERE DO SUBTRACTION, ABS (Y) > ABS (X) */
- z.re_fraction[i] = c - small_fraction[i - med];
- c = 0;
-
- /* BORROW GENERATED HERE */
- if (z.re_fraction[i] < 0)
- {
- c = -1;
- z.re_fraction[i] += BASE;
- }
- }
-
- for (; i >= med; i--)
- {
- c = big_fraction[i] + c - small_fraction[i - med];
- if (c >= 0)
- {
- /* NO BORROW GENERATED HERE */
- z.re_fraction[i] = c;
- c = 0;
- }
- else
- {
- /* BORROW GENERATED HERE */
- z.re_fraction[i] = c + BASE;
- c = -1;
- }
- }
-
- for (; i >= 0; i--)
- {
- c = big_fraction[i] + c;
-
- if (c >= 0)
- {
- z.re_fraction[i] = c;
- i--;
-
- /* NO CARRY POSSIBLE HERE */
- for (; i >= 0; i--)
- z.re_fraction[i] = big_fraction[i];
-
- break;
- }
-
- z.re_fraction[i] = c + BASE;
- c = -1;
- }
- }
-
- mp_normalize (ref z);
-
- return z;
- }
-
- private Number multiply_real (Number y)
- {
- /* x*0 = 0*y = 0 */
- if (re_sign == 0 || y.re_sign == 0)
- return new Number.integer (0);
-
- var z = new Number.integer (0);
- z.re_sign = re_sign * y.re_sign;
- z.re_exponent = re_exponent + y.re_exponent;
-
- var r = new Number.integer (0);
-
- /* PERFORM MULTIPLICATION */
- var c = 8;
- for (var i = 0; i < T; i++)
- {
- var xi = re_fraction[i];
-
- /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
- if (xi == 0)
- continue;
-
- /* Computing MIN */
- for (var j = 0; j < int.min (T, T + 3 - i); j++)
- r.re_fraction[i+j+1] += xi * y.re_fraction[j];
- c--;
- if (c > 0)
- continue;
-
- /* CHECK FOR LEGAL BASE B DIGIT */
- if (xi < 0 || xi >= BASE)
- {
- mperr ("*** ILLEGAL BASE B DIGIT IN CALL TO MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- return new Number.integer (0);
- }
-
- /* PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
- * FASTER THAN DOING IT EVERY TIME.
- */
- for (var j = T + 3; j >= 0; j--)
- {
- int ri = r.re_fraction[j] + c;
- if (ri < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN MULTIPLY, B TOO LARGE ***");
- return new Number.integer (0);
- }
- c = ri / BASE;
- r.re_fraction[j] = ri - BASE * c;
- }
- if (c != 0)
- {
- mperr ("*** ILLEGAL BASE B DIGIT IN CALL TO MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- return new Number.integer (0);
- }
- c = 8;
- }
-
- if (c != 8)
- {
- c = 0;
- for (var i = T + 3; i >= 0; i--)
- {
- int ri = r.re_fraction[i] + c;
- if (ri < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN MULTIPLY, B TOO LARGE ***");
- return new Number.integer (0);
- }
- c = ri / BASE;
- r.re_fraction[i] = ri - BASE * c;
- }
-
- if (c != 0)
- {
- mperr ("*** ILLEGAL BASE B DIGIT IN CALL TO MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- return new Number.integer (0);
- }
- }
-
- /* Clear complex part */
- z.im_sign = 0;
- z.im_exponent = 0;
- for (var i = 0; i < z.im_fraction.length; i++)
- z.im_fraction[i] = 0;
-
- /* NORMALIZE AND ROUND RESULT */
- // FIXME: Use stack variable because of mp_normalize brokeness
- for (var i = 0; i < SIZE; i++)
- z.re_fraction[i] = r.re_fraction[i];
- mp_normalize (ref z);
-
+ var tmp = z.re_num;
+ tmp.root (re_num, (ulong) n, Round.NEAREST);
+ z.re_num = tmp;
return z;
}
- private Number multiply_integer_real (int64 y)
+ private Number ln_real ()
{
- /* x*0 = 0*y = 0 */
- if (is_zero () || y == 0)
- return new Number.integer (0);
-
- /* x*1 = x, x*-1 = -x */
- // FIXME: Why is this not working? mp_ext is using this function to do a normalization
- /*if (y == 1 || y == -1)
- {
- if (y < 0)
- z = invert_sign ();
- else
- z = this;
- return z;
- }*/
-
- /* Copy x as z may also refer to x */
- var z = new Number.integer (0);
-
- if (y < 0)
- {
- y = -y;
- z.re_sign = -re_sign;
- }
- else
- z.re_sign = re_sign;
- z.re_exponent = re_exponent + 4;
-
- /* FORM PRODUCT IN ACCUMULATOR */
- int64 c = 0;
-
- /* IF y*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
- * DOUBLE-PRECISION MULTIPLICATION.
- */
-
- /* Computing MAX */
- if (y >= int.max (BASE << 3, 32767 / BASE))
- {
- /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
- var j1 = y / BASE;
- var j2 = y - j1 * BASE;
-
- /* FORM PRODUCT */
- for (var i = T + 3; i >= 0; i--)
- {
- var c1 = c / BASE;
- var c2 = c - BASE * c1;
- var ix = 0;
- if (i > 3)
- ix = re_fraction[i - 4];
-
- var t = j2 * ix + c2;
- var is = t / BASE;
- c = j1 * ix + c1 + is;
- z.re_fraction[i] = (int) (t - BASE * is);
- }
- }
- else
- {
- int64 ri = 0;
- for (var i = T + 3; i >= 4; i--)
- {
- ri = y * re_fraction[i - 4] + c;
- c = ri / BASE;
- z.re_fraction[i] = (int) (ri - BASE * c);
- }
-
- /* CHECK FOR INTEGER OVERFLOW */
- if (ri < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN multiply_integer, B TOO LARGE ***");
- return new Number.integer (0);
- }
-
- /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */
- for (var i = 3; i >= 0; i--)
- {
- var t = c;
- c = t / BASE;
- z.re_fraction[i] = (int) (t - BASE * c);
- }
- }
-
- /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
- while (c != 0)
- {
- for (var i = T + 3; i >= 1; i--)
- z.re_fraction[i] = z.re_fraction[i - 1];
- var t = c;
- c = t / BASE;
- z.re_fraction[0] = (int) (t - BASE * c);
- z.re_exponent++;
- }
-
- if (c < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN multiply_integer, B TOO LARGE ***");
- return new Number.integer (0);
- }
-
- z.im_sign = 0;
- z.im_exponent = 0;
- for (var i = 0; i < z.im_fraction.length; i++)
- z.im_fraction[i] = 0;
- mp_normalize (ref z);
-
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.log (re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
return z;
}
@@ -2586,215 +1293,19 @@ public class Number : Object
/* 1/0 invalid */
if (is_zero ())
{
- mperr (_("Reciprocal of zero is undefined"));
+ error = _("Reciprocal of zero is undefined");
return new Number.integer (0);
}
- /* Start by approximating value using floating point */
- var t1 = copy ();
- t1.re_exponent = 0;
- t1 = new Number.double (1.0 / t1.to_double ());
- t1.re_exponent -= re_exponent;
-
- var t = 3;
- var it0 = t;
- Number t2;
- while (true)
- {
- /* t1 = t1 - (t1 * ((x * t1) - 1)) (2*t1 - t1^2*x) */
- t2 = multiply (t1);
- t2 = t2.add (new Number.integer (-1));
- t2 = t1.multiply (t2);
- t1 = t1.subtract (t2);
- if (t >= T)
- break;
-
- /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
- * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
- */
- var ts3 = t;
- var ts2 = 0;
- t = T;
- do
- {
- ts2 = t;
- t = (t + it0) / 2;
- } while (t > ts3);
- t = int.min (ts2, T);
- }
-
- /* RETURN IF NEWTON ITERATION WAS CONVERGING */
- if (t2.re_sign != 0 && (t1.re_exponent - t2.re_exponent) << 1 < T - it0)
- {
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
- * OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
- */
- mperr ("*** ERROR OCCURRED IN RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
- }
-
- return t1;
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.set_unsigned_integer (1, Round.NEAREST);
+ tmp.divide (tmp, re_num, Round.NEAREST);
+ var z = copy ();
+ z.re_num.clear ();
+ z.re_num = tmp;
+ return z;
}
- private Number divide_integer_real (int64 y)
- {
- /* x/0 */
- if (y == 0)
- {
- /* Translators: Error displayed attempted to divide by zero */
- mperr (_("Division by zero is undefined"));
- return new Number.integer (0);
- }
-
- /* 0/y = 0 */
- if (is_zero ())
- return new Number.integer (0);
-
- /* Division by -1 or 1 just changes re_sign */
- if (y == 1 || y == -1)
- {
- if (y < 0)
- return invert_sign ();
- else
- return this;
- }
-
- var z = new Number.integer (0);
- if (y < 0)
- {
- y = -y;
- z.re_sign = -re_sign;
- }
- else
- z.re_sign = re_sign;
- z.re_exponent = re_exponent;
-
- int64 c = 0;
- int64 i = 0;
-
- /* IF y*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
- * LONG DIVISION. ASSUME AT LEAST 16-BIT WORD.
- */
-
- /* Computing MAX */
- var b2 = int.max (BASE << 3, 32767 / BASE);
- if (y < b2)
- {
- /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
- int64 r1 = 0;
- do
- {
- c = BASE * c;
- if (i < T)
- c += re_fraction[i];
- i++;
- r1 = c / y;
- if (r1 < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***");
- return new Number.integer (0);
- }
- } while (r1 == 0);
-
- /* ADJUST re_exponent AND GET T+4 DIGITS IN QUOTIENT */
- z.re_exponent += (int) (1 - i);
- z.re_fraction[0] = (int) r1;
- c = BASE * (c - y * r1);
- int64 kh = 1;
- if (i < T)
- {
- kh = T + 1 - i;
- for (var k = 1; k < kh; k++)
- {
- c += re_fraction[i];
- z.re_fraction[k] = (int) (c / y);
- c = BASE * (c - y * z.re_fraction[k]);
- i++;
- }
- if (c < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***");
- return new Number.integer (0);
- }
- }
-
- for (var k = kh; k < T + 4; k++)
- {
- z.re_fraction[k] = (int) (c / y);
- c = BASE * (c - y * z.re_fraction[k]);
- }
- if (c < 0)
- {
- mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***");
- return new Number.integer (0);
- }
-
- mp_normalize (ref z);
- return z;
- }
-
- /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
- var j1 = y / BASE;
- var j2 = y - j1 * BASE;
-
- /* LOOK FOR FIRST NONZERO DIGIT */
- var c2 = 0;
- do
- {
- c = BASE * c + c2;
- c2 = i < T ? re_fraction[i] : 0;
- i++;
- } while (c < j1 || (c == j1 && c2 < j2));
-
- /* COMPUTE T+4 QUOTIENT DIGITS */
- z.re_exponent += (int) (1 - i);
- i--;
-
- /* MAIN LOOP FOR LARGE ABS (y) CASE */
- for (var k = 1; k <= T + 4; k++)
- {
- /* GET APPROXIMATE QUOTIENT FIRST */
- var ir = c / (j1 + 1);
-
- /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */
- var iq = c - ir * j1;
- if (iq >= b2)
- {
- /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */
- ir++;
- iq -= j1;
- }
-
- iq = iq * BASE - ir * j2;
- if (iq < 0)
- {
- /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
- ir--;
- iq += y;
- }
-
- if (i < T)
- iq += re_fraction[i];
- i++;
- var iqj = iq / y;
-
- /* r (K) = QUOTIENT, C = REMAINDER */
- z.re_fraction[k - 1] = (int) (iqj + ir);
- c = iq - y * iqj;
-
- if (c < 0)
- {
- /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
- mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***");
- return new Number.integer (0);
- }
- }
-
- mp_normalize (ref z);
-
- /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
- mperr ("*** INTEGER OVERFLOW IN DIVIDE_INTEGER, B TOO LARGE ***");
- return new Number.integer (0);
- }
private Number from_radians (AngleUnit unit)
{
@@ -2829,145 +1340,23 @@ public class Number : Object
}
}
- /* z = sin (x) -1 >= x >= 1, do_sin = 1
- * z = cos (x) -1 >= x >= 1, do_sin = 0
- */
- private Number sin1 (bool do_sin)
- {
- /* sin (0) = 0, cos (0) = 1 */
- if (is_zero ())
- {
- if (do_sin)
- return new Number.integer (0);
- else
- return new Number.integer (1);
- }
-
- var t2 = multiply (this);
- if (t2.compare (new Number.integer (1)) > 0)
- mperr ("*** ABS (X) > 1 IN CALL TO SIN1 ***");
-
- Number t1;
- int i;
- Number z;
- if (do_sin)
- {
- t1 = this;
- z = t1;
- i = 2;
- }
- else
- {
- t1 = new Number.integer (1);
- z = new Number.integer (0);
- i = 1;
- }
-
- /* Taylor series */
- /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
- var b2 = 2 * int.max (BASE, 64);
- do
- {
- if (T + t1.re_exponent <= 0)
- break;
-
- /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
- * DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
- */
- t1 = t2.multiply (t1);
- if (i > b2)
- {
- t1 = t1.divide_integer (-i);
- t1 = t1.divide_integer (i + 1);
- }
- else
- t1 = t1.divide_integer (-i * (i + 1));
- z = t1.add (z);
-
- i += 2;
- } while (t1.re_sign != 0);
-
- if (!do_sin)
- z = z.add (new Number.integer (1));
-
- return z;
- }
-
private Number sin_real (AngleUnit unit)
{
- /* sin (0) = 0 */
- if (is_zero ())
- return new Number.integer (0);
-
var x_radians = to_radians (unit);
-
- var xs = x_radians.re_sign;
- x_radians = x_radians.abs ();
-
- /* USE SIN1 IF ABS (X) <= 1 */
- Number z;
- if (x_radians.compare (new Number.integer (1)) <= 0)
- z = x_radians.sin1 (true);
- /* FIND ABS (X) MODULO 2PI */
- else
- {
- z = new Number.pi ().divide_integer (4);
- x_radians = x_radians.divide (z);
- x_radians = x_radians.divide_integer (8);
- x_radians = x_radians.fractional_component ();
-
- /* SUBTRACT 1/2, SAVE re_sign AND TAKE ABS */
- x_radians = x_radians.add (new Number.fraction (-1, 2));
- xs = -xs * x_radians.re_sign;
- if (xs == 0)
- return new Number.integer (0);
-
- x_radians.re_sign = 1;
- x_radians = x_radians.multiply_integer (4);
-
- /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */
- if (x_radians.re_exponent > 0)
- x_radians = x_radians.add (new Number.integer (-2));
-
- if (x_radians.is_zero ())
- return new Number.integer (0);
-
- x_radians.re_sign = 1;
- x_radians = x_radians.multiply_integer (2);
-
- /* NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE
- * POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT
- */
- if (x_radians.re_exponent > 0)
- {
- x_radians = x_radians.add (new Number.integer (-2));
- x_radians = x_radians.multiply (z);
- z = x_radians.sin1 (false);
- }
- else
- {
- x_radians = x_radians.multiply (z);
- z = x_radians.sin1 (true);
- }
- }
-
- z.re_sign = xs;
+ var z = new Number.integer (0);
+ var tmp = z.re_num;
+ tmp.sin (x_radians.re_num, Round.NEAREST);
+ z.re_num = tmp;
return z;
}
private Number cos_real (AngleUnit unit)
{
- /* cos (0) = 1 */
- if (is_zero ())
- return new Number.integer (1);
-
- /* Use power series if |x| <= 1 */
- var z = to_radians (unit).abs ();
- if (z.compare (new Number.integer (1)) <= 0)
- return z.sin1 (false);
- else
- /* cos (x) = sin (π/2 - |x|) */
- return new Number.pi ().divide_integer (2).subtract (z).sin (AngleUnit.RADIANS);
+ var x_radians = to_radians (unit);
+ var tmp = MPFloat.init2 ((Precision) precision);
+ tmp.cos (x_radians.re_num, Round.NEAREST);
+ var z = new Number.mpfloat (tmp);
+ return z;
}
private Number bitwise (Number y, BitwiseFunc bitwise_operator, int wordlen)
@@ -2981,7 +1370,7 @@ public class Number : Object
offset_out = offset1 > offset2 ? offset1 : offset2;
if (offset_out > 0 && (offset_out < offset1 || offset_out < offset2))
{
- mperr ("Overflow. Try a bigger word size");
+ error = ("Overflow. Try a bigger word size");
return new Number.integer (0);
}
@@ -3031,29 +1420,6 @@ public class Number : Object
// FIXME: Re-add overflow and underflow detection
-static string? mp_error = null;
-
-/* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
- * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
- */
-public void mperr (string text)
-{
- mp_error = text;
-}
-
-/* Returns error string or null if no error */
-// FIXME: Global variable
-public string mp_get_error ()
-{
- return mp_error;
-}
-
-/* Clear any current error */
-public void mp_clear_error ()
-{
- mp_error = null;
-}
-
/* Sets z from a string representation in 'text'. */
public Number? mp_set_from_string (string str, int default_base = 10)
{
@@ -3102,6 +1468,7 @@ public Number? mp_set_from_string (string str, int default_base = 10)
/* Convert integer part */
var z = new Number.integer (0);
+
while (str.get_next_char (ref index, out c))
{
var i = char_val (c, number_base);
@@ -3233,73 +1600,6 @@ private Number? set_from_sexagesimal (string str)
return null;
}
-/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
- * GREATEST COMMON DIVISOR OF K AND L.
- * SAVE INPUT PARAMETERS IN LOCAL VARIABLES
- */
-public void mp_gcd (ref int64 k, ref int64 l)
-{
- var i = k.abs ();
- var j = l.abs ();
- if (j == 0)
- {
- /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */
- k = 1;
- l = 0;
- if (i == 0)
- k = 0;
- return;
- }
-
- /* EUCLIDEAN ALGORITHM LOOP */
- do
- {
- i %= j;
- if (i == 0)
- {
- k = k / j;
- l = l / j;
- return;
- }
- j %= i;
- } while (j != 0);
-
- /* HERE J IS THE GCD OF K AND L */
- k = k / i;
- l = l / i;
-}
-
-// FIXME: Is r.re_fraction large enough? It seems to be in practise but it may be T+4 instead of T
-// FIXME: There is some sort of stack corruption/use of unitialised variables here. Some functions are
-// using stack variables as x otherwise there are corruption errors. e.g. "Cos (45) - 1/Sqrt (2) = -0"
-// (try in scientific mode)
-public void mp_normalize (ref Number x)
-{
- int start_index;
-
- /* Find first non-zero digit */
- for (start_index = 0; start_index < SIZE && x.re_fraction[start_index] == 0; start_index++);
-
- /* Mark as zero */
- if (start_index >= SIZE)
- {
- x.re_sign = 0;
- x.re_exponent = 0;
- return;
- }
-
- /* Shift left so first digit is non-zero */
- if (start_index > 0)
- {
- x.re_exponent -= start_index;
- var i = 0;
- for (; (i + start_index) < SIZE; i++)
- x.re_fraction[i] = x.re_fraction[i + start_index];
- for (; i < SIZE; i++)
- x.re_fraction[i] = 0;
- }
-}
-
/* Returns true if x is cannot be represented in a binary word of length 'wordlen' */
public bool mp_is_overflow (Number x, int wordlen)
{
diff --git a/src/serializer.vala b/src/serializer.vala
index 3568b33..7cb78f2 100644
--- a/src/serializer.vala
+++ b/src/serializer.vala
@@ -32,6 +32,9 @@ public class Serializer : Object
private unichar tsep; /* Locale specific thousands separator. */
private int tsep_count; /* Number of digits between separator. */
+ /* is set when an error (for example precision error while converting) occurs */
+ public string? error { get; set; default = null; }
+
public Serializer (DisplayFormat format, int number_base, int trailing_digits)
{
var radix_string = Posix.nl_langinfo (Posix.NLItem.RADIXCHAR);
@@ -318,7 +321,18 @@ public class Serializer : Object
var t3 = temp.subtract (t2);
var d = t3.to_integer ();
- string.prepend_c (d < 16 ? digits[d] : '?');
+
+ if (d < 16 && d >= 0)
+ {
+ string.prepend_c (digits[d]);
+ }
+ else
+ {
+ string.prepend_c ('?');
+ error = _("Precision error");
+ string.assign ("0");
+ break;
+ }
n_digits++;
temp = t;
diff --git a/src/test-equation.vala b/src/test-equation.vala
index 48044ca..423d309 100644
--- a/src/test-equation.vala
+++ b/src/test-equation.vala
@@ -1,13 +1,13 @@
/*
* Copyright (C) 2008-2012 Robert Ancell.
- *
+ *
* This program is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License as published by the Free Software
* Foundation, either version 2 of the License, or (at your option) any later
* version. See http://www.gnu.org/copyleft/gpl.html the full text of the
* license.
*/
-
+
private int number_base = 10;
private int wordlen = 32;
private AngleUnit angle_units = AngleUnit.DEGREES;
@@ -20,7 +20,7 @@ private int pass_count = 0;
private string error_code_to_string (ErrorCode error)
{
if (error == ErrorCode.MP)
- return "ErrorCode.MP(\"%s\")".printf (mp_get_error ());
+ return "ErrorCode.MP(\"%s\")".printf (Number.error);
return mp_error_code_to_string (error);
}
@@ -115,7 +115,7 @@ private class TestEquation : Equation
{
if (!enable_conversions)
return null;
-
+
return UnitManager.get_default ().convert_by_symbol (x, x_units, z_units);
}
}
@@ -133,19 +133,19 @@ private void test_conversions ()
test ("100 gradians in degrees", "90", 0);
/* Length */
- test ("1 meter in mm", "1000", 0);
+ test ("1 meter in mm", "1000", 0);
test ("1m in mm", "1000", 0);
test ("1 inch in cm", "2.54", 0);
/* Area */
test ("1m² in mm²", "1000000", 0);
-
+
/* Volume */
- test ("1m³ in mm³", "1000000000", 0);
+ test ("1m³ in mm³", "1000000000", 0);
/* Weight */
test ("1 kg in pounds", "2.204622622", 0);
-
+
/* Duration */
test ("1 minute in seconds", "60", 0);
test ("1s in ms", "1000", 0);
@@ -223,7 +223,7 @@ private void test_equations ()
test ("١٢٣٤٥٦٧٨٩٠", "1234567890", 0);
test ("۱۲۳۴۵۶۷۸۹۰", "1234567890", 0);
-/*
+/*
//test ("2A", "2000000000000000", 0);
test ("2T", "2000000000000", 0);
test ("2G", "2000000000", 0);
@@ -263,7 +263,7 @@ private void test_equations ()
test ("2y", "6", 0);
test ("y2", "", ErrorCode.INVALID);
test ("y 2", "", ErrorCode.INVALID);
- test ("2z", "", ErrorCode.UNKNOWN_VARIABLE);
+ test ("2z", "", ErrorCode.UNKNOWN_VARIABLE);
test ("z2", "", ErrorCode.UNKNOWN_VARIABLE);
test ("z 2", "", ErrorCode.UNKNOWN_VARIABLE);
test ("z(2)", "", ErrorCode.UNKNOWN_FUNCTION);
@@ -414,7 +414,7 @@ private void test_equations ()
test ("-21 mod 9", "6", 0);
test ("-21 mod -9", "−3", 0);
- test ("sgn 0", "0", 0);
+ test ("sgn 0", "0", 0);
test ("sgn 3", "1", 0);
test ("sgn −3", "−1", 0);
test ("⌊3⌋", "3", 0);
@@ -530,16 +530,16 @@ private void test_equations ()
test ("i", "i", 0);
test ("−i", "−i", 0);
test ("2i", "2i", 0);
- test ("1+i", "1+i", 0);
+ test ("1+i", "1+i", 0);
test ("i+1", "1+i", 0);
- test ("1−i", "1−i", 0);
+ test ("1−i", "1−i", 0);
test ("i−1", "−1+i", 0);
test ("i×i", "−1", 0);
test ("i÷i", "1", 0);
test ("1÷i", "−i", 0);
test ("|i|", "1", 0);
test ("|3+4i|", "5", 0);
- test ("arg 0", "", ErrorCode.MP);
+ test ("arg 0", "", ErrorCode.MP);
test ("arg 1", "0", 0);
test ("arg (1+i)", "45", 0);
test ("arg i", "90", 0);
@@ -565,9 +565,9 @@ private void test_equations ()
test ("1 and 1", "1", 0);
test ("3 and 5", "1", 0);
- test ("0 or 0", "0", 0);
+ test ("0 or 0", "0", 0);
test ("1 or 0", "1", 0);
- test ("0 or 1", "1", 0);
+ test ("0 or 1", "1", 0);
test ("1 or 1", "1", 0);
test ("3 or 5", "7", 0);
diff --git a/src/test-number.vala b/src/test-number.vala
index b94535b..cb346a4 100644
--- a/src/test-number.vala
+++ b/src/test-number.vala
@@ -72,21 +72,6 @@ private void test_fraction ()
pass ();
}
-private void test_float ()
-{
- for (var a = -10.0f; a <= 10.0f; a += 0.5f)
- {
- var z = new Number.float (a);
- if (z.to_float () != a)
- {
- fail ("Number.float (%f).to_float () -> %f, expected %f".printf (a, z.to_float (), a));
- return;
- }
- }
-
- pass ();
-}
-
private void test_double ()
{
for (var a = -10.0; a <= 10.0; a += 0.5)
@@ -1089,7 +1074,6 @@ static int main (string[] args)
test_integer ();
test_unsigned_integer ();
test_fraction ();
- test_float ();
test_double ();
test_complex ();
test_polar ();
diff --git a/src/unit.vala b/src/unit.vala
index c6eb866..a051588 100644
--- a/src/unit.vala
+++ b/src/unit.vala
@@ -13,7 +13,7 @@ private UnitManager? default_unit_manager = null;
public class UnitManager : Object
{
private List<UnitCategory> categories;
-
+
public UnitManager ()
{
categories = new List<UnitCategory> ();
@@ -313,7 +313,7 @@ private class UnitSolveEquation : Equation
base (function);
this.x = x;
}
-
+
public override bool variable_is_defined (string name)
{
return true;
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
