← Back to branch summary

~ubuntu-branches/ubuntu/trusty/librep/trusty

« back to all changes in this revision

Viewing changes to src/numbers.c

  • Committer: Bazaar Package Importer
  • Author(s): Christian Marillat
  • Date: 2005-01-14 14:18:11 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050114141811-k2x3wczuc17qai2v
Tags: 0.17-7
Build with -Oo for amd64

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/numbers.c
1
1
/* numbers.c -- Implement the tower of numeric types
2
2
   Copyright (C) 1993, 1994, 2000 John Harper <john@dcs.warwick.ac.uk>
3
 
   $Id: numbers.c,v 1.56 2001/10/21 02:48:27 jsh Exp $
 
3
   $Id: numbers.c,v 1.63 2002/10/07 07:42:09 jsh Exp $
4
4
 
5
5
   This file is part of librep.
6
6
 
43
43
#include <stdlib.h>
44
44
#include <assert.h>
45
45
#include <math.h>
 
46
#include <float.h>
46
47
#include <ctype.h>
47
48
#include <limits.h>
48
49
#include <errno.h>
1236
1237
            : rep_signal_missing_arg (1));
1237
1238
}
1238
1239
 
 
1240
static inline repv
 
1241
number_foldv (int argc, repv *argv, repv (*op) (repv, repv))
 
1242
{
 
1243
    repv sum;
 
1244
    int i;
 
1245
 
 
1246
    if (argc < 1)
 
1247
        return rep_signal_missing_arg (1);
 
1248
    if (!rep_NUMERICP (argv[0]))
 
1249
        return rep_signal_arg_error (argv[0], 1);
 
1250
 
 
1251
    sum = argv[0];
 
1252
    for (i = 1; i < argc; i++)
 
1253
    {
 
1254
        if (!rep_NUMERICP (argv[i]))
 
1255
            return rep_signal_arg_error (argv[i], i + 1);
 
1256
 
 
1257
        sum = op (sum, argv[i]);
 
1258
    }
 
1259
 
 
1260
    return sum;
 
1261
}
 
1262
 
1239
1263
repv
1240
1264
rep_integer_foldl (repv args, repv (*op)(repv, repv))
1241
1265
{
1259
1283
            : rep_signal_missing_arg (1));
1260
1284
}
1261
1285
 
 
1286
static inline repv
 
1287
integer_foldv (int argc, repv *argv, repv (*op) (repv, repv))
 
1288
{
 
1289
    repv sum;
 
1290
    int i;
 
1291
 
 
1292
    if (argc < 1)
 
1293
        return rep_signal_missing_arg (1);
 
1294
    if (!rep_INTEGERP (argv[0]))
 
1295
        return rep_signal_arg_error (argv[0], 1);
 
1296
 
 
1297
    sum = argv[0];
 
1298
    for (i = 1; i < argc; i++)
 
1299
    {
 
1300
        if (!rep_INTEGERP (argv[i]))
 
1301
            return rep_signal_arg_error (argv[i], i + 1);
 
1302
 
 
1303
        sum = op (sum, argv[i]);
 
1304
    }
 
1305
 
 
1306
    return sum;
 
1307
}
 
1308
 
1262
1309
repv
1263
1310
rep_foldl (repv args, repv (*op)(repv, repv))
1264
1311
{
1279
1326
    return rep_signal_missing_arg (1);
1280
1327
}
1281
1328
 
 
1329
static inline repv
 
1330
foldv (int argc, repv *argv, repv (*op) (repv, repv))
 
1331
{
 
1332
    repv sum;
 
1333
    int i;
 
1334
 
 
1335
    if (argc < 1)
 
1336
        return rep_signal_missing_arg (1);
 
1337
 
 
1338
    sum = argv[0];
 
1339
    for (i = 1; i < argc; i++)
 
1340
    {
 
1341
        sum = op (sum, argv[i]);
 
1342
    }
 
1343
 
 
1344
    return sum;
 
1345
}
 
1346
 
1282
1347
repv
1283
1348
rep_number_add (repv x, repv y)
1284
1349
{
1724
1789
    return out;
1725
1790
}
1726
1791
 
1727
 
DEFUN("+", Fplus, Splus, (repv args), rep_SubrN) /*
 
1792
DEFUN("+", Fplus, Splus, (int argc, repv *argv), rep_SubrV) /*
1728
1793
::doc:rep.lang.math#+::
1729
1794
+ NUMBERS...
1730
1795
 
1731
1796
Adds all NUMBERS together. If no arguments are given returns 0.
1732
1797
::end:: */
1733
1798
{
1734
 
    if (args == Qnil)
 
1799
    if (argc == 0)
1735
1800
        return rep_MAKE_INT (0);
1736
1801
    else
1737
 
        return rep_number_foldl (args, rep_number_add);
 
1802
        return number_foldv (argc, argv, rep_number_add);
1738
1803
}
1739
1804
 
