ratio implementation

ratio implementation

[Bill Schottstaedt]

There is an implementation of ratios for Guile (based on CVS version
of 27-Jul-03) at ccrma-ftp.stanford.edu:/pub/Lisp/gratio.tar.gz.
Rather than send a huge diff, I placed the new versions of the changed
files (all from the libguile directory) in the tarball along with
*.diff showing the changes.  I added numerator, denominator,
rationalize and ratio?, and at the C level scm_make_ratio and
scm_i_ratio2real (should it be scm_ratio2dbl?).  "ratio?" is needed
because "rational?" returns #t if passed a real -- there has to be
some way to distiguish a ratio from a real. An alternative would be to
make "rational?"  rational.

I don't know how the FSF/GPL copyright stuff works, but I did this
work on my own time, did not look at any other implementation, and
hereby donate the code to you.  I'd be happy to "sign the papers".

I use longs for the numerator and denominator, so if bignums are
encountered, I fallback on the old method using scm_divide.  I didn't
try to provide the standard set of C-level type conversions -- not
sure these are needed in this case. I decided to follow Common Lisp
and reduce ratios, and return an integer if the ratio denominator is 1
(or if the numerator is 0).  One ugly change is the SCM_NUMP macro in
numbers.h -- the original of this macro strikes me as a kludge; the
new version is worse.  A questionable change (I followed Dybvig's book
p 126): (/ 17) -> 1/17, (/ 2 3) -> 2/3.  I notice that r5rs seems to
imply that (inexact->exact .3) should return 3/10 (see the rationalize
example which is assuming this) -- should I implement this?  What is
"reasonably close" in this case?  Another gray area involves functions
like round -- in r5rs (round 7/2) returns 4, whereas Guile returns 4.0
(I left it this way because (round 7) currently returns 7.0 rather
than 7) -- I decided to make minimal changes, but handling of
exact/inexact distinctions in Guile could use some work (I am willing
do this, if others approve).  Also, should random accept ratios?
And currently (format #f "~F" 2/3) hangs, but so does
(format #f "~B" 1.5) (in Guile 1.6.4 you get an error, but
in the CVS Guile it hangs in mutex_lock).


The changed files are:

smob.c
hash.c
eq.c
goops.c
gc-card.c
objects.c
numbers.c (nearly all the changes are in this file)
eval.c

numbers.h
tags.h
objects.h

Also in the tarball: ratio-tests.scm.  




_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

ratio implementation

[Han-Wen Nienhuys]

bil@ccrma.Stanford.EDU writes:
> There is an implementation of ratios for Guile (based on CVS version
> of 27-Jul-03) at ccrma-ftp.stanford.edu:/pub/Lisp/gratio.tar.gz.
> Rather than send a huge diff, I placed the new versions of the changed
> files (all from the libguile directory) in the tarball along with
> *.diff showing the changes.  I added numerator, denominator,
> rationalize and ratio?, and at the C level scm_make_ratio and
> scm_i_ratio2real (should it be scm_ratio2dbl?).  "ratio?" is needed
> because "rational?" returns #t if passed a real -- there has to be
> some way to distiguish a ratio from a real. An alternative would be to
> make "rational?"  rational.
> 
> I don't know how the FSF/GPL copyright stuff works, but I did this
> work on my own time, did not look at any other implementation, and
> hereby donate the code to you.  I'd be happy to "sign the papers".


Cool!

(If that goes through, then we can scrap our own support for Rational
numbers in LilyPond!)


-- 

Han-Wen Nienhuys   |   hanwen@cs.uu.nl   |   http://www.xs4all.nl/~hanwen 


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:
>
> I use longs for the numerator and denominator, so if bignums are
> encountered, I fallback on the old method using scm_divide.

It'd be nice to allow bignum numerators and denominators.  It's rather
easy to make a calculation where the denominator starts to blow up.

The gmp mpq functions might help, but there's no particular need to
use them, most of the same things could be done directly just with a
pair of scheme integers.

> I notice that r5rs seems to
> imply that (inexact->exact .3) should return 3/10 (see the rationalize
> example which is assuming this)

Isn't that the effect of the rationalize function, rather than
inexact->exact?


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Bill Schottstaedt]

Kevin Ryde wrote:

>>I notice that r5rs seems to
>>imply that (inexact->exact .3) should return 3/10 (see the rationalize
>>example which is assuming this)
>>    
>>
>
>Isn't that the effect of the rationalize function, rather than
>inexact->exact?
>  
>
The example I was referring to was:

(rationalize
  (inexact->exact .3) 1/10)            ==> 1/3    ; exact


which only makes sense if (inexact->exact .3) does not
return 0, but 3/10.  Then rationalize of that with an
error of 1/10 finds the simpler ratio 1/3. In Guile
currently (I'm typing from work, so I can't actually
check this...), the (inexact->exact .3) returns 0,
so rationalize would also return 0.





_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:
>
> which only makes sense if (inexact->exact .3) does not
> return 0, but 3/10.

I might have imagined it returning the binary fraction corresponding
to a float 0.3, eg. 0x13333333333333/0x40000000000000 for a double.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

> There is an implementation of ratios for Guile (based on CVS version
> of 27-Jul-03) at ccrma-ftp.stanford.edu:/pub/Lisp/gratio.tar.gz.
> Rather than send a huge diff, I placed the new versions of the changed
> files (all from the libguile directory) in the tarball along with
> *.diff showing the changes.

