Skip to content

Commit

Permalink
feat(stdlib): Reimplement Number.sin, Number.cos, Number.tan (#…
Browse files Browse the repository at this point in the history
  • Loading branch information
spotandjake authored Sep 21, 2024
1 parent e99dcba commit f97c011
Show file tree
Hide file tree
Showing 15 changed files with 1,035 additions and 0 deletions.
108 changes: 108 additions & 0 deletions compiler/test/stdlib/number.test.gr
Original file line number Diff line number Diff line change
Expand Up @@ -949,3 +949,111 @@ assert Number.linearMap(
0.5
) ==
150

// Number.sin - 0 to pi/2
assert Number.sin(0) == 0
assert Number.isClose(Number.sin(Number.pi / 6), 1/2)
assert Number.isClose(Number.sin(Number.pi / 4), Number.sqrt(2) / 2)
assert Number.isClose(Number.sin(Number.pi / 3), Number.sqrt(3) / 2)
assert Number.isClose(Number.sin(Number.pi / 2), 1)
// Number.sin - pi/2 to 2pi
assert Number.isClose(Number.sin(2 * Number.pi / 3), Number.sqrt(3) / 2)
assert Number.isClose(Number.sin(3 * Number.pi / 4), Number.sqrt(2) / 2)
assert Number.isClose(Number.sin(5 * Number.pi / 6), 1/2)
// Note: This has an absolute error of 1e-15 because `Number.pi` is not exact
assert Number.isClose(Number.sin(Number.pi), 0, absoluteTolerance=1e-15)
// Number.sin - 2pi to 3pi/2
assert Number.isClose(Number.sin(7 * Number.pi / 6), -1/2)
assert Number.isClose(Number.sin(5 * Number.pi / 4), Number.sqrt(2) / -2)
assert Number.isClose(Number.sin(4 * Number.pi / 3), Number.sqrt(3) / -2)
assert Number.isClose(Number.sin(3 * Number.pi / 2), -1)
// Number.sin - 3pi/2 to 0
assert Number.isClose(Number.sin(5 * Number.pi / 3), Number.sqrt(3) / -2)
assert Number.isClose(Number.sin(7 * Number.pi / 4), Number.sqrt(2) / -2)
assert Number.isClose(Number.sin(11 * Number.pi / 6), -1/2)
// Note: This has an absolute error of 1e-15 because `Number.pi` is not exact
assert Number.isClose(Number.sin(2 * Number.pi), 0, absoluteTolerance=1e-15)
// Number.sin - special cases
assert Number.sin(1/2) == Number.sin(0.5)
assert Number.sin(1/4) == Number.sin(0.25)
assert Number.isClose(Number.sin(18_446_744_073_709_551_615), 0.0235985099) // 2^64-1
assert Number.isClose(Number.sin(-18_446_744_073_709_551_615), -0.0235985099) // -2^64+1
assert Number.isClose(Number.sin(18_446_744_073_709_551_616), 0.0235985099) // 2^64
assert Number.isClose(Number.sin(-18_446_744_073_709_551_616), -0.0235985099) // -2^64
assert Number.isClose( // Note: We lose a lot of precision here do to ieee754 representation
Number.sin(1.7976931348623157e+308),
0.0049619,
absoluteTolerance=1e7
) // Max F64
assert Number.isClose(
Number.sin(-1.7976931348623157e+308),
0.00496,
absoluteTolerance=1e7
) // Max F64
assert Number.isNaN(Number.sin(Infinity))
assert Number.isNaN(Number.sin(-Infinity))
assert Number.isNaN(Number.sin(NaN))

// Number.cos - 0 to pi/2
assert Number.cos(0) == 1
assert Number.isClose(Number.cos(Number.pi / 6), Number.sqrt(3) / 2)
assert Number.isClose(Number.cos(Number.pi / 4), Number.sqrt(2) / 2)
assert Number.isClose(Number.cos(Number.pi / 3), 1/2)
// Note: This has an absolute error of 1e-15 because `Number.pi` is not exact
assert Number.isClose(Number.cos(Number.pi / 2), 0, absoluteTolerance=1e-15)
// Number.cos - pi/2 to 2pi
assert Number.isClose(Number.cos(2 * Number.pi / 3), -1/2)
assert Number.isClose(Number.cos(3 * Number.pi / 4), Number.sqrt(2) / -2)
assert Number.isClose(Number.cos(5 * Number.pi / 6), Number.sqrt(3) / -2)
assert Number.isClose(Number.cos(Number.pi), -1)
// Number.cos - 2pi to 3pi/2
assert Number.isClose(Number.cos(7 * Number.pi / 6), Number.sqrt(3) / -2)
assert Number.isClose(Number.cos(5 * Number.pi / 4), Number.sqrt(2) / -2)
assert Number.isClose(Number.cos(4 * Number.pi / 3), -1/2)
// Note: This has an absolute error of 1e-15 because `Number.pi` is not exact
assert Number.isClose(Number.cos(3 * Number.pi / 2), 0, absoluteTolerance=1e-15)
// Number.cos - 3pi/2 to 0
assert Number.isClose(Number.cos(5 * Number.pi / 3), 1/2)
assert Number.isClose(Number.cos(7 * Number.pi / 4), Number.sqrt(2) / 2)
assert Number.isClose(Number.cos(11 * Number.pi / 6), Number.sqrt(3) / 2)
assert Number.isClose(Number.cos(2 * Number.pi), 1)
// Number.cos - special cases
assert Number.cos(1/2) == Number.cos(0.5)
assert Number.cos(1/4) == Number.cos(0.25)
assert Number.isClose(Number.cos(18_446_744_073_709_551_615), -0.99972151638) // 2^64-1
assert Number.isClose(Number.cos(-18_446_744_073_709_551_615), -0.99972151638) // -2^64+1
assert Number.isClose(Number.cos(18_446_744_073_709_551_616), -0.99972151638) // 2^64
assert Number.isClose(Number.cos(-18_446_744_073_709_551_616), -0.99972151638) // -2^64
assert Number.isClose(Number.cos(1.7976931348623157e+308), -0.99998768942) // Max F64
assert Number.isClose(Number.cos(-1.7976931348623157e+308), -0.99998768942) // Max F64
assert Number.isNaN(Number.cos(Infinity))
assert Number.isNaN(Number.cos(-Infinity))
assert Number.isNaN(Number.cos(NaN))

// Number.tan - base cases
assert Number.tan(0) == 0
assert Number.isClose(Number.tan(Number.pi / 6), 1 / Number.sqrt(3))
assert Number.isClose(Number.tan(Number.pi / 4), 1)
assert Number.isClose(Number.tan(Number.pi / 3), Number.sqrt(3))
// Note: one might expect this to produce infinity but instead we produce 16331239353195370 because pi can not be represented accurately in iee754, This logic follows c
assert Number.isClose(Number.tan(Number.pi / 2), 16331239353195370)
// Number.tan - special cases
assert Number.tan(1/2) == Number.tan(0.5)
assert Number.tan(1/4) == Number.tan(0.25)
assert Number.isClose(Number.tan(18_446_744_073_709_551_615), -0.02360508353) // 2^64-1
assert Number.isClose(Number.tan(-18_446_744_073_709_551_615), 0.02360508353) // -2^64+1
assert Number.isClose(Number.tan(18_446_744_073_709_551_616), -0.02360508353) // 2^64
assert Number.isClose(Number.tan(-18_446_744_073_709_551_616), 0.02360508353) // -2^64
assert Number.isClose( // Note: We lose a lot of precision here do to ieee754 representation
Number.tan(1.7976931348623157e+308),
-0.00496201587,
absoluteTolerance=1e7
) // Max F64
assert Number.isClose(
Number.tan(-1.7976931348623157e+308),
-0.00496201587,
absoluteTolerance=1e7
) // Max F64
assert Number.isNaN(Number.tan(Infinity))
assert Number.isNaN(Number.tan(-Infinity))
assert Number.isNaN(Number.tan(NaN))
68 changes: 68 additions & 0 deletions stdlib/number.gr
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,8 @@ from "runtime/atoi/parse" include Parse as Atoi
from "runtime/atof/parse" include Parse as Atof
from "runtime/unsafe/tags" include Tags
from "runtime/exception" include Exception
from "runtime/math/trig" include Trig
use Trig.{ sin, cos, tan }

use Atoi.{ type ParseIntError }

Expand Down Expand Up @@ -1072,3 +1074,69 @@ provide let linearMap = (inputRange, outputRange, current) => {
clamp(outputRange, mapped)
}
}