1740
 
DEFUN("-", Fminus, Sminus, (repv args), rep_SubrN) /*
 
1805
DEFUN("-", Fminus, Sminus, (int argc, repv *argv), rep_SubrV) /*
1741
1806
::doc:rep.lang.math#-::
1742
1807
- NUMBER [NUMBERS...]
1743
1808
 
1745
1810
NUMBERS
1746
1811
::end:: */
1747
1812
{
1748
 
    if (args == Qnil)
 
1813
    if (argc == 0)
1749
1814
        return rep_signal_missing_arg (1);
1750
 
    else if (!rep_CONSP (rep_CDR (args)))
1751
 
        return rep_number_neg (rep_CAR (args));
 
1815
    else if (argc == 1)
 
1816
        return rep_number_neg (argv[0]);
1752
1817
    else
1753
 
        return rep_number_foldl (args, rep_number_sub);
 
1818
        return number_foldv (argc, argv, rep_number_sub);
1754
1819
}
1755
1820
 
1756
 
DEFUN("*", Fproduct, Sproduct, (repv args), rep_SubrN) /*
 
1821
DEFUN("*", Fproduct, Sproduct, (int argc, repv *argv), rep_SubrV) /*
1757
1822
::doc:rep.lang.math#*::
1758
1823
* NUMBERS...
1759
1824
 
1760
1825
Multiplies all NUMBERS together. If no numbers are given returns 1.
1761
1826
::end:: */
1762
1827
{
1763
 
    if (args == Qnil)
 
1828
    if (argc == 0)
1764
1829
        return rep_MAKE_INT (1);
1765
1830
    else
1766
 
        return rep_number_foldl (args, rep_number_mul);
 
1831
        return number_foldv (argc, argv, rep_number_mul);
1767
1832
}
1768
1833
 
1769
 
DEFUN("/", Fdivide, Sdivide, (repv args), rep_SubrN) /*
 
1834
DEFUN("/", Fdivide, Sdivide, (int argc, repv *argv), rep_SubrV) /*
1770
1835
::doc:rep.lang.math#/::
1771
1836
/ NUMBERS...
1772
1837
 
1773
1838
Divides NUMBERS (in left-to-right order).
1774
1839
::end:: */
1775
1840
{
1776
 
    if (args == Qnil)
 
1841
    if (argc == 0)
1777
1842
        return rep_signal_missing_arg (1);
1778
 
    else if (!rep_CONSP (rep_CDR (args)))
1779
 
        return rep_number_div (rep_MAKE_INT (1), rep_CAR (args));
1780
 
    return rep_number_foldl (args, rep_number_div);
 
1843
    else if (argc == 1)
 
1844
        return rep_number_div (rep_MAKE_INT (1), argv[0]);
 
1845
    else
 
1846
        return number_foldv (argc, argv, rep_number_div);
1781
1847
}
1782
1848
 
1783
1849
DEFUN("remainder", Fremainder, Sremainder, (repv n1, repv n2), rep_Subr2) /*
1920
1986
    return rep_number_lognot (num);
1921
1987
}
1922
1988
 
1923
 
DEFUN("logior", Flogior, Slogior, (repv args), rep_SubrN) /*
 
1989
DEFUN("logior", Flogior, Slogior, (int argc, repv *argv), rep_SubrV) /*
1924
1990
::doc:rep.lang.math#logior::
1925
1991
logior NUMBERS...
1926
1992
 
1927
1993
Returns the bitwise logical `inclusive-or' of its arguments.
1928
1994
::end:: */
1929
1995
{
1930
 
    if (args == Qnil)
 
1996
    if (argc == 0)
1931
1997
        return rep_MAKE_INT (0);
1932
1998
    else
1933
 
        return rep_number_foldl (args, rep_number_logior);
 
1999
        return number_foldv (argc, argv, rep_number_logior);
1934
2000
}
1935
2001
 
