~vcs-imports/gawk/master

« back to all changes in this revision

Viewing changes to mpfr.c

  • Committer: Arnold D. Robbins
  • Date: 2013-01-11 12:30:08 UTC
  • Revision ID: git-v1:0bd22a097fcde68cf8586e8737ac7ad8f4286669
Make mpfr and, or, xor, accept >= 2 arguments.

Show diffs side-by-side

added added

removed removed

Lines of Context:
48
48
static mpfr_exp_t min_exp = MPFR_EMIN_DEFAULT;
49
49
static mpfr_exp_t max_exp = MPFR_EMAX_DEFAULT;
50
50
 
51
 
/* temporaries used in bit ops */
52
 
static NODE *_tz1;
53
 
static NODE *_tz2;
54
 
static mpz_t _mpz1;
55
 
static mpz_t _mpz2;
56
 
static mpz_ptr mpz1;
57
 
static mpz_ptr mpz2;
58
 
 
59
 
static NODE *get_bit_ops(const char *op);
60
 
#define free_bit_ops()  (DEREF(_tz1), DEREF(_tz2))
61
 
 
62
51
/* temporary MPFR floats used to hold converted GMP integer operands */
63
52
static mpfr_t _mpf_t1;
64
53
static mpfr_t _mpf_t2;
93
82
        mpz_init(MFNR);
94
83
        do_ieee_fmt = false;
95
84
 
96
 
        mpz_init(_mpz1);
97
 
        mpz_init(_mpz2);
98
85
        mpfr_init2(_mpf_t1, PRECISION_MIN);
99
86
        mpfr_init2(_mpf_t2, PRECISION_MIN);
100
87
        mpz_init(mpzval);
837
824
        return r;
838
825
}
839
826
 
840
 
 
841
 
/*
842
 
 * get_bit_ops --- get the numeric operands of a binary function.
843
 
 *      Returns a copy of the operand if either is inf or nan. Otherwise
844
 
 *      each operand is converted to an integer if necessary, and
845
 
 *      the results are placed in the variables mpz1 and mpz2.
846
 
 */
847
 
 
848
 
static NODE *
849
 
get_bit_ops(const char *op)
 
827
/* get_intval --- get the (converted) integral operand of a binary function. */
 
828
 
 
829
static mpz_ptr
 
