[racket] Math library ready for testing
Neil Toronto
for wider testing. Almost everything ready for use is documented, the
tests keep passing, and everything *seems* to work.
We (the Racket development team) need people to use the heck out of it,
to find the majority of the remaining logic and design errors before the
next release. This is a big chunk of new code - ~20k lines in ~100
nontrivial files - so the more eyes, the better.
If all you have time for is checking documentation, we'll take it. The
latest docs are here:
http://pre.racket-lang.org/docs/html/math/
The nightly builds are here:
http://pre.racket-lang.org/installers/
For Windows and Mac, the nightlies should include two new libraries used
by `math/bigfloat': libgmp and libmpfr. Linux users almost certainly
already have them; many packages (e.g. gcc) depend on libmpfr.
If you prefer to build from source, clone the git repository here:
If you're building on Windows or Mac, you can download pre-built libgmp
and libmpfr from here:
http://download.racket-lang.org/libs/10/
************
Once you have a pre-release build, do (require math), and away you go!
Well, it's probably not obvious where to start. The following are the
modules that `math' re-exports, and for some, what needs testing. Please
feel free to skim. In decreasing order of attention they need:
(require math/array)
NumPy/SciPy/Repa-like arrays, Typed-Racket-style. Provides a
covariant `Array' type and several mutable subtypes, array literals,
functional constructors, pointwise operations with automatic
broadcasting in SciPy/NumPy or permissive mode, slicing, indexing
tricks, folds, easy arbitrary transformations, and parallel FFT.
This one needs a lot of usability testing. In particular...
We're trying something that's been done only once before, in
Haskell's Repa: by default, arrays' elements are *non-strict*,
meaning that they're recomputed when referenced. Cool things about
this: pretty much everything can in principle be parallelized, and
most intermediate arrays use almost no memory. But it requires more
thought from users about when arrays should be made strict (computed
and stored all at once). If it totally sucks, we can change it. If
it's a worthwhile imposition, we'll leave it. If it's a wash, we
could parameterize the default behavior.
Also, I've been considering allowing negative axis numbers and row
indexes. Need feedback on how helpful it would be vs. confusing and
error-hiding.
Lastly, I'm looking for a dual of `array-axis-reduce' (a kind of
fold) that makes it easy to write `list-array->array', the inverse of
the existing `array->list-array'. If you're into functional design
puzzles, give this one a shot.
(require math/bigfloat)
Floating-point numbers with arbitrarily large precision and
large exponents. Also elementary and special functions, whose results
are proved to be correct in many theses and dissertations. This is a
Racket interface to MPFR (www.mpfr.org).
This module needs platform-specific testing.
I think we've verified that `math/bigfloat' works on all the
supported 64-bit platforms, but not whether it works on supported
32-bit platforms. It should. ;)
There's an error in the current nightly build, which causes an
infinite loop when libmpfr isn't installed and a bigfloat function is
applied. I just pushed a fix; it should fail properly tomorrow.
(require math/distributions)
Probability distributions: discrete, integer-valued and real-valued.
Distribution objects can compute pdfs, cdfs and inverse cdfs,
optionally in log space, and for cdfs and inverse cdfs, optionally
for upper tails. They can also generate samples.
Design ideas are welcome. More distributions are planned. Let me know
which you need, and I'll concentrate on those first.
Watch out, R. Our gamma cdf is more accurate.
(require math/special-functions)
Non-elementary functions like gamma, erf, zeta, Lambert W, and
incomplete gamma and beta integrals. Most of these should be fairly
accurate; i.e. they compute answers with apparently < 5 ulps error
for the majority of their domains. But floating-point domains are
huge, so the more use, the better.
(require math/number-theory)
Number theory! Chinese Remainder solutions, coprimality and primality
testing, factorization, congruence arithmetic parameterized on a
current modulus, integer nth roots, multiplicative functions, common
number sequences (Bernoulli, Eulerian, tangent, etc.), combinatorics,
primitive roots, and more.
Please thank Jens Axel Søgaard, who wrote this module, for his
excellent work.
(require math/statistics)
Functions to compute summary values for collections of data, and to
manage weighted samples. Every function for which it makes sense
accepts weighted values, including the O(1)-space running statistics
update.
I'm in the middle of documenting this. I'll get around to documenting
correlation, kth-value order statistics (via quickselect), and
counting/binning at least by Dec. 18.
(require math/flonum)
Re-exports `racket/flonum' along with many more floating-point
goodies. If you're a floating-point nut, you'll love it. If you're
not, at least check out `flsum'. If you're doing statistics, look
into log-space arithmetic (`lg+' and friends).
(require math/base)
Re-exports `racket/math'; also exports some more constants like
gamma.0 (Euler-Mascheroni), a bit more float-complex support, inverse
hyperbolic functions, large random numbers, error measurement.
Lastly, there's `math/matrix', which is currently not re-exported by
`math' because I haven't reviewed it yet. This is more from Jens Axel,
and if it's like his other fine work, it's correct and well-tested. I'll
probably get to it by Christmas.
To sum up, there are 8 modules that need pounding on, and we need YOU to
do some pounding. If you're not terribly busy, please pick one that
looks interesting/fun/useful, and try it out.
Thanks!
Neil ⊥
____________________
Racket Users list:
http://lists.racket-lang.org/users
Laurent
Neil Toronto
I should make that list and some others on my to-do list public, like
the "planned special functions" list. Where's a good place for that?
Neil ⊥
On 12/07/2012 01:58 PM, Ray Racine wrote:
> Amazing. Tons of hard work. Thank you. I'm going to use it asap.
> Grab a beer and if you're making a list, please add Dirichlet
> Distribution to the Phase 2 list.
>
>
> On Fri, Dec 7, 2012 at 3:33 PM, Neil Toronto <neil.t...@gmail.com
> <mailto:neil.t...@gmail.com>> wrote:
>
> I've just pushed the last commits that make the new math library
> ready for wider testing. Almost everything ready for use is
> documented, the tests keep passing, and everything *seems* to work.
>
> We (the Racket development team) need people to use the heck out of
> it, to find the majority of the remaining logic and design errors
> before the next release. This is a big chunk of new code - ~20k
> lines in ~100 nontrivial files - so the more eyes, the better.
>
> If all you have time for is checking documentation, we'll take it.
> The latest docs are here:
>
> http://pre.racket-lang.org/__docs/html/math/
> <http://pre.racket-lang.org/docs/html/math/>
>
> The nightly builds are here:
>
> http://pre.racket-lang.org/__installers/
> <http://pre.racket-lang.org/installers/>
>
> For Windows and Mac, the nightlies should include two new libraries
> used by `math/bigfloat': libgmp and libmpfr. Linux users almost
> certainly already have them; many packages (e.g. gcc) depend on libmpfr.
>
> If you prefer to build from source, clone the git repository here:
>
> https://github.com/plt/racket
>
> If you're building on Windows or Mac, you can download pre-built
> libgmp and libmpfr from here:
>
> http://download.racket-lang.__org/libs/10/
> Racket interface to MPFR (www.mpfr.org <http://www.mpfr.org>).
> http://lists.racket-lang.org/__users
> <http://lists.racket-lang.org/users>
Asumu Takikawa
> I should make that list and some others on my to-do list public,
> like the "planned special functions" list. Where's a good place for
> that?
https://github.com/plt/racket/wiki
Cheers,
Asumu
Ray Racine
Neil Toronto
> On 2012-12-07 14:31:36 -0700, Neil Toronto wrote:
>> I should make that list and some others on my to-do list public,
>> like the "planned special functions" list. Where's a good place for
>> that?
>
> The github wiki would be good. We already have some stuff there:
> https://github.com/plt/racket/wiki
Excellent, thanks. I've just started "Math Library Features" here:
https://github.com/plt/racket/wiki/Math-Library-Features
Neil ⊥
Hendrik Boom
>
> Also, I've been considering allowing negative axis numbers and row
> indexes. Need feedback on how helpful it would be vs. confusing and
> error-hiding.
Torrix, a library for Algol 68. The convention that they used there was
that matrices had infinitely many rows and columns, indexed by the
integers, but that all but a finite number of entries were zero. This
operations on what the user thought of as 7x7 matrices acted as if they
were infinite-extent matrices instead. The difference came in the ease
which which they could handle matrices of different sizes (even with
different index origins) in a semantically simple way. That apparently
was appreciated by its users, though I might have considered some of
these machinations to be usefully reported as type violations.
-- hendrik
Neil Toronto
Oh, Typed Racket would never let us get away with that, except with the
FlArray and FCArray subtypes. Unless every array somehow knew what its
"default" or "zero" value was, or every (Array A) was actually an (Array
(U Zero A)). The latter approach is truly icky, and I don't like the
former much, either.
FWIW, the `flomap' type from `images/flomap' does this. With color
components premultiplied by alpha, it makes all the image operations
Just Work near the borders without treating them as special cases.
Anyway, I should have been more clear. Negative indexes would count
backwards from the size of whatever they're indexing. For a
5-dimensional array, axis number -1 would be the same as axis number 4.
For an array with 10 rows, row -1 would be the same as row 9. It can be
useful to write code that says, "Give me the last row; I don't care how
many there are," or "fold over the second-to-last axis."
I'm more sold on doing this for axis numbers, but still waffling. It's
an irreversible decision. For row indexes, it matches perfectly with
permissive broadcasting rules, but those aren't the default rules.
Neil ⊥
Michael Wilber
your work, Neil.
During winter break once I'm finished with grad school apps, I'll
convert my final Signals and Systems Matlab project to racket/math to
get a feel for the library.
One note off-hand: I really don't like the (array #[#[#[...]]]) syntax.
In Matlab, a 2D array is just:
[1 2; 3 4]
That's 10 characters. In Python, the equivalent array is (including whitespace):
[[1, 2], [3, 4]]
But in Racket, the equivalent array is:
(array #[#[1 1] #[3 4]])
This syntax, especially the visual weight of all the hash marks #[#[,
makes it harder for me to visually scan the line and see the actual
array values. To see this, note that I miscopied one of the numbers
when I wrote the racket array. How many seconds did it take you to find
my mistake? Now imagine doing that at 3AM when drunk and you get an idea
of my reason for saying this. This is the kind of syntax that will trip
up undergraduates who are trying to debug their programs.
Another syntax might be something akin to vectors: #[[1 2] [3 4]]
Or just (array (1 2) (3 4))
Or maybe using a : or something to separate rows for 2D arrays?
(array 1 2 : 3 4)
(array 1 0 0 0 0 0
: 0 1 0 0 0 0
: 0 0 1 0 0 0
: 0 0 0 1 0 0
: 0 0 0 0 1 0
: 0 0 0 0 0 1)
All of these feel much less "racket-y," but they would eliminate all the
visual cruft.
Robby Findler
Danny Yoo
Or maybe using a : or something to separate rows for 2D arrays?
(array 1 2 : 3 4)
(array 1 0 0 0 0 0
: 0 1 0 0 0 0
: 0 0 1 0 0 0
: 0 0 0 1 0 0
: 0 0 0 0 1 0
: 0 0 0 0 0 1)
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1)
Neil Toronto
This works right now, and is pretty close:
(make-mutable-array
#(6 6)
#(1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1))
It gets a little uglier when you want to override TR's generalization
for polymorphic, mutable types. For example, you'd have to apply (inst
make-mutable-array Byte) to get a (Mutable-Array Byte) instead of a
(Mutable-Array Integer). No problem in untyped Racket, of course.
There's no shape check on the data, though, so it'd be really easy to
make an array with values in the wrong places.
As far as using (array [[...]]) - no hashes - it used to work like that.
The `array' macro would check the shape of the parens to determine
whether (e ...) was meant to be a row like [list 1 2 3] or an element
like (list 1 2 3). The problem is, that's a syntax property, which
macros rarely preserve, and Scribble couldn't render them properly. In
general, it would be too easy for array literals computed or transformed
at syntax time to lose their shape and meaning.
Using commas, dots and semicolons is out of the question. Colons would
clash with TR's use in types. Other punctuation is sinfully ugly. I'm
afraid the only real option for something with a little less visual
clutter is a reader extension or a custom language. I considered it, but
in the end I rejected #lang math. I just didn't have the hubris for it. :D
Neil Toronto
> This is such an exciting time to be a Racket programmer! Thanks for all
> your work, Neil.
You're welcome!
You're also right, it is an exciting time. I wouldn't have even
attempted the distributions and special functions without Typed Racket's
optimizer and recent JIT improvements. Or I might have, using nothing
but unsafe operations, and driven myself insane. :D
So I'm passing some of that thanks on to Matthew, Sam and Vincent.
(Oh, and bigfloats were critical to making the functions accurate. Also
plot. So thanks again to Matthew, and to Eli, for the FFI, to the French
floating-point geeks for MPFR, and again to myself for plot, and Matthew
again for racket/draw. Go us. Go us all!)
> During winter break once I'm finished with grad school apps, I'll
> convert my final Signals and Systems Matlab project to racket/math to
> get a feel for the library.
Oooh neet. I'm looking forward to your feedback. Be brutally honest. And
stay sober! Or at least have a designated coder.
Justin R. Slepak
---
Justin Slepak
PhD student, Computer Science dept.
Noel Welsh
I've just pushed the last commits that make the new math library ready for wider testing. Almost everything ready for use is documented, the tests keep passing, and everything *seems* to work.
We (the Racket development team) need people to use the heck out of it, to find the majority of the remaining logic and design errors before the next release. This is a big chunk of new code - ~20k lines in ~100 nontrivial files - so the more eyes, the better.
Jens Axel Søgaard
>> Or maybe using a : or something to separate rows for 2D arrays?
>> (array 1 2 : 3 4)
The matrix library uses arrays to represent the matrices.
An 2x3 matrix can be constructed as:
(matrix/dim 2 3
1 2 3
4 5 6)
(list->matrix '[[1 2 3] [4 5 6]])
(for/matrix: : Number 2 3 ([i (in-naturals)]) i)
These expressions will all return the same 2D array.
Given that the shape of a matrix is fixed, the problems
with square parentheses in constructor notation disappear.
That is, one could let matrix be a macro, and it could
accept both kinds of parentheses:
(matrix [[1 2 3] [4 5 6]])
(matrix ((1 2 3) (4 5 6)))
See more matrix examples here:
https://github.com/plt/racket/blob/master/collects/math/tests/matrix-tests.rkt
--
Jens Axel Søgaard
Martin Neal
Neil Toronto
> The matrix library uses arrays to represent the matrices.
> An 2x3 matrix can be constructed as:
>
> (matrix/dim 2 3
> 1 2 3
> 4 5 6)
>
> (list->matrix '[[1 2 3] [4 5 6]])
>
> (for/matrix: : Number 2 3 ([i (in-naturals)]) i)
>
> These expressions will all return the same 2D array.
>
> Given that the shape of a matrix is fixed, the problems
> with square parentheses in constructor notation disappear.
> That is, one could let matrix be a macro, and it could
> accept both kinds of parentheses:
>
> (matrix [[1 2 3] [4 5 6]])
> (matrix ((1 2 3) (4 5 6)))
I wouldn't mind having a special constructor like that for matrices.
BTW, I've just pushed the documentation for `math/statistics' v1.0, so
it's ready for testing. That also means I'll be reviewing `math/matrix'
soon. Woo hoo!
Neil ⊥
Jens Axel Søgaard
> First - I love the library. There's so much functionality that it's going
> to take quite some time to explore it all
>
> I'd love to see an `in-primes` function added to math/number-theory.
Excellent idea. I have added it to the todo-list.
You make an excellent point. The performance of prime? needs to
be improved for small primes.
Currently prime? works as follows:
For a small (less than 10000) number n prime? looks up primality of
with a simple (vector-ref small-primes n).
For a large number (greater than 10000) a pseudo prime test is used.
Originally the limit for small primes was a million, but that limit is too
large for a general purpose library.
The first change will be to replace the vector with a bit-vector
in order to store larger table without using too much space.
The second will be to experiment with "medium numbers"
and see whether trial division is faster than the pseudo
prime test, when the small prime limit is raised to larger number
again.
Here is how I plan to use bit-vectors:
To test whether a small number n is prime,
step one is to check whether it is divisible
by 2, 3 or 5. If not then if the number is small
then its primality is looked up in a bit-vector.
Below 60 there are only 16 numbers not divisible by 2, 3 or 5.
The number of bytes needed to store primality results for the
first million numbers will therefore be
1000000 * 16/60 * 1/8 ~ 33333.3
How much space would be reasonable to use in the library?
Robby Findler
but my apriori guess would be to say that 32k isn't too bad, but that
it also doesn't seem necessary in this case. Martin Neal's micro
benchmark already saw lots of improvements with just a different
algorithm so perhaps that's enough?
Robby
Pierpaolo Bernardi
<jens...@soegaard.net> wrote:
> Below 60 there are only 16 numbers not divisible by 2, 3 or 5.
> The number of bytes needed to store primality results for the
> first million numbers will therefore be
> 1000000 * 16/60 * 1/8 ~ 33333.3
>
> How much space would be reasonable to use in the library?
Would make sense to have a user-tuneable parameter (with a small
default) for this?
So a prime intensive computation could do (precompute-primes-up-to
(expt 10 9) ).
What do Mathematica, Maple, etc, do in this case?
P.
Jens Axel Søgaard
> On Tue, Dec 11, 2012 at 12:01 PM, Jens Axel Søgaard
> <jens...@soegaard.net> wrote:
>
>> Below 60 there are only 16 numbers not divisible by 2, 3 or 5.
>> The number of bytes needed to store primality results for the
>> first million numbers will therefore be
>> 1000000 * 16/60 * 1/8 ~ 33333.3
>>
>> How much space would be reasonable to use in the library?
>
> Would make sense to have a user-tuneable parameter (with a small
> default) for this?
> So a prime intensive computation could do (precompute-primes-up-to
> (expt 10 9) ).
That would make sense.
> What do Mathematica, Maple, etc, do in this case?
Here is what Maxima does (line 628 and below):
https://github.com/andrejv/maxima/blob/master/src/ifactor.lisp
They keep a table of the primes below 10000. Small numbers
are simply looked up in the table.
Medium numbers are divided into ranges. In each range
Miller-Rabin is used with pre-selected moduli. This makes
the test deterministic.
For large numbers a pseudo test is used.
--
Jens Axel Søgaard
Stephen Bloch
For a small (less than 10000) number n prime? looks up primality of
with a simple (vector-ref small-primes n).
For a large number (greater than 10000) a pseudo prime test is used.
Originally the limit for small primes was a million, but that limit is too
large for a general purpose library.
Jens Axel Søgaard
Sam Tobin-Hochstadt
<jens...@soegaard.net> wrote:
>
> The second will be to experiment with "medium numbers"
> and see whether trial division is faster than the pseudo
> prime test, when the small prime limit is raised to larger number
> again.
One thing that might be worth considering, esp. in the context of
`in-primes`, might be a stateful version of `prime?` that keeps
previously-found primes around to improve the performance of trial
division.
Sam
Jens Axel Søgaard
> 2012/12/11 Jens Axel Søgaard <jens...@soegaard.net>:
>>
>> https://github.com/andrejv/maxima/blob/master/src/ifactor.lisp
>
> tester is called new-prime? and the old one prime?.
> The old one seems to be faster on the simple
> "sum primes below a million test".
>
> (time (for/sum ([n (in-range 1000000)]
> #:when (new-prime? n))
> n))
>
> ; cpu time: 31264 real time: 31427 gc time: 340
> ; 37550402023
>
> (time (for/sum ([n (in-range 1000000)]
> #:when (prime? n))
> n))
>
> ; cpu time: 8944 real time: 9079 gc time: 257
> ; 37550402023
>
> --
> Jens Axel Søgaard
Sam Tobin-Hochstadt
means that we can't easily just take their code and add it to Racket.
Sam
Jens Axel Søgaard
> Would it perhaps make more sense for small-primes to contain primes
> themselves, in increasing order so one can be found by binary search, rather
> than booleans? The O(1) behavior would be replaced by O(log(limit)), but
> perhaps you would save enough memory to put the limit higher.
I think there are too many primes.
Since
> (require math)
> (nth-prime 78498)
1000003
there are 78497 primes below a million. On a 64 bit machine
that requires 8*78497 = 627976 bytes.
Stephen Bloch
On Dec 11, 2012, at 9:03 AM, Jens Axel Søgaard wrote:
> 2012/12/11 Stephen Bloch <bl...@adelphi.edu>:
>
>> Would it perhaps make more sense for small-primes to contain primes
>> themselves, in increasing order so one can be found by binary search, rather
>> than booleans? The O(1) behavior would be replaced by O(log(limit)), but
>> perhaps you would save enough memory to put the limit higher.
>
> I think there are too many primes.
>
> Since
>
>> (require math)
>> (nth-prime 78498)
> 1000003
>
> there are 78497 primes below a million. On a 64 bit machine
> that requires 8*78497 = 627976 bytes.
How many primes are below ten million? A hundred million? At some point storing the primes will take less memory than storing primality flags, but that point may be above the size of tables we can realistically store today.
Wait: it's conceivable that there is no such crossing point. As the numbers get big, it takes O(log(limit)) bits to store numbers less than limit. The number of primes less than limit is Theta(limit/log(limit)), so storing them all takes Theta(limit) space, asymptotically the same as storing the flags.
Stephen Bloch
sbl...@adelphi.edu
Phil Bewig
Martin Neal
Phil Bewig
Pierpaolo Bernardi
> How many primes are below ten million? A hundred million? At some point storing the primes will take less memory than storing primality flags, but that point may be above the size of tables we can realistically store today.
(printf "~a: ~a~n" i (boolean-vector-count (boolean-sieve (expt 10 i)))))
3: 168
4: 1229
5: 9592
6: 78498
7: 664579
8: 5761455
(some values may be higher than the right number, as my boolean-sieve
function rounds up the limit to the next 8 multiple)
P.
Jens Axel Søgaard
> I can't type. The correct count of primes less than a million is 78498. And
> there is something wrong if (nth-prime 78498) returns 100003, as the 78498th
> prime is 999983.
That's because nth-prime counts from 0 and Mathematica counts from 1.
--
Jens Axel Søgaard
Jens Axel Søgaard
I have just pushed an improvement for prime? on
small numbers.
Before:
>> (time (for/sum ([n (in-range 1000000)]
>> #:when (prime? n))
>> n))
>>
>> ; cpu time: 8944 real time: 9079 gc time: 257
>> ; 37550402023
Now:
cpu time: 1511 real time: 1541 gc time: 87
37550402023
https://github.com/plt/racket/commit/e5016951d06b87104317d4c0b499ddb769d24e20
Hendrik Boom
> On Tue, Dec 11, 2012 at 5:11 PM, Stephen Bloch <sbl...@adelphi.edu> wrote:
>
>
> > How many primes are below ten million? A hundred million? At some
> > point storing the primes will take less memory than storing
> > primality flags, but that point may be above the size of tables we
> > can realistically store today.
might find that the increasing sparseness of primes is eaten up by the
increasingly larger number of digits used to express the primes. I
don't know how mich gzipping the lise would help.
-- hendrik
Pierpaolo Bernardi
vector of the differences between consecutive primes. Storing a byte
per prime goes a long way.
Let's see... there are 50847534 primes under 10^9, the max difference
between a consecutive pair of them is 282. We can store that value /2,
so the max number we need to store is 141, which fits in a byte.
So this table would need 48MB in this representation, and (/ (expt 10
9) 8) ==> 119M bytes as a bit-vector.
Just another possibility...
Cheers
P.
Pierpaolo Bernardi
> I've just pushed the last commits that make the new math library ready for
> wider testing. Almost everything ready for use is documented, the tests keep
> passing, and everything *seems* to work.
(max-math-threads) → Positive-Integer
(max-math-threads num) → void?
num : Positive-Integer
The maximum number of threads a parallelized math function will use.
The default value is (max 1 (processor-count))."
Isn't (max 1 (processor-count)) the same as (processor-count) ?
Or does Racket runs on machines with less than 1 processors? 8^)
J. Ian Johnson
-Ian
----- Original Message -----
From: "Pierpaolo Bernardi" <olop...@gmail.com>
To: "Neil Toronto" <neil.t...@gmail.com>
Cc: us...@racket-lang.org
Sent: Wednesday, December 12, 2012 9:17:57 AM GMT -05:00 US/Canada Eastern
Subject: Re: [racket] Math library ready for testing
Robby Findler
Robby
Neil Toronto
Neil ⊥
Pierpaolo Bernardi
one, which I assume there will be a bunch in the new docs?
Better a private mail to Neil?
P.
Neil Toronto
If you're not sure whether some behavior is an error, email the list.
Even if you're wrong, you're not wrong: if you misunderstood something,
it's probably the documentation's fault. Also, other people probably
misunderstood it in the same way, so resolving it in public is good.
Otherwise, submit a bug report. Here, the word "bug" means "bug the
person responsible," which ensures the *error* will get fixed. :D
Neil ⊥
Harry Spier
> The default value is (max 1 (processor-count))."
Pierpaolo Bernardi
====
flulp-error
"For non-rational arguments such as +nan.0, flulp-error returns 0.0 if
(eqv? x r); otherwise it returns +inf.0."
However one example is:
> (flulp-error +inf.0 +nan.0)
+nan.0
Which does not agree with the doc.
====
"2.3.2 Flonum Constants
-max.0 : Flonum
-min.0 : Flonum
+min.0 : Flonum
+max.0 : Flonum
The rational flonums with maximum and minimum magnitude.
Example:
> (list -max.0 -min.0 +min.0 +max.0)
'(-1.7976931348623157e+308
-4.9406564584125e-324
4.9406564584125e-324
1.7976931348623157e+308)"
The minimum magnitude rational flonums are 0.0 and -0.0. Perhaps
"minimum non-zero magnitude"?
====
gamma-inc
"The following identities should hold:
• (gamma-inc k x) = 0"
I'm no expert, but this identity smells 8^)
====
pairwise-coprime?
"Returns #t if the integers a b ... are pairwise coprime, meaning that
each adjacent pair of integers is coprime."
Either are coprime pairwise or only adjacent pairs. (I cannot check
what the implementation does, at the moment. I guess the
implementation is correct and the doc is badly worded).
====
"4.5 Multiplicative Functions
The functions in this section are multiplicative. In number theory, a
multiplicative function is a function f such that (f a b) = (* (f a)
(f b)) for all coprime natural numbers a and b."
Probably should be: (f (* a b)) = (* (f a) (f b)) ?
====
moebius-mu
"Returns:
• 1 if n is a product of an even number of primes
• -1 if n is a product of an odd number of primes
• 0 if n has a multiple prime factor"
Doesn't look right.
(moebius-mu 4) ==> ??
should be 1 according to first rule.
should be 0 according to third rule.
OK. checking wikipedia:
"μ(n) is defined for all positive integers n and has its values in
{−1, 0, 1} depending on the factorization of n into prime factors. It
is defined as follows:
μ(n) = 1 if n is a square-free positive integer with an even number
of prime factors.
μ(n) = −1 if n is a square-free positive integer with an odd number
of prime factors.
μ(n) = 0 if n is not square-free."
it's the square-free condition that's missing.
====
prime-omega
"Counting multiplicities the number of prime factors of n is returned."
Doesn't parse (maybe it's a bug in my English parser, though).
"Note: The function prime-omega is multiplicative."
Since it's in a section named Multiplicative Functions, the note
appears redundant.
====
farey-sequence
The examples show always an increasing sequence. If this is guaranted,
it is well worth to put it in the spec.
If it's not guaranted, then please ignore this comment.
====
multinomial
"(multinomial n ks) → Natural
n : Integer
ks : (Listof Integer)
...
> (multinomial 5 3 2)
10"
The example does not agree with the spec.
====
triangle-number?, triangle-number
I think it's customary to call these numbers 'triangular'. Also for
consistency with the other polygonal number function names.
====
"4.11 The group Zn and Primitive Roots
...
A generator the group Un is called a primitive root modolo n. Note
that g is a primitive root if and only if order(g)=phi(n), where phi
is Eulers totient. A group with a generator is called cyclic."
Ehm... why not write simply order(g)=totient(n), since this function
has been calles 'totient' in this library.
(also: "modolo")
====
"(unit-group-orders n) → (Listf Positive-Integer)
n : Integer"
Typo
====
Unless I missed it, there's no evaluator of polynomials? seems a
useful and common operation, and surely there must be one or more
implementations somewhere in the internals of the lib?
====
The prefix for bigfloat operations sometimes is "bigfloat-", sometimes
"bf-", sometimes "bf", with no evident (to me) pattern.
My two neurons protest vehemently.
====
Unless I missed it, there's no discussion of the result precision of
mixed precision bf operations. Maybe it should be mentioned?
====
bflog2, bflog10, bfexp2, bfexp10
Unless I have missed them, there are no fl and generic equivalent of these.
These are commonly used functions. Maybe should be provided?
====
bfagm
Wins the prize for the most fortranesque name.
====
Arrays
why math/arrays and not data/arrays ?
====
diagonal-array
"returns a square array..."
hmmm... is 'square' appropriate for more than 2-dimensional arrays?
====
Great work!
====
Cheers
P.
Jens Axel Søgaard
2012/12/17 Pierpaolo Bernardi <olop...@gmail.com>:
> Some comments after scanning the manual.
> ====
> gamma-inc
>
> "The following identities should hold:
>
> • (gamma-inc k x) = 0"
>
> I'm no expert, but this identity smells 8^)
I think it should be (gamma-inc k x) >= 0 ?
>
> ====
> pairwise-coprime?
>
> "Returns #t if the integers a b ... are pairwise coprime, meaning that
> each adjacent pair of integers is coprime."
>
> Either are coprime pairwise or only adjacent pairs. (I cannot check
> what the implementation does, at the moment. I guess the
> implementation is correct and the doc is badly worded).
The word "adjacent" should be removed from the docs.
The implementation checks each pair.
> ====
> "4.5 Multiplicative Functions
>
> The functions in this section are multiplicative. In number theory, a
> multiplicative function is a function f such that (f a b) = (* (f a)
> (f b)) for all coprime natural numbers a and b."
>
> Probably should be: (f (* a b)) = (* (f a) (f b)) ?
Yes.
> ====
> moebius-mu
>
> "Returns:
> • 1 if n is a product of an even number of primes
> • -1 if n is a product of an odd number of primes
> • 0 if n has a multiple prime factor"
>
> Doesn't look right.
>
> (moebius-mu 4) ==> ??
>
> should be 1 according to first rule.
> should be 0 according to third rule.
>
> OK. checking wikipedia:
>
> "μ(n) is defined for all positive integers n and has its values in
> {−1, 0, 1} depending on the factorization of n into prime factors. It
> is defined as follows:
> μ(n) = 1 if n is a square-free positive integer with an even number
> of prime factors.
> μ(n) = −1 if n is a square-free positive integer with an odd number
> of prime factors.
> μ(n) = 0 if n is not square-free."
>
> it's the square-free condition that's missing.
Would it be enough to add "distinct" ?
> • 1 if n is a product of an even number of distinct primes
> • -1 if n is a product of an odd number of distinct primes
> • 0 if n has a multiple prime factor"
Or is square-free easier to understand?
> ====
> prime-omega
>
> "Counting multiplicities the number of prime factors of n is returned."
>
> Doesn't parse (maybe it's a bug in my English parser, though).
I'll let a native speaker find a better wording.
> "Note: The function prime-omega is multiplicative."
>
> Since it's in a section named Multiplicative Functions, the note
> appears redundant.
True.
>
> ====
> farey-sequence
>
> The examples show always an increasing sequence. If this is guaranted,
> it is well worth to put it in the spec.
> If it's not guaranted, then please ignore this comment.
It is guaranted.
> ====
> multinomial
>
> "(multinomial n ks) → Natural
> n : Integer
> ks : (Listof Integer)
> ...
>> (multinomial 5 3 2)
> 10"
>
> The example does not agree with the spec.
A proper formula would be easier.
5! / ( 3! * 2! ) = 120 / ( 6*2 ) = 10
I read the spec to give this:
> (apply / (factorial 5) (map factorial (list 2 3)))
10
> ====
> triangle-number?, triangle-number
>
> I think it's customary to call these numbers 'triangular'. Also for
> consistency with the other polygonal number function names.
I think both are used, but I am not attached to the name.
> ====
> "4.11 The group Zn and Primitive Roots
> ...
> A generator the group Un is called a primitive root modolo n. Note
> that g is a primitive root if and only if order(g)=phi(n), where phi
> is Eulers totient. A group with a generator is called cyclic."
>
> Ehm... why not write simply order(g)=totient(n), since this function
> has been calles 'totient' in this library.
Agree.
> (also: "modolo")
> ====
> "(unit-group-orders n) → (Listf Positive-Integer)
> n : Integer"
>
> Typo
Thanks again for the careful proof reading.
/Jens Axel
Pierpaolo Bernardi
<jens...@soegaard.net> wrote:
>> ====
>> moebius-mu
> Would it be enough to add "distinct" ?
IMO, yes.
>> ====
>> multinomial
>>
>> "(multinomial n ks) → Natural
>> n : Integer
>> ks : (Listof Integer)
>> ...
>>> (multinomial 5 3 2)
>> 10"
>>
>> The example does not agree with the spec.
>
> A proper formula would be easier.
>
> 5! / ( 3! * 2! ) = 120 / ( 6*2 ) = 10
>
> I read the spec to give this:
>
>> (apply / (factorial 5) (map factorial (list 2 3)))
> 10
Here there is a misunderstanding. I was not complaining about the formulas.
The spec says it is a 2-arg function, an integer and a list of
integers. In the example it takes three integers.
====
> Thanks again for the careful proof reading.
Thanks for the excellent work!
P.
Neil Toronto
> Some comments after scanning the manual.
>
> ====
> flulp-error
>
> "For non-rational arguments such as +nan.0, flulp-error returns 0.0 if
> (eqv? x r); otherwise it returns +inf.0."
>
> However one example is:
>
>> (flulp-error +inf.0 +nan.0)
> +nan.0
>
> Which does not agree with the doc.
Whoops! I'll fix this one.
> "2.3.2 Flonum Constants
>
> -max.0 : Flonum
> -min.0 : Flonum
> +min.0 : Flonum
> +max.0 : Flonum
>
> The rational flonums with maximum and minimum magnitude.
>
> Example:
>
>> (list -max.0 -min.0 +min.0 +max.0)
> '(-1.7976931348623157e+308
> -4.9406564584125e-324
> 4.9406564584125e-324
> 1.7976931348623157e+308)"
>
> The minimum magnitude rational flonums are 0.0 and -0.0. Perhaps
> "minimum non-zero magnitude"?
Yes. Or "The nonzero, rational flonums with maximum and minimum
magnitude." I'll fix it.
> gamma-inc
>
> "The following identities should hold:
>
> • (gamma-inc k x) = 0"
>
> I'm no expert, but this identity smells 8^)
Looks like another case of eyes just sliding over an obvious mistake
because the brain thinks it already knows what's there. It should be
(gamma-inc k 0) = 0. I'll fix it.
> multinomial
>
> "(multinomial n ks) → Natural
> n : Integer
> ks : (Listof Integer)
> ...
>> (multinomial 5 3 2)
> 10"
>
> The example does not agree with the spec.
I didn't update the docs after changing all the multinomial function
argument types. I'll fix this.
> Unless I missed it, there's no evaluator of polynomials? seems a
> useful and common operation, and surely there must be one or more
> implementations somewhere in the internals of the lib?
There are, for Chebyshev polynomials. When I considered making
`math/polynomial', I discovered that I needed at least two kinds:
Chebyshev basis and power basis (e.g. a + bx + cx^2 + ...). Not only
that, but I needed functions for flonum, real, complex, and bigfloat
polynomials, for both kinds.
I'd like to provide a general type for polynomials with any orthogonal
basis, parameterized on type, with polymorphic operations that dispatch
to flonum/real/complex/bigfloat implementations. (And more. For example,
it's sensible to use polynomials as coefficients. Using n-dimensional
monomial coefficients gives you supersparse n-dimensional polynomials.)
That'll require support for generics in Typed Racket, which doesn't exist.
At very least, if `math/polynomial' doesn't have all that nifty
parameterization, I'd like it to be designed to not break user programs
when we upgrade it.
It's a lot of work either way, so I decided to hold back for Math v1.0.
> The prefix for bigfloat operations sometimes is "bigfloat-", sometimes
> "bf-", sometimes "bf", with no evident (to me) pattern.
Those that have "bigfloat" in the name are struct accessors, the
bigfloat predicate, and conversion functions.
Values with the "bf-" prefix don't operate on bigfloats: they're
parameters or constants.
The docs should explain this. Might be enough to group them better. Will
fix.
> My two neurons protest vehemently.
Classic!
> Unless I missed it, there's no discussion of the result precision of
> mixed precision bf operations. Maybe it should be mentioned?
The docs for `bf-precision' cover it:
A parameter that determines the precision of new bigfloat values.
With a few noted exceptions, bigfloat functions conceptually compute in
infinite precision, round the answers to (bf-precision) bits using
(bf-rounding-mode), and return the rounded answers. (In practice, they
employ finite algorithms that have been painstakingly proved to be
equivalent to the aforementioned infinite process.)
**
It's been bugging me that this isn't discussed earlier in the docs. I'll
try to do something about it. Advice welcome.
> bflog2, bflog10, bfexp2, bfexp10
>
> Unless I have missed them, there are no fl and generic equivalent of these.
> These are commonly used functions. Maybe should be provided?
I'm not sure why MPFR provides these. (bfexp2 x) is equivalent to
(bfexpt (bf 2) x), and similarly for `bfexp10'. If they had `bflogb'
(log with base), the former two would be similarly redundant.
I'm planning `fllogb' instead of `fllog2' and `fllog10'. Just haven't
gotten around to it.
> bfagm
>
> Wins the prize for the most fortranesque name.
Woo! Do I get a toaster? Luggage? An all-expenses-paid vacation to
Wisconsin?
> why math/arrays and not data/arrays ?
Good question.
> diagonal-array
>
> "returns a square array..."
>
> hmmm... is 'square' appropriate for more than 2-dimensional arrays?
Nope. Will fix.
> Great work!
Thanks, and thanks for the feedback!
Neil ⊥