1936
 
DEFUN("logxor", Flogxor, Slogxor, (repv args), rep_SubrN) /*
 
2002
DEFUN("logxor", Flogxor, Slogxor, (int argc, repv *argv), rep_SubrV) /*
1937
2003
::doc:rep.lang.math#logxor::
1938
2004
logxor NUMBERS...
1939
2005
 
1940
2006
Returns the bitwise logical `exclusive-or' of its arguments.
1941
2007
::end:: */
1942
2008
{
1943
 
    return rep_number_foldl (args, rep_number_logxor);
 
2009
    return number_foldv (argc, argv, rep_number_logxor);
1944
2010
}
1945
2011
 
1946
 
DEFUN("logand", Flogand, Slogand, (repv args), rep_SubrN) /*
 
2012
DEFUN("logand", Flogand, Slogand, (int argc, repv *argv), rep_SubrV) /*
1947
2013
::doc:rep.lang.math#logand::
1948
2014
logand NUMBERS...
1949
2015
 
1950
2016
Returns the bitwise logical `and' of its arguments.
1951
2017
::end:: */
1952
2018
{
1953
 
    return rep_number_foldl (args, rep_number_logand);
 
2019
    return number_foldv (argc, argv, rep_number_logand);
1954
2020
}
1955
2021
 
1956
2022
DEFUN("eql", Feql, Seql, (repv arg1, repv arg2), rep_Subr2) /*
2491
2557
    return out;
2492
2558
}
2493
2559
 
2494
 
DEFUN("gcd", Fgcd, Sgcd, (repv args), rep_SubrN) /*
 
2560
DEFUN("gcd", Fgcd, Sgcd, (int argc, repv *argv), rep_SubrV) /*
2495
2561
::doc:rep.lang.math#gcd::
2496
2562
gcd ...
2497
2563
 
2499
2565
is always non-negative. Returns 0 with arguments.
2500
2566
::end:: */
2501
2567
{
2502
 
    if (args == Qnil)
 
2568
    if (argc == 0)
2503
2569
        return rep_MAKE_INT (0);
2504
 
    else if (rep_CONSP (args) && rep_CDR (args) == Qnil)
 
2570
    else if (argc == 1)
2505
2571
    {
2506
 
        rep_DECLARE1 (rep_CAR (args), rep_INTEGERP);
2507
 
        return rep_integer_gcd (rep_CAR (args), rep_CAR (args));
 
2572
        rep_DECLARE1 (argv[0], rep_INTEGERP);
 
2573
        return rep_integer_gcd (argv[0], argv[0]);
2508
2574
    }
2509
2575
    else
2510
 
        return rep_integer_foldl (args, rep_integer_gcd);
 
2576
        return integer_foldv (argc, argv, rep_integer_gcd);
2511
2577
}
2512
2578
 
2513
2579
DEFUN("numberp", Fnumberp, Snumberp, (repv arg), rep_Subr1) /*
2582
2648
        return rep_make_float (rep_get_float (arg), rep_TRUE);
2583
2649
}
2584
2650
 
 
2651
static void
 
2652
rationalize (repv arg, double *numerator, double *denominator)
 
2653
{
 
2654
    double x, y;
 
2655
    int expt;
 
2656
 
 
2657
    /* X/Y always equals the input value. Tactic is to iteratively
 
2658
       multiply both X and Y by 2 until X is an integer. We bound
 
2659
       the number of iterations to the size of the mantissa
 
2660
       by starting with the normalized value... */
 
2661
 
 
2662
    x = frexp (rep_get_float (arg), &expt);
 
2663
    y = pow (2.0, -expt);
 
2664
 
 
2665
    while (x - floor (x) > DBL_EPSILON)
 