830
get_intval(NODE *t1, int argnum, const char *op)
850
831
{
851
 
        _tz2 = POP_SCALAR();
852
 
        _tz1 = POP_SCALAR();
853
 
 
854
 
        if (do_lint) {
855
 
                if ((_tz1->flags & (NUMCUR|NUMBER)) == 0)
856
 
                        lintwarn(_("%s: received non-numeric first argument"), op);
857
 
                if ((_tz2->flags & (NUMCUR|NUMBER)) == 0)
858
 
                        lintwarn(_("%s: received non-numeric second argument"), op);
859
 
        }
860
 
 
861
 
        force_number(_tz1);
862
 
        force_number(_tz2);
863
 
 
864
 
        if (is_mpg_float(_tz1)) {
865
 
                mpfr_ptr left = _tz1->mpg_numbr;
 
832
        mpz_ptr pz;
 
833
 
 
834
        if (do_lint && (t1->flags & (NUMCUR|NUMBER)) == 0)
 
835
                lintwarn(_("%s: received non-numeric argument #%d"), op, argnum);
 
836
 
 
837
        (void) force_number(t1);
 
838
 
 
839
        if (is_mpg_float(t1)) {
 
840
                mpfr_ptr left = t1->mpg_numbr;
866
841
                if (! mpfr_number_p(left)) {
867
842
                        /* inf or NaN */
868
 
                        NODE *res;
869
 
                        res = mpg_float();
870
 
                        mpfr_set(res->mpg_numbr, _tz1->mpg_numbr, ROUND_MODE);
871
 
                        return res;
 
843
                        if (do_lint)
 
844
                                lintwarn("%s",
 
845
                mpg_fmt(_("%s: argument #%d has invalid value %Rg, using 0"),
 
846
                                        op, argnum, left)
 
847
                                );
 
848
 
 
849
                        emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
 
850
                        mpz_init(pz);
 
851
                        return pz;      /* should be freed */
872
852
                }
873
853
 
874
854
                if (do_lint) {
875
855
                        if (mpfr_sgn(left) < 0)
876
856
                                lintwarn("%s",
877
 
                        mpg_fmt(_("%s(%Rg, ..): negative values will give strange results"),
878
 
                                                op, left)
879
 
                                        );
 
857
                mpg_fmt(_("%s: argument #%d negative value %Rg will give strange results"),
 
858
                                        op, argnum, left)
 
859
                                );
 
860
 
880
861
                        if (! mpfr_integer_p(left))
881
862
                                lintwarn("%s",
882
 
                        mpg_fmt(_("%s(%Rg, ..): fractional values will be truncated"),
883
 
                                                op, left)
884
 
                                );
885
 
                }
886
 
                
887
 
                mpfr_get_z(_mpz1, left, MPFR_RNDZ);     /* float to integer conversion */
888
 
                mpz1 = _mpz1;
889
 
        } else {
890
 
                /* (_tz1->flags & MPZN) != 0 */ 
891
 
                mpz1 = _tz1->mpg_i;
892
 
                if (do_lint) {
893
 
                        if (mpz_sgn(mpz1) < 0)
894
 
                                lintwarn("%s",
895
 
                        mpg_fmt(_("%s(%Zd, ..): negative values will give strange results"),
896
 
                                                op, mpz1)
897
 
                                        );
898
 
                }
899
 
        }
900
 
 
901
 
        if (is_mpg_float(_tz2)) {
902
 
                mpfr_ptr right = _tz2->mpg_numbr;
903
 
                if (! mpfr_number_p(right)) {
904
 
                        /* inf or NaN */
905
 
                        NODE *res;
906
 
                        res = mpg_float();
907
 
                        mpfr_set(res->mpg_numbr, _tz2->mpg_numbr, ROUND_MODE);
908
 
                        return res;
909
 
                }
910
 
 
911
 
                if (do_lint) {
912
 
                        if (mpfr_sgn(right) < 0)
913
 
                                lintwarn("%s",
914
 
                        mpg_fmt(_("%s(.., %Rg): negative values will give strange results"),
915
 
                                                op, right)
916
 
                                        );
917
 
                        if (! mpfr_integer_p(right))
918
 
                                lintwarn("%s",
919
 
                        mpg_fmt(_("%s(.., %Rg): fractional values will be truncated"),
920
 
                                                op, right)
921
 
                                );
922
 
                }
923
 
 
924
 
                mpfr_get_z(_mpz2, right, MPFR_RNDZ);    /* float to integer conversion */
925
 
                mpz2 = _mpz2;
926
 
        } else {
927
 
                /* (_tz2->flags & MPZN) != 0 */ 
928
 
                mpz2 = _tz2->mpg_i;
929
 
                if (do_lint) {
930
 
                        if (mpz_sgn(mpz2) < 0)
931
 
                                lintwarn("%s",
932
 
                        mpg_fmt(_("%s(.., %Zd): negative values will give strange results"),
933
 
                                                op, mpz2)
934
 
                                        );
935
 
                }
936
 
        }
937
 
 
938
 
        return NULL;
939
 
}
 
863
                mpg_fmt(_("%s: argument #%d fractional value %Rg will be truncated"),
 
864
                                        op, argnum, left)
 
865
                                );
 
866
                }
 
867
 
 
868
                emalloc(pz, mpz_ptr, sizeof (mpz_t), "get_intval");
 
869
                mpz_init(pz);
 
870
                mpfr_get_z(pz, left, MPFR_RNDZ);        /* float to integer conversion */
 
871
                return pz;      /* should be freed */
 
872
        } 
 
873
        /* (t1->flags & MPZN) != 0 */ 
 
874
        pz = t1->mpg_i;
 
875
        if (do_lint) {
 
876
                if (mpz_sgn(pz) < 0)
 
877
                        lintwarn("%s",
 
878
        mpg_fmt(_("%s: argument #%d negative value %Zd will give strange results"),
 
879
                                        op, argnum, pz)
 
880
                                );
 
881
        }
 
882
        return pz;      /* must not be freed */
 
883
}
 
884
 
 
885
 
 
886
/* free_intval --- free the converted integer value returned by get_intval() */
 
887
 
 
888
static inline void
 
889
free_intval(NODE *t, mpz_ptr pz)
 
890
{
 
891
        if ((t->flags & MPZN) == 0) {
 
892
                mpz_clear(pz);
 
893
                efree(pz);
 
894
        }
 
895
}
 
896
 
940
897
 
941
898
/* do_mpfr_lshift --- perform a << operation */
942
899
 