Good work, thanks!  However, I don't think we should restrict us to
longs as the numerator/denominator, we should use fixnums and bignums.
Your patch is a very good basis for this and it should be quite
straightforward to make it use SCM integers as the
numerator/denominator.

That is, instead of using "long", I'd say we simply use "SCM" and
instead of "+", "==", etc, we use use "scm_sum", "scm_num_eq_p", etc.
Also, instead of mallocing scm_t_ratio, we can use double cells, which
should be more efficient.

Any takers?

(Also I didn't really check whether your rationals behave like R5RS
demands it.  Did you?  We should be sure to follow R5RS.)

> I added numerator, denominator, rationalize and ratio?, and at the C
> level scm_make_ratio and scm_i_ratio2real (should it be
> scm_ratio2dbl?).  "ratio?" is needed because "rational?" returns #t
> if passed a real -- there has to be some way to distiguish a ratio
> from a real.

A ratio is exact, while a real is not.  Thus, I'd say that ratio? is
not really necessary.  There should be no reason to add anything eyond
R5RS for ratios.

> An alternative would be to make "rational?"  rational.

That would violate R5RS, wouldn't it?

> I don't know how the FSF/GPL copyright stuff works, but I did this
> work on my own time, did not look at any other implementation, and
> hereby donate the code to you.  I'd be happy to "sign the papers".

Thanks!  I'll get back to you about the papers when necessary.

> I decided to make minimal changes, but handling of exact/inexact
> distinctions in Guile could use some work (I am willing do this, if
> others approve).

That would be great!

> And currently (format #f "~F" 2/3) hangs, but so does (format #f
> "~B" 1.5) (in Guile 1.6.4 you get an error, but in the CVS Guile it
> hangs in mutex_lock).

This should be fixed in CVS.

-- 
GPG: D5D4E405 - 2F9B BCCC 8527 692A 04E3  331E FAF8 226A D5D4 E405


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Rob Browning]

Marius Vollmer <mvo@zagadka.de> writes:

> (Also I didn't really check whether your rationals behave like R5RS
> demands it.  Did you?  We should be sure to follow R5RS.)

And ideally, it'd be nice to have tests in "make check" to verify as
much as possible.  Makes working on the code a lot less
nerve-wracking.

-- 
Rob Browning
rlb @defaultvalue.org and @debian.org; previously @cs.utexas.edu
GPG starting 2002-11-03 = 14DD 432F AE39 534D B592  F9A0 25C8 D377 8C7E 73A4


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Dirk Herrmann]

Marius Vollmer wrote:

 >Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:
 >
 >>There is an implementation of ratios for Guile (based on CVS version
 >>of 27-Jul-03) at ccrma-ftp.stanford.edu:/pub/Lisp/gratio.tar.gz.
 >>Rather than send a huge diff, I placed the new versions of the changed
 >>files (all from the libguile directory) in the tarball along with
 >>*.diff showing the changes.
 >
 >
 >Good work, thanks!  However, I don't think we should restrict us to
 >longs as the numerator/denominator, we should use fixnums and bignums.
 >Your patch is a very good basis for this and it should be quite
 >straightforward to make it use SCM integers as the
 >numerator/denominator.
 >
 >That is, instead of using "long", I'd say we simply use "SCM" and
 >instead of "+", "==", etc, we use use "scm_sum", "scm_num_eq_p", etc.
 >Also, instead of mallocing scm_t_ratio, we can use double cells, which
 >should be more efficient.
 >
 >Any takers?

Hello together,

the last couple of days I have been looking closely at Guile's 
typecodes, and
I figured that we could improve the treating of numbers in a general way.  I
have to admit that I have not looked into Rob's code, thus my 
observations and
suggestions are based on the current CVS solution only.

Currently, fixnums are realized as immediates, while all other number types
(bignums, reals and complexs) are realized as smobs.  On the other hand, 
there
are at the moment two unused tc7 codes.

My suggestion is to use one of the free tc7 codes to represent a 'number'
type.  First, it would simplify parts of the smob code and free three smob
numbers.  Second, it would make type dispatch for numbers a little faster at
some places.  Third, in contrast to smobs which use 16 type code bits, a
smaller number of typecode bits may be required, making it possible to store
more information in the SCM_CELL_TYPE mask.

In addition to the tc7 bits, additional bits can be used for subtyping the
suggested scm_tc7_number type.  The following subtypes need to be provided:

 - bignums, holding a pointer to a gmp structure
 - doubles, holding a 64 bit double precision number
 - fraction, holding two SCM values for nominator and denominator
 - complex, holding two SCM values for real and imaginary part

The following subtypes can be provided additionally to support certain
features:

 - fixnums, holding an inum SCM value.  This type can be provided to allow
   inexact fixnum values.
 - floats, holding a 32 bit single precision number
 - further floating point representations like mpfr numbers
 - fixed point numbers
 - ...

For example, the four bits above the tc7 bits could be defined as follows
(still using only 11 bits for type information, compared to 16 bits with the
current solution):

 xsss