/**
* Computes the sine of a number (in radians).
*
* @param radians: The input in radians
* @returns The computed sine
*
* @example Number.sin(0) == 0
*
* @since v0.7.0
*/
@unsafe
provide let sin = (radians: Number) => {
use WasmF64.{ (==) }
let xval = coerceNumberToWasmF64(radians)
let value = sin(xval)
return if (!isFloat(radians) && value == WasmF64.trunc(value)) {
WasmI32.toGrain(reducedInteger(WasmI64.truncF64S(value))): Number
} else {
WasmI32.toGrain(newFloat64(value)): Number
}
}

/**
* Computes the cosine of a number (in radians).
*
* @param radians: The input in radians
* @returns The computed cosine
*
* @example Number.cos(0) == 1
*
* @since v0.7.0
*/
@unsafe
provide let cos = (radians: Number) => {
use WasmF64.{ (==) }
let xval = coerceNumberToWasmF64(radians)
let value = cos(xval)
return if (!isFloat(radians) && value == WasmF64.trunc(value)) {
WasmI32.toGrain(reducedInteger(WasmI64.truncF64S(value))): Number
} else {
WasmI32.toGrain(newFloat64(value)): Number
}
}

/**
* Computes the tangent of a number (in radians).
*
* @param radians: The input in radians
* @returns The computed tangent
*
* @example Number.tan(0) == 0
*
* @since v0.7.0
*/
@unsafe
provide let tan = (radians: Number) => {
use WasmF64.{ (==) }
let xval = coerceNumberToWasmF64(radians)
let value = tan(xval)
return if (!isFloat(radians) && value == WasmF64.trunc(value)) {
WasmI32.toGrain(reducedInteger(WasmI64.truncF64S(value))): Number
} else {
WasmI32.toGrain(newFloat64(value)): Number
}
}
93 changes: 93 additions & 0 deletions stdlib/number.md
Original file line number Diff line number Diff line change
Expand Up @@ -1566,3 +1566,96 @@ Throws:
* When `outputRange` is not finite
* When `outputRange` includes NaN