943
900
NODE *
944
901
do_mpfr_lshift(int nargs)
945
902
{
946
 
        NODE *res;
 
903
        NODE *t1, *t2, *res;
947
904
        unsigned long shift;
948
 
 
949
 
        if ((res = get_bit_ops("lshift")) == NULL) {
950
 
 
951
 
                /*
952
 
                 * mpz_get_ui: If op is too big to fit an unsigned long then just
953
 
                 * the least significant bits that do fit are returned.
954
 
                 * The sign of op is ignored, only the absolute value is used.
955
 
                 */
956
 
 
957
 
                shift = mpz_get_ui(mpz2);       /* GMP integer => unsigned long conversion */
958
 
                res = mpg_integer();
959
 
                mpz_mul_2exp(res->mpg_i, mpz1, shift);          /* res = mpz1 * 2^shift */
960
 
        }
961
 
        free_bit_ops();
 
905
        mpz_ptr pz1, pz2;
 
906
 
 
907
        t2 = POP_SCALAR();
 
908
        t1 = POP_SCALAR();
 
909
 
 
910
        pz1 = get_intval(t1, 1, "lshift");
 
911
        pz2 = get_intval(t2, 2, "lshift");
 
912
 
 
913
        /*
 
914
         * mpz_get_ui: If op is too big to fit an unsigned long then just
 
915
         * the least significant bits that do fit are returned.
 
916
         * The sign of op is ignored, only the absolute value is used.
 
917
         */
 
918
 
 
919
        shift = mpz_get_ui(pz2);        /* GMP integer => unsigned long conversion */
 
920
        res = mpg_integer();
 
921
        mpz_mul_2exp(res->mpg_i, pz1, shift);           /* res = pz1 * 2^shift */
 
922
 
 
923
        free_intval(t1, pz1);
 
924
        free_intval(t2, pz2);
 
925
        DEREF(t2);
 
926
        DEREF(t1);
962
927
        return res;
963
928
}
964
929
 
965
930
/* do_mpfr_rshift --- perform a >> operation */
966
931
 
967
932
NODE *
968
 
do_mpfr_rhift(int nargs)
 
933
do_mpfr_rshift(int nargs)
969
934
{
970
 
        NODE *res;
 
935
        NODE *t1, *t2, *res;
971
936
        unsigned long shift;
972
 
 
973
 
        if ((res = get_bit_ops("rshift")) == NULL) {
974
 
                /*
975
 
                 * mpz_get_ui: If op is too big to fit an unsigned long then just
976
 
                 * the least significant bits that do fit are returned.
977
 
                 * The sign of op is ignored, only the absolute value is used.
978
 
                 */
979
 
 
980
 
                shift = mpz_get_ui(mpz2);       /* GMP integer => unsigned long conversion */
981
 
                res = mpg_integer();
982
 
                mpz_fdiv_q_2exp(res->mpg_i, mpz1, shift);       /* res = mpz1 / 2^shift, round towards −inf */
983
 
        }
984
 
        free_bit_ops();
 
937
        mpz_ptr pz1, pz2;
 
938
 
 
939
        t2 = POP_SCALAR();
 
940
        t1 = POP_SCALAR();
 
941
 
 
942
        pz1 = get_intval(t1, 1, "rshift");
 
943
        pz2 = get_intval(t2, 2, "rshift");
 
944
 
 
945
        /* N.B: See do_mpfp_lshift. */
 
946
        shift = mpz_get_ui(pz2);        /* GMP integer => unsigned long conversion */
 
947
        res = mpg_integer();
 
948
        mpz_fdiv_q_2exp(res->mpg_i, pz1, shift);        /* res = pz1 / 2^shift, round towards −inf */
 
949
 
 
950
        free_intval(t1, pz1);
 
951
        free_intval(t2, pz2);
 
952
        DEREF(t2);
 
953
        DEREF(t1);
985
954
        return res;
986
955
}
987
956
 
 
957
 
988
958
/* do_mpfr_and --- perform an & operation */
989
959
 