with x indicating the exactness, and sss indicating the subtype.  The 
subtypes
could, for example, be assigned as follows (certainly, only the mandatory
types listed above would need to be supported):

 000 - Fixnum (inum in cell word #2), optional to allow inexact fixnums
 001 - Bignum (gmp struct ptr in cell word #2), typically exact, optionally
       inexact to allow inexact bignums.
 010 - fraction (all combinations of fixnums and bignums stored as SCM
       values in double cell words #2 #3), typically exact, optionally 
inexact
       if inexact fixnums or inexact bignums are to be supported.
 011 - fixed point number (implementation not clear yet, just listed for
       completeness), exact or inexact
 100 - float (32 bit float in cell word #2), always inexact
 101 - double (64 bit float in double cell words #2 #3), always inexact
 110 - mpfr (gmp struct ptr in cell word #2), always inexact
 111 - complex, stored as SCM values in double cell words #2 #3, possible
       combinations are (double double) - always inexact, (float, float) -
       always inexact, (fixedpoint fixedpoint) - inexact or exact,
       (mpfr mpfr) - always inexact, or all combinations of fixnums, bignums
       and fractions - typically exact, optionally inexact

Since fractions and complex numbers are represented in terms of pairs of 
other
numbers, operations on fractions and complex numbers can be delegated to
operations that work on the remaining, more basic representations.

The limitations on the possible combinations for real and imaginary part of
complex numbers has the background that it does not seem to make sense to me
to represent real and imaginary part with different precision.

Even if no inexact fixnums, bignums and rationals are to be supported it 
does
still make sense to provide the exactness bit for faster operation of the
exact? predicate.


I had been planning, as a first step, to simply allocate that scm_tc7_number
typecode and turn bignums, doubles and complex numbers into subtypes of that
type.  I had not originally intended to add rationals in the same run.  
Maybe
I could help to prepare their addition.  However, there are still a 
number of
changes with respect to separating memoization and execution in the 
evaluator,
so it may take some time until I start to work on it.

Best regards
Dirk Herrmann



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

Dirk Herrmann <dirk@dirk-herrmanns-seiten.de> writes:
>
> Since fractions and complex numbers are represented in terms of pairs
> of other
> numbers, operations on fractions and complex numbers can be delegated to
> operations that work on the remaining, more basic representations.

For complex numbers where real and imag are "double"s it'd be nice to
make use of the c99 libm "complex double" functions where available.
That can probably come about naturally anyway, but just something to
bear in mind.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Bill Schottstaedt]

 > That is, instead of using "long", I'd say we simply use "SCM" and
 > instead of "+", "==", etc, we use use "scm_sum", "scm_num_eq_p", etc.
 > Also, instead of mallocing scm_t_ratio, we can use double cells, which
 > should be more efficient.
 >
 > Any takers?

Looks straightforward, but I may not get around to it right away.
Should I wait for Dirk's type-related changes?  Also, I haven't
looked closely at double cells, though I assume (from a glance
at the scm_t_double stuff) that I can include my "reduce" flag
as well as the numerator and denominator.  The point being to
reduce calls on gcd as much as possible (this being Guile's integer
divide, there's some desire to keep it from crawling).  But, it's
an optional optimization.


 > (Also I didn't really check whether your rationals behave like R5RS
 > demands it.  Did you?  We should be sure to follow R5RS.)

Yes, I think it follows R5RS.  See test-ratios.scm for some tests.
I could make these fit into 'make check'.


 > A ratio is exact, while a real is not.  Thus, I'd say that ratio? is
 > not really necessary.  There should be no reason to add anything eyond
 > R5RS for ratios.

Ok -- "exact?" used to return #f for 3/4 -- hence the confusion.


 >> An alternative would be to make "rational?"  rational.
 > That would violate R5RS, wouldn't it?

I don't think so.  You are allowed to make rational? be the same
as real?, but you can also distinguish them.  My druthers would
be for rational? to return #t for a ratio, #f for a float.
Why have it at all if it is just a synonym for "real?"?




_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Rob Browning]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

> Also, I haven't looked closely at double cells, though I assume
> (from a glance at the scm_t_double stuff) that I can include my
> "reduce" flag as well as the numerator and denominator.

A while back, at Marius' suggestion, I reworked my original gmp bignum
implementation to allocate the mpz_t's using a double cell (which
definitely did help allocation) -- so you might want to check out the
bignum usage in numbers.c for an example.  Using a double cell meant I
didn't need to do anything to free the mpz_t memory in gc-card.c
(i.e. I just needed to mpz_clear the mpz_t).

-- 
Rob Browning
rlb @defaultvalue.org and @debian.org; previously @cs.utexas.edu
GPG starting 2002-11-03 = 14DD 432F AE39 534D B592  F9A0 25C8 D377 8C7E 73A4


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Dirk Herrmann]

Bill Schottstaedt wrote:

> Looks straightforward, but I may not get around to it right away.
> Should I wait for Dirk's type-related changes?

Just for your information: I have submitted a patch which does the first 
part: Introducing scm_tc7_number and turning bignums, reals and complex 
numbers from smobs into subtypes of scm_tc7_number. The patch, however, 
is quite simplistic and does not yet optimize the way the subtypes of 
scm_tc7_number are allocated. It should nevertheless provide the 
type-system infrastructure for easy integration of the fraction type.

> > A ratio is exact, while a real is not.  Thus, I'd say that ratio? is
> > not really necessary.  There should be no reason to add anything eyond
> > R5RS for ratios.
>
> Ok -- "exact?" used to return #f for 3/4 -- hence the confusion.
>
>
> >> An alternative would be to make "rational?"  rational.
> > That would violate R5RS, wouldn't it?
>
> I don't think so.  You are allowed to make rational? be the same
> as real?, but you can also distinguish them.  My druthers would
> be for rational? to return #t for a ratio, #f for a float.
> Why have it at all if it is just a synonym for "real?"? 

Just a side note: We should in guile distinguish between the internal 
representation of values, and their mathematical meaning. For example, 
fixnums and bignums are two representation of integer types. They 
should, however, also report #t when queried about being rational, real 
or complex, since they are elements of all these sets. When introducing 
exact rationals, which are represented internally as a pair of two exact 
integer values, we should not name them rationals, but, for example, 
fractions.

The reason I mention this is that I want to avoid a macro SCM_RATIONAL_P 
to be introduced: The meaning of such a macro would not be clear, since 
it might either be expected to return the equivalent of rational? or 
just return true if an object was detected with an internal 
representation of a fraction.

A nice set of macros to check for numbers of a specific internal 
representation could be (names are just given as examples):

SCM_FIXNUM_P (currently named SCM_INUMP)
SCM_BIGNUM_P (currently named SCM_BIGP)
SCM_FRACTNUM_P
SCM_FLOATNUM_P (currently named SCM_REALP, which is confusing)
SCM_FLOATCPLXNUM_P (currently named SCM_COMPLEXP)

while a nice set of macros for testing for a mathematical meaning could 
be (names are probably acceptable):

SCM_EXACT_P
SCM_NUMBER_P
SCM_INTEGER_P
SCM_RATIONAL_P
SCM_REAL_P
SCM_COMPLEX_P

In the long term, we should rethink our current set of confusing names. 
For the short term, we should at least avoid to introduce more confusing 
terms, especially when introducing the new internal representation of 
fractions.

Best regards
Dirk Herrmann



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

>  > That is, instead of using "long", I'd say we simply use "SCM" and
>  > instead of "+", "==", etc, we use use "scm_sum", "scm_num_eq_p", etc.
>  > Also, instead of mallocing scm_t_ratio, we can use double cells, which
>  > should be more efficient.
>  >
>  > Any takers?
>
> Looks straightforward, but I may not get around to it right away.
> Should I wait for Dirk's type-related changes?

No, I don't think that is needed.  Dirk, what do you think?

> Also, I haven't looked closely at double cells, though I assume
> (from a glance at the scm_t_double stuff) that I can include my
> "reduce" flag as well as the numerator and denominator.

Yes, double cells are just like ordinary cells, they only have four
slots instead of two.  One (half) slot is taken for type information,
leaving three for the fraction stuff.  There are unused bits in the
first slot, and you could put your 'reduce' flag there as well.

>  > (Also I didn't really check whether your rationals behave like R5RS
>  > demands it.  Did you?  We should be sure to follow R5RS.)
>
> Yes, I think it follows R5RS.  See test-ratios.scm for some tests.
> I could make these fit into 'make check'.

That would be good.

>  > A ratio is exact, while a real is not.  Thus, I'd say that ratio? is
>  > not really necessary.  There should be no reason to add anything eyond
>  > R5RS for ratios.
>
> Ok -- "exact?" used to return #f for 3/4 -- hence the confusion.

Yes, but only because we didn't have exact rationals, right?

>  >> An alternative would be to make "rational?"  rational.
>  > That would violate R5RS, wouldn't it?
>
> I don't think so.  You are allowed to make rational? be the same
> as real?, but you can also distinguish them.  My druthers would
> be for rational? to return #t for a ratio, #f for a float.
> Why have it at all if it is just a synonym for "real?"?

Hmmm, I see.  My take would be that our floats are always inexact
rationals.  An example for a real that wouldn't also be a rational
would be an exact representation of pi, I'd say.  But I also agree
that this way lays madness and that you'd probably need a full blown
symbolic algebra system for this to actually be useful.

-- 
GPG: D5D4E405 - 2F9B BCCC 8527 692A 04E3  331E FAF8 226A D5D4 E405


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Dirk Herrmann <dirk@dirk-herrmanns-seiten.de> writes:

> When introducing exact rationals, which are represented internally
> as a pair of two exact integer values, we should not name them
> rationals, but, for example, fractions.
>
> The reason I mention this is that I want to avoid a macro
> SCM_RATIONAL_P to be introduced: The meaning of such a macro would not
> be clear, since it might either be expected to return the equivalent
> of rational? or just return true if an object was detected with an
> internal representation of a fraction.

Yes, very good point.  Bill, please name your data type "fraction".

-- 
GPG: D5D4E405 - 2F9B BCCC 8527 692A 04E3  331E FAF8 226A D5D4 E405


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Bill Schottstaedt]

I made the requested changes (using the Oct-9 CVS Guile):

    ftp://ccrma-ftp.stanford.edu/pub/Lisp/gratio-1.tar.gz

I haven't yet translated the test suite to Guile's form.


I notice that Guile doesn't support "exact complex", but
it appears that some other Schemes do -- is this a desired
extension?

guile> (exact? 2/3+i)
#f



Under scm_ash, the doc string:

	    "Formally, the function returns an integer equivalent to\n"
	    "@code{(inexact->exact (floor (* @var{n} (expt 2 @var{cnt}))))}.\n"

is incorrect:

guile> (ash .1 2)
0.4

but ash is not handling negative shifts correctly:

guile> (ash .1 -2)

Backtrace:
In standard input:
    4: 0* [ash 0.1 -2]

standard input:4:1: In procedure quotient in expression (ash 0.1 -2):
standard input:4:1: Wrong type argument in position 1: 0.1
ABORT: (wrong-type-arg)

The scm_quotient function assumes its arguments are ints/bignums,
so it can't be used directly in ash.


Also: (sqrt 9/49) -> 3/7?  Seems like it should parallel:
guile> (expt 2 -2)
1/4



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

> I made the requested changes (using the Oct-9 CVS Guile):
>
>     ftp://ccrma-ftp.stanford.edu/pub/Lisp/gratio-1.tar.gz

Excellent!  I'm going to look it over and apply it once we have the
papers from you.  Thanks a lot!

(More comments later.)


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

> I notice that Guile doesn't support "exact complex", but
> it appears that some other Schemes do -- is this a desired
> extension?
>
> guile> (exact? 2/3+i)
> #f

Yes, it would be welcome.

> Under scm_ash, the doc string:
>
> 	    "Formally, the function returns an integer equivalent to\n"
> 	    "@code{(inexact->exact (floor (* @var{n} (expt 2 @var{cnt}))))}.\n"
>
> is incorrect:
> [...]
> The scm_quotient function assumes its arguments are ints/bignums,
> so it can't be used directly in ash.

Right!  Do you have a patch? :-) I will record this in our bug
database nevertheless.

> Also: (sqrt 9/49) -> 3/7?  Seems like it should parallel:
> guile> (expt 2 -2)
> 1/4

Hmm.  We would then also need

  (sqrt 4) => 2 (exact)

I'm not sure whether exact square roots are imortant: they will only
be useful when both the numerator and denominator are squares of
integers and those pairs will be rare, I'd say.  We could use GMP for
computing bignum square roots...


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

> I made the requested changes (using the Oct-9 CVS Guile):
>
>     ftp://ccrma-ftp.stanford.edu/pub/Lisp/gratio-1.tar.gz

With your new files, I now get

  guile> (inexact->exact 123456789123456789.0)
  0

This might have to do with your clever continued fraction
approximation procedure.

Also:

  guile> (define pi (* 2 (acos 0)))
  guile> (- (inexact->exact pi) pi)
  3.31628058347633e-10

Shouldn't we be able to do better than this?  Double precision is good
to about 2e-16 and I would expect the fraction returned from
inexact->exact to be accurate within that margin.  (I don't know
inexact->really why I expect this, it just seems to make sense.)

What about this approach: find the integer correspondig to the bits of
the mantissa of the double number (a bignum) and then correct for the
scaling and the exponent either by multiplying with a power of two, or
building a fraction with that power of two.

I don't have the details ready yet, but I think I'll try to come up
with code for this... (when it isn't already in GMP).


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

[-- Attachment #1: Type: text/plain, Size: 819 bytes --]

Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:
>
>   (sqrt 4) => 2 (exact)
>
> I'm not sure whether exact square roots are imortant: they will only
> be useful when both the numerator and denominator are squares of
> integers and those pairs will be rare, I'd say.

I'd wondered if some sort of isqrt or isqrt+remainder would be of more
value to those wanting integer roots.  (Though I'm aware r5rs says
exact results for exact perfect squares is desirable.)

> We could use GMP for
> computing bignum square roots...

Rob mentioned this to me at one stage, I actually started on an
attempt at it.  Code below, which might sort of work, maybe.

One concern though is to ensure the root calculated is the same for a
"double" input (which uses sqrt()) as compared to some bignum root
(using gmp + conversions).


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: sqrt.c --]
[-- Type: text/x-csrc, Size: 3815 bytes --]


/* sqrt and mpz_sqrt+mpz_get_d don't round the same

   */


#include <float.h>


#ifndef DBL_MANT_DIG
#define DBL_MANT_DIG  53  /* by default assume IEEE double */
#endif


SCM scm_sqrt (SCM x);
SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0,
            (SCM x),
	    "Return the square root of @var{x}.\n"
	    "\n"
	    "If @var{x} is exact and a perfect square then the return is\n"
	    "exact, otherwise the return is inexact (either real or\n"
	    "complex).")
#define FUNC_NAME s_scm_sqrt
{
  double root_d;
  SCM    root;
  int    neg;

  if (SCM_INUMP (x))
    {
      long x_l = SCM_INUM (x);
      neg = (x_l < 0);
      x_l = (neg ? -x_l : x_l);

      /* inums which fit the mantissa of a double are handled with sqrt().
         On 32-bit systems inums always fit, but on 64-bit systems we need
         to check.  */
#if SCM_I_FIXNUM_BIT > DBL_MANT_DIG
      if ((unsigned long) x_l < (1UL << DBL_MANT_DIG))
#endif
        {
          root_d = sqrt ((double) x_l);
          if (! neg && floor (root_d) == root_d)
            {
              /* perfect square, return as inum */
              long root_l = (long) root_d;
              return SCM_MAKINUM (neg ? -root_l : root_l);
            }
          goto inexact;
        }

      /* inums which don't fit the mantissa of a double are handled with the
         bignum code, to ensure full accuracy in the result  */
      root = scm_i_mkbig ();
      mpz_set_si (SCM_I_BIG_MPZ (root), x_l);
      x = root;
      goto big;
    }
  else if (SCM_BIGP (x))
    {
      int  scale;

      root = scm_i_mkbig ();
    big:
      if (mpz_perfect_square_p (SCM_I_BIG_MPZ (x)))
        {
          /* Perfect square, return as exact.
             Enhancement: mpz_perfect_square_p already calculates the root
             when proving squareness, it'd be nice to have some sort of
             "square root if perfect square" function.  */
          mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (x));
          scm_remember_upto_here_1 (x);
          return scm_i_normbig (root);
        }

      neg = (mpz_sgn (SCM_I_BIG_MPZ (x)) < 0);
      mpz_abs (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (x));
      scm_remember_upto_here_1 (x);

      /* We need DBL_MANT_DIG bits of root, which means twice that of
         operand.  If we've got less, then pad with zeros, effectively
         generating fractional bits for the root.  If we've got more, then
         truncate so as to save work in mpz_sqrt.  */
      scale = ((long) mpz_sizeinbase (SCM_I_BIG_MPZ (root), 2)
               - 1 - 2*DBL_MANT_DIG) / 2;
      if (scale > 0)
        mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root),
                         2 * scale);
      else
        mpz_mul_2exp (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root),
                      - 2 * scale);

      mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root));
      root_d = ldexp (mpz_get_d (SCM_I_BIG_MPZ (root)), scale);
      scm_remember_upto_here_1 (root);
      goto inexact;
    }
  else if (SCM_REALP (x))
    {
      double x_d = SCM_REAL_VALUE (x);
      neg = (x_d < 0.0);
      root_d = sqrt (fabs (x_d));

    inexact:
      if (neg)
        return scm_make_complex (0.0, root_d);
      else
        return scm_make_real (root_d);
    }
  else if (SCM_COMPLEXP (x))
    {
      /* sqrt(x) = (make-polar (sqrt (magnitude x)) (/ (angle x) 2))
         ENHANCE-ME: Use csqrt() when available.  */
      double re = SCM_COMPLEX_REAL (x);
      double im = SCM_COMPLEX_IMAG (x);
      double mag = sqrt (hypot (re, im));
      double ang = atan2 (re, im) / 2.0;
      double s, c;
#if HAVE_SINCOS
      sincos (ang, &s, &c);
#else
      s = sin (ang);
      c = cos (ang);
#endif
      return scm_make_complex (mag * c, mag * s);
    }
  else
    SCM_WRONG_TYPE_ARG (SCM_ARG1, x);
}
#undef FUNC_NAME