2666
    {
 
2667
        x = x * 2.0;
 
2668
        y = y * 2.0;
 
2669
    }
 
2670
 
 
2671
    if (numerator != NULL)
 
2672
        *numerator = x;
 
2673
    if (denominator != NULL)
 
2674
        *denominator = y;
 
2675
}
 
2676
 
2585
2677
DEFUN("inexact->exact", Finexact_to_exact,
2586
2678
      Sinexact_to_exact, (repv arg), rep_Subr1) /*
2587
2679
::doc:rep.lang.math#inexact->exact::
2592
2684
::end:: */
2593
2685
{
2594
2686
    rep_DECLARE1(arg, rep_NUMERICP);
 
2687
 
2595
2688
    if (rep_INTP (arg) || !rep_NUMBER_FLOAT_P (arg))
2596
2689
        return arg;
2597
 
    else
2598
 
    {
2599
 
        /* XXX is there a way to decide if it's rationalizable? */
2600
 
        double d = floor (rep_get_float (arg));
2601
 
        if (d >= rep_LISP_MIN_INT && d <= rep_LISP_MAX_INT)
2602
 
            return rep_MAKE_INT ((long) d);
2603
 
        else
2604
 
        {
2605
 
            rep_number_z *z = make_number (rep_NUMBER_BIGNUM);
 
2690
 
2606
2691
#ifdef HAVE_GMP
2607
 
            mpz_init_set_d (z->z, d);
 
2692
    else
 
2693
    {
 
2694
        rep_number_q *q;
 
2695
 
 
2696
        q = make_number (rep_NUMBER_RATIONAL);
 
2697
        mpq_init (q->q);
 
2698
        mpq_set_d (q->q, rep_get_float (arg));
 
2699
 
 
2700
        return maybe_demote (rep_VAL (q));
 
2701
    }
2608
2702
#else
2609
 
            if (d >= BIGNUM_MAX)
2610
 
                z->z = BIGNUM_MAX;
2611
 
            else if (d <= BIGNUM_MIN)
2612
 
                z->z = BIGNUM_MIN;
2613
 
            else
2614
 
                z->z = (rep_long_long) d;
 
2703
    else
 
2704
    {
 
2705
        double x, y;
 
2706
        rep_number_z *z;
 
2707
 
 
2708
        rationalize (arg, &x, &y);
 
2709
        z = make_number (rep_NUMBER_BIGNUM);
 
2710
        z->z = x / y;
 
2711
 
 
2712
        return maybe_demote (rep_VAL (z));
 
2713
    }
2615
2714
#endif
2616
 
            return rep_VAL (z);
2617
 
        }
2618
 
    }
2619
2715
}
2620
2716
 
2621
 
DEFUN("numerator", Fnumerator, Snumerator, (repv x), rep_Subr1) /*
 
2717
DEFUN("numerator", Fnumerator, Snumerator, (repv arg), rep_Subr1) /*
2622
2718
::doc:rep.lang.math#numerator::
2623
2719
numerator X
2624
2720
 
2625
2721
Return the numerator of rational number X.
2626
2722
::end:: */
2627
2723
{
2628
 
    rep_DECLARE1(x, rep_NUMERICP);
2629
 
    if (rep_INTP (x) || rep_NUMBER_BIGNUM_P (x))
2630
 
        return x;
 
2724
    rep_bool inexact = rep_FALSE;
 
2725
    double x;
 
2726
 
 
2727
    rep_DECLARE1(arg, rep_NUMERICP);
 
2728
 
2631
2729
#ifdef HAVE_GMP
2632
 
    else if (rep_NUMBER_RATIONAL_P (x))
 
2730
    if (rep_NUMBER_RATIONAL_P (arg))
2633
2731
    {
2634
2732
        rep_number_z *z = make_number (rep_NUMBER_BIGNUM);
2635
 
        mpz_init_set (z->z, mpq_numref (rep_NUMBER(x,q)));
 
2733
        mpz_init_set (z->z, mpq_numref (rep_NUMBER(arg,q)));
2636
2734
        return maybe_demote (rep_VAL (z));
2637
2735
    }
2638
2736
#endif
2639
 
    else
2640
 
        return rep_signal_arg_error (x, 1);
 
2737
 
 
2738
    if (rep_NUMBER_INEXACT_P (arg))
 
2739
        inexact = rep_TRUE;
 
2740
 
 
2741
    rationalize (arg, &x, NULL);
 
2742
 
 
2743
    return rep_make_float (x, inexact);
2641
2744
}
2642
2745
 
