Well, since I voted for the basic library with sin, cos, pi, gcd, etc, I felt obligated to implement a basic math library with sin, cos, pi, gcd, closely related fucntions, and everything I needed or thought I'd need to implement those.
Implementations guaranteed not to be optimally efficient:
Edit: Fixed a mistake in bringing out-of-period numbers into period in sin and cos
(= pi 3.141592653589793238)
(= 2*pi (* 2 pi))
(def signum (n)
(if (> n 0)
1
(< n 0)
-1
0))
(def floor (x)
(if (< x 0)
(- (trunc x) 1)
(trunc x)))
(def ceil (x)
(if (< x 0)
(trunc x)
(+ (trunc x) 1)))
(defmemo fac (n)
(if (is n 0)
1
(* n (fac (- n 1)))))
;Existing mod only accepts integers. Left intact as "modulo"
;This mod copies the behavior of Ruby's mod
(def mod (dividend divisor)
"Returns the remainder from dividing dividend by divisor.
Adds divisor once more if they are of different signs so that the result is always of the same sign as divisor."
(if (is divisor 0)
(error "modulo undefined for 0")
(isnt (signum dividend) (signum divisor))
(+ divisor (- dividend (* divisor (trunc (/ dividend divisor)))))
(- dividend (* divisor (trunc (/ dividend divisor))))))
(def sin (x)
"Returns the sine of x in radians."
(let x (let red (mod x 2*pi)
(if (> (abs red) pi)
(- red (* (signum red) 2*pi))
red))
;Taylor polynomial; 0.0 is to cast to float
(- (+ 0.0 x (/ (expt x 5) (fac 5)) (/ (expt x 9) (fac 9)))
(/ (expt x 3) (fac 3)) (/ (expt x 7) (fac 7)))))
(def cos (x)
"Returns the cosine of x in radians."
(let x (let red (mod x 2*pi)
(if (> (abs red) pi)
(- red (* (signum red) 2*pi))
red))
;Taylor polynomial
(- (+ 1.0 (/ (expt x 4) (fac 4)) (/ (expt x 8) (fac 8)))
(/ (expt x 2) (fac 2)) (/ (expt x 6) (fac 6)))))
(def tan (x)
"Returns the tangent of x in radians."
;Lazy definition
(/ (sin x) (cos x)))
(def int? (x)
"Returns whether x is an integer"
(is (mod x 1.0) 0.0))
(defmemo prime (n)
"Returns the nth prime. 2 is the 0th prime."
(if (< n 0)
nil
(is n 0)
2
(let prev-primes (map prime (range 0 (- n 1)))
((afn (i)
(if (no (breakable:each p prev-primes ;Each always returns nil, so a break returns t
(if (int? (/ i p))
(break t))))
i
(self (+ i 1))))
(+ 1 (last prev-primes))))))
(def prime-factorization (n)
"Returns a list each prime up to the greatest prime in n paired with the power of that prime in n.
E.g.: (prime-factorization 20) returns ((2 2) (3 0) (5 1)).
Use (reduce * (map [apply expt _] ...)) to change a prime factorization into the number."
(rev:accum keep
(let p-ord 0
(while (> n 1)
(with (p (prime p-ord)
pow 0)
(until (isnt (mod n p) 0)
(++ pow)
(zap [/ _ p] n))
(keep (list p pow)))
(++ p-ord)))))
(def gcd (x y)
"Returns the greatest common divisor of x and y."
(reduce * (map [apply expt _]
(map (fn (a b)
(list (car a) (min (cadr a) (cadr b))))
(prime-factorization x) (prime-factorization y)))))
(def lcm (x y)
"Returns the least common multiple of x and y."
(/ (* x y) (gcd x y)))
Nice. I only did the gcd function, so you've got an amazing head start on me. However, the version of gcd that I did was Euclid's algorithm. Probably faster than prime factorization; and I included a case for only one number, and a list of greater than two.
Since I couldn't commit to github for some reason, here's the code:
;;Int (predicate)
(def int (n)
"Returns true if n is an integer; i.e. has no decimal part."
(is (trunc n) n)) ; return true if the number is the same as itself without the decimal part
;;Greatest Common Denominator
(def gcd l
"returns the greatest common denominator, or divisor, of a list of numbers. Numbers should be integers,
or the result will default to 0."
(with (a (car l) c (cdr l))
(if (len> c 1)
(gcd a (gcd (car c) (cdr c))) ; handle lists of more than two numbers
(no a) 0
(no c) (abs a)
(let b (car (flat c))
(if (or (~int a) (~int b)) 0 ; skip non-integers
(is a b) a ; return common divisor
(if (> a b)
(gcd (- a b) b)
(gcd (- b a) a)))))))
Update:
Figured out what I was doing wrong with github. Math library now started in lib/math.arc. Currently it only contains my gcd function and the 'int test, but feel free to add any other math functions.
In the future, we may need to make it a math folder, with individual files for trig functions, hyperbolic functions, calculus, matrix and vector function, quaternions, etc. But for now one file is more than enough to contain my lowly gcd :)
My original int? was identical to yours, but that's why I changed it. Although, now that I think about it, since prime-factorization discards the results of operations that cast n to a rational, that's unnecessary. Still, there will need to be separate functions for testing if a num's type is int (which yours does), and testing if a num lacks decimal places (which mine does).
(And yes, Euclid's is definitely faster than prime factorization -- I originally looked at better algorithms for prime finding and gcd but I decided to settle for the "pure" approach -- sieve instead of some Fermat's or another fast primality test, prime factorization based instead of Euclid's.)
Ah, I didn't notice that. Apparently I should test things more carefully next time :)
I don't know if you have github access, but if you do, you can feel free to add your own code to the math.arc lib. I probably won't have a chance to add it myself in the next week or so. If I can, though, I will definitely try to do so.
In doing so, I noticed one bug in Darmani's original implementation of 'floor and 'ceil. 'floor would return incorrect results on negative integers (e.g. (floor -1) => -2), and 'ceil on positive integers (e.g. (ceil 1) => 2). This has been corrected on Anarki.
I also used mzscheme's 'sin, 'cos, and 'tan instead of Darmani's, not because of speed issues, but because of decreased precision in those functions. In order to get maximum precision it would be necessary to calculate the Taylor series an extra couple of terms, which I didn't feel like doing at the time.
I didn't commit 'signum, 'mod, 'prime, or 'prime-factorization, because I wasn't sure if they were needed except for computing 'sin, 'cos, and 'gcd... but feel free to commit them if you want.
1) Isn't pushing those math functions straight from scheme sort of cheating? I mean, maybe I'm just wrong, but wouldn't the solution be more long term if we avoided scheme and implemented the math in arc?
2) Shouldn't fac be tail recursive? Or is it, and I just can't tell? Or are you just expecting that no one will try and compute that large of a factorial
3) If some one did compute that large of a factorial, is there some way for arc to handle arbitrarily sized integers?
1) No, you should implement in the math in the underlying machine instructions, which are guaranteed to be as precise and as fast as the manufacturer can make it. The underlying machine instructions are fortunately possible to access in standard C libraries, and the standard C library functions are wrapped by mzscheme, which we then import in arc.
2) It should be, and it isn't.
(defmemo fac (n)
((afn (n a)
(if (> n 1)
(self (- n 1) (* a n))
a))
n 1))
3) Yes, arc-on-mzscheme handles this automagically. arc2c does not (I think it'll overflow)
Implementing numerically stable and accurate transcendental functions is rather difficult. If you're going down that road, please don't just use Taylor series, but look up good algorithms that others have implemented. One source is http://developer.intel.com/technology/itj/q41999/pdf/transen...
That said, I don't see much value in re-implementing math libraries in Arc, given that Arc is almost certainly going to be running on a platform that already has good native math libraries.
I figured that being close to machine instructions was a good thing, but I thought that we should do that via some other method, not necessarily scheme, which may or may not remain the base of arc in the future.
That being said, if you think that pulling from scheme is a good idea, why don't we just pull all of the other math functions from there as well?
Actually I think it might be better if we had a spec which says "A Good Arc Implementation (TM) includes the following functions when you (require "lib/math.arc"): ...." Then the programmer doesn't even have to care about "scheme functions" or "java functions" or "c functions" or "machine language functions" or "SKI functions" - the implementation imports it by whatever means it wants.
Maybe also spec that the implementation can reserve the plain '$ for implementation-specific stuff.