[-- Attachment #3: Type: text/plain, Size: 142 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:
>
> 	    "Formally, the function returns an integer equivalent to\n"
> 	    "@code{(inexact->exact (floor (* @var{n} (expt 2 @var{cnt}))))}.\n"

I think I removed words like that from manual, on the basis that they
weren't as clear as they could be.

> guile> (ash .1 2)
> 0.4

I guess that's only a side effect of using scm_product.  It really
shouldn't be doing that of course, it ought to be bit shifting not
multiplying.

No doubt ash could work on flonums in a sensible way, but maybe
there's no particular need for that, if other bitwise logical
functions (ie. everything under "Bitwise Operations" in the manual)
only take integers.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:
>
> I don't have the details ready yet, but I think I'll try to come up
> with code for this... (when it isn't already in GMP).

gmp has an mpq_set_d, which does what you imagine, extract the bits to
make a binary fraction.  The actual implementation is uglified by the
details of gmp internals though.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Bill Schottstaedt]

 > With your new files, I now get
 >  guile> (inexact->exact 123456789123456789.0)
 >  0

I noticed this, but wasn't sure how to proceed;
in the previous guile, you'd get (in a sense) equally
bogus results:

guile> (inexact->exact 17452826108659293487.0)
17452826108659294208

guile> (= (truncate 17452826108659293487.3) 17452826108659293487.0)
#f

I wasn't even sure whether you wanted inexact->exact to
be changed -- I left the old code in place.

I wrote a gmp version of the continued fraction code about
10 years ago -- I think I can still find it.

 >  guile> (define pi (* 2 (acos 0)))
 >  guile> (- (inexact->exact pi) pi)
 >  3.31628058347633e-10
 > Shouldn't we be able to do better than this?

I think so; but in the current version, I'm setting the
minimum error to 1/INT_MAX, which looks like it's in the
ballpark of e-10.



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Bill Schottstaedt]