2643
 
DEFUN("denominator", Fdenominator, Sdenominator, (repv x), rep_Subr1) /*
 
2746
DEFUN("denominator", Fdenominator, Sdenominator, (repv arg), rep_Subr1) /*
2644
2747
::doc:rep.lang.math#denominator::
2645
2748
denominator X
2646
2749
 
2647
2750
Return the denominator of rational number X.
2648
2751
::end:: */
2649
2752
{
2650
 
    rep_DECLARE1(x, rep_NUMERICP);
2651
 
    if (rep_INTP (x) || rep_NUMBER_BIGNUM_P (x))
2652
 
        return rep_MAKE_INT (1);
 
2753
    rep_bool inexact = rep_FALSE;
 
2754
    double y;
 
2755
 
 
2756
    rep_DECLARE1(arg, rep_NUMERICP);
 
2757
 
2653
2758
#ifdef HAVE_GMP
2654
 
    else if (rep_NUMBER_RATIONAL_P (x))
 
2759
    if (rep_NUMBER_RATIONAL_P (arg))
2655
2760
    {
2656
2761
        rep_number_z *z = make_number (rep_NUMBER_BIGNUM);
2657
 
        mpz_init_set (z->z, mpq_denref (rep_NUMBER(x,q)));
 
2762
        mpz_init_set (z->z, mpq_denref (rep_NUMBER(arg,q)));
2658
2763
        return maybe_demote (rep_VAL (z));
2659
2764
    }
2660
2765
#endif
2661
 
    else
2662
 
        return rep_signal_arg_error (x, 1);
 
2766
 
 
2767
    if (rep_NUMBER_INEXACT_P (arg))
 
2768
        inexact = rep_TRUE;
 
2769
 
 
2770
    rationalize (arg, NULL, &y);
 
2771
 
 
2772
    return rep_make_float (y, inexact);
2663
2773
}
2664
2774
 
2665
 
DEFUN("max", Fmax, Smax, (repv args), rep_SubrN) /*
 
2775
DEFUN("max", Fmax, Smax, (int argc, repv *argv), rep_SubrV) /*
2666
2776
::doc:rep.lang.math#max::
2667
2777
max ARGS...
2668
2778
 
2671
2781
result to be inexact.
2672
2782
::end:: */
2673
2783
{
2674
 
    return rep_foldl (args, rep_number_max);
 
2784
    return foldv (argc, argv, rep_number_max);
2675
2785
}
2676
2786
 
2677
 
DEFUN("min", Fmin, Smin, (repv args), rep_SubrN) /*
 
2787
DEFUN("min", Fmin, Smin, (int argc, repv *argv), rep_SubrV) /*
2678
2788
::doc:rep.lang.math#min::
2679
2789
min ARGS...
2680
2790
 
2683
2793
result to be inexact.
2684
2794
::end:: */
2685
2795
{
2686
 
    return rep_foldl (args, rep_number_min);
 
2796
    return foldv (argc, argv, rep_number_min);
2687
2797
}
2688
2798
 
2689
2799
DEFUN("string->number", Fstring_to_number,
2797
2907
 
2798
2908
/* Random number generation */
2799
2909
 
2800
 
#if defined (HAVE_GMP) && defined (HAVE_GMP_RANDINIT)
 
2910
#if defined (HAVE_GMP) && defined (HAVE_GMP_RANDINIT) && __GNU_MP__ >= 4
2801
2911
 
2802
2912
static gmp_randstate_t random_state;
2803
2913
 
2911
3021
 
2912
3022
    if (arg == Qt)
2913
3023
    {
2914
 
        random_seed (time (0));
 
3024
        u_long seed = time (0);
 
3025
        seed = (seed << 8) | (rep_getpid () & 0xff);
 
3026
        random_seed (seed);
2915
3027
        return Qt;
2916
3028
    }
2917
3029