From 896e93057b109e59b1953e0f7d725cbebb717ffc Mon Sep 17 00:00:00 2001 From: "Walton, Izaak (VEN)" Date: Tue, 10 Oct 2023 14:13:51 -0700 Subject: [PATCH] Move pi to trig and e to exp, fix atan2 --- library/big-float/impl-default.lisp | 35 +++++++++------- library/big-float/impl-sbcl.lisp | 55 ++++++++++++------------- library/math/dual.lisp | 8 +++- library/math/elementary.lisp | 62 ++++++++++++++--------------- 4 files changed, 86 insertions(+), 74 deletions(-) diff --git a/library/big-float/impl-default.lisp b/library/big-float/impl-default.lisp index 35d3207c2..a7cffd5c9 100644 --- a/library/big-float/impl-default.lisp +++ b/library/big-float/impl-default.lisp @@ -121,7 +121,9 @@ (declare ilog2-abs (Integer -> Integer)) (define (ilog2-abs x) - (- (bit-length (abs x)) 1)) + (- (bit-length (abs x)) 1))) + +(coalton-toplevel (define-type Big-Float (BFValue Dyadic) @@ -540,10 +542,13 @@ returns the nth SeriesSplit, return the series evaluated to the Nth element." (define (make-recip-results n start f) "Reciprocal of `make-result`." (recip-result - (eval-series start (fn (m) (f (toInteger m))) n))) + (eval-series start (fn (m) (f (toInteger m))) n)))) + +(coalton-toplevel ;; XXX: We can pre-compute some values here - (define (bf-pi _) + (declare bf-pi (Unit -> Big-Float)) + (define (bf-pi) "Return the value of pi to the currently set precision." ;; Chudnovsky algorithm (progn @@ -675,8 +680,11 @@ returns the nth SeriesSplit, return the series evaluated to the Nth element." (define (asin x) (atan (* x (reciprocal-sqrt (- 1 (* x x)))))) (define (acos x) - (- (scale (bf-pi) -1) (asin x)))) + (- (scale (bf-pi) -1) (asin x))) + (define pi (BFconst bf-pi)))) +(coalton-toplevel + (define-instance (Exponentiable Big-Float) (define (exp x) (let prec = (get-precision)) @@ -760,9 +768,15 @@ returns the nth SeriesSplit, return the series evaluated to the Nth element." ((Tuple x (BFConst f)) (pow x (f))) (_ (exp (* y (ln x)))))) (define (log b x) - (/ (ln x) (ln b)))) + (/ (ln x) (ln b))) + (define ee (BFConst bf-ee))) + + (declare bf-ee (Unit -> Big-Float)) + (define (bf-ee) + "Return the value of ee = exp(1) to the currently set precision." + (exp 1)) - (define-instance (Radical Big-Float) + (define-instance (Radical Big-Float) (define (sqrt x) (reciprocal (reciprocal-sqrt x))) (define (nth-root n x) @@ -783,14 +797,7 @@ returns the nth SeriesSplit, return the series evaluated to the Nth element." (define (polar z) (Tuple (magnitude z) (phase z)))) - (declare bf-ee (Unit -> Big-Float)) - (define (bf-ee) - "Return the value of ee = exp(1) to the currently set precision." - (exp 1)) - - (define-instance (Elementary Big-Float) - (define pi (BFConst bf-pi)) - (define ee (BFConst bf-ee))) + (define-instance (Elementary Big-Float)) (define (big-float->string x) "Returns a Big-Float in scientific notation." diff --git a/library/big-float/impl-sbcl.lisp b/library/big-float/impl-sbcl.lisp index c3794160d..eaf64709e 100644 --- a/library/big-float/impl-sbcl.lisp +++ b/library/big-float/impl-sbcl.lisp @@ -255,6 +255,12 @@ (coalton-library/math/complex::%define-standard-complex-instances Big-Float) (coalton-toplevel + + (define (bf-pi _) + "Return the value of pi to the currently set precision." + (lisp Big-Float () + (cl:values (sb-mpfr:const-pi)))) + ;; Trig (define-instance (Trigonometric Big-Float) (define (sin x) @@ -274,8 +280,18 @@ (cl:values (sb-mpfr:acos x)))) (define (atan x) (lisp Big-Float (x) - (cl:values (sb-mpfr:atan x))))) - + (cl:values (sb-mpfr:atan x)))) + ;; BUG: Pi is calculated just once, so if we change precision, + ;; it will *NOT* get updated. + (define pi + "Return the value of pi to the currently set precision." + (bf-pi))) + + (define (bf-ee _) + "Return the value of ee = exp(1) to the currently set precision." + (lisp Big-Float () + (cl:values (sb-mpfr:exp (sb-mpfr:coerce 1 'sb-mpfr:mpfr-float))))) + ;; Exp/Log (define-instance (Exponentiable Big-Float) (define (exp x) @@ -299,7 +315,12 @@ (True (/ (ln x) (ln n))))) (define (ln x) (lisp Big-Float (x) - (cl:values (sb-mpfr:log x))))) + (cl:values (sb-mpfr:log x)))) + ;; BUG: EE is calculated just once, so if we change precision, + ;; it will *NOT* get updated. + (define ee + "Return the value of ee = exp(1) to the currently set precision." + (bf-ee))) (define-instance (Radical Big-Float) (define (sqrt x) @@ -310,34 +331,14 @@ (define-instance (Polar Big-Float) (define (phase z) - (let x = (real-part z)) - (let y = (imag-part z)) - (match (Tuple (<=> x 0) (<=> y 0)) - ((Tuple (GT) _) (atan (/ y x))) - ((Tuple (LT) (LT)) (- (atan (/ y x)) (bf-pi))) - ((Tuple (LT) _) (+ (atan (/ y x)) (bf-pi))) - ((Tuple (EQ) (GT)) (/ (bf-pi) 2)) - ((Tuple (EQ) (LT)) (/ (bf-pi) -2)) - ((Tuple (EQ) (EQ)) 0))) + (atan2 (imag-part z) (real-part z))) (define (polar z) (Tuple (magnitude z) (phase z)))) - ;; Elementary - (define (bf-pi _) - "Return the value of pi to the currently set precision." - (lisp Big-Float () - (cl:values (sb-mpfr:const-pi)))) - (define (bf-ee _) - "Return the value of ee = exp(1) to the currently set precision." - (lisp Big-Float () - (cl:values (sb-mpfr:exp (sb-mpfr:coerce 1 'sb-mpfr:mpfr-float))))) - ;; BUG: These are calculated just once, so if we change precision, - ;; these will *NOT* get updated. - (define-instance (Elementary Big-Float) - (define pi (bf-pi)) - (define ee (bf-ee))) -) ; COALTON-TOPLEVEL + (define-instance (Elementary Big-Float))) + +;COALTON-TOPLEVEL #+sb-package-locks (sb-ext:lock-package "COALTON-LIBRARY/BIG-FLOAT") diff --git a/library/math/dual.lisp b/library/math/dual.lisp index ddb082a7e..d1902d9bb 100644 --- a/library/math/dual.lisp +++ b/library/math/dual.lisp @@ -150,7 +150,9 @@ Note: `Eq`, and `Ord` and `Hash` only make use of the primal component." (define (atan (Dual p1 d1)) (Dual (atan p1) - (/ d1 (+ 1 (sq p1)))))) + (/ d1 (+ 1 (sq p1))))) + + (define pi (Dual pi 0))) (define-instance ((Num :t) (Exponentiable :t) (Reciprocable :t) => (Exponentiable (Dual :t))) (define (exp (Dual p1 d1)) @@ -165,7 +167,9 @@ Note: `Eq`, and `Ord` and `Hash` only make use of the primal component." (exp (* dual2 (ln dual1)))) (define (log dual1 dual2) - (/ (ln dual2) (ln dual1)))) + (/ (ln dual2) (ln dual1))) + + (define ee (Dual ee 0))) (define-instance ((Num :t) (Radical :t) (Reciprocable :t) (Exponentiable :t) => (Radical (Dual :t))) (define (nth-root n (Dual p1 d1)) diff --git a/library/math/elementary.lisp b/library/math/elementary.lisp index a98194ba5..0cadf1e8b 100644 --- a/library/math/elementary.lisp +++ b/library/math/elementary.lisp @@ -48,13 +48,26 @@ (tan (:a -> :a)) (asin (:a -> :a)) (acos (:a -> :a)) - (atan (:a -> :a))) + (atan (:a -> :a)) + (pi (:a))) (declare sincos (Trigonometric :a => :a -> (Tuple :a :a))) (define (sincos x) "Computes the sine and cosine of X." (Tuple (sin x) (cos x))) + (declare atan2 ((Ord :f) (Trigonometric :f) (Reciprocable :f) => :f -> :f -> :f)) + (define (atan2 y x) + "Computes the two-argument arctangent of y and x, which is roughly the same +as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." + (match (Tuple (<=> x 0) (<=> y 0)) + ((Tuple (GT) _) (atan (/ y x))) + ((Tuple (LT) (LT)) (- (atan (/ y x)) pi)) + ((Tuple (LT) _) (+ (atan (/ y x)) pi)) + ((Tuple (EQ) (GT)) (/ pi 2)) + ((Tuple (EQ) (LT)) (/ pi -2)) + ((Tuple (EQ) (EQ)) 0))) + (define-class (Exponentiable :a) "Exponential maps obeying: @@ -66,7 +79,8 @@ (exp (:a -> :a)) (pow (:a -> :a -> :a)) (ln (:a -> :a)) - (log (:a -> :a -> :a))) + (log (:a -> :a -> :a)) + (ee (:a))) (define-class (Radical :a) "Obeys: @@ -106,9 +120,7 @@ ((Reciprocable :a) (Polar :a) (Trigonometric :a) (Exponentiable :a) (Radical :a) => Elementary :a) - "Numbers that can be can be passed to elementary functions." - (ee :a) - (pi :a)) + "`Elementary` is a marker class, providing `Reciprocable`, `Polar`, `Trigonometric`, `Exponentiable`, and `Radical`.") ;; See http://clhs.lisp.se/Body/f_sinh_.htm @@ -136,18 +148,6 @@ (define (atanh x) (/ (- (ln (+ 1 x)) (ln (- 1 x))) (fromInt 2))) - (declare atan2 ((Ord :f) (Elementary :f) => :f -> :f -> :f)) - (define (atan2 y x) - "Computes the two-argument arctangent of y and x, which is roughly the same -as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." - (match (Tuple (<=> x 0) (<=> y 0)) - ((Tuple (GT) _) (atan (/ y x))) - ((Tuple (LT) (LT)) (- (atan (/ y x)) pi)) - ((Tuple (LT) _) (+ (atan (/ y x)) pi)) - ((Tuple (EQ) (GT)) (/ pi 2)) - ((Tuple (EQ) (LT)) (/ pi -2)) - ((Tuple (EQ) (EQ)) 0))) - (define (canonical-nth-root n x) "By definition implementation of `nth-root` for reals" (pow x (reciprocal (fromInt n)))) @@ -248,7 +248,10 @@ as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." (lisp ,coalton-type (x) (#+(not ccl) cl:progn #+ccl ff:with-float-traps-masked #+ccl cl:t - (cl:atan x))))))) + (cl:atan x)))))) + (define pi + (lisp ,coalton-type () + (cl:coerce cl:pi ',underlying-type)))) (define-instance (Polar ,coalton-type) (define (phase x) @@ -323,7 +326,10 @@ as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." (lisp ,coalton-type (x) (cl:log x))) ((< x 0) nan) - (True negative-infinity)))) + (True negative-infinity))) + (define ee + (lisp ,coalton-type () + (cl:exp (cl:coerce 1 ',underlying-type))))) (define-instance (Radical ,coalton-type) (define (sqrt x) @@ -336,13 +342,7 @@ as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." (define (nth-root n x) (canonical-nth-root n x))) - (define-instance (Elementary ,coalton-type) - (define ee - (lisp ,coalton-type () - (cl:exp (cl:coerce 1 ',underlying-type)))) - (define pi - (lisp ,coalton-type () - (cl:coerce cl:pi ',underlying-type)))))) + (define-instance (Elementary ,coalton-type)))) (%define-real-float-elementary Single-Float cl:single-float) (%define-real-float-elementary Double-Float cl:double-float) @@ -377,7 +377,8 @@ as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." (define (pow x y) (exp (* y (ln x)))) (define (log b x) - (/ (ln x) (ln b)))) + (/ (ln x) (ln b))) + (define ee (Complex ee 0))) (define-instance ((Elementary :a) => Radical (Complex :a)) (define (sqrt z) @@ -431,7 +432,8 @@ as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." ;; atan = (- i/2 (ln (i - z)/(i+z)) (* (complex 0 (/ -1 2)) (ln (/ (- ii z) - (+ ii z)))))) + (+ ii z))))) + (define pi (Complex pi 0))) ;; This doesn't have much mathematical meaning (define-instance ((Elementary :a) => Polar (Complex :a)) @@ -448,9 +450,7 @@ as (atan (/ y x)) when defined and accounting for the quadrant of the (x,y)." (* 2 (atan (/ y (+ r x)))))) (Tuple r p))) - (define-instance ((Elementary :a) => Elementary (Complex :a)) - (define ee (complex ee 0)) - (define pi (complex pi 0)))) + (define-instance ((Elementary :a) => Elementary (Complex :a)))) #+sb-package-locks (sb-ext:lock-package "COALTON-LIBRARY/MATH/ELEMENTARY")