If Guile wants to support "bignum floats" so to speak, I
think it could use gmp floats -- this strikes me as a big
task; I don't think its current dependence on doubles
and double->long conversions can be pushed very far.

I'll add a check in scm_rationalize for an error whose
inverse can't fit in a long and switch over to gmp in
that case, but I think the 0 error case should
use the simple version -- if someone wants pi to 10000000
decimals, he should ask for it explicitly.



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Kevin Ryde <user42@zip.com.au> writes:

> Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:
>>
>> I don't have the details ready yet, but I think I'll try to come up
>> with code for this... (when it isn't already in GMP).
>
> gmp has an mpq_set_d, which does what you imagine, extract the bits to
> make a binary fraction.  The actual implementation is uglified by the
> details of gmp internals though.

Yes, that's what I was looking for.  I now have in scm_inexact_to_exact

  else if (SCM_REALP (z))
    {
      mpq_t frac;
      SCM q;
      mpq_init (frac);
      mpq_set_d (frac, SCM_REAL_VALUE (z));
      q = scm_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
			  scm_i_mpz2num (mpq_denref (frac)));
      mpq_clear (frac);
      return q;
    }

using the helper

    static SCM_C_INLINE_KEYWORD SCM
    scm_i_mpz2num (mpz_t b)
    {
      /* convert a mpz number to a SCM number. */
      if (mpz_fits_slong_p (b))
        {
          long val = mpz_get_si (b);
          if (SCM_FIXABLE (val))
            return SCM_MAKINUM (val);
        }

      {
        SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
        mpz_init_set (SCM_I_BIG_MPZ (z), b);
        return z;
      }
    }