### Number.**sin**

<details disabled>
<summary tabindex="-1">Added in <code>next</code></summary>
No other changes yet.
</details>

```grain
sin : (radians: Number) => Number
```

Computes the sine of a number (in radians).

Parameters:

|param|type|description|
|-----|----|-----------|
|`radians`|`Number`|The input in radians|

Returns:

|type|description|
|----|-----------|
|`Number`|The computed sine|

Examples:

```grain
Number.sin(0) == 0
```

### Number.**cos**

<details disabled>
<summary tabindex="-1">Added in <code>next</code></summary>
No other changes yet.
</details>

```grain
cos : (radians: Number) => Number
```

Computes the cosine of a number (in radians).

Parameters:

|param|type|description|
|-----|----|-----------|
|`radians`|`Number`|The input in radians|

Returns:

|type|description|
|----|-----------|
|`Number`|The computed cosine|

Examples:

```grain
Number.cos(0) == 1
```

### Number.**tan**

<details disabled>
<summary tabindex="-1">Added in <code>next</code></summary>
No other changes yet.
</details>

```grain
tan : (radians: Number) => Number
```

Computes the tangent of a number (in radians).

Parameters:

|param|type|description|
|-----|----|-----------|
|`radians`|`Number`|The input in radians|

Returns:

|type|description|
|----|-----------|
|`Number`|The computed tangent|

Examples:

```grain
Number.tan(0) == 0
```

70 changes: 70 additions & 0 deletions stdlib/runtime/math/kernel/cos.gr
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
module Cosine

from "runtime/unsafe/wasmf64" include WasmF64

/*
* Source: https://git.musl-libc.org/cgit/musl/tree/src/math/__cos.c
*
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
*
* Algorithm
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
* 3. cos(x) is approximated by a polynomial of degree 14 on
* [0,pi/4]
* 4 14
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
* where the remez error is
*
* | 2 4 6 8 10 12 14 | -58
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
* | |
*
* 4 6 8 10 12 14
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
* cos(x) ~ 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* ~ cos(x) - x*y,
* a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
* For better accuracy, rearrange to
* cos(x+y) ~ w + (tmp + (r-x*y))
* where w = 1 - x*x/2 and tmp is a tiny correction term
* (1 - x*x/2 == w + tmp exactly in infinite precision).
* The exactness of w + tmp in infinite precision depends on w
* and tmp having the same precision as x. If they have extra
* precision due to compiler bugs, then the extra precision is
* only good provided it is retained in all terms of the final
* expression for cos(). Retention happens in all cases tested
* under FreeBSD, so don't pessimize things by forcibly clipping
* any extra precision in w.
*/
@unsafe
provide let cos = (x: WasmF64, y: WasmF64) => {
use WasmF64.{ (+), (-), (*) }
let c1 = 4.16666666666666019037e-02W /* 0x3FA55555, 0x5555554C */
let c2 = -1.38888888888741095749e-03W /* 0xBF56C16C, 0x16C15177 */
let c3 = 2.48015872894767294178e-05W /* 0x3EFA01A0, 0x19CB1590 */
let c4 = -2.75573143513906633035e-07W /* 0xBE927E4F, 0x809C52AD */
let c5 = 2.08757232129817482790e-09W /* 0x3E21EE9E, 0xBDB4B1C4 */
let c6 = -1.13596475577881948265e-11W /* 0xBDA8FAE9, 0xBE8838D4 */

let z = x * x
let w = z * z
let r = z * (c1 + z * (c2 + z * c3)) + w * w * (c4 + z * (c5 + z * c6))
let hz = 0.5W * z
let w = 1.0W - hz
w + (1.0W - w - hz + (z * r - x * y))
}
14 changes: 14 additions & 0 deletions stdlib/runtime/math/kernel/cos.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
---
title: Cosine
---

## Values

Functions and constants included in the Cosine module.

### Cosine.**cos**

```grain
cos : (x: WasmF64, y: WasmF64) => WasmF64
```

Loading

0 comments on commit f97c011

Please sign in to comment.