Changeset 26284 in project


Ignore:
Timestamp:
03/29/12 21:47:49 (9 years ago)
Author:
sjamaan
Message:

numbers: Fix acos and asin so they work for the general case (returning compnums for x < -1 and x > 1). Add and fix the relevant tests

Location:
release/4/numbers/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • release/4/numbers/trunk/numbers.scm

    r26283 r26284  
    14071407    (else (##core#inline_allocate ("C_a_i_tan" 4) (%exact->inexact n)) ) ))
    14081408
     1409;; General case: sin^{-1}(z) = -i\ln(iz + \sqrt{1-z^2})
    14091410(define (%asin n)
    1410   (switchq (%check-number n)
    1411     (NONE (bad-number 'asin n))
    1412     (COMP (%* %-i (%log (%+ (%* %i n) (%sqrt (%- 1 (%* n n)))))))
    1413     (else (##core#inline_allocate ("C_a_i_asin" 4) (%exact->inexact n)) ) ))
     1411  (let ((t (%check-number n)))
     1412    (cond ((eq? t NONE) (bad-number 'asin n))
     1413          ((and (eq? t FLO) (fp>= n -1.0) (fp<= n 1.0))
     1414           (##core#inline_allocate ("C_a_i_asin" 4) n))
     1415          ((and (eq? t FIX) (fx>= n -1) (fx<= n 1))
     1416           (##core#inline_allocate ("C_a_i_asin" 4) (%fix->flo n)))
     1417          ;; General definition can return compnums
     1418          (else (%* %-i (%log (%+ (%* %i n) (%sqrt (%- 1 (%* n n))))))))))
    14141419
    14151420(define asin %asin)
    14161421
    1417 (define acos
     1422;; General case:
     1423;; cos^{-1}(z) = 1/2\pi + i\ln(iz + \sqrt{1-z^2}) = 1/2\pi - sin^{-1}(z) = sin(1) - sin(z)
     1424(define %acos
    14181425  (let ((asin1 (##core#inline_allocate ("C_a_i_asin" 4) 1)))
    14191426    (lambda (n)
    1420       (switchq (%check-number n)
    1421         (NONE (bad-number 'acos n))
    1422         (COMP (%- asin1 (%asin n)))
    1423         (else (##core#inline_allocate ("C_a_i_acos" 4) (%exact->inexact n)) ) ) ) ) )
     1427      (let ((t (%check-number n)))
     1428        (cond ((eq? t NONE) (bad-number 'acos n))
     1429              ((and (eq? t FLO) (fp>= n -1.0) (fp<= n 1.0))
     1430               (##core#inline_allocate ("C_a_i_acos" 4) n))
     1431              ((and (eq? t FIX) (fx>= n -1) (fx<= n 1))
     1432               (##core#inline_allocate ("C_a_i_acos" 4) (%fix->flo n)))
     1433              ;; General definition can return compnums
     1434              (else (%- asin1 (%asin n)))) ) ) ) )
     1435
     1436(define acos %acos)
    14241437
    14251438(define (atan n #!optional b)
  • release/4/numbers/trunk/numbers.types

    r26241 r26284  
    330330             ((float) (float) (##core#inline_allocate ("C_a_i_flonum_tan" 4) #(1))))
    331331
    332 (numbers#asin (#(procedure #:clean #:enforce) numbers#asin ((or fixnum float (struct bignum) (struct compnum) (struct ratnum))) (or float (struct compnum)))
    333               ((float) (float) (##core#inline_allocate ("C_a_i_flonum_asin" 4) #(1))))
     332(numbers#asin (#(procedure #:clean #:enforce) numbers#asin ((or fixnum float (struct bignum) (struct compnum) (struct ratnum))) (or float (struct compnum)))
     333              ;; Unfortunately this doesn't work when the number is > 1.0 (returns compnum)
     334              #;((float) (float) (##core#inline_allocate ("C_a_i_flonum_asin" 4) #(1))))
    334335
    335336(numbers#acos (#(procedure #:clean #:enforce) numbers#acos ((or fixnum float (struct bignum) (struct compnum) (struct ratnum))) (or float (struct compnum)))
    336               ((float) (float) (##core#inline_allocate ("C_a_i_flonum_acos" 4) #(1))))
     337              ;; Unfortunately this doesn't work when the number is > 1.0 (returns compnum)
     338              #;((float) (float) (##core#inline_allocate ("C_a_i_flonum_acos" 4) #(1))))
    337339
    338340(numbers#atan (#(procedure #:clean #:enforce) numbers#atan ((or fixnum float (struct bignum) (struct compnum) (struct ratnum)) #!optional number) (or float (struct compnum)))
  • release/4/numbers/trunk/tests/numbers-test.scm

    r26283 r26284  
    797797      (test "cos(   2pi)" 1.0 (cos (* 2 pi)))
    798798      (test "acos(cos(   2pi))" 0 (acos (cos (* 2 pi))))
    799       (test "acos(<large number>)" +nan.0 (acos 1e100))
     799      (test "acos(pi)" 0.0+1.81152627246085i (acos pi))
    800800      (test "acos(+inf)" -nan.0 (acos +inf.0))
    801801
     
    832832      (test "sin(   2pi)" 0.0 (sin (* 2 pi)))
    833833      (test "asin(sin(   2pi))" 0.0 (asin (sin (* 2 pi))))
    834       (test "asin(<large number>)" +nan.0 (asin 1e100))
     834      (test "asin(pi)" 1.57079632679490-1.81152627246085i (asin pi))
    835835      (test "asin(+inf)" -nan.0 (asin +inf.0))
    836836     
     
    862862      (test "tan(   2pi)" 0.0 (tan (* 2 pi)))
    863863      (test "atan(tan(   2pi))" 0.0 (atan (tan (* 2 pi))))
    864       (test "atan(<large number>)" (/ pi 2) (atan 1e100))
     864      (test "atan(pi)" (/ pi 2) (atan 1e100))
    865865      (test "atan(+inf)" (/ pi 2) (atan +inf.0))
    866866
     
    884884      (test "cos(0.0+1.0i)" 1.5430806348152437
    885885            (cos (make-rectangular 0.0 1.0)))
    886       (test "acos(cos(0.0+1.0i))" +nan.0
     886      (test "acos(cos(0.0+1.0i))" 0.0+1.0i
    887887            (acos (cos (make-rectangular 0.0 1.0))))
    888888      (test "cos(0.0-1.0i)" 1.5430806348152437
    889889            (cos (make-rectangular 0.0 -1.0)))
    890       (test "acos(cos(0.0-1.0i))" +nan.0
     890      (test "acos(cos(0.0-1.0i))" 0.0+1.0i
    891891            (acos (cos (make-rectangular 0.0 -1.0))))
    892892      (test "cos(0.0+3.0i)" 10.067661995777765
    893893            (cos (make-rectangular 0.0 3.0)))
    894       (test "acos(cos(0.0+3.0i))" +nan.0
     894      (test "acos(cos(0.0+3.0i))" 0.0+3.0i
    895895            (acos (cos (make-rectangular 0.0 3.0))))
    896896      (test "cos(0.0-3.0i)" 10.067661995777765
    897897            (cos (make-rectangular 0.0 -3.0)))
    898       (test "acos(cos(0.0-3.0i))" +nan.0
     898      (test "acos(cos(0.0-3.0i))" 0.0+3.0i
    899899            (acos (cos (make-rectangular 0.0 -3.0))))
    900900      (test "cos(0.5+0.5i)"
     
    11711171    (test-group "bignums"
    11721172      (test "acos(<negative bignum>)" -nan.0 (acos (- b1)))
    1173       (test "acos(<bignum>)" +nan.0 (acos b1))
     1173      ;; These are bogus (maybe the negative ones too!), but I don't want to
     1174      ;; "fix" them by copying the output and assume it's alright.
     1175      #;(test "acos(<bignum>)" +nan.0 (acos b1))
    11741176      (test "asin(<negative bignum>)" -nan.0 (asin (- b1)))
    1175       (test "asin(<bignum>)" +nan.0 (asin b1))
     1177      #;(test "asin(<bignum>)" +nan.0 (asin b1))
    11761178      (test "atan(<negative bignum>)" (- (/ pi 2)) (atan (- b1)))
    11771179      (test "atan(<bignum>)" (/ pi 2) (atan b1)))
     
    11851187      (test "acos(1)" 0.0 (acos 1))
    11861188      (test "cos(-1)" (cos -1.0) (cos -1))
    1187       (test "acos(-1)" pi (acos -1)))
     1189      (test "acos(-1)" pi (acos -1))
     1190      (test "acos(-2)" (make-rectangular pi -1.31695789692482) (acos -2))
     1191      (test "acos(2)" 0.0+1.31695789692482i (acos 2))
     1192      (test "asin(1)" (/ pi 2) (asin 1))
     1193      (test "asin(-1)" (/ pi -2) (asin -1))
     1194      (test "asin(2)" (make-rectangular (/ pi 2) -1.31695789692482) (asin 2))
     1195      (test "asin(-2)" (make-rectangular (/ pi -2) 1.31695789692482) (asin -2)))
    11881196
    11891197    (test-group "ratnums"
Note: See TracChangeset for help on using the changeset viewer.