Kevin, does this look OK from a GMP point of view?  I.e., no memory
leaks, etc?


That code gives perfect accuracy but that can be strange as well:

    guile> (exact->inexact 1/1000)
    0.001
    guile> (inexact->exact (exact->inexact 1/1000))
    1152921504606847/1152921504606846976
    guile> (rationalize (exact->inexact 1/1000) 1e-16)
    1/1000

The above behavior of inexact->exact better fits R5RS, I'd say, since
R5RS calls for the _closest_ representable exact number, not for the
simplest, as with rationalize.  So I'm inclined to use it intead of
rationalize.  Ok?


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Bill Schottstaedt <bil@ccrma.Stanford.EDU> writes:

>  > With your new files, I now get
>  >  guile> (inexact->exact 123456789123456789.0)
>  >  0
>
> I noticed this, but wasn't sure how to proceed;
> in the previous guile, you'd get (in a sense) equally
> bogus results:
>
> guile> (inexact->exact 17452826108659293487.0)
> 17452826108659294208

Hmm, I'd say that this is not nearly as bogus.  It is not very
accurate (since doubles do not have that much precision), but much
better than zero.  (Did I use 'accurate' and 'precise' correctly? :)

> I wasn't even sure whether you wanted inexact->exact to
> be changed -- I left the old code in place.

