1/*
2 * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
3 * Copyright (C) 2008-2012 Robert Ancell
4 *
5 * This program is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU General Public License as published by the Free Software
7 * Foundation, either version 3 of the License, or (at your option) any later
8 * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
9 * license.
10 */
11
12/*  This maths library is based on the MP multi-precision floating-point
13 *  arithmetic package originally written in FORTRAN by Richard Brent,
14 *  Computer Centre, Australian National University in the 1970's.
15 *
16 *  It has been converted from FORTRAN into C using the freely available
17 *  f2c translator, available via netlib on research.att.com.
18 *
19 *  The subsequently converted C code has then been tidied up, mainly to
20 *  remove any dependencies on the libI77 and libF77 support libraries.
21 *
22 *  FOR A GENERAL DESCRIPTION OF THE PHILOSOPHY AND DESIGN OF MP,
23 *  SEE - R. P. BRENT, A FORTRAN MULTIPLE-PRECISION ARITHMETIC
24 *  SOME ADDITIONAL DETAILS ARE GIVEN IN THE SAME ISSUE, 71-81.
25 *  FOR DETAILS OF THE IMPLEMENTATION, CALLING SEQUENCES ETC. SEE
26 *  THE MP USERS GUIDE.
27 */
28
29using MPC;
30
31private delegate int BitwiseFunc (int v1, int v2);
32
33public enum AngleUnit
34{
35    RADIANS,
36    DEGREES,
37    GRADIANS
38}
39
40/* Object for a high precision floating point number representation */
41public class Number : Object
42{
43    /* real and imaginary part of a Number */
44    private Complex num = Complex (precision);
45
46    public static MPFR.Precision precision { get; set; default = 1000; }
47
48    /* Stores the error msg if an error occurs during calculation. Otherwise should be null */
49    public static string? error { get; set; default = null; }
50
51    public Number.integer (int64 real, int64 imag = 0)
52    {
53        num.set_signed_integer ((long) real, (long) imag);
54    }
55
56    public Number.unsigned_integer (uint64 real, uint64 imag = 0)
57    {
58        num.set_unsigned_integer ((ulong) real, (ulong) imag);
59    }
60
61    public Number.fraction (int64 numerator, int64 denominator)
62    {
63        if (denominator < 0)
64        {
65            numerator = -numerator;
66            denominator = -denominator;
67        }
68
69        this.integer (numerator);
70        if (denominator != 1)
71        {
72            num.divide_unsigned_integer (num, (long) denominator);
73        }
74    }
75
76    /* Helper constructor. Creates new Number from already existing MPFR.Real. */
77    public Number.mpreal (MPFR.Real real, MPFR.Real? imag = null)
78    {
79        num.set_mpreal (real, imag);
80    }
81
82    public Number.double (double real, double imag = 0)
83    {
84        num.set_double (real, imag);
85    }
86
87    public Number.complex (Number r, Number i)
88    {
89        num.set_mpreal (r.num.get_real ().val, i.num.get_real ().val);
90    }
91
92    public Number.polar (Number r, Number theta, AngleUnit unit = AngleUnit.RADIANS)
93    {
94        var x = theta.cos (unit);
95        var y = theta.sin (unit);
96        this.complex (x.multiply (r), y.multiply (r));
97    }
98
99    public Number.eulers ()
100    {
101        num.get_real ().val.set_unsigned_integer (1);
102        /* e^1, since mpfr doesn't have a function to return e */
103        num.get_real ().val.exp (num.get_real ().val);
104        num.get_imag ().val.set_zero ();
105    }
106
107    public Number.i ()
108    {
109        num.set_signed_integer (0, 1);
110    }
111
112    public Number.pi ()
113    {
114        num.get_real ().val.const_pi ();
115        num.get_imag ().val.set_zero ();
116    }
117
118    public Number.tau ()
119    {
120        num.get_real ().val.const_tau ();
121        num.get_imag ().val.set_zero ();
122    }
123
124    /* Sets z to be a uniform random number in the range [0, 1] */
125    public Number.random ()
126    {
127        this.double (Random.next_double ());
128    }
129
130    public int64 to_integer ()
131    {
132        return num.get_real ().val.get_signed_integer ();
133    }
134
135    public uint64 to_unsigned_integer ()
136    {
137        return num.get_real ().val.get_unsigned_integer ();
138    }
139
140    public float to_float ()
141    {
142        return num.get_real ().val.get_float (MPFR.Round.NEAREST);
143    }
144
145    public double to_double ()
146    {
147        return num.get_real ().val.get_double (MPFR.Round.NEAREST);
148    }
149
150    /* Return true if the value is x == 0 */
151    public bool is_zero ()
152    {
153        return num.is_zero ();
154    }
155
156    /* Return true if x < 0 */
157    public bool is_negative ()
158    {
159        return num.get_real ().val.sgn () < 0;
160    }
161
162    /* Return true if x is integer */
163    public bool is_integer ()
164    {
165        if (is_complex ())
166            return false;
167
168        return num.get_real ().val.is_integer () != 0;
169    }
170
171    /* Return true if x is a positive integer */
172    public bool is_positive_integer ()
173    {
174        if (is_complex ())
175            return false;
176        else
177            return num.get_real ().val.sgn () >= 0 && is_integer ();
178    }
179
180    /* Return true if x is a natural number (an integer ≥ 0) */
181    public bool is_natural ()
182    {
183        if (is_complex ())
184            return false;
185        else
186            return num.get_real ().val.sgn () > 0 && is_integer ();
187    }
188
189    /* Return true if x has an imaginary component */
190    public bool is_complex ()
191    {
192        return !num.get_imag ().val.is_zero ();
193    }
194
195    /* Return error if overflow or underflow */
196    public static void check_flags ()
197    {
198        if (MPFR.mpfr_is_underflow () != 0)
199        {
200            /* Translators: Error displayed when underflow error occured */
201            error = _("Underflow error");
202        }
203        else if (MPFR.mpfr_is_overflow () != 0)
204        {
205            /* Translators: Error displayed when overflow error occured */
206            error = _("Overflow error");
207        }
208    }
209
210    /* Return true if x == y */
211    public bool equals (Number y)
212    {
213        return num.is_equal (y.num);
214    }
215
216    /* Returns:
217     *  0 if x == y
218     * <0 if x < y
219     * >0 if x > y
220     */
221    public int compare (Number y)
222    {
223        return num.get_real ().val.cmp (y.num.get_real ().val);
224    }
225
226    /* Sets z = sgn (x) */
227    public Number sgn ()
228    {
229        var z = new Number.integer (num.get_real ().val.sgn ());
230        return z;
231    }
232
233    /* Sets z = −x */
234    public Number invert_sign ()
235    {
236        var z = new Number ();
237        z.num.neg (num);
238        return z;
239    }
240
241    /* Sets z = |x| */
242    public Number abs ()
243    {
244      var z = new Number ();
245      z.num.get_imag ().val.set_zero ();
246      MPC.abs (z.num.get_real ().val, num);
247      return z;
248    }
249
250    /* Sets z = Arg (x) */
251    public Number arg (AngleUnit unit = AngleUnit.RADIANS)
252    {
253        if (is_zero ())
254        {
255            /* Translators: Error display when attempting to take argument of zero */
256            error = _("Argument not defined for zero");
257            return new Number.integer (0);
258        }
259        var z = new Number ();
260        z.num.get_imag ().val.set_zero ();
261        MPC.arg (z.num.get_real ().val, num);
262        mpc_from_radians (z.num, z.num, unit);
263        // MPC returns -π for the argument of negative real numbers if
264        // their imaginary part is -0 (which it is in the numbers
265        // created by test-equation), we want +π for all real negative
266        // numbers
267        if (!is_complex () && is_negative ())
268            z.num.get_real ().val.abs (z.num.get_real ().val);
269
270        return z;
271    }
272
273    /* Sets z = ‾̅x */
274    public Number conjugate ()
275    {
276        var z = new Number ();
277        z.num.conj (num);
278        return z;
279    }
280
281    /* Sets z = Re (x) */
282    public Number real_component ()
283    {
284        var z = new Number ();
285        z.num.set_mpreal (num.get_real ().val);
286        return z;
287    }
288
289    /* Sets z = Im (x) */
290    public Number imaginary_component ()
291    {
292        /* Copy imaginary component to real component */
293        var z = new Number ();
294        z.num.set_mpreal (num.get_imag ().val);
295        return z;
296    }
297
298    public Number integer_component ()
299    {
300        var z = new Number ();
301        z.num.get_imag ().val.set_zero ();
302        z.num.get_real ().val.trunc (num.get_real ().val);
303        return z;
304    }
305
306    /* Sets z = x mod 1 */
307    public Number fractional_component ()
308    {
309        var z = new Number ();
310        z.num.get_imag ().val.set_zero ();
311        z.num.get_real ().val.frac (num.get_real ().val);
312        return z;
313    }
314
315    /* Sets z = {x} */
316    public Number fractional_part ()
317    {
318        return subtract (floor ());
319    }
320
321    /* Sets z = ⌊x⌋ */
322    public Number floor ()
323    {
324        var z = new Number ();
325        z.num.get_imag ().val.set_zero ();
326        z.num.get_real ().val.floor (num.get_real ().val);
327        return z;
328    }
329
330    /* Sets z = ⌈x⌉ */
331    public Number ceiling ()
332    {
333        var z = new Number ();
334        z.num.get_imag ().val.set_zero ();
335        z.num.get_real ().val.ceil (num.get_real ().val);
336        return z;
337    }
338
339    /* Sets z = [x] */
340    public Number round ()
341    {
342        var z = new Number ();
343        z.num.get_imag ().val.set_zero ();
344        z.num.get_real ().val.round (num.get_real ().val);
345        return z;
346    }
347
348    /* Sets z = 1 ÷ x */
349    public Number reciprocal ()
350    {
351        var z = new Number ();
352        z.num.set_signed_integer (1);
353        z.num.mpreal_divide (z.num.get_real ().val, num);
354        return z;
355    }
356
357    /* Sets z = e^x */
358    public Number epowy ()
359    {
360        var z = new Number ();
361        z.num.exp (num);
362        return z;
363    }
364
365    /* Sets z = x^y */
366    public Number xpowy (Number y)
367    {
368        /* 0^-n invalid */
369        if (is_zero () && y.is_negative ())
370        {
371            /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */
372            error = _("The power of zero is undefined for a negative exponent");
373            return new Number.integer (0);
374        }
375
376        /* 0^0 is indeterminate */
377        if (is_zero () && y.is_zero ())
378        {
379            /* Translators: Error displayed when attempted to raise 0 to power of zero */
380            error = _("Zero raised to zero is undefined");
381            return new Number.integer (0);
382        }
383        if (!is_complex () && !y.is_complex () && !y.is_integer ())
384        {
385            var reciprocal = y.reciprocal ();
386            if (reciprocal.is_integer ())
387                return root (reciprocal.to_integer ());
388        }
389
390        var z = new Number ();
391        z.num.power (num, y.num);
392        return z;
393    }
394
395    /* Sets z = x^y */
396    public Number xpowy_integer (int64 n)
397    {
398        /* 0^-n invalid */
399        if (is_zero () && n < 0)
400        {
401            /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */
402            error = _("The power of zero is undefined for a negative exponent");
403            return new Number.integer (0);
404        }
405
406        /* 0^0 is indeterminate */
407        if (is_zero () && n == 0)
408        {
409            /* Translators: Error displayed when attempted to raise 0 to power of zero */
410            error = _("Zero raised to zero is undefined");
411            return new Number.integer (0);
412        }
413        var z = new Number ();
414        z.num.power_integer (num, (long) n);
415        return z;
416    }
417
418    /* Sets z = n√x */
419    public Number root (int64 n)
420    {
421        uint64 p;
422        var z = new Number ();
423        if (n < 0)
424        {
425            z.num.unsigned_integer_divide (1, num);
426            if (n == int64.MIN)
427                p = (uint64) int64.MAX + 1;
428            else
429                p = -n;
430        } else if (n > 0) {
431            z.num.@set (num);
432            p = n;
433        } else {
434            error = _("The zeroth root of a number is undefined");
435            return new Number.integer (0);
436        }
437
438        if (!is_complex () && (!is_negative () || (p & 1) == 1))
439        {
440            z.num.get_real ().val.root (z.num.get_real ().val, (ulong) p);
441            z.num.get_imag().val.set_zero();
442        } else {
443            var tmp = MPFR.Real (precision);
444            tmp.set_unsigned_integer ((ulong) p);
445            tmp.unsigned_integer_divide (1, tmp);
446            z.num.power_mpreal (z.num, tmp);
447        }
448        return z;
449    }
450
451    /* Sets z = √x */
452    public Number sqrt ()
453    {
454        return root(2);
455    }
456
457    /* Sets z = ln x */
458    public Number ln ()
459    {
460        /* ln (0) undefined */
461        if (is_zero ())
462        {
463            /* Translators: Error displayed when attempting to take logarithm of zero */
464            error = _("Logarithm of zero is undefined");
465            return new Number.integer (0);
466        }
467
468        /* ln (-x) complex */
469        /* FIXME: Make complex numbers optional */
470        /*if (is_negative ())
471        {
472            // Translators: Error displayed attempted to take logarithm of negative value
473            mperr (_("Logarithm of negative values is undefined"));
474            return new Number.integer (0);
475        }*/
476
477        var z = new Number ();
478        z.num.log (num);
479        // MPC returns -π for the imaginary part of the log of
480        // negative real numbers if their imaginary part is -0 (which
481        // it is in the numbers created by test-equation), we want +π
482        if (!is_complex () && is_negative ())
483            z.num.get_imag ().val.abs (z.num.get_imag ().val);
484
485        return z;
486    }
487
488    /* Sets z = log_n x */
489    public Number logarithm (int64 n)
490    {
491        /* log (0) undefined */
492        if (is_zero ())
493        {
494            /* Translators: Error displayed when attempting to take logarithm of zero */
495            error = _("Logarithm of zero is undefined");
496            return new Number.integer (0);
497        }
498
499        /* logn (x) = ln (x) / ln (n) */
500        var t1 = new Number.integer (n);
501        return ln ().divide (t1.ln ());
502    }
503
504    /* Sets z = x! */
505    public Number factorial ()
506    {
507        /* 0! == 1 */
508        if (is_zero ())
509            return new Number.integer (1);
510        if (!is_natural ())
511        {
512
513             /* Factorial Not defined for Complex or for negative numbers */
514            if (is_negative () || is_complex ())
515            {
516                /* Translators: Error displayed when attempted take the factorial of a negative or complex number */
517                error = _("Factorial is only defined for non-negative real numbers");
518                return new Number.integer (0);
519            }
520
521            var tmp = add (new Number.integer (1));
522            var tmp2 = MPFR.Real (precision);
523
524            /* Factorial(x) = Gamma(x+1) - This is the formula used to calculate Factorial.*/
525            tmp2.gamma (tmp.num.get_real ().val);
526
527            return new Number.mpreal (tmp2);
528        }
529
530        /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
531        var value = to_integer ();
532        var z = this;
533        for (var i = 2; i < value; i++)
534            z = z.multiply_integer (i);
535
536        return z;
537    }
538
539    /* Sets z = x + y */
540    public Number add (Number y)
541    {
542        var z = new Number ();
543        z.num.add (num, y.num);
544        return z;
545    }
546
547    /* Sets z = x − y */
548    public Number subtract (Number y)
549    {
550        var z = new Number ();
551        z.num.subtract (num, y.num);
552        return z;
553    }
554
555    /* Sets z = x × y */
556    public Number multiply (Number y)
557    {
558        var z = new Number ();
559        z.num.multiply (num, y.num);
560        return z;
561    }
562
563    /* Sets z = x × y */
564    public Number multiply_integer (int64 y)
565    {
566        var z = new Number ();
567        z.num.multiply_signed_integer (num, (long) y);
568        return z;
569    }
570
571    /* Sets z = x ÷ y */
572    public Number divide (Number y)
573    {
574        if (y.is_zero ())
575        {
576            /* Translators: Error displayed attempted to divide by zero */
577            error = _("Division by zero is undefined");
578            return new Number.integer (0);
579        }
580
581        var z = new Number ();
582        z.num.divide (num, y.num);
583        return z;
584    }
585
586    /* Sets z = x ÷ y */
587    public Number divide_integer (int64 y)
588    {
589        return divide (new Number.integer (y));
590    }
591
592    /* Sets z = x mod y */
593    public Number modulus_divide (Number y)
594    {
595        if (!is_integer () || !y.is_integer ())
596        {
597            /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
598            error = _("Modulus division is only defined for integers");
599            return new Number.integer (0);
600        }
601
602        var t1 = divide (y).floor ();
603        var t2 = t1.multiply (y);
604        var z = subtract (t2);
605
606        t1 = new Number.integer (0);
607        if ((y.compare (t1) < 0 && z.compare (t1) > 0) || (y.compare (t1) > 0 && z.compare (t1) < 0))
608            z = z.add (y);
609
610        return z;
611    }
612
613    /* Sets z = x ^ y mod p */
614    public Number modular_exponentiation (Number exp, Number mod)
615    {
616        var base_value = copy ();
617        if (exp.is_negative ())
618            base_value = base_value.reciprocal ();
619        var exp_value = exp.abs ();
620        var ans = new Number.integer (1);
621        var two = new Number.integer (2);
622        while (!exp_value.is_zero ())
623        {
624            bool is_even = exp_value.modulus_divide (two).is_zero ();
625            if (!is_even)
626            {
627                ans = ans.multiply (base_value);
628                ans = ans.modulus_divide (mod);
629            }
630            base_value = base_value.multiply (base_value);
631            base_value = base_value.modulus_divide (mod);
632            exp_value = exp_value.divide_integer (2).floor ();
633        }
634        return ans.modulus_divide (mod);
635    }
636
637    /* Sets z = sin x */
638    public Number sin (AngleUnit unit = AngleUnit.RADIANS)
639    {
640        var z = new Number ();
641        if (is_complex ())
642          z.num.@set (num);
643        else
644          mpc_to_radians (z.num, num, unit);
645        z.num.sin (z.num);
646        return z;
647    }
648
649    /* Sets z = cos x */
650    public Number cos (AngleUnit unit = AngleUnit.RADIANS)
651    {
652        var z = new Number ();
653        if (is_complex ())
654          z.num.@set (num);
655        else
656          mpc_to_radians (z.num, num, unit);
657        z.num.cos (z.num);
658        return z;
659    }
660
661    /* Sets z = tan x */
662    public Number tan (AngleUnit unit = AngleUnit.RADIANS)
663    {
664        /* Check for undefined values */
665        var x_radians = to_radians (unit);
666        var check = x_radians.subtract (new Number.pi ().divide_integer (2)).divide (new Number.pi ());
667
668        if (check.is_integer ())
669        {
670            /* Translators: Error displayed when tangent value is undefined */
671            error = _("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)");
672            return new Number.integer (0);
673        }
674
675        var z = new Number ();
676        if (is_complex ())
677          z.num.@set (num);
678        else
679          mpc_to_radians (z.num, num, unit);
680        z.num.tan (z.num);
681        return z;
682    }
683
684    /* Sets z = sin⁻¹ x */
685    public Number asin (AngleUnit unit = AngleUnit.RADIANS)
686    {
687        if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0)
688        {
689            /* Translators: Error displayed when inverse sine value is undefined */
690            error = _("Inverse sine is undefined for values outside [-1, 1]");
691            return new Number.integer (0);
692        }
693
694        var z = new Number ();
695        z.num.asin (num);
696        if (!z.is_complex ())
697          mpc_from_radians (z.num, z.num, unit);
698        return z;
699    }
700
701    /* Sets z = cos⁻¹ x */
702    public Number acos (AngleUnit unit = AngleUnit.RADIANS)
703    {
704        if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0)
705        {
706            /* Translators: Error displayed when inverse cosine value is undefined */
707            error = _("Inverse cosine is undefined for values outside [-1, 1]");
708            return new Number.integer (0);
709        }
710
711        var z = new Number ();
712        z.num.acos (num);
713        if (!z.is_complex ())
714          mpc_from_radians (z.num, z.num, unit);
715        return z;
716    }
717
718    /* Sets z = tan⁻¹ x */
719    public Number atan (AngleUnit unit = AngleUnit.RADIANS)
720    {
721        /* Check x != i and x != -i */
722        if (equals (new Number.integer (0,1)) || equals (new Number.integer(0,-1)))
723        {
724            /* Translators: Error displayed when trying to calculate undefined atan(i) or atan (-i) */
725            error = _("Arctangent function is undefined for values i and -i");
726            return new Number.integer (0);
727        }
728        var z = new Number ();
729        z.num.atan (num);
730        if (!z.is_complex ())
731          mpc_from_radians (z.num, z.num, unit);
732        return z;
733    }
734
735    /* Sets z = sinh x */
736    public Number sinh ()
737    {
738        var z = new Number ();
739        z.num.sinh (num);
740        return z;
741    }
742
743    /* Sets z = cosh x */
744    public Number cosh ()
745    {
746        var z = new Number ();
747        z.num.cosh (num);
748        return z;
749    }
750
751    /* Sets z = tanh x */
752    public Number tanh ()
753    {
754        var z = new Number ();
755        z.num.tanh (num);
756        return z;
757    }
758
759    /* Sets z = sinh⁻¹ x */
760    public Number asinh ()
761    {
762        var z = new Number ();
763        z.num.asinh (num);
764        return z;
765    }
766
767    /* Sets z = cosh⁻¹ x */
768    public Number acosh ()
769    {
770        /* Check x >= 1 */
771        var t = new Number.integer (1);
772        if (compare (t) < 0)
773        {
774            /* Translators: Error displayed when inverse hyperbolic cosine value is undefined */
775            error = _("Inverse hyperbolic cosine is undefined for values less than one");
776            return new Number.integer (0);
777        }
778
779        var z = new Number ();
780        z.num.acosh (num);
781        return z;
782    }
783
784    /* Sets z = tanh⁻¹ x */
785    public Number atanh ()
786    {
787        /* Check -1 <= x <= 1 */
788        if (compare (new Number.integer (1)) >= 0 || compare (new Number.integer (-1)) <= 0)
789        {
790            /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */
791            error = _("Inverse hyperbolic tangent is undefined for values outside [-1, 1]");
792            return new Number.integer (0);
793        }
794
795        var z = new Number ();
796        z.num.atanh (num);
797        return z;
798    }
799
800    /* Sets z = boolean AND for each bit in x and z */
801    public Number and (Number y)
802    {
803        if (!
804        is_positive_integer () || !y.is_positive_integer ())
805        {
806            /* Translators: Error displayed when boolean AND attempted on non-integer values */
807            error = _("Boolean AND is only defined for positive integers");
808        }
809
810        return bitwise (y, (v1, v2) => { return v1 & v2; }, 0);
811    }
812
813    /* Sets z = boolean OR for each bit in x and z */
814    public Number or (Number y)
815    {
816        if (!is_positive_integer () || !y.is_positive_integer ())
817        {
818            /* Translators: Error displayed when boolean OR attempted on non-integer values */
819            error = _("Boolean OR is only defined for positive integers");
820        }
821
822        return bitwise (y, (v1, v2) => { return v1 | v2; }, 0);
823    }
824
825    /* Sets z = boolean XOR for each bit in x and z */
826    public Number xor (Number y)
827    {
828        if (!is_positive_integer () || !y.is_positive_integer ())
829        {
830            /* Translators: Error displayed when boolean XOR attempted on non-integer values */
831            error = _("Boolean XOR is only defined for positive integers");
832        }
833
834        return bitwise (y, (v1, v2) => { return v1 ^ v2; }, 0);
835    }
836
837    /* Sets z = boolean NOT for each bit in x and z for word of length 'wordlen' */
838    public Number not (int wordlen)
839    {
840        if (!is_positive_integer ())
841        {
842            /* Translators: Error displayed when boolean XOR attempted on non-integer values */
843            error = _("Boolean NOT is only defined for positive integers");
844        }
845
846        return bitwise (new Number.integer (0), (v1, v2) => { return v1 ^ 0xF; }, wordlen);
847    }
848
849    /* Sets z = x masked to 'wordlen' bits */
850    public Number mask (Number x, int wordlen)
851    {
852        /* Convert to a hexadecimal string and use last characters */
853        var text = x.to_hex_string ();
854        var len = text.length;
855        var offset = wordlen / 4;
856        offset = len > offset ? (int) len - offset: 0;
857        return mp_set_from_string (text.substring (offset), 16);
858    }
859
860    /* Sets z = x shifted by 'count' bits.  Positive shift increases the value, negative decreases */
861    public Number shift (int64 count)
862    {
863        if (!is_integer ())
864        {
865            /* Translators: Error displayed when bit shift attempted on non-integer values */
866            error = _("Shift is only possible on integer values");
867            return new Number.integer (0);
868        }
869
870        if (count >= 0)
871        {
872            var multiplier = 1;
873            for (var i = 0; i < count; i++)
874                multiplier *= 2;
875            return multiply_integer (multiplier);
876        }
877        else
878        {
879            var multiplier = 1;
880            for (var i = 0; i < -count; i++)
881                multiplier *= 2;
882            return divide_integer (multiplier).floor ();
883        }
884    }
885
886    /* Sets z to be the ones complement of x for word of length 'wordlen' */
887    public Number ones_complement (int wordlen)
888    {
889        return bitwise (new Number.integer (0), (v1, v2) => { return v1 ^ v2; }, wordlen).not (wordlen);
890    }
891
892    /* Sets z to be the twos complement of x for word of length 'wordlen' */
893    public Number twos_complement (int wordlen)
894    {
895        return ones_complement (wordlen).add (new Number.integer (1));
896    }
897
898    /* Returns a list of all prime factors in x as Numbers */
899    public List<Number?> factorize ()
900    {
901        var factors = new List<Number?> ();
902
903        var value = abs ();
904
905        if (value.is_zero ())
906        {
907            factors.append (value);
908            return factors;
909        }
910
911        if (value.equals (new Number.integer (1)))
912        {
913            factors.append (this);
914            return factors;
915        }
916
917        // if value < 2^64-1, call for factorize_uint64 function which deals in integers
918
919        uint64 num = 1;
920        num = num << 63;
921        num += (num - 1);
922        var int_max = new Number.unsigned_integer (num);
923
924        if (value.compare (int_max) <= 0)
925        {
926            var factors_int64 = factorize_uint64 (value.to_unsigned_integer ());
927            if (is_negative ())
928                factors_int64.data = factors_int64.data.invert_sign ();
929            return factors_int64;
930        }
931
932        var divisor = new Number.integer (2);
933        while (true)
934        {
935            var tmp = value.divide (divisor);
936            if (tmp.is_integer ())
937            {
938                value = tmp;
939                factors.append (divisor);
940            }
941            else
942                break;
943        }
944
945        divisor = new Number.integer (3);
946        var root = value.sqrt ();
947        while (divisor.compare (root) <= 0)
948        {
949            var tmp = value.divide (divisor);
950            if (tmp.is_integer ())
951            {
952                value = tmp;
953                root = value.sqrt ();
954                factors.append (divisor);
955            }
956            else
957            {
958                tmp = divisor.add (new Number.integer (2));
959                divisor = tmp;
960            }
961        }
962
963        if (value.compare (new Number.integer (1)) > 0)
964            factors.append (value);
965
966        if (is_negative ())
967            factors.data = factors.data.invert_sign ();
968
969        return factors;
970    }
971
972    public List<Number?> factorize_uint64 (uint64 n)
973    {
974        var factors = new List<Number?> ();
975        while (n % 2 == 0)
976        {
977            n /= 2;
978            factors.append (new Number.unsigned_integer (2));
979        }
980
981        for (uint64 divisor = 3; divisor <= n / divisor; divisor += 2)
982        {
983            while (n % divisor == 0)
984            {
985                n /= divisor;
986                factors.append (new Number.unsigned_integer (divisor));
987            }
988        }
989
990        if (n > 1)
991            factors.append (new Number.unsigned_integer (n));
992        return factors;
993    }
994
995    private Number copy ()
996    {
997        var z = new Number ();
998        z.num.@set (num);
999        return z;
1000    }
1001
1002    private static void mpc_from_radians (Complex res, Complex op, AngleUnit unit)
1003    {
1004        int i;
1005
1006        switch (unit)
1007        {
1008            default:
1009            case AngleUnit.RADIANS:
1010                if (res != op)
1011                    res.@set (op);
1012                return;
1013
1014            case AngleUnit.DEGREES:
1015                i = 180;
1016                break;
1017
1018            case AngleUnit.GRADIANS:
1019                i=200;
1020                break;
1021
1022        }
1023        var scale = MPFR.Real (precision);
1024        scale.const_pi ();
1025        scale.signed_integer_divide (i, scale);
1026        res.multiply_mpreal (op, scale);
1027    }
1028
1029    private static void mpc_to_radians (Complex res, Complex op, AngleUnit unit)
1030    {
1031        int i;
1032
1033        switch (unit)
1034        {
1035            default:
1036            case AngleUnit.RADIANS:
1037                if (res != op)
1038                    res.@set (op);
1039                return;
1040
1041            case AngleUnit.DEGREES:
1042                i = 180;
1043                break;
1044
1045            case AngleUnit.GRADIANS:
1046                i=200;
1047                break;
1048        }
1049        var scale = MPFR.Real (precision);
1050        scale.const_pi ();
1051        scale.divide_signed_integer (scale, i);
1052        res.multiply_mpreal (op, scale);
1053    }
1054
1055    /* Convert x to radians */
1056    private Number to_radians (AngleUnit unit)
1057    {
1058        var z = new Number ();
1059        mpc_to_radians (z.num, num, unit);
1060        return z;
1061    }
1062
1063    private Number bitwise (Number y, BitwiseFunc bitwise_operator, int wordlen)
1064    {
1065        var text1 = to_hex_string ();
1066        var text2 = y.to_hex_string ();
1067        var offset1 = text1.length - 1;
1068        var offset2 = text2.length - 1;
1069        var offset_out = wordlen / 4 - 1;
1070        if (offset_out <= 0)
1071            offset_out = offset1 > offset2 ? offset1 : offset2;
1072        if (offset_out > 0 && (offset_out < offset1 || offset_out < offset2))
1073        {
1074            error = ("Overflow. Try a bigger word size");
1075            return new Number.integer (0);
1076        }
1077
1078        var text_out = new char[offset_out + 2];
1079
1080        /* Perform bitwise operator on each character from right to left */
1081        for (text_out[offset_out+1] = '\0'; offset_out >= 0; offset_out--)
1082        {
1083            int v1 = 0, v2 = 0;
1084            const char digits[] = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F' };
1085
1086            if (offset1 >= 0)
1087            {
1088                v1 = hex_to_int (text1[offset1]);
1089                offset1--;
1090            }
1091            if (offset2 >= 0)
1092            {
1093                v2 = hex_to_int (text2[offset2]);
1094                offset2--;
1095            }
1096            text_out[offset_out] = digits[bitwise_operator (v1, v2)];
1097        }
1098
1099        return mp_set_from_string ((string) text_out, 16);
1100    }
1101
1102    private int hex_to_int (char digit)
1103    {
1104        if (digit >= '0' && digit <= '9')
1105            return digit - '0';
1106        if (digit >= 'A' && digit <= 'F')
1107            return digit - 'A' + 10;
1108        if (digit >= 'a' && digit <= 'f')
1109            return digit - 'a' + 10;
1110        return 0;
1111    }
1112
1113    private string to_hex_string ()
1114    {
1115        var serializer = new Serializer (DisplayFormat.FIXED, 16, 0);
1116        return serializer.to_string (this);
1117    }
1118}
1119
1120private static int parse_literal_prefix (string str, ref int prefix_len)
1121{
1122    var new_base = 0;
1123
1124    if (str.length < 3 || str[0] != '0')
1125        return new_base;
1126
1127    var prefix = str[1].tolower ();
1128
1129    if (prefix == 'b')
1130        new_base = 2;
1131    else if (prefix == 'o')
1132        new_base = 8;
1133    else if (prefix == 'x')
1134        new_base = 16;
1135
1136    if (new_base != 0)
1137        prefix_len = 2;
1138
1139    return new_base;
1140}
1141
1142// FIXME: Should all be in the class
1143
1144// FIXME: Re-add overflow and underflow detection
1145
1146/* Sets z from a string representation in 'text'. */
1147public Number? mp_set_from_string (string str, int default_base = 10)
1148{
1149    if (str.index_of_char ('°') >= 0)
1150        return set_from_sexagesimal (str);
1151
1152    /* Find the base */
1153    const unichar base_digits[] = {'₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'};
1154    var index = 0;
1155    var base_prefix = 0;
1156    unichar c;
1157    while (str.get_next_char (ref index, out c));
1158    var end = index;
1159    var number_base = 0;
1160    var literal_base = 0;
1161    var base_multiplier = 1;
1162    while (str.get_prev_char (ref index, out c))
1163    {
1164        var value = -1;
1165        for (var i = 0; i < base_digits.length; i++)
1166        {
1167            if (c == base_digits[i])
1168            {
1169                value = i;
1170                break;
1171            }
1172        }
1173        if (value < 0)
1174            break;
1175
1176        end = index;
1177        number_base += value * base_multiplier;
1178        base_multiplier *= 10;
1179    }
1180
1181    literal_base = parse_literal_prefix (str, ref base_prefix);
1182
1183    if (number_base != 0 && literal_base != 0 && literal_base != number_base)
1184        return null;
1185
1186    if (number_base == 0)
1187        number_base = (literal_base != 0) ? literal_base : default_base;
1188
1189    /* Check if this has a sign */
1190    var negate = false;
1191    index = base_prefix;
1192    str.get_next_char (ref index, out c);
1193    if (c == '+')
1194        negate = false;
1195    else if (c == '-' || c == '−')
1196        negate = true;
1197    else
1198        str.get_prev_char (ref index, out c);
1199
1200    /* Convert integer part */
1201    var z = new Number.integer (0);
1202
1203    while (str.get_next_char (ref index, out c))
1204    {
1205        var i = char_val (c, number_base);
1206        if (i > number_base)
1207            return null;
1208        if (i < 0)
1209        {
1210            str.get_prev_char (ref index, out c);
1211            break;
1212        }
1213
1214        z = z.multiply_integer (number_base).add (new Number.integer (i));
1215    }
1216
1217    /* Look for fraction characters, e.g. ⅚ */
1218    const unichar fractions[] = {'½', '⅓', '⅔', '¼', '¾', '⅕', '⅖', '⅗', '⅘', '⅙', '⅚', '⅛', '⅜', '⅝', '⅞'};
1219    const int numerators[]    = { 1,   1,   2,   1,   3,   1,   2,   3,   4,   1,   5,   1,   3,   5,   7};
1220    const int denominators[]  = { 2,   3,   3,   4,   4,   5,   5,   5,   5,   6,   6,   8,   8,   8,   8};
1221    var has_fraction = false;
1222    if (str.get_next_char (ref index, out c))
1223    {
1224        for (var i = 0; i < fractions.length; i++)
1225        {
1226            if (c == fractions[i])
1227            {
1228                var fraction = new Number.fraction (numerators[i], denominators[i]);
1229                z = z.add (fraction);
1230
1231                /* Must end with fraction */
1232                if (!str.get_next_char (ref index, out c))
1233                    return z;
1234                else
1235                    return null;
1236            }
1237        }
1238
1239        /* Check for decimal point */
1240        if (c == '.')
1241            has_fraction = true;
1242        else
1243            str.get_prev_char (ref index, out c);
1244    }
1245
1246    /* Convert fractional part */
1247    if (has_fraction)
1248    {
1249        var numerator = new Number.integer (0);
1250        var denominator = new Number.integer (1);
1251
1252        while (str.get_next_char (ref index, out c))
1253        {
1254            var i = char_val (c, number_base);
1255            if (i < 0)
1256            {
1257                str.get_prev_char (ref index, out c);
1258                break;
1259            }
1260
1261            denominator = denominator.multiply_integer (number_base);
1262            numerator = numerator.multiply_integer (number_base);
1263            numerator = numerator.add (new Number.integer (i));
1264        }
1265
1266        numerator = numerator.divide (denominator);
1267        z = z.add (numerator);
1268    }
1269
1270    if (index != end)
1271        return null;
1272
1273    if (negate)
1274        z = z.invert_sign ();
1275
1276    return z;
1277}
1278
1279private int char_val (unichar c, int number_base)
1280{
1281    if (!c.isxdigit ())
1282        return -1;
1283
1284    var value = c.xdigit_value ();
1285
1286    if (value >= number_base)
1287        return -1;
1288
1289    return value;
1290}
1291
1292private Number? set_from_sexagesimal (string str)
1293{
1294    var degree_index = str.index_of_char ('°');
1295    if (degree_index < 0)
1296        return null;
1297    var degrees = mp_set_from_string (str.substring (0, degree_index));
1298    if (degrees == null)
1299        return null;
1300    var minute_start = degree_index;
1301    unichar c;
1302    str.get_next_char (ref minute_start, out c);
1303
1304    if (str[minute_start] == '\0')
1305        return degrees;
1306    var minute_index = str.index_of_char ('\'', minute_start);
1307    if (minute_index < 0)
1308        return null;
1309    var minutes = mp_set_from_string (str.substring (minute_start, minute_index - minute_start));
1310    if (minutes == null)
1311        return null;
1312    degrees = degrees.add (minutes.divide_integer (60));
1313    var second_start = minute_index;
1314    str.get_next_char (ref second_start, out c);
1315
1316    if (str[second_start] == '\0')
1317        return degrees;
1318    var second_index = str.index_of_char ('"', second_start);
1319    if (second_index < 0)
1320        return null;
1321    var seconds = mp_set_from_string (str.substring (second_start, second_index - second_start));
1322    if (seconds == null)
1323        return null;
1324    degrees = degrees.add (seconds.divide_integer (3600));
1325    str.get_next_char (ref second_index, out c);
1326
1327    /* Skip over second marker and expect no more characters */
1328    if (str[second_index] == '\0')
1329        return degrees;
1330    else
1331        return null;
1332}
1333
1334/* Returns true if x is cannot be represented in a binary word of length 'wordlen' */
1335public bool mp_is_overflow (Number x, int wordlen)
1336{
1337    var t2 = new Number.integer (2).xpowy_integer (wordlen);
1338    return t2.compare (x) > 0;
1339}
1340