990
960
NODE *
991
961
do_mpfr_and(int nargs)
992
962
{
993
 
        NODE *res;
994
 
 
995
 
        if ((res = get_bit_ops("and")) == NULL) {
996
 
                res = mpg_integer();
997
 
                mpz_and(res->mpg_i, mpz1, mpz2);
 
963
        NODE *t1, *t2, *res;
 
964
        mpz_ptr pz1, pz2;
 
965
        int i;
 
966
 
 
967
        if (nargs < 2)
 
968
                fatal(_("and: called with less than two arguments"));
 
969
 
 
970
        t2 = POP_SCALAR();
 
971
        pz2 = get_intval(t2, nargs, "and");
 
972
 
 
973
        res = mpg_integer();
 
974
        for (i = 1; i < nargs; i++) {
 
975
                t1 = POP_SCALAR();
 
976
                pz1 = get_intval(t1, nargs - i, "and");
 
977
                mpz_and(res->mpg_i, pz1, pz2);
 
978
                free_intval(t1, pz1);
 
979
                DEREF(t1);
 
980
                if (i == 1) {
 
981
                        free_intval(t2, pz2);
 
982
                        DEREF(t2);
 
983
                }
 
984
                pz2 = res->mpg_i;
998
985
        }
999
 
        free_bit_ops();
1000
986
        return res;
1001
987
}
1002
988
 
 
989
 
1003
990
/* do_mpfr_or --- perform an | operation */
1004
991
 
1005
992
NODE *
1006
993
do_mpfr_or(int nargs)
1007
994
{
1008
 
        NODE *res;
1009
 
 
1010
 
        if ((res = get_bit_ops("or")) == NULL) {
1011
 
                res = mpg_integer();
1012
 
                mpz_ior(res->mpg_i, mpz1, mpz2);
1013
 
        }
1014
 
        free_bit_ops();
 
995
        NODE *t1, *t2, *res;
 
996
        mpz_ptr pz1, pz2;
 
997
        int i;
 
998
 
 
999
        if (nargs < 2)
 
1000
                fatal(_("or: called with less than two arguments"));
 
1001
 
 
1002
        t2 = POP_SCALAR();
 
1003
        pz2 = get_intval(t2, nargs, "or");
 
1004
 
 
1005
        res = mpg_integer();
 
1006
        for (i = 1; i < nargs; i++) {
 
1007
                t1 = POP_SCALAR();
 
1008
                pz1 = get_intval(t1, nargs - i, "or");
 
1009
                mpz_ior(res->mpg_i, pz1, pz2);
 
1010
                free_intval(t1, pz1);
 
1011
                DEREF(t1);
 
1012
                if (i == 1) {
 
1013
                        free_intval(t2, pz2);
 
1014
                        DEREF(t2);
 
1015
                }
 
1016
                pz2 = res->mpg_i;
 
1017
        }
 
1018
        return res;
 
1019
}
 
1020
 
 
1021
/* do_mpfr_xor --- perform an ^ operation */
 
1022
 
 
1023
NODE *
 
1024
do_mpfr_xor(int nargs)
 
1025
{
 
1026
        NODE *t1, *t2, *res;
 
1027
        mpz_ptr pz1, pz2;
 
1028
        int i;
 
1029
 
 
1030
        if (nargs < 2)
 
1031
                fatal(_("xor: called with less than two arguments"));
 
1032
 
 
1033
        t2 = POP_SCALAR();
 
1034
        pz2 = get_intval(t2, nargs, "xor");
 
1035
 
 
1036
        res = mpg_integer();
 
1037
        for (i = 1; i < nargs; i++) {
 
1038
                t1 = POP_SCALAR();
 
1039
                pz1 = get_intval(t1, nargs - i, "xor");
 
1040
                mpz_xor(res->mpg_i, pz1, pz2);
 
1041
                free_intval(t1, pz1);
 
1042
                DEREF(t1);
 
1043
                if (i == 1) {
 
1044
                        free_intval(t2, pz2);
 
1045
                        DEREF(t2);
 
1046
                }
 
1047
                pz2 = res->mpg_i;
 
1048
        }
1015
1049
        return res;
1016
1050
}
1017
1051
 
1047
1081
        return r;
1048
1082
}
1049
1083
 
1050
 
/* do_mpfr_xor --- perform an ^ operation */
1051
 
 
1052
 
NODE *
1053
 
do_mpfr_xor(int nargs)
1054
 
{
1055
 
        NODE *res;
1056
 
 
1057
 
        if ((res = get_bit_ops("xor")) == NULL) {
1058
 
                res = mpg_integer();
1059
 
                mpz_xor(res->mpg_i, mpz1, mpz2);
1060
 
        }
1061
 
        free_bit_ops();
1062
 
        return res;
1063
 
}
1064
 
 
1065
1084
 
1066
1085
static bool firstrand = true;
1067
1086
static gmp_randstate_t state;