Yes, it needs to be changed.  People who have assumed that
inexact->exact always returns an integer need to think again (I'm
afraid I'm one of them...).

> I wrote a gmp version of the continued fraction code about 10 years
> ago -- I think I can still find it.

Hmm, what about a SCM version of the code?  I'm about to write such a
thing, but I first needed a proper scm_floor...

>  >  guile> (define pi (* 2 (acos 0)))
>  >  guile> (- (inexact->exact pi) pi)
>  >  3.31628058347633e-10
>  > Shouldn't we be able to do better than this?
>
> I think so; but in the current version, I'm setting the
> minimum error to 1/INT_MAX, which looks like it's in the
> ballpark of e-10.

The SCM version of rationalize can work to arbitrary error margins
(hopefully), by using generic SCM arithmetic instead of the longs and
doubles that your version uses.  That will make rationalize a bit
slower, but hey.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Bill Schottstaedt]

 > The above behavior of inexact->exact better fits R5RS, I'd say, since
 > R5RS calls for the _closest_ representable exact number, not for the
 > simplest, as with rationalize.  So I'm inclined to use it intead of
 > rationalize.  Ok?

I didn't notice the distinction.  I'd go with your version.

 > Hmm, what about a SCM version of the code?  I'm about to write such a
 > thing, but I first needed a proper scm_floor...

Yes, that's a better approach.



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Rob Browning]

Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:

> Yes, that's what I was looking for.  I now have in scm_inexact_to_exact
>
>   else if (SCM_REALP (z))
>     {
>       mpq_t frac;
>       SCM q;
>       mpq_init (frac);
>       mpq_set_d (frac, SCM_REAL_VALUE (z));

Do you need a scm_remember_upto_here1(z) here, or is z already
protected?

>       q = scm_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
> 			  scm_i_mpz2num (mpq_denref (frac)));

Also, can scm_make_ratio throw?  If so, then I'm wondering if frac
might be a potential memory leak, in which case, one solution would be
to just allocate frac via scm_i_mkbig().

>       mpq_clear (frac);
>       return q;
>     }

-- 
Rob Browning
rlb @defaultvalue.org and @debian.org; previously @cs.utexas.edu
GPG starting 2002-11-03 = 14DD 432F AE39 534D B592  F9A0 25C8 D377 8C7E 73A4


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Rob Browning]

Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:

>> I wasn't even sure whether you wanted inexact->exact to
>> be changed -- I left the old code in place.
>
> Yes, it needs to be changed.  People who have assumed that
> inexact->exact always returns an integer need to think again (I'm
> afraid I'm one of them...).

Yep, I'm sure there are plenty of people in that situation.  I suspect
it's been used heavily in places where a (C) function require an
integer argument, i.e. (foo (inexact->exact (round x))).  Though as
long as (inexact->exact (round x)) is guaranteed to return an integer,
perhaps that's sufficient.

-- 
Rob Browning
rlb @defaultvalue.org and @debian.org; previously @cs.utexas.edu
GPG starting 2002-11-03 = 14DD 432F AE39 534D B592  F9A0 25C8 D377 8C7E 73A4


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Marius Vollmer]

Rob Browning <rlb@defaultvalue.org> writes:

> Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:
>
>> [...] People who have assumed that
>> inexact->exact always returns an integer need to think again (I'm
>> afraid I'm one of them...).
>
> Yep, I'm sure there are plenty of people in that situation.  I suspect
> it's been used heavily in places where a (C) function require an
> integer argument, i.e. (foo (inexact->exact (round x))).  Though as
> long as (inexact->exact (round x)) is guaranteed to return an integer,
> perhaps that's sufficient.

(inexact->exact (round x)) is OK, but just (inexact->exact x) and
assuming that inexact->exact implicitely rounds to an integer will
break.


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: ratio implementation

[Kevin Ryde]

Marius Vollmer <marius.vollmer@uni-dortmund.de> writes:
>
>       mpq_init (frac);
>       mpq_set_d (frac, SCM_REAL_VALUE (z));
>       q = scm_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
> 			  scm_i_mpz2num (mpq_denref (frac)));
>       mpq_clear (frac);

Yep, that's about it.

I guess the num and den from the mpq_t could be copied directly into
bignums and used as mpz_t's from then on, rather than initializing new
mpz_t's.

Just have to think a bit whether that would be compatible going
forward.  We already allow arbitrary mpz calculations on the two
parts, just have to decide if that can include clearing :).


_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

ash using shifts (was: ratio implementation)

[Kevin Ryde]

[-- Attachment #1: Type: text/plain, Size: 256 bytes --]

Speaking of ash, the bit of code below is what I'm looking at for
shifts instead of mul or div.  Untested as yet, but hopefully not too
far wrong.  Anyone who wants to take it forward should feel free ...
or kick it into touch and do something different.


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: ash.c --]
[-- Type: text/x-csrc, Size: 3395 bytes --]

SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
            (SCM n, SCM cnt),
	    "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
	    "if @var{cnt} is negative.  This is an ``arithmetic'' shift.\n"
	    "\n"
	    "This is effectively a multiplication by @m{2^{cnt},\n"
	    "2^@var{cnt}}, and when @var{cnt} is negative it's a division,\n"
	    "rounded towards negative infinity.  (Note that this is not the\n"
	    "same rounding as @code{quotient} does.)\n"
	    "\n"
	    "With @var{n} viewed as an infinite precision twos complement,\n"
	    "@code{ash} means a left shift introducing zero bits, or a right\n"
	    "shift dropping bits.\n"
	    "\n"
	    "@lisp\n"
	    "(number->string (ash #b1 3) 2)     @result{} \"1000\"\n"
	    "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
	    "\n"
	    ";; -23 is bits ...11101001, -6 is bits ...111010\n"
	    "(ash -23 -2) @result{} -6\n"
	    "@end lisp")
#define FUNC_NAME s_scm_ash
{
  long  bits_to_shift;
  SCM_VALIDATE_INUM (2, cnt);
  bits_to_shift = SCM_INUM (cnt);

  if (SCM_INUMP (y)) {
    long in = SCM_INUM (n);

    if (bits_to_shift > 0)
      {
        /* Left shift of more than SCM_I_FIXNUM_BIT-1 will certainly
           overflow a non-zero fixnum.  For smaller shifts we check the bits
           going into positions above SCM_I_FIXNUM_BIT-1.  If they're all 0s
           for in>=0, or all 1s for in<0 then there's no overflow.  Those
           bits are "in >> (SCM_I_FIXNUM_BIT-1 - bits_to_shift".  */

        if (in == 0)
          return n;

        /* FIXME: This relies on signed right shifts being arithmetic, which
           is not guaranteed by C99.  */
        if (bits_to_shift < SCM_I_FIXNUM_BIT-1
            && ((unsigned long) ((in >> (SCM_I_FIXNUM_BIT-1 - bits_to_shift))
                                 + 1 <= 1)))
          {
            return SCM_MKINUM (in << bits_to_shift);
          }
        else
          {
            SCM result = scm_i_long2big (z);
            mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
                          bits_to_shift);
            return result;
          }
      }
    else
      {
        bits_to_shift = -bits_to_shift;
        if (bits_to_shift >= LONG_BIT)
          return (in >= ? 0 : -1);
        else
          {
            /* FIXME: This relies on signed right shifts being arithmetic,
               which is not guaranteed by C99.  */
            return SCM_MKINUM (in >> bits_to_shift);
          }
      }

  } else if (SCM_BIGP (n)) {
    SCM result;

    if (bits_to_shift == 0)
      return n;

    result = scm_i_mkbig ();
    if (bits_to_shift >= 0)
      {
        mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
                      bits_to_shift);
        return result;
      }
    else
      {
        /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
           we have to allocate a bignum even if the result is going to be a
           fixnum.  We could detect the case of bits_to_shift being so big
           as to leave us with only 0 or -1, and avoid allocating a bignum,
           but that doesn't seem worth worrying about.  */
        mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
                         -bits_to_shift);
        return scm_i_normbig (result);
      }

  } else {
    SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
  }
}
#undef FUNC_NAME

[-- Attachment #3: Type: text/plain, Size: 142 bytes --]

_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel

Re: fractions.test

[Bill Schottstaedt]

Only the bignum cases after the comment about clisp (around
line 346 in my version) are from clisp; I wrote the rest by hand,
although I looked at the plt-scheme test suite for ideas.
I also looked at Dybvig's book and borrowed some of his
examples (they are marked "from Dybvig").  I'm sorry I
haven't translated that file to "the Guile Way" -- no time.



_______________________________________________
Guile-devel mailing list
Guile-devel@gnu.org
http://mail.gnu.org/mailman/listinfo/guile-devel