/*---------------------------------------------------------------------------
  Copyright 2017-2022, Microsoft Research, Daan Leijen.

  This is free software; you can redistribute it and/or modify it under the
  terms of the Apache License, Version 2.0. A copy of the License can be
  found in the LICENSE file at the root of this distribution.
---------------------------------------------------------------------------*/

/* 128-bit double-double floating point numbers.

The `:ddouble` type implements [double double][ddwiki] 128-bit floating point numbers
as a pair of IEEE `:float64` values. This extends the precision to 31 decimal digits
(versus 15 for `:float64`), but keeps the same range as
a `:float64` with a maximum  value of about 1.8·10^308^ (`dd-max`). Because
`:float64`s usually have hardware support, a `:ddouble` is usually much faster
than arbitrary precision floating point numbers.

## Precision {-}

If you add two regular 64-bit `:float64` values you can quickly notice the imprecision
due to decimal fractions that cannot be represented precisely. For example:
```
> 0.1 + 0.2
0.30000000000000004
```

This happens because the decimal `0.1` and `0.2` cannot be represented exactly
by a `:float64` which is encoded in base 2. For example, if we show `0.1` to 20 digits, we get:
```
> 0.1.show(20)
"0.10000000000000000555"
```

When we convert the constant `0.1` to a `:ddouble` we can see with even more precision how it is
approximated:

```
> 0.1.ddouble
0.1000000000000000055511151231258
```

However, if we convert the number `0.1` directly to a `:ddouble` (instead of going through
a `:float64` constant), we can represent `0.1` more precisely:

```
> ddouble("0.1")   // for convenience; using `ddouble-exp(1,-1)` is more efficient
0.1
```

This is possible because a `:ddouble` uses two `:float64`s to represent numbers; in this
case the first (_hi_) float64 is (the inexact) `0.1` while the second one (_lo_) is the
correction to that number. The `show-sum`  (or `show-hex`) function shows this internal pair of numbers:

```
> ddouble("0.1").show-sum(20)
"0.10000000000000000555 + -5.55111512312578301027e-18"

> ddouble("0.1") + ddouble("0.2")
0.3
```

Generally, a `:ddouble` _d_ is represented by a pair of `:float64`s, _hi_ and _lo_,
such that _d_ equals _hi_+_lo_, where \|_lo_\| ≤ 0.5·ulp(_hi_).
Note that despite the extra precision, underflow/overflow and machine precision are
not as well-defined as with a regular `:float64` [@lin;@shewchuk;@Dekker].
Nevertheless, using a `:ddouble` often prevents instability in practice over
`:float64` computations due to  rounding errors or combining very large and
small numbers.

Take the "thin triangle" problem for example [@Goldberg:float;@Kahan:triangle].
The challenge is to compute
the area of a very thin triangle with sides $a$, $b$, and $c$. The short sides
 $b$ and $c$ are just 3 units in the last place (ulp) shorter than
the longest side $a$. Using the value $s = (a + b + c)/2$, the area of
a triangle is $A = \sqrt{s(s-a)(s-b)(s-c)}$. This is troublesome if $s$ is
close to $a$ since $s-a$ magnifies any rounding error. Let's take $a=7$,
and $b = c = 7/2 + 3\cdot2^{-111}$; according to [@gustafson:posit],
“If the units are in light-years, then the short sides are only longer than
half the long side by $1/200$th the diameter of a proton. Yet that pops the triangle
up to about the width of a doorway at the thickest point”. Using various
128-bit numbers we get:
````
exact          : 3.14784204874900425235885265494550774498...e-16
128-bit ddouble: 3.147842048749004252358852654945\([47]{color:#F88}\)e-16
128-bit posit  : 3.147842048749004252358852654945507744\([39]{color:#F88}\)e-16
128-bit ieee   : 3.\([63481490842332134725920516158057682878]{color:#F88}\)e-16
````
For this kind of example, a `:ddouble` has better precision than a
regular 128-bit IEEE float since it can combine very large and
small values. Note that Kahan [@Kahan:triangle] shows how to rewrite the area equation
to avoid magnifying rounding errors -- in that case the result for IEEE 128-bit floats becomes:
````
128-bit ieee x : 3.147842048749004252358852654945507\([92210]{color:#F88}\)e-16
````

The implementation is based closely on the excellent [QD] C++ library [@Hida:qd;@Hida:qdlib],
and assumes proper 64-bit IEEE `:float64`s with correct rounding.
Integers can be represented precisely up to 30 decimal digits (and a bit more...
up to 2^106^ - 2).

## References {-}

~ Bibliography { caption:"0" }

~~ BibItem { #dekker; bibitem-label:"[1]"; searchterm:"Dekker+A+Floating-Point+Technique+for+Extending+the+Available+Precision" }
T.\ Dekker.
_A Floating-Point Technique for Extending the Available Precision_.
Numerische\ Mathematik, vol. 18 (3), June 1971, 224--242.
[pdf](https://gdz.sub.uni-goettingen.de/dms/resolveppn/?PPN=GDZPPN001170007).
~~

~~ BibItem { #Goldberg:float; bibitem-label:"[2]"}
David Goldberg.
_What Every Computer Scientist Should Know About Floating-point Arithmetic_.
ACM Computing Survey. vol. 23 (1), March 1991. doi: [10.1145/103162.103163](https://doi.org/10.1145/103162.103163).
[pdf](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.22.6768&rep=rep1&type=pdf).
~~

~~ BibItem { #Gustafson:posit; bibitem-label:"[3]"; searchterm:"Gustafson+John+Beating+floating+point+at+its+own+game+posit+arithmetic"; }
John\ L. Gustafson and Isaac\ Yonemoto.
_Beating Floating Point at its Own Game: Posit Arithmetic_. 2017.
[pdf](https://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf).
~~

~~ BibItem { #Hida:qd; bibitem-label:"[4]"; searchterm:"Hida+Quad-Double+Arithmetic:+Algorithms,+Implementation,+and+Application" }
Yozo Hida,\ Xiaoye S.\ Li, and David\ H.\ Bailey.
_Quad-Double Arithmetic: Algorithms, Implementation, and Application_.
Lawrence Berkeley National Laboratory Technical Report LBNL-46996. 2000. [pdf](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.4.5769).
~~

~~ BibItem { #Hida:qdlib; bibitem-label:"[5]"; searchterm:"Hida+Library+for+double-double+and+quad-double+arithmetic" }
Yozo Hida,\ Xiaoye\ S.\ Li, and David\ H.\ Bailey.
_Library for double-double and quad-double arithmetic._
(2007). [pdf](https://www.jaist.ac.jp/~s1410018/papers/qd.pdf).
~~

~~ Bibitem { #lin; bibitem-label:"[6]"; searchterm:'Linnainmaa+"Software+for+Doubled-Precision+Floating-Point+Computations"' }
Seppo\ Linnainmaa.
_Software for Doubled-Precision Floating-Point Computations_.
ACM Transactions on Mathematical Software (TOMS), vol. 7 (3), Sept. 1981, 272--283.
~~

~~ Bibitem { #Kahan:triangle; bibitem-label:"[7]"}
William\ Kahan.
_Miscalculating Area and Angles of a Needle-like Triangle_.
Lecture notes, 2004, [pdf](https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf).
~~

~~ BibItem { #shewchuk; bibitem-label:"[8]"; }
Jonathan\ Richard\ Shewchuk.
_Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates_.
Discrete & Computational Geometry, vol. 18, 305--363, 1997. [pdf](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
~~

~

[ddwiki]: https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic
[qd]:     https://crd-legacy.lbl.gov/~dhbailey/mpdist

\/
*/
module std/num/ddoublestd/num/ddouble

import std/core/undivstd/core/undiv
import std/num/float64std/num/float64
import std/num/decimalstd/num/decimal
import std/text/parsestd/text/parse

/* The `:ddouble` type implements [double double][ddwiki] 128-bit floating point numbers
as a pair of IEEE `:float64` values. This extends the precision to 31 decimal digits
(versus 15 for `:float64`), but keeps the same range as
a `:float64` with a maximum  value of about 1.8·10^308^. Because
`:float64`s usually have hardware support, a `:ddouble` is usually much faster
than arbitrary precision floating point numbers.

Internally a `:ddouble` _d_ is represented as a pair of `:float64`s, _hi_ and _lo_,
such that the number represented by _d_ is _hi_+_lo_, where \|_lo_\| ≤ 0.5·ulp(_hi_).
*/
abstract value struct ddoublestd/num/ddouble/ddouble: V
  histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 : float64std/core/types/float64: V
  lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64 : float64std/core/types/float64: V

// Create a `:ddouble` from a `:float64`.
pub fun float64/ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble( dd: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V
  Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(dd: float64,0.0literal: float64
hex64= 0x0p+0
) val maxprecisestd/num/ddouble/maxprecise: int : intstd/core/types/int: V = 9007199254740991literal: int
dec = 9007199254740991
hex64= 0x001FFFFFFFFFFFFF
bit64= 0b0000000000011111111111111111111111111111111111111111111111111111
val minprecisestd/num/ddouble/minprecise: int : intstd/core/types/int: V = ~std/core/int/(~): (i : int) -> intmaxprecisestd/num/ddouble/maxprecise: int fun is-precisestd/num/ddouble/is-precise: (i : int) -> bool(ii: int : intstd/core/types/int: V )result: -> total bool : boolstd/core/types/bool: V (ii: int >=std/core/int/(>=): (x : int, y : int) -> bool minprecisestd/num/ddouble/minprecise: int &&std/core/types/(&&): (x : bool, y : bool) -> bool ii: int <=std/core/int/(<=): (x : int, y : int) -> bool maxprecisestd/num/ddouble/maxprecise: int) // Create a `:ddouble` from an `:int`. // A `:ddouble` can represent integers precisely up to 30 digits. // If an integer is passed that is out of range an infinity is returned. pub fun int/ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble( ii: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V ddouble-int-expstd/num/ddouble/ddouble-int-exp: (i : int, e : int) -> ddouble(ii: int,0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
) fun ddouble-int-expstd/num/ddouble/ddouble-int-exp: (i : int, e : int) -> ddouble( ii: int : intstd/core/types/int: V, ee: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if ii: int.is-precisestd/num/ddouble/is-precise: (i : int) -> bool then small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( ii: int, ee: int ) else // ensure we convert at least 31 (and at most 42) digits precisely val pp: int = ii: int.count-digitsstd/core/int/count-digits: (i : int) -> int val xx: int = ii: int val pxpx: int = pp: int -std/core/int/(-): (x : int, y : int) -> int 14literal: int
dec = 14
hex8 = 0x0E
bit8 = 0b00001110
val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)hihi: int,yy: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = xx: int.cdivmod-exp10std/core/int/cdivmod-exp10: (i : int, n : int) -> (int, int)(pxpx: int) val pypy: int = pxpx: int -std/core/int/(-): (x : int, y : int) -> int 14literal: int
dec = 14
hex8 = 0x0E
bit8 = 0b00001110
if pypy: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then // if 28 or less digits we can convert with two doubles // trace("ddouble-i: " + i.show + ", hi:" + hi.show + ", lo: " + y.show) small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble(hihi: int, pxpx: int +std/core/int/(+): (x : int, y : int) -> int ee: int) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( yy: int, ee: int) else // otherwise use 2 initial 14 digits and a trailing value val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)midmid: int,zz: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = yy: int.cdivmod-exp10std/core/int/cdivmod-exp10: (i : int, n : int) -> (int, int)(pypy: int) val pzpz: int = pypy: int -std/core/int/(-): (x : int, y : int) -> int 14literal: int
dec = 14
hex8 = 0x0E
bit8 = 0b00001110
val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)lolo: int,ploplo: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = if pzpz: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)zz: int,0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) else (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)zz: int.cdiv-exp10std/core/int/cdiv-exp10: (i : int, n : int) -> int(pzpz: int), pzpz: int)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) //trace("ddouble-i: " + i.show + ", hi:" + hi.show + ", mid: " + mid.show + ", lo: " + lo.show) //trace(" px: " + px.show + ", py : " + py.show + ", plo: " + plo.show + ", e: " + e.show ) small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( hihi: int, pxpx: int +std/core/int/(+): (x : int, y : int) -> int ee: int) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( midmid: int, pypy: int +std/core/int/(+): (x : int, y : int) -> int ee: int) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( lolo: int, ploplo: int +std/core/int/(+): (x : int, y : int) -> int ee: int)
) fun small-expstd/num/ddouble/small-exp: (i : int, e : int) -> ddouble( ii: int : intstd/core/types/int: V, ee: int : intstd/core/types/int: V )result: -> total ddouble val dddd: ddouble = ii: int.float64std/num/float64/float64: (i : int) -> float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble if ee: int.is-zerostd/core/int/is-zero: (x : int) -> bool then dddd: ddouble else dddd: ddouble.mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble(ee: int) // &pi; pub val dd-pistd/num/ddouble/dd-pi: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p+1literal: float64
hex64= 0x1.921fb54442d18p1
, 0x1.1A62633145C07p-53literal: float64
hex64= 0x1.1a62633145c07p-53
) // 3.1415926535897932384626433832795028841971693993751 // pub val dd-pi-qd = Ddouble(3.1415926535897931160, 1.2246467991473529607e-16) // QD's pi constant is one bit off... // 2&pi; pub val dd-twopistd/num/ddouble/dd-twopi: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p+2literal: float64
hex64= 0x1.921fb54442d18p2
, 0x1.1A62633145C07p-52literal: float64
hex64= 0x1.1a62633145c07p-52
) // &pi;/2 pub val dd-pi2std/num/ddouble/dd-pi2: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p+0literal: float64
hex64= 0x1.921fb54442d18p0
, 0x1.1A62633145C07p-54literal: float64
hex64= 0x1.1a62633145c07p-54
) // &pi;/4 pub val dd-pi4std/num/ddouble/dd-pi4: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.921FB54442D18p-1literal: float64
hex64= 0x1.921fb54442d18p-1
, 0x1.1A62633145C07p-55literal: float64
hex64= 0x1.1a62633145c07p-55
) // QD: 7.853981633974482790e-01,3.061616997868383018e-17) // &pi;/16 val dd-pi16std/num/ddouble/dd-pi16: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0x1.921FB54442D18p+1literal: float64
hex64= 0x1.921fb54442d18p1
, 0x1.1A62633145C07p-53literal: float64
hex64= 0x1.1a62633145c07p-53
) // 3&pi;/4 pub val dd-pi34std/num/ddouble/dd-pi34: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.2D97C7F3321D2p+1literal: float64
hex64= 0x1.2d97c7f3321d2p1
, 0x1.A79394C9E8A0Bp-54literal: float64
hex64= 0x1.a79394c9e8a0bp-54
) // QD: 2.356194490192344837e+00,9.1848509936051484375e-17) // QD is again one bit off // The _e_ constant. pub val dd-estd/num/ddouble/dd-e: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.5BF0A8B145769p+1literal: float64
hex64= 0x1.5bf0a8b145769p1
, 0x1.4D57EE2B1013Ap-53literal: float64
hex64= 0x1.4d57ee2b1013ap-53
) // 2.7182818284590452353602874713526624977572470937 // The natural logarithm of 2 pub val dd-ln2std/num/ddouble/dd-ln2: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.62E42FEFA39EFp-1literal: float64
hex64= 0x1.62e42fefa39efp-1
, 0x1.ABC9E3B39803Fp-56literal: float64
hex64= 0x1.abc9e3b39803fp-56
) // 0.69314718055994530941723212145817656807550013436026 // The natural logarithm of 10 pub val dd-ln10std/num/ddouble/dd-ln10: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.26BB1BBB55516p+1literal: float64
hex64= 0x1.26bb1bbb55516p1
, -0x1.F48AD494EA3E9p-53literal: float64
hex64= -0x1.f48ad494ea3e9p-53
) // 2.3025850929940456840179914546843642076011014886288 // The base-2 logarithm of _e_. pub val dd-log2estd/num/ddouble/dd-log2e: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.71547652B82FEp+0literal: float64
hex64= 0x1.71547652b82fep0
, 0x1.777D0FFDA0D23p-56literal: float64
hex64= 0x1.777d0ffda0d23p-56
) // 1.442695040888963407359924681001892137426645954153 // The base-10 logarithm of _e_. pub val dd-log10estd/num/ddouble/dd-log10e: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.BCB7B1526E50Ep-2literal: float64
hex64= 0x1.bcb7b1526e50ep-2
, 0x1.95355BAAAFAD3p-57literal: float64
hex64= 0x1.95355baaafad3p-57
) // 0.43429448190325182765112891891660508229439700580367 // The square-root of 2 pub val dd-sqrt2std/num/ddouble/dd-sqrt2: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.6A09E667F3BCDp+0literal: float64
hex64= 0x1.6a09e667f3bcdp0
, -0x1.BDD3413B26456p-54literal: float64
hex64= -0x1.bdd3413b26456p-54
) // 1.4142135623730950488016887242096980785696718753769 // `1.0 / sqrt(2.0)` pub val dd-sqrt12std/num/ddouble/dd-sqrt12: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.6A09E667F3BCDp-1literal: float64
hex64= 0x1.6a09e667f3bcdp-1
, -0x1.BDD3413B26456p-55literal: float64
hex64= -0x1.bdd3413b26456p-55
) // 0.70710678118654752440084436210484903928483593768847 // [Euler's constant](https://en.wikipedia.org/wiki/Euler%E2%80%93Mascheroni_constant) pub val dd-eulerstd/num/ddouble/dd-euler: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.2788CFC6FB619p-1literal: float64
hex64= 0x1.2788cfc6fb619p-1
, -0x1.6CB90701FBFABp-58literal: float64
hex64= -0x1.6cb90701fbfabp-58
) //0.57721566490153286060651209008240243104215933593992 // The 'machine epsilon': this is not well-defined for a `:ddouble` in general since // the difference between 1.0 and the next representable `:ddouble` value is `dd-true-min`. // Instead, we take the square of `flt-epsilon`, i.e. 2^-104^. pub val dd-epsilonstd/num/ddouble/dd-epsilon: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1p-104literal: float64
hex64= 0x1p-104
,0.0literal: float64
hex64= 0x0p+0
) // 8*dd-epsilon val dd-epsilon8std/num/ddouble/dd-epsilon8: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(0x1.0p-101literal: float64
hex64= 0x1p-101
, 0x0.0p+0literal: float64
hex64= 0x0p+0
) // The maximum representable `:ddouble` pub val dd-maxstd/num/ddouble/dd-max: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0x1.FFFFFFFFFFFFFp+1023literal: float64
hex64= 0x1.fffffffffffffp1023
, 0x1.FFFFFFFFFFFFFp+969literal: float64
hex64= 0x1.fffffffffffffp969
) // The smallest positive `:ddouble` that is still normalized pub val dd-minstd/num/ddouble/dd-min: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( flt-minstd/num/float64/flt-min: float64, 0.0literal: float64
hex64= 0x0p+0
) // normalized = 2^-1022 // The smallest positive `:ddouble` (which is subnormal). pub val dd-true-minstd/num/ddouble/dd-true-min: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( flt-true-minstd/num/float64/flt-true-min: float64, 0.0literal: float64
hex64= 0x0p+0
) // Not-A-Number pub val dd-nanstd/num/ddouble/dd-nan: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(nanstd/num/float64/nan: float64,0.0literal: float64
hex64= 0x0p+0
) // Positive infinity pub val dd-posinfstd/num/ddouble/dd-posinf: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(posinfstd/num/float64/posinf: float64,0.0literal: float64
hex64= 0x0p+0
) // Negative infinity pub val dd-neginfstd/num/ddouble/dd-neginf: ddouble = Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(neginfstd/num/float64/neginf: float64,0.0literal: float64
hex64= 0x0p+0
) // maximal precision in decimal digits of a `:ddouble`. pub val dd-max-precstd/num/ddouble/dd-max-prec: int = 31literal: int
dec = 31
hex8 = 0x1F
bit8 = 0b00011111
val dd-default-precstd/num/ddouble/dd-default-prec: int = -31literal: int
dec = -31
hex8 = 0xE1
bit8 = 0b11100001
/*------------------------------------------------------ Compare etc. ------------------------------------------------------*/ // Zero constant pub val zerostd/num/ddouble/zero: ddouble : ddoublestd/num/ddouble/ddouble: V = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(0.0literal: float64
hex64= 0x0p+0
) // One pub val onestd/num/ddouble/one: ddouble : ddoublestd/num/ddouble/ddouble: V = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(1.0literal: float64
hex64= 0x1p0
) // Ten (`10.ddouble`) pub val tenstd/num/ddouble/ten: ddouble : ddoublestd/num/ddouble/ddouble: V = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(10.0literal: float64
hex64= 0x1.4p3
) val twostd/num/ddouble/two: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(2.0literal: float64
hex64= 0x1p1
) // Is this `:ddouble` equal to is-zero pub fun is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-zerostd/num/float64/is-zero: (d : float64) -> bool // Is this `:ddouble` negative? pub fun is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-negstd/num/float64/is-neg: (d : float64) -> bool // Is this `:ddouble` positive? pub fun is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-posstd/num/float64/is-pos: (d : float64) -> bool // Return the sign of the `:ddouble`. pub fun is-signstd/num/ddouble/is-sign: (x : ddouble) -> order( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total order : orderstd/core/types/order: V if xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64==std/num/float64/(==): (x : float64, y : float64) -> bool0.0literal: float64
hex64= 0x0p+0
then Eqstd/core/types/Eq: order elif xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 <std/num/float64/(<): (x : float64, y : float64) -> bool 0.0literal: float64
hex64= 0x0p+0
then Ltstd/core/types/Lt: order else Gtstd/core/types/Gt: order
// Is this `:ddouble` not-a-number? pub fun is-nanstd/num/ddouble/is-nan: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-nanstd/num/float64/is-nan: (d : float64) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-nanstd/num/float64/is-nan: (d : float64) -> bool // Is this a finite `:ddouble`? (i.e. not `is-nan` or `is-inf`) pub fun is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool &&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool // Is this an infinite value. pub fun is-infstd/num/ddouble/is-inf: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-infstd/num/float64/is-inf: (d : float64) -> bool // Does `x` equal positive infinity? pub fun is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-posinfstd/num/float64/is-posinf: (d : float64) -> bool // Does `x` equal negative infinity? pub fun is-neginfstd/num/ddouble/is-neginf: (x : ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-neginfstd/num/float64/is-neginf: (d : float64) -> bool // Compare two `:ddouble` values. pub fun cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total order : orderstd/core/types/order: V match cmpstd/num/float64/cmp: (x : float64, y : float64) -> order(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) Eqstd/core/types/Eq: order -> cmpstd/num/float64/cmp: (x : float64, y : float64) -> order(xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64,yy: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64) ordord: order -> ordord: order pub fun (>)std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool (xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) ==std/core/order/(==): (x : order, y : order) -> bool Gtstd/core/types/Gt: order } pub fun (>=)std/num/ddouble/(>=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) !=std/core/order/(!=): (x : order, y : order) -> bool Ltstd/core/types/Lt: order } pub fun (==)std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) ==std/core/order/(==): (x : order, y : order) -> bool Eqstd/core/types/Eq: order } pub fun (!=)std/num/ddouble/(!=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) !=std/core/order/(!=): (x : order, y : order) -> bool Eqstd/core/types/Eq: order } pub fun (<)std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool (xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) ==std/core/order/(==): (x : order, y : order) -> bool Ltstd/core/types/Lt: order } pub fun (<=)std/num/ddouble/(<=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V)result: -> total bool : boolstd/core/types/bool: V { cmpstd/num/ddouble/cmp: (x : ddouble, y : ddouble) -> order(xx: ddouble,yy: ddouble) !=std/core/order/(!=): (x : order, y : order) -> bool Gtstd/core/types/Gt: order } // The minimum of `x` and `y`. pub fun minstd/num/ddouble/min: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble <=std/num/ddouble/(<=): (x : ddouble, y : ddouble) -> bool yy: ddouble then xx: ddouble else yy: ddouble // The maximum of `x` and `y` pub fun maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble >=std/num/ddouble/(>=): (x : ddouble, y : ddouble) -> bool yy: ddouble then xx: ddouble else yy: ddouble /*------------------------------------------------------ Addition ------------------------------------------------------*/ value struct edoublestd/num/ddouble/edouble: V numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64 : float64std/core/types/float64: V errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 : float64std/core/types/float64: V // often called `twosum` in literature (see [@shewchuk]) fun sumstd/num/ddouble/sum: (x : float64, y : float64) -> edouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V val zz: float64 = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64 val diffdiff: float64 = zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64 val errerr: float64 = (xx: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64)) +std/num/float64/(+): (x : float64, y : float64) -> float64 (yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64) Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64) fun dsumstd/num/ddouble/dsum: (x : float64, y : float64) -> ddouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val zz: float64 = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64 val diffdiff: float64 = zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64 val errerr: float64 = (xx: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64)) +std/num/float64/(+): (x : float64, y : float64) -> float64 (yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 diffdiff: float64) Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64) // As `sum` but with `x.abs >= y.abs` fun quicksumstd/num/ddouble/quicksum: (x : float64, y : float64) -> edouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V val zz: float64 = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64 val errerr: float64 = yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64) Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64) fun dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if !std/core/types/bool/(!): (b : bool) -> boolxx: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool returnreturn: ddouble ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(xx: float64)std/core/types/Unit: () val zz: float64 = xx: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 yy: float64 val errerr: float64 = yy: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 (zz: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: float64) Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(zz: float64, if zz: float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool then errerr: float64 else zz: float64) // Add two `:ddouble`s pub fun (+)std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val z1z1: edouble = sumstd/num/ddouble/sum: (x : float64, y : float64) -> edouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) val lolo: edouble = sumstd/num/ddouble/sum: (x : float64, y : float64) -> edouble(xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64,yy: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64) val e1e1: float64 = z1z1: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 lolo: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64 val z2z2: edouble = quicksumstd/num/ddouble/quicksum: (x : float64, y : float64) -> edouble(z1z1: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64,e1e1: float64) val e2e2: float64 = z2z2: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 lolo: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(z2z2: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64,e2e2: float64) // Create a `:ddouble` as the sum of two `:float64`'s. pub fun ddoublestd/num/ddouble/ddouble: (x : float64, y : float64) -> ddouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if yy: float64.is-zerostd/num/float64/is-zero: (d : float64) -> bool then ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(xx: float64) else dsumstd/num/ddouble/dsum: (x : float64, y : float64) -> ddouble(xx: float64,yy: float64) // Negate a `:ddouble`. pub fun (~)std/num/ddouble/(~): (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(~std/num/float64/(~): (f : float64) -> float64xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,~std/num/float64/(~): (f : float64) -> float64xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64) // Subtract two values. pub fun (-)std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (~std/num/ddouble/(~): (x : ddouble) -> ddoubleyy: ddouble) // Return the absolute value. pub fun absstd/num/ddouble/abs: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddoublexx: ddouble else xx: ddouble // Increment by one. pub fun incstd/num/ddouble/inc: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble :ddoublestd/num/ddouble/ddouble: V xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble // Decrement by one. pub fun decstd/num/ddouble/dec: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble :ddoublestd/num/ddouble/ddouble: V xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble /*------------------------------------------------------ Multiply and divide ------------------------------------------------------*/ /* val splitter = 0x1.0000002p+27 // 134217729.0 = 2^27 + 1 val splitbound = 0x1.0p996 // 6.696928794914171e+299 = 2^996 // Note, QD seems one bit off val two28 = 0x1.0p28 // 268435456.0 = 2^28 fun split( x : float64 ) : (float64,float64) if x > splitbound || x < ~splitbound val y = x * 0x1.0p-28 // 3.7252902984619140625e-09 = 2^-28 val t = y * splitter val hi = t - (t - y) val lo = y - hi (hi*two28, lo*two28) else val t = x * splitter val hi = t - (t - x) val lo = x - hi (hi,lo) // often called `twoproduct` in literature (see [@shewchuk]) fun prod( x : float64, y : float64 ) : edouble val z = x*y val (xhi,xlo) = split(x) val (yhi,ylo) = split(y) val err = ((xhi*yhi - z) + (xhi*ylo + xlo*yhi)) + (xlo*ylo) Edouble(z,if z.is-finite then err else z) fun prodsqr( x : float64 ) : edouble val z = x*x val (h,l) = split(x) val err = ((h*h - z) + (2.0*h*l)) + (l*l) Edouble(z,if z.is-finite then err else z) */ // often called `twoproduct` in literature (see [@shewchuk]) fun prodstd/num/ddouble/prod: (x : float64, y : float64) -> edouble( xx: float64 : float64std/core/types/float64: V, yy: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V val zz: float64 = xx: float64*std/num/float64/(*): (x : float64, y : float64) -> float64yy: float64 val errerr: float64 = fmaddstd/num/float64/fmadd: (x : float64, y : float64, z : float64) -> float64(xx: float64,yy: float64,~std/num/float64/(~): (f : float64) -> float64zz: float64) Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64,errerr: float64) fun prodsqrstd/num/ddouble/prodsqr: (x : float64) -> edouble( xx: float64 : float64std/core/types/float64: V )result: -> total edouble : edoublestd/num/ddouble/edouble: V val zz: float64 = xx: float64*std/num/float64/(*): (x : float64, y : float64) -> float64xx: float64 val errerr: float64 = fmaddstd/num/float64/fmadd: (x : float64, y : float64, z : float64) -> float64(xx: float64,xx: float64,~std/num/float64/(~): (f : float64) -> float64zz: float64) Edoublestd/num/ddouble/Edouble: (num : float64, err : float64) -> edouble(zz: float64,errerr: float64) // Multiply two `:ddouble`s pub fun (*)std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val zz: edouble = prodstd/num/ddouble/prod: (x : float64, y : float64) -> edouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64, yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) val ee: float64 = zz: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 (xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64yy: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(zz: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64, ee: float64) // Multiply `x` with itself. pub fun sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val zz: edouble = prodsqrstd/num/ddouble/prodsqr: (x : float64) -> edouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) val ee: float64 = (zz: edouble.errstd/num/ddouble/edouble/err: (edouble : edouble) -> float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 (2.0literal: float64
hex64= 0x1p1
*std/num/float64/(*): (x : float64, y : float64) -> float64xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)) +std/num/float64/(+): (x : float64, y : float64) -> float64 (xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64*std/num/float64/(*): (x : float64, y : float64) -> float64xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64) dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(zz: edouble.numstd/num/ddouble/edouble/num: (edouble : edouble) -> float64,ee: float64
) // Divide two `:ddouble`s pub fun (/)std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val q1q1: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 /std/num/float64/(/): (x : float64, y : float64) -> float64 yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) if !std/core/types/bool/(!): (b : bool) -> boolq1q1: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool !std/core/types/bool/(!): (b : bool) -> boolyy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-finitestd/num/float64/is-finite: (d : float64) -> bool returnreturn: ddouble q1q1: ddouble val r1r1: ddouble = xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (yy: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble q1q1: ddouble) val q2q2: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(r1r1: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 /std/num/float64/(/): (x : float64, y : float64) -> float64 yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) val r2r2: ddouble = r1r1: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (yy: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble q2q2: ddouble) val q3q3: ddouble = ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(r2r2: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 /std/num/float64/(/): (x : float64, y : float64) -> float64 yy: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) val qq: ddouble = dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(q1q1: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,q2q2: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) qq: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble q3q3: ddouble // Remainder of two `:ddouble`s pub fun (%)std/num/ddouble/(%): (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val nn: ddouble = (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble yy: ddouble).roundstd/num/ddouble/round: (x : ddouble) -> ddouble xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (nn: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubleyy: ddouble) // Division and remainder of two `:ddouble`s pub fun divremstd/num/ddouble/divrem: (x : ddouble, y : ddouble) -> (ddouble, ddouble)( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total (ddouble, ddouble) : (std/core/types/tuple2: (V, V) -> Vddoublestd/num/ddouble/ddouble: V,ddoublestd/num/ddouble/ddouble: V) val nn: ddouble = (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble yy: ddouble).roundstd/num/ddouble/round: (x : ddouble) -> ddouble (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)nn: ddouble, xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (nn: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubleyy: ddouble))std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) /*------------------------------------------------------ Rounding ------------------------------------------------------*/ // Convert a `:ddouble` to a `:float64` (losing precision) pub fun float64std/num/ddouble/float64: (x : ddouble) -> float64( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total float64 : float64std/core/types/float64: V xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 // Convert a `:ddouble` to the nearest integer (rounding to the nearest even number in case of a tie) pub fun intstd/num/ddouble/int: (x : ddouble, nonfin : ? int) -> int( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nonfinnonfin: ? int : intstd/core/types/int: V = 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
)result: -> total int : intstd/core/types/int: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then nonfinnonfin: int else xx: ddouble.roundstd/num/ddouble/round: (x : ddouble) -> ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.intstd/num/decimal/int: (x : decimal, rnd : ? round) -> int
// Round a `:ddouble` to the nearest integer, rounding to the nearest even number in case of a tie. pub fun roundstd/num/ddouble/round: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val rr: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.roundstd/num/float64/round: (d : float64) -> float64 val diffdiff: float64 = rr: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 if diffdiff: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool 0.0literal: float64
hex64= 0x0p+0
then dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(rr: float64,xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.roundstd/num/float64/round: (d : float64) -> float64) elif diffdiff: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool 0.5literal: float64
hex64= 0x1p-1
&&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-negstd/num/float64/is-neg: (d : float64) -> bool then ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(rr: float64 -std/num/float64/(-): (x : float64, y : float64) -> float64 1.0literal: float64
hex64= 0x1p0
) elif diffdiff: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool -0.5literal: float64
hex64= -0x1p-1
&&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.is-posstd/num/float64/is-pos: (d : float64) -> bool then ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(rr: float64 +std/num/float64/(+): (x : float64, y : float64) -> float64 1.0literal: float64
hex64= 0x1p0
) else ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble(rr: float64
) // Round to negative infinity. pub fun floorstd/num/ddouble/floor: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val rr: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.floorstd/num/float64/floor: (d : float64) -> float64 if rr: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 then dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(rr: float64,xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.floorstd/num/float64/floor: (d : float64) -> float64) else Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(rr: float64,0.0literal: float64
hex64= 0x0p+0
) // Round to positive infinity. pub fun ceilingstd/num/ddouble/ceiling: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val rr: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.ceilingstd/num/float64/ceiling: (d : float64) -> float64 if rr: float64 ==std/num/float64/(==): (x : float64, y : float64) -> bool xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 then dquicksumstd/num/ddouble/dquicksum: (x : float64, y : float64) -> ddouble(rr: float64,xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.ceilingstd/num/float64/ceiling: (d : float64) -> float64) else Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(rr: float64,0.0literal: float64
hex64= 0x0p+0
) // Round towards zero. pub fun truncatestd/num/ddouble/truncate: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then xx: ddouble.ceilingstd/num/ddouble/ceiling: (x : ddouble) -> ddouble else xx: ddouble.floorstd/num/ddouble/floor: (x : ddouble) -> ddouble // The fraction of `x` such that `x.truncate + x.fraction == x`. pub fun fractionstd/num/ddouble/fraction: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.truncatestd/num/ddouble/truncate: (x : ddouble) -> ddouble // The _floored_ fraction of `x`. This is always positive, such that `x.floor + x.ffraction == x`. pub fun ffractionstd/num/ddouble/ffraction: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.floorstd/num/ddouble/floor: (x : ddouble) -> ddouble // Round a `:ddouble` to a specified precision. // Uses `round` if the precision is smaller or equal to zero. pub fun round-to-precstd/num/ddouble/round-to-prec: (x : ddouble, prec : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if precprec: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then xx: ddouble.roundstd/num/ddouble/round: (x : ddouble) -> ddouble elif precprec: int >std/core/int/(>): (x : int, y : int) -> bool dd-max-precstd/num/ddouble/dd-max-prec: int then xx: ddouble else val pp: ddouble : ddoublestd/num/ddouble/ddouble: V = powi10std/num/ddouble/powi10: (exp : int) -> ddouble(precprec: int) (xx: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble pp: ddouble).roundstd/num/ddouble/round: (x : ddouble) -> ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble
pp: ddouble /*------------------------------------------------------ Exponents (needed for showing ddouble's) ------------------------------------------------------*/ // 'Load exponent': returns `x`&middot;2^`exp`^. pub fun ldexpstd/num/ddouble/ldexp: (x : ddouble, exp : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(ldexpstd/num/float64/ldexp: (x : float64, e : int) -> float64(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,expexp: int),ldexpstd/num/float64/ldexp: (x : float64, e : int) -> float64(xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64,expexp: int)) fun npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, accacc: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if nn: int <=std/core/int/(<=): (x : int, y : int) -> bool 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then accacc: ddouble elif nn: int.is-evenstd/core/int/is-even: (i : int) -> bool then npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble( sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(xx: ddouble), accacc: ddouble, pretend-decreasingstd/core/undiv/pretend-decreasing: (x : int) -> int(nn: int/std/core/int/(/): (x : int, y : int) -> int2literal: int
dec = 2
hex8 = 0x02
bit8 = 0b00000010
) ) else npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble( xx: ddouble, xx: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubleaccacc: ddouble, pretend-decreasingstd/core/undiv/pretend-decreasing: (x : int) -> int(nn: int.decstd/core/int/dec: (i : int) -> int)
) fun npwrstd/num/ddouble/npwr: (x : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if nn: int==std/core/int/(==): (x : int, y : int) -> bool0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then (if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then dd-nanstd/num/ddouble/dd-nan: ddouble else onestd/num/ddouble/one: ddouble) elif nn: int==std/core/int/(==): (x : int, y : int) -> bool1literal: int
dec = 1
hex8 = 0x01
bit8 = 0b00000001
then xx: ddouble else xx: ddouble.npwr-accstd/num/ddouble/npwr-acc: (x : ddouble, acc : ddouble, n : int) -> ddouble(onestd/num/ddouble/one: ddouble,nn: int
) // Return `x` to the power of `n`. fun powistd/num/ddouble/powi: (x : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val pp: ddouble = npwrstd/num/ddouble/npwr: (x : ddouble, n : int) -> ddouble(xx: ddouble,nn: int.absstd/core/int/abs: (i : int) -> int) if nn: int.is-negstd/core/int/is-neg: (i : int) -> bool then (onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble pp: ddouble) else pp: ddouble // Return 10 to the power of `exp`. fun powi10std/num/ddouble/powi10: (exp : int) -> ddouble( expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V powistd/num/ddouble/powi: (x : ddouble, n : int) -> ddouble(tenstd/num/ddouble/ten: ddouble,expexp: int) // Return 10 to the power of `exp`. pub fun exp10std/num/ddouble/exp10: (exp : ddouble) -> ddouble( expexp: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V powstd/num/ddouble/pow: (x : ddouble, y : ddouble) -> ddouble(tenstd/num/ddouble/ten: ddouble,expexp: ddouble) // Return 2 to the power of `exp`. pub fun exp2std/num/ddouble/exp2: (exp : ddouble) -> ddouble( expexp: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V powstd/num/ddouble/pow: (x : ddouble, y : ddouble) -> ddouble(twostd/num/ddouble/two: ddouble,expexp: ddouble) // Create a `:ddouble` `x` such that `x` equals `d`&middot;10^`e`^. pub fun float64/ddouble-expstd/num/ddouble/float64/ddouble-exp: (d : float64, e : int) -> ddouble( dd: float64 : float64std/core/types/float64: V, ee: int : intstd/core/types/int: V )result: -> total ddouble dd: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble.mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble(ee: int) // Create a `:ddouble` `x` such that `x` equals `i`&middot;10^`e`^. pub fun int/ddouble-expstd/num/ddouble/int/ddouble-exp: (i : int, exp : int) -> ddouble( ii: int : intstd/core/types/int: V, expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V ii: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble.mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble(expexp: int) fun mul-exp10std/num/ddouble/mul-exp10: (x : ddouble, exp : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, expexp: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if expexp: int.is-zerostd/core/int/is-zero: (x : int) -> bool then xx: ddouble else xx: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble powi10std/num/ddouble/powi10: (exp : int) -> ddouble(expexp: int) /*------------------------------------------------------ Show ------------------------------------------------------*/ // Decode a `:ddouble` `d` into two doubles `(hi,lo)` such that `d` equals `hi`+`lo`, // where `lo` &le; 0.5&middot;ulp(`hi`). pub fun decodestd/num/ddouble/decode: (d : ddouble) -> (float64, float64)( dd: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total (float64, float64) : (std/core/types/tuple2: (V, V) -> Vfloat64std/core/types/float64: V,float64std/core/types/float64: V) (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)dd: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64,dd: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) // Encode a `:ddouble` `d` from two doubles `(hi,lo)` such that `d` equals `hi`+`lo`. pub fun encodestd/num/ddouble/encode: (hi : float64, lo : float64) -> ddouble( hihi: float64 : float64std/core/types/float64: V, lolo: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V hihi: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble lolo: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble // Convert a `:ddouble` to a `:decimal` up to a given precision `prec` (= `-1`). // A negative precision converts precisely. Returns 0 for non-finite `:ddouble`'s. pub fun decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = -1literal: int
dec = -1
hex8 = 0xFF
bit8 = 0b11111111
)result: -> total decimal : decimalstd/num/decimal/decimal: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
.decimalstd/num/decimal/int/decimal: (i : int, exp : ? int) -> decimal else xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.decimalstd/num/decimal/float64/decimal: (d : float64, max-prec : ? int) -> decimal(precprec: int) +std/num/decimal/(+): (x : decimal, y : decimal) -> decimal xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.decimalstd/num/decimal/float64/decimal: (d : float64, max-prec : ? int) -> decimal(precprec: int
) /* Show a `:ddouble` `x` precisely as the sum of two `:float64`s in [hexadecimal notation](https://books.google.com/books?id=FgMsCwAAQBAJ&pg=PA41). Use this if you need to guarantee that you can parse back `:ddouble`s exactly, i.e. `x == x.show-hex.ddouble`. ``` > 0.1.ddouble.show-hex "0x1.999999999999Ap-4 + 0x0.0p+0" > "0.1".ddouble.show-hex "0x1.999999999999Ap-4 + -0x1.999999999999Ap-58" > dd-pi.show-hex "0x1.921FB54442D18p+1 + 0x1.1A62633145C07p-53" > dd-max.show-hex "0x1.FFFFFFFFFFFFFp+1023 + 0x1.FFFFFFFFFFFFFp+969" ``` . */ pub fun show-hexstd/num/ddouble/show-hex: (x : ddouble, width : ? int, use-capitals : ? bool, pre : ? string) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, widthwidth: ? int : intstd/core/types/int: V = 1literal: int
dec = 1
hex8 = 0x01
bit8 = 0b00000001
, use-capitalsuse-capitals: ? bool : boolstd/core/types/bool: V = Truestd/core/types/True: bool, prepre: ? string : stringstd/core/types/string: V = "0x"literal: string
count= 2
)result: -> total string : stringstd/core/types/string: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string else xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.show-hexstd/num/float64/show-hex: (d : float64, width : ? int, use-capitals : ? bool, pre : ? string) -> string(widthwidth: int,use-capitalsuse-capitals: bool,prepre: string) ++std/core/types/(++): (x : string, y : string) -> string " + "literal: string
count= 3
++std/core/types/(++): (x : string, y : string) -> string xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.show-hexstd/num/float64/show-hex: (d : float64, width : ? int, use-capitals : ? bool, pre : ? string) -> string(widthwidth: int,use-capitalsuse-capitals: bool,prepre: string
) // Show a `:ddouble` `x` with a given precision `prec` (=`-31`). // The precision specifies the number of digits after the dot (in either scientific of fixed-point notation). // If the precision is negative, _at most_ `prec` digits are displayed, while for a positive // precision, exactly `prec` digits behind the dot are displayed. // This uses `show-fixed` when the exponent of `x` in scientific notation is larger than -5 // and smaller than the precision (or 15 in case of a negative precision), otherwise it uses `show-exp`. pub fun showstd/num/ddouble/show: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = dd-default-precstd/num/ddouble/dd-default-prec: int )result: -> total string : stringstd/core/types/string: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string else xx: ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.showstd/num/decimal/show: (d : decimal, prec : ? int) -> string(precprec: int) /* Show a ddouble `x` with a given precision `prec` (=`-31`) in scientific notation. The precision specifies the number of digits after the dot, i.e. the number of significant digits is `prec+1`. If the precision is negative, _at most_ `prec` digits are displayed, and if it is positive exactly `prec` digits are used. ``` > 1.1.ddouble.show-exp "1.1000000000000000888178419700125" > 1.1.ddouble.show-exp(-100) "1.100000000000000088817841970012523233890533447265625" > 1.1.ddouble.show-exp(5) "1.10000" > 1.1.ddouble.show-exp(-5) "1.1" ``` . */ pub fun show-expstd/num/ddouble/show-exp: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = dd-default-precstd/num/ddouble/dd-default-prec: int )result: -> total string : stringstd/core/types/string: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string else xx: ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.show-expstd/num/decimal/show-exp: (d : decimal, prec : ? int) -> string(precprec: int) /* Show a ddouble `x` with a given precision `prec` (=`-31`) in fixed-point notation. The precision specifies the number of digits after the dot. If the precision is negative, _at most_ `prec` digits after the dot are displayed, while for a positive precision, exactly `prec` digits are used. ``` > 0.1.ddouble.show-fixed "0.1000000000000000055511151231258" > 0.1.ddouble.show-fixed(-100) "0.1000000000000000055511151231257827021181583404541015625" > 0.1.ddouble.show-fixed(5) "0.10000" > 0.1.ddouble.show-fixed(-5) "0.1" ``` . */ pub fun show-fixedstd/num/ddouble/show-fixed: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = dd-default-precstd/num/ddouble/dd-default-prec: int )result: -> total string : stringstd/core/types/string: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string else xx: ddouble.decimalstd/num/ddouble/decimal: (x : ddouble, prec : ? int) -> decimal.show-fixedstd/num/decimal/show-fixed: (d : decimal, prec : ? int) -> string(precprec: int) // Show a `:ddouble` as the sum of `:float64`s with an optional precision. // Note: use `show-hex` for reliable round-trip parsing. pub fun show-sumstd/num/ddouble/show-sum: (x : ddouble, prec : ? int) -> string( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, precprec: ? int : intstd/core/types/int: V = -17literal: int
dec = -17
hex8 = 0xEF
bit8 = 0b11101111
)result: -> total string : stringstd/core/types/string: V if !std/core/types/bool/(!): (b : bool) -> boolxx: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string else xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string(precprec: int) ++std/core/types/(++): (x : string, y : string) -> string " + "literal: string
count= 3
++std/core/types/(++): (x : string, y : string) -> string xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64.showstd/num/float64/show: (d : float64, precision : ? int) -> string(precprec: int
) /*------------------------------------------------------ Parsing ------------------------------------------------------*/ // Parse a floating point number with up to 31 digits precision. // Return `dd-nan` if the string is an invalid number. pub fun string/ddoublestd/num/ddouble/string/ddouble: (s : string) -> ddouble( ss: string : stringstd/core/types/string: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V ss: string.parse-ddoublestd/num/ddouble/parse-ddouble: (s : string) -> maybe<ddouble>.defaultstd/core/maybe/default: (m : maybe<ddouble>, nothing : ddouble) -> ddouble(dd-nanstd/num/ddouble/dd-nan: ddouble) /* // Parse a floating point number with up to 31 digits precision. // Returns `Nothing` on invalid input. pub fun parse-ddouble( s : string ) : maybe<ddouble> val t = s.trim match(t.find(rx-ddouble)) Just(cap) -> val whole = cap.groups[1] val frac = cap.groups[2].trim-right("0") val exp = cap.groups[3].parse-int.default(0) val w = (whole + frac).parse-int.default(0) val e = exp - frac.count val x = ddouble-int-exp(w,e) // trace("parse: s: " + s + "\n w: " + w.show + ", e: " + e.show + "\n wx: " + w.ddouble.show-sum + ", x: " + x.show + "\n " + x.show-sum) Just(x) Nothing -> match(t.find(rx-special)) Just(cap) -> val special = if !cap.groups[1].is-empty then dd-nan elif cap.groups[2] == "-" then dd-neginf else dd-posinf Just(special) Nothing -> match(t.find(rx-float64-sum)) Just(cap) -> match cap.groups[1].parse-float64 Just(hi) -> match cap.groups[2].parse-float64 Just(lo) -> Just( hi.ddouble + lo.ddouble ) Nothing -> Nothing Nothing -> Nothing Nothing -> Nothing val rx-ddouble = regex(r"^([\-\+]?\d+)(?:\.(\d*))?(?:[eE]([\-\+]?\d+))?$") val rx-float64-sum = regex(r"^(\S*\d) *\+ *(\S+)$") val rx-special = regex(r"^(nan)|([\+\-])?inf(?:inity)?$",ignorecase=True) */ pub fun parse-ddoublestd/num/ddouble/parse-ddouble: (s : string) -> maybe<ddouble>( ss: string : stringstd/core/types/string: V )result: -> total maybe<ddouble> : maybestd/core/types/maybe: V -> V<ddoublestd/num/ddouble/ddouble: V> ss: string.trimstd/core/string/trim: (s : string) -> string.to-lowerstd/core/string/to-lower: (s : string) -> string.slicestd/core/sslice/slice: (s : string) -> sslice.parse-eofstd/text/parse/parse-eof: (input : sslice, p : () -> parse ddouble) -> parse-error<ddouble>(pddoublestd/num/ddouble/pddouble: () -> parse ddouble).maybestd/text/parse/maybe: (perr : parse-error<ddouble>) -> maybe<ddouble> pub fun pddoublestd/num/ddouble/pddouble: () -> parse ddouble()result: -> parse ddouble : parsestd/text/parse/parse: (E, V) -> V ddoublestd/num/ddouble/ddouble: V (pddouble-sumstd/num/ddouble/pddouble-sum: () -> parse ddouble ||std/text/parse/(||): (p1 : parser<total,ddouble>, p2 : parser<total,ddouble>) -> parse ddouble pddouble-normalstd/num/ddouble/pddouble-normal: () -> parse ddouble) fun pddouble-sumstd/num/ddouble/pddouble-sum: () -> parse ddouble()result: -> parse ddouble : parsestd/text/parse/parse: (E, V) -> V ddoublestd/num/ddouble/ddouble: V val hihi: float64 = pdoublestd/num/float64/pdouble: () -> parse float64() if !std/core/types/bool/(!): (b : bool) -> parse boolhihi: float64.is-finitestd/num/float64/is-finite: (d : float64) -> parse bool returnreturn: ddouble hihi: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> parse ddouble // NaN or +-Infinity pstringstd/text/parse/pstring: (s : string) -> parse string(" + "literal: string
count= 3
) val lolo: float64 = pdoublestd/num/float64/pdouble: () -> parse float64() returnreturn: ddouble (hihi: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> parse ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> parse ddouble lolo: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> parse ddouble
) fun pddouble-normalstd/num/ddouble/pddouble-normal: () -> parse ddouble()result: -> parse ddouble : parsestd/text/parse/parse: (E, V) -> V ddoublestd/num/ddouble/ddouble: V val negneg: bool = signstd/text/parse/sign: () -> parse bool() val wholewhole: string = digitsstd/text/parse/digits: () -> parse string() val fracfrac: string = optionalstd/text/parse/optional: (default : string, p : parser<total,string>) -> parse string(""literal: string
count= 0
, { charstd/text/parse/char: (c : char) -> parse char('.'literal: char
unicode= u002E
); digitsstd/text/parse/digits: () -> parse string() } ).trim-rightstd/core/sslice/trim-right: (s : string, sub : string) -> parse string("0"literal: string
count= 1
) val expexp: int = optionalstd/text/parse/optional: (default : int, p : parser<total,int>) -> parse int(0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
, { charstd/text/parse/char: (c : char) -> parse char('e'literal: char
unicode= u0065
); pintstd/text/parse/pint: () -> parse int() } ) val ww: int = ((if negneg: bool then "-"literal: string
count= 1
else ""literal: string
count= 0
) ++std/core/types/(++): (x : string, y : string) -> parse string wholewhole: string ++std/core/types/(++): (x : string, y : string) -> parse string fracfrac: string).parse-intstd/core/int/parse-int: (s : string, hex : ? bool) -> parse maybe<int>.defaultstd/core/maybe/default: (m : maybe<int>, nothing : int) -> parse int(0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
) val ee: int = expexp: int -std/core/int/(-): (x : int, y : int) -> parse int fracfrac: string.countstd/core/string/count: (s : string) -> parse int val xx: ddouble = ddouble-int-expstd/num/ddouble/ddouble-int-exp: (i : int, e : int) -> parse ddouble(ww: int,ee: int) // trace("parse-normal: w: " + w.show + ", e: " ++ e.show ++ "\n wx: " ++ w.ddouble.show-sum ++ ", x: " ++ x.show ++ "\n " ++ x.show-sum)
xx: ddouble /*------------------------------------------------------ Advanced operations ------------------------------------------------------*/ /*------------------------------------------------------ Roots ------------------------------------------------------*/ // The square root of a non-negative `:ddouble` `x`. // For negative `x`, `dd-nan` is returned. pub fun sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V // From the QD library: // Use Karp's strategy: if a is an approximation to sqrt(x), then // sqrt(x) = x*a + (x - (x*a)^2) * a/2 if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble elif xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then dd-nanstd/num/ddouble/dd-nan: ddouble else val aa: float64 = 1.0literal: float64
hex64= 0x1p0
/std/num/float64/(/): (x : float64, y : float64) -> float64 sqrtstd/num/float64/sqrt: (d : float64) -> float64(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) val t1t1: float64 = xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 aa: float64 val t2t2: float64 = (xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(t1t1: float64.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble)).histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 aa: float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 0.5literal: float64
hex64= 0x1p-1
dsumstd/num/ddouble/dsum: (x : float64, y : float64) -> ddouble(t1t1: float64,t2t2: float64
) // The `n`-th root of a `:ddouble` number `x`. // `n` must be positive, and if `n` is even, then // `x` must not be negative. pub fun nrootstd/num/ddouble/nroot: (x : ddouble, n : int) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, nn: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if nn: int==std/core/int/(==): (x : int, y : int) -> bool2literal: int
dec = 2
hex8 = 0x02
bit8 = 0b00000010
then xx: ddouble.sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble elif nn: int==std/core/int/(==): (x : int, y : int) -> bool1literal: int
dec = 1
hex8 = 0x01
bit8 = 0b00000001
then xx: ddouble elif (nn: int<=std/core/int/(<=): (x : int, y : int) -> bool0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
||std/core/types/(||): (x : bool, y : bool) -> bool (nn: int.is-evenstd/core/int/is-even: (i : int) -> bool &&std/core/types/(&&): (x : bool, y : bool) -> bool xx: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool)) then dd-nanstd/num/ddouble/dd-nan: ddouble elif xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else val rr: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble val a0a0: ddouble = expstd/num/float64/exp: (p : float64) -> float64( ~std/num/float64/(~): (f : float64) -> float64lnstd/num/float64/ln: (f : float64) -> float64(rr: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64) /std/num/float64/(/): (x : float64, y : float64) -> float64 nn: int.float64std/num/float64/float64: (i : int) -> float64 ).ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble val a1a1: ddouble = a0a0: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ((a0a0: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble rr: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble powistd/num/ddouble/powi: (x : ddouble, n : int) -> ddouble(a0a0: ddouble,nn: int))) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble nn: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble) (if xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.is-negstd/num/float64/is-neg: (d : float64) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble else onestd/num/ddouble/one: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble a1a1: ddouble
// Multiply `x` by a `:float64` `p` where `p` must be a power of 2. fun mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, pp: float64 : float64std/core/types/float64: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 pp: float64, xx: ddouble.lostd/num/ddouble/ddouble/lo: (ddouble : ddouble) -> float64 *std/num/float64/(*): (x : float64, y : float64) -> float64 pp: float64) val one-halfstd/num/ddouble/one-half: ddouble : ddoublestd/num/ddouble/ddouble: V = 0.5literal: float64
hex64= 0x1p-1
.ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble
fun halfstd/num/ddouble/half: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> ddouble( xx: ddouble, 0.5literal: float64
hex64= 0x1p-1
) fun twicestd/num/ddouble/twice: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> ddouble( xx: ddouble, 2.0literal: float64
hex64= 0x1p1
) /*------------------------------------------------------ Exponents and logarithms ------------------------------------------------------*/ // Return _e_ (`dd-e`) to the power of `x`. pub fun expstd/num/ddouble/exp: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V /* From the QD library: Strategy: We first reduce the size of x by noting that exp(kr + m * ln(2)) = 2^m * exp(r)^k where m and k are integers. By choosing m appropriately we can make |kr| <= ln(2) / 2 = 0.347. Then exp(r) is evaluated using the familiar Taylor series. Reducing the argument substantially speeds up the convergence. */ val kk: float64 : float64std/core/types/float64: V = 512.0literal: float64
hex64= 0x1p9
val inv-kinv-k: float64 = 1.0literal: float64
hex64= 0x1p0
/std/num/float64/(/): (x : float64, y : float64) -> float64 kk: float64 if xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 <=std/num/float64/(<=): (x : float64, y : float64) -> bool -709.0literal: float64
hex64= -0x1.628p9
then zerostd/num/ddouble/zero: ddouble elif xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 >=std/num/float64/(>=): (x : float64, y : float64) -> bool 709.0literal: float64
hex64= 0x1.628p9
then dd-posinfstd/num/ddouble/dd-posinf: ddouble elif xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then onestd/num/ddouble/one: ddouble elif xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolonestd/num/ddouble/one: ddouble then dd-estd/num/ddouble/dd-e: ddouble else val mm: int : intstd/core/types/int: V = floorstd/num/float64/floor: (d : float64) -> (local<$4252>) float64( (xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> (local<$4252>) float64 /std/num/float64/(/): (x : float64, y : float64) -> (local<$4252>) float64 dd-ln2std/num/ddouble/dd-ln2: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> (local<$4252>) float64) +std/num/float64/(+): (x : float64, y : float64) -> (local<$4252>) float64 0.5literal: float64
hex64= 0x1p-1
).intstd/num/float64/int: (f : float64) -> (local<$4252>) int val rr: ddouble = mul-pwr2std/num/ddouble/mul-pwr2: (x : ddouble, p : float64) -> (local<$4252>) ddouble( xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> (local<$4252>) ddouble (mm: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> (local<$4252>) ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> (local<$4252>) ddouble dd-ln2std/num/ddouble/dd-ln2: ddouble), inv-kinv-k: float64 ) val pp: ddouble = rr: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> (local<$4252>) ddouble val tt: ddouble = rr: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$4252>) ddouble pp: ddouble.halfstd/num/ddouble/half: (x : ddouble) -> (local<$4252>) ddouble var accacc: local-var<$4252,ddouble> := exp-approxstd/num/ddouble/exp-approx: (p : ddouble, t : ddouble, r : ddouble, eps : float64, fs : list<ddouble>, s : ? ddouble) -> (local<$4252>) ddouble( pp: ddouble, tt: ddouble, rr: ddouble, inv-kinv-k: float64 *std/num/float64/(*): (x : float64, y : float64) -> (local<$4252>) float64 dd-epsilonstd/num/ddouble/dd-epsilon: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> (local<$4252>) float64, exp-factorsstd/num/ddouble/exp-factors: list<ddouble> ) repeatstd/core/repeat: (n : int, action : () -> (local<$4252>) ()) -> (local<$4252>) ()(9literal: int
dec = 9
hex8 = 0x09
bit8 = 0b00001001
) accacc: local-var<$4252,ddouble> :=std/core/types/local-set: (v : local-var<$4252,ddouble>, assigned : ddouble) -> (local<$4252>) () accacc: ddouble.twicestd/num/ddouble/twice: (x : ddouble) -> (local<$4252>) ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$4252>) ddouble accacc: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> (local<$4252>) ddouble accacc: local-var<$4252,ddouble> :=std/core/types/local-set: (v : local-var<$4252,ddouble>, assigned : ddouble) -> (local<$4252>) () accacc: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$4252>) ddouble onestd/num/ddouble/one: ddouble ldexpstd/num/ddouble/ldexp: (x : ddouble, exp : int) -> (local<$4252>) ddouble( accacc: ddouble, mm: int
) val exp-factorsstd/num/ddouble/exp-factors: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> = [std/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a> Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(1.66666666666666657e-01literal: float64
hex64= 0x1.5555555555555p-3
, 9.25185853854297066e-18literal: float64
hex64= 0x1.5555555555555p-57
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(4.16666666666666644e-02literal: float64
hex64= 0x1.5555555555555p-5
, 2.31296463463574266e-18literal: float64
hex64= 0x1.5555555555555p-59
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(8.33333333333333322e-03literal: float64
hex64= 0x1.1111111111111p-7
, 1.15648231731787138e-19literal: float64
hex64= 0x1.1111111111111p-63
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(1.38888888888888894e-03literal: float64
hex64= 0x1.6c16c16c16c17p-10
, -5.30054395437357706e-20literal: float64
hex64= -0x1.f49f49f49f49fp-65
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(1.98412698412698413e-04literal: float64
hex64= 0x1.a01a01a01a01ap-13
, 1.72095582934207053e-22literal: float64
hex64= 0x1.a01a01a01a01ap-73
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble(2.48015873015873016e-05literal: float64
hex64= 0x1.a01a01a01a01ap-16
, 2.15119478667758816e-23literal: float64
hex64= 0x1.a01a01a01a01ap-76
)
]std/core/types/Nil: forall<a> list<a> fun exp-approxstd/num/ddouble/exp-approx: (p : ddouble, t : ddouble, r : ddouble, eps : float64, fs : list<ddouble>, s : ? ddouble) -> ddouble( pp: ddouble : ddoublestd/num/ddouble/ddouble: V, tt: ddouble : ddoublestd/num/ddouble/ddouble: V, rr: ddouble : ddoublestd/num/ddouble/ddouble: V, epseps: float64 : float64std/core/types/float64: V, fsfs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V>, ss: ? ddouble : ddoublestd/num/ddouble/ddouble: V = zerostd/num/ddouble/zero: ddouble )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V match fsfs: list<ddouble> Nilstd/core/types/Nil: forall<a> list<a> -> ss: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble tt: ddouble Consstd/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a>(ff: ddouble,fs1fs1: list<ddouble>) -> val s1s1: ddouble = ss: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble tt: ddouble val p1p1: ddouble = pp: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble rr: ddouble val t1t1: ddouble = p1p1: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble ff: ddouble if t1t1: ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64.absstd/num/float64/abs: (f : float64) -> float64 <=std/num/float64/(<=): (x : float64, y : float64) -> bool epseps: float64 then ss: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble tt: ddouble else exp-approxstd/num/ddouble/exp-approx: (p : ddouble, t : ddouble, r : ddouble, eps : float64, fs : list<ddouble>, s : ? ddouble) -> ddouble( p1p1: ddouble, t1t1: ddouble, rr: ddouble, epseps: float64, fs1fs1: list<ddouble>, s1s1: ddouble ) // The natural logarithm (in base _e_) of `x`. pub fun lnstd/num/ddouble/ln: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V /* From QD Library: Strategy. The Taylor series for log converges much more slowly than that of exp, due to the lack of the factorial term in the denominator. Hence this routine instead tries to determine the root of the function f(x) = exp(x) - a using Newton iteration. The iteration is given by x' = x - f(x)/f'(x) = x - (1 - a * exp(-x)) = x + a * exp(-x) - 1. Only one iteration is needed, since Newton's iteration approximately doubles the number of digits per iteration. */ if xx: ddouble <=std/num/ddouble/(<=): (x : ddouble, y : ddouble) -> bool zerostd/num/ddouble/zero: ddouble then (if xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolzerostd/num/ddouble/zero: ddouble then dd-neginfstd/num/ddouble/dd-neginf: ddouble else dd-nanstd/num/ddouble/dd-nan: ddouble) elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then zerostd/num/ddouble/zero: ddouble elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool dd-estd/num/ddouble/dd-e: ddouble then onestd/num/ddouble/one: ddouble elif xx: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then xx: ddouble else val a0a0: ddouble = lnstd/num/float64/ln: (f : float64) -> float64(xx: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64).ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble // approximate val a1a1: ddouble = a0a0: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (xx: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble expstd/num/ddouble/exp: (x : ddouble) -> ddouble(~std/num/ddouble/(~): (x : ddouble) -> ddoublea0a0: ddouble) -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble) a1a1: ddouble // Return the logarithm in some base `b` of a `:ddouble` `x` pub fun logstd/num/ddouble/log: (x : ddouble, base : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, basebase: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(basebase: ddouble) // The logarithm in base 10 of `x`. pub fun log10std/num/ddouble/log10: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble dd-ln10std/num/ddouble/dd-ln10: ddouble // The logarithm in base 2 of `x`. pub fun log2std/num/ddouble/log2: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble dd-ln2std/num/ddouble/dd-ln2: ddouble // `x` to the power of `y` both as `:ddouble` pub fun powstd/num/ddouble/pow: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V expstd/num/ddouble/exp: (x : ddouble) -> ddouble(yy: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble)) // Return `ln(1.0 + x)`. // Avoids potential imprecision for small `x` where adding `1.0` explicitly // may lead to rounding errors. pub fun ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then xx: ddouble else val yy: ddouble = onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble xx: ddouble val zz: ddouble = yy: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble if zz: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then xx: ddouble else lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(yy: ddouble) *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble zz: ddouble) // Return `exp(x - 1.0)`. // Avoids rounding errors for values of `x` very close to `1.0`. pub fun expm1std/num/ddouble/expm1: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then xx: ddouble else val yy: ddouble = expstd/num/ddouble/exp: (x : ddouble) -> ddouble(xx: ddouble) if yy: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolonestd/num/ddouble/one: ddouble then xx: ddouble else val ymym: ddouble = yy: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble if ymym: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble then ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble else ymym: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(yy: ddouble)) fun log2p1std/num/ddouble/log2p1: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V dd-log2estd/num/ddouble/dd-log2e: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble(xx: ddouble) fun exp2m1std/num/ddouble/exp2m1: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V expm1std/num/ddouble/expm1: (x : ddouble) -> ddouble(dd-ln2std/num/ddouble/dd-ln2: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble xx: ddouble) // Returns `ln(exp(x) + exp(y))`. // Avoids overlow/underflow errors. pub fun lnaddexpstd/num/ddouble/lnaddexp: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolyy: ddouble then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble dd-ln2std/num/ddouble/dd-ln2: ddouble else val zz: ddouble = xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble yy: ddouble if zz: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble(expstd/num/ddouble/exp: (x : ddouble) -> ddouble(~std/num/ddouble/(~): (x : ddouble) -> ddoublezz: ddouble)) else yy: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ln1pstd/num/ddouble/ln1p: (x : ddouble) -> ddouble(expstd/num/ddouble/exp: (x : ddouble) -> ddouble(zz: ddouble)) // Returns `log2( exp2(x) + exp2(y) )`. // Avoids overlow/underflow errors. pub fun logaddexp2std/num/ddouble/logaddexp2: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> boolyy: ddouble then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble else val zz: ddouble = xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble yy: ddouble if zz: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble log2p1std/num/ddouble/log2p1: (x : ddouble) -> ddouble(exp2std/num/ddouble/exp2: (exp : ddouble) -> ddouble(~std/num/ddouble/(~): (x : ddouble) -> ddoublezz: ddouble)) else yy: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble log2p1std/num/ddouble/log2p1: (x : ddouble) -> ddouble(exp2std/num/ddouble/exp2: (exp : ddouble) -> ddouble(zz: ddouble)) // Return if two `:ddouble`s are nearly equal with respect to some `epsilon` (=`8*dd-epsilon`). // The epsilon is the nearest difference for numbers around 1.0. The routine automatically // scales the epsilon for larger and smaller numbers, and for numbers close to zero. pub fun nearly-eqstd/num/ddouble/nearly-eq: (x : ddouble, y : ddouble, epsilon : ? ddouble) -> bool( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V, epsilonepsilon: ? ddouble : ddoublestd/num/ddouble/ddouble: V = dd-epsilon8std/num/ddouble/dd-epsilon8: ddouble )result: -> total bool : boolstd/core/types/bool: V if xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool yy: ddouble returnreturn: bool Truestd/core/types/True: bool val diffdiff: ddouble = (xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble yy: ddouble).absstd/num/ddouble/abs: (x : ddouble) -> ddouble if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool yy: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool ||std/core/types/(||): (x : bool, y : bool) -> bool diffdiff: ddouble <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool dd-minstd/num/ddouble/dd-min: ddouble then // very close to zero, scale the epsilon for subnormal numbers (twostd/num/ddouble/two: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublediffdiff: ddouble <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool (epsilonepsilon: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble dd-minstd/num/ddouble/dd-min: ddouble)) else val sumsum: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble yy: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble ((twostd/num/ddouble/two: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublediffdiff: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble (if sumsum: ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool dd-maxstd/num/ddouble/dd-max: ddouble then dd-maxstd/num/ddouble/dd-max: ddouble else sumsum: ddouble)) <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool epsilonepsilon: ddouble) // Return if two `:ddouble`s are nearly equal with respect to an `epsilon` of `10*dd-epsilon`. // See also `nearly-eq` which takes an explicit `epsilon`. pub fun (~=)std/num/ddouble/(~=): (x : ddouble, y : ddouble) -> bool(xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total bool : boolstd/core/types/bool: V nearly-eqstd/num/ddouble/nearly-eq: (x : ddouble, y : ddouble, epsilon : ? ddouble) -> bool(xx: ddouble,yy: ddouble) /*------------------------------------------------------ Summing ------------------------------------------------------*/ // Return the sum of a list of doubles. // Uses [Kahan-Babu&scaron;kan-Neumaier summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements) // to minimize rounding errors. This // is more precise as Kahan summation and about as fast.\ // `[1.0e3,1.0e97,1.0e3,-1.0e97].sum == 2000.0`\ // A. Neumaier, _Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen_. // Math. Mechanik, 54:39--51, 1974. pub fun list/sumstd/num/ddouble/list/sum: (xs : list<ddouble>) -> ddouble( xsxs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V var totaltotal: local-var<$1697,ddouble> : ddoublestd/num/ddouble/ddouble: V := zerostd/num/ddouble/zero: ddouble var compcomp: local-var<$1697,ddouble> : ddoublestd/num/ddouble/ddouble: V := zerostd/num/ddouble/zero: ddouble xsxs: list<ddouble>.foreachstd/core/list/foreach: (xs : list<ddouble>, action : (ddouble) -> (local<$1697>) ()) -> (local<$1697>) () fnfn: (x : ddouble) -> (local<$1697>) ()(xx: ddouble) val tt: ddouble = totaltotal: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble xx: ddouble val cc: ddouble = if totaltotal: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> (local<$1697>) ddouble >=std/num/ddouble/(>=): (x : ddouble, y : ddouble) -> (local<$1697>) bool xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> (local<$1697>) ddouble then (totaltotal: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble tt: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble xx: ddouble else (xx: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble tt: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble totaltotal: ddouble compcomp: local-var<$1697,ddouble> :=std/core/types/local-set: (v : local-var<$1697,ddouble>, assigned : ddouble) -> (local<$1697>) () compcomp: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble cc: ddouble totaltotal: local-var<$1697,ddouble> :=std/core/types/local-set: (v : local-var<$1697,ddouble>, assigned : ddouble) -> (local<$1697>) () tt: ddouble totaltotal: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> (local<$1697>) ddouble compcomp: ddouble; // The hypotenuse of `x` and `y`: `sqrt(x*x + y*y)`. // Prevents overflow for large numbers. pub fun hypotstd/num/ddouble/hypot: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val xxxx: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble) val yyyy: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(yy: ddouble) val lolo: ddouble = minstd/num/ddouble/min: (x : ddouble, y : ddouble) -> ddouble(xxxx: ddouble,yyyy: ddouble) val hihi: ddouble = maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(xxxx: ddouble,yyyy: ddouble) if hihi: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else val zz: ddouble = lolo: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble hihi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble zz: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublezz: ddouble ) // The square root of the sum of the squares of three doubles. // Prevents overflow for large numbers. pub fun xyz/hypotstd/num/ddouble/xyz/hypot: (x : ddouble, y : ddouble, z : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V, zz: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val xxxx: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble) val yyyy: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(yy: ddouble) val zzzz: ddouble = absstd/num/ddouble/abs: (x : ddouble) -> ddouble(zz: ddouble) val hihi: ddouble = maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(xxxx: ddouble,yyyy: ddouble),zzzz: ddouble) if hihi: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else hihi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(xxxx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(yyyy: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(zzzz: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) ) // The square root of the sum of squares of a list of doubles. // Prevents overflow for large numbers and uses Kahan-Babu&scaron;kan-Neumaier summation // for precision. pub fun list/hypotstd/num/ddouble/list/hypot: (xs : list<ddouble>) -> ddouble( xsxs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val hihi: ddouble = xsxs: list<ddouble>.abs-maxstd/num/ddouble/list/abs-max: (xs : list<ddouble>) -> ddouble if hihi: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else hihi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble xsxs: list<ddouble>.mapstd/core/list/map: (xs : list<ddouble>, f : (ddouble) -> ddouble) -> list<ddouble>(fnfn: (x : ddouble) -> ddouble(xx: ddouble){ sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble(xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble hihi: ddouble) }).sumstd/num/ddouble/list/sum: (xs : list<ddouble>) -> ddouble.sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble // The maximum of the absolute values. fun abs-maxstd/num/ddouble/abs-max: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble),absstd/num/ddouble/abs: (x : ddouble) -> ddouble(yy: ddouble)) // The maximum of a list of absolute values. fun list/abs-maxstd/num/ddouble/list/abs-max: (xs : list<ddouble>) -> ddouble( xsxs: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V xsxs: list<ddouble>.foldrstd/core/list/foldr: (xs : list<ddouble>, z : ddouble, f : (ddouble, ddouble) -> ddouble) -> ddouble(ddouble/zerostd/num/ddouble/zero: ddouble,fnfn: (x : ddouble, m : ddouble) -> ddouble(xx: ddouble,mm: ddouble){ maxstd/num/ddouble/max: (x : ddouble, y : ddouble) -> ddouble(absstd/num/ddouble/abs: (x : ddouble) -> ddouble(xx: ddouble),mm: ddouble) }) /*------------------------------------------------------ Trigonometry ------------------------------------------------------*/ // The sine function of a given angle in radians. pub fun sinstd/num/ddouble/sin: (rad : ddouble) -> ddouble( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V radrad: ddouble.sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble).fststd/core/types/tuple2/fst: (tuple2 : (ddouble, ddouble)) -> ddouble // The cosine function of a given angle in radians. pub fun cosstd/num/ddouble/cos: (rad : ddouble) -> ddouble( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V radrad: ddouble.sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble).sndstd/core/types/tuple2/snd: (tuple2 : (ddouble, ddouble)) -> ddouble // The tangent of a given angle in radians. pub fun tanstd/num/ddouble/tan: (rad : ddouble) -> ddouble( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)ss: ddouble,cc: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = radrad: ddouble.sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble) ss: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble cc: ddouble // Calculate sine and cosine on an angle in radians. pub fun sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble)( radrad: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total (ddouble, ddouble) : (std/core/types/tuple2: (V, V) -> Vddoublestd/num/ddouble/ddouble: V,ddoublestd/num/ddouble/ddouble: V) // quick approximation for small values if radrad: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.absstd/num/float64/abs: (f : float64) -> float64 <std/num/float64/(<): (x : float64, y : float64) -> bool 1.0e-11literal: float64
hex64= 0x1.5fd7fe1796495p-37
then // return (rad, one - 0.5.ddouble*sqr(rad)) val ss: ddouble = radrad: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (radrad: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble 3literal: int
dec = 3
hex8 = 0x03
bit8 = 0b00000011
.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)) val cc: ddouble = (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble).sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble returnreturn: (ddouble, ddouble) (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)ss: ddouble,cc: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) /* To compute sin(x), we choose integers a, b so that rad = s + a * (pi/2) + b * (pi/16) and |s| <= pi/32. A precomputed sin(pi/16) and Chebyshev series, we can compute very precisely. */ // find a, b, and s. /* val z = (rad / dd-twopi).round val r = rad - (dd-twopi * z) val qa : float64 = floor( (r.hi / dd-pi2.hi) + 0.5) val a = qa.truncate.int val t = r - (dd-pi2 * qa.ddouble) val qb : float64 = floor((t.hi / dd-pi16.hi) + 0.5) val b = qb.truncate.int val s = t - (dd-pi16 * qb.ddouble) */ val x1x1: ddouble = radrad: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble dd-twopistd/num/ddouble/dd-twopi: ddouble val x3x3: ddouble = x1x1: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble x1x1: ddouble.roundstd/num/ddouble/round: (x : ddouble) -> ddouble // s = x - a*pi/2 - b*pi/16 val x32x32: ddouble = x3x3: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble x3x3: ddouble val x34x34: ddouble = x32x32: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble x32x32: ddouble // 4.ddouble * x3 val aa: int = x34x34: ddouble.intstd/num/ddouble/int: (x : ddouble, nonfin : ? int) -> int val bb: int = (8literal: int
dec = 8
hex8 = 0x08
bit8 = 0b00001000
.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (x34x34: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble aa: int.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)).intstd/num/ddouble/int: (x : ddouble, nonfin : ? int) -> int val ss: ddouble = dd-pistd/num/ddouble/dd-pi: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble (x32x32: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble ((8literal: int
dec = 8
hex8 = 0x08
bit8 = 0b00001000
*std/core/int/(*): (int, int) -> intaa: int +std/core/int/(+): (x : int, y : int) -> int bb: int).ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble 16literal: int
dec = 16
hex8 = 0x10
bit8 = 0b00010000
.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)) val s2s2: ddouble = ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble // Use the Chebyshev series for best precision. val sinssins: ddouble = ss: ddouble *std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble ch-factorsstd/num/ddouble/ch-factors: list<ddouble>.foldlstd/core/list/foldl: (xs : list<ddouble>, z : ddouble, f : (ddouble, ddouble) -> ddouble) -> ddouble(ddouble/zerostd/num/ddouble/zero: ddouble) fnfn: (acc : ddouble, f : ddouble) -> ddouble(accacc: ddouble,ff: ddouble) { ff: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble accacc: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoubles2s2: ddouble } val cosscoss: ddouble = (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble).sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble // sinb = sin(b*pi/16) val sinbsinb: ddouble = if bb: int >=std/core/int/(>=): (x : int, y : int) -> bool 0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then sin16std/num/ddouble/sin16: (i : int) -> ddouble(bb: int) else ~std/num/ddouble/(~): (x : ddouble) -> ddoublesin16std/num/ddouble/sin16: (i : int) -> ddouble(~std/core/int/(~): (i : int) -> intbb: int) val cosbcosb: ddouble = sin16std/num/ddouble/sin16: (i : int) -> ddouble(8literal: int
dec = 8
hex8 = 0x08
bit8 = 0b00001000
-std/core/int/(-): (x : int, y : int) -> int bb: int.absstd/core/int/abs: (i : int) -> int) if aa: int==std/core/int/(==): (x : int, y : int) -> bool0literal: int
dec = 0
hex8 = 0x00
bit8 = 0b00000000
then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble, cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) elif aa: int==std/core/int/(==): (x : int, y : int) -> bool1literal: int
dec = 1
hex8 = 0x01
bit8 = 0b00000001
then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble, ~std/num/ddouble/(~): (x : ddouble) -> ddoublecosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) elif aa: int==std/core/int/(==): (x : int, y : int) -> bool -1literal: int
dec = -1
hex8 = 0xFF
bit8 = 0b11111111
then (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble, cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) else (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)~std/num/ddouble/(~): (x : ddouble) -> ddoublesinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble, sinssins: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublesinbsinb: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble cosscoss: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddoublecosbcosb: ddouble
)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) // |a| == 2 val ch-factorsstd/num/ddouble/ch-factors: list<ddouble> : liststd/core/types/list: V -> V<ddoublestd/num/ddouble/ddouble: V> = [std/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a> Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 1.6056491947130061858e-10literal: float64
hex64= 0x1.6116039167de1p-33
, 6.1925234565562595936e-27literal: float64
hex64= 0x1.ea9f4c1702653p-88
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( -2.5052108052208301379e-8literal: float64
hex64= -0x1.ae64561f4804fp-26
, -3.6598195022865791876e-25literal: float64
hex64= -0x1.c5104f49cf934p-82
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.0000027557319223964441402literal: float64
hex64= 0x1.71de3a556b36bp-19
, -2.0315661398415506513e-22literal: float64
hex64= -0x1.eb33da7509008p-73
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( -0.00019841269841269841253literal: float64
hex64= -0x1.a01a01a01a01ap-13
, 6.8577289081075077176e-21literal: float64
hex64= 0x1.0313e2634850bp-67
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.0083333333333333332177literal: float64
hex64= 0x1.1111111111111p-7
, 1.1563735775184918100e-19literal: float64
hex64= 0x1.110a7e6657814p-63
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( -0.16666666666666665741literal: float64
hex64= -0x1.5555555555555p-3
, -9.2518585321663028924e-18literal: float64
hex64= -0x1.5555555162e8fp-57
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 1.0000000000000000000literal: float64
hex64= 0x1p0
, -6.0239567712403467977e-31literal: float64
hex64= -0x1.86fa265ab5ed7p-101
)
]std/core/types/Nil: forall<a> list<a> // Return sin(i*pi/16) for 0 <= i <= 8 fun sin16std/num/ddouble/sin16: (i : int) -> ddouble( ii: int : intstd/core/types/int: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V sin16-tablestd/num/ddouble/sin16-table: vector<ddouble>.atstd/core/vector/at: (v : vector<ddouble>, index : int) -> maybe<ddouble>(ii: int).defaultstd/core/maybe/default: (m : maybe<ddouble>, nothing : ddouble) -> ddouble(dd-nanstd/num/ddouble/dd-nan: ddouble) val sin16-tablestd/num/ddouble/sin16-table: vector<ddouble> : vectorstd/core/types/vector: V -> V<ddoublestd/num/ddouble/ddouble: V> = [std/core/types/Cons: forall<a> (head : a, tail : list<a>) -> list<a> ddouble/zerostd/num/ddouble/zero: ddouble, Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.19509032201612827584literal: float64
hex64= 0x1.8f8b83c69a60bp-3
, -7.9910790684617312634e-18literal: float64
hex64= -0x1.26d19b9ff8d82p-57
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.38268343236508978178literal: float64
hex64= 0x1.87de2a6aea963p-2
, -1.0050772696461587612e-17literal: float64
hex64= -0x1.72cedd3d5a61p-57
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.55557023301960217765literal: float64
hex64= 0x1.1c73b39ae68c8p-1
, 4.7094109405616768214e-17literal: float64
hex64= 0x1.b25dd267f66p-55
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.70710678118654757274literal: float64
hex64= 0x1.6a09e667f3bcdp-1
, -4.8336466567264561092e-17literal: float64
hex64= -0x1.bdd3413b26455p-55
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.83146961230254523567literal: float64
hex64= 0x1.a9b66290ea1a3p-1
, 1.4073856984728100930e-18literal: float64
hex64= 0x1.9f630e8b6dafp-60
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.92387953251128673848literal: float64
hex64= 0x1.d906bcf328d46p-1
, 1.7645047084336683223e-17literal: float64
hex64= 0x1.457e610231ac4p-56
), Ddoublestd/num/ddouble/Ddouble: (hi : float64, lo : float64) -> ddouble( 0.98078528040323043058literal: float64
hex64= 0x1.f6297cff75cbp-1
, 1.8546939997825014970e-17literal: float64
hex64= 0x1.562172a361fd6p-56
), onestd/num/ddouble/one: ddouble ]std/core/types/Nil: forall<a> list<a>.vectorstd/core/vector/list/vector: (xs : list<ddouble>) -> vector<ddouble>
// Return `x` with the sign of `y`. pub fun with-sign-ofstd/num/ddouble/with-sign-of: (x : ddouble, y : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V, yy: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if yy: ddouble.is-negstd/num/ddouble/is-neg: (x : ddouble) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddouble(xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble) else xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble // The arc-sine of `x`. Returns the angle in radians. pub fun asinstd/num/ddouble/asin: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val aa: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble if aa: ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble elif aa: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-pi2std/num/ddouble/dd-pi2: ddouble.with-sign-ofstd/num/ddouble/with-sign-of: (x : ddouble, y : ddouble) -> ddouble(xx: ddouble) else atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble( xx: ddouble, sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble) ) // The arc-cosine of `x`. Returns the angle in radians. pub fun acosstd/num/ddouble/acos: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V val aa: ddouble = xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble if aa: ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble elif aa: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then (if xx: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else dd-pistd/num/ddouble/dd-pi: ddouble) else atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble( sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble), xx: ddouble ) // The arc-tangent of `x`. Returns the angle in radians. pub fun atanstd/num/ddouble/atan: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble(xx: ddouble,onestd/num/ddouble/one: ddouble) // The arc-tangent of a point (`x`,`y`). Returns the angle with respect to the x-axis in radians between -&pi; and &pi;. pub fun atan2std/num/ddouble/atan2: (y : ddouble, x : ddouble) -> ddouble( yy: ddouble : ddoublestd/num/ddouble/ddouble: V, xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then if yy: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble /* dd-nan */ else dd-pi2std/num/ddouble/dd-pi2: ddouble.with-sign-ofstd/num/ddouble/with-sign-of: (x : ddouble, y : ddouble) -> ddouble(yy: ddouble) elif yy: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then if xx: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble else dd-pistd/num/ddouble/dd-pi: ddouble elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool yy: ddouble then if yy: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then dd-pi4std/num/ddouble/dd-pi4: ddouble else ~std/num/ddouble/(~): (x : ddouble) -> ddoubledd-pi34std/num/ddouble/dd-pi34: ddouble elif xx: ddouble ==std/num/ddouble/(==): (x : ddouble, y : ddouble) -> bool ~std/num/ddouble/(~): (x : ddouble) -> ddoubleyy: ddouble then if yy: ddouble.is-posstd/num/ddouble/is-pos: (x : ddouble) -> bool then dd-pi34std/num/ddouble/dd-pi34: ddouble else ~std/num/ddouble/(~): (x : ddouble) -> ddoubledd-pi4std/num/ddouble/dd-pi4: ddouble else val rr: ddouble = sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble yy: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble ) val xrxr: ddouble = xx: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble rr: ddouble val yryr: ddouble = yy: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble rr: ddouble val zz: ddouble = atan2std/num/float64/atan2: (x : float64, y : float64) -> float64( yy: ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64, xx: ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 ).ddoublestd/num/ddouble/float64/ddouble: (d : float64) -> ddouble // approximation val (std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b)sinzsinz: ddouble,coszcosz: ddouble)std/core/types/Tuple2: forall<a,b> (fst : a, snd : b) -> (a, b) = sincosstd/num/ddouble/sincos: (rad : ddouble) -> (ddouble, ddouble)(zz: ddouble) if xrxr: ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64.absstd/num/float64/abs: (f : float64) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool yryr: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.histd/num/ddouble/ddouble/hi: (ddouble : ddouble) -> float64 then zz: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ((yryr: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble sinzsinz: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble coszcosz: ddouble) else zz: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble ((xrxr: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble coszcosz: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble sinzsinz: ddouble) // The hyperbolic sine of `x`. pub fun sinhstd/num/ddouble/sinh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble elif xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool 0.05literal: float64
hex64= 0x1.999999999999ap-5
then val exex: ddouble = xx: ddouble.expstd/num/ddouble/exp: (x : ddouble) -> ddouble if !std/core/types/bool/(!): (b : bool) -> boolexex: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then exex: ddouble else ( exex: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble exex: ddouble)).halfstd/num/ddouble/half: (x : ddouble) -> ddouble else // small x so use Taylor series to avoid cancellation val x2x2: ddouble = xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble xx: ddouble*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (x2x2: ddouble/std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble6literal: int
dec = 6
hex8 = 0x06
bit8 = 0b00000110
.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (x2x2: ddouble/std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble20literal: int
dec = 20
hex8 = 0x14
bit8 = 0b00010100
.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)*std/num/ddouble/(*): (x : ddouble, y : ddouble) -> ddouble(onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (x2x2: ddouble/std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble42literal: int
dec = 42
hex8 = 0x2A
bit8 = 0b00101010
.ddoublestd/num/ddouble/int/ddouble: (i : int) -> ddouble)))
) // The hyperbolic cosine of `x`. pub fun coshstd/num/ddouble/cosh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then onestd/num/ddouble/one: ddouble elif xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool 0.05literal: float64
hex64= 0x1.999999999999ap-5
then val exex: ddouble = xx: ddouble.expstd/num/ddouble/exp: (x : ddouble) -> ddouble if !std/core/types/bool/(!): (b : bool) -> boolexex: ddouble.is-finitestd/num/ddouble/is-finite: (x : ddouble) -> bool then exex: ddouble else (exex: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble exex: ddouble)).halfstd/num/ddouble/half: (x : ddouble) -> ddouble else val ss: ddouble = xx: ddouble.sinhstd/num/ddouble/sinh: (x : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble
) // The hyperbolic tangent of `x`. pub fun tanhstd/num/ddouble/tanh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then zerostd/num/ddouble/zero: ddouble elif xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble.float64std/num/ddouble/float64: (x : ddouble) -> float64 >std/num/float64/(>): (x : float64, y : float64) -> bool 0.05literal: float64
hex64= 0x1.999999999999ap-5
then val exex: ddouble = xx: ddouble.expstd/num/ddouble/exp: (x : ddouble) -> ddouble val iexiex: ddouble = onestd/num/ddouble/one: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble exex: ddouble if exex: ddouble.is-zerostd/num/ddouble/is-zero: (x : ddouble) -> bool then ~std/num/ddouble/(~): (x : ddouble) -> ddoubleonestd/num/ddouble/one: ddouble elif exex: ddouble.is-posinfstd/num/ddouble/is-posinf: (x : ddouble) -> bool then onestd/num/ddouble/one: ddouble else ((exex: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble iexiex: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble (exex: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble iexiex: ddouble)) else val ss: ddouble = xx: ddouble.sinhstd/num/ddouble/sinh: (x : ddouble) -> ddouble val cc: ddouble = sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble( onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble ss: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble ) returnreturn: ddouble (ss: ddouble /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble cc: ddouble
) // The area hyperbolic sine of `x`. pub fun asinhstd/num/ddouble/asinh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V lnstd/num/ddouble/ln: (x : ddouble) -> ddouble( xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble)) // The area hyperbolic cosine of `x`. pub fun acoshstd/num/ddouble/acosh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble <std/num/ddouble/(<): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble else lnstd/num/ddouble/ln: (x : ddouble) -> ddouble(xx: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble sqrtstd/num/ddouble/sqrt: (x : ddouble) -> ddouble(xx: ddouble.sqrstd/num/ddouble/sqr: (x : ddouble) -> ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble onestd/num/ddouble/one: ddouble)) // The area hyperbolic tangent of `x`. pub fun atanhstd/num/ddouble/atanh: (x : ddouble) -> ddouble( xx: ddouble : ddoublestd/num/ddouble/ddouble: V )result: -> total ddouble : ddoublestd/num/ddouble/ddouble: V if xx: ddouble.absstd/num/ddouble/abs: (x : ddouble) -> ddouble >std/num/ddouble/(>): (x : ddouble, y : ddouble) -> bool onestd/num/ddouble/one: ddouble then dd-nanstd/num/ddouble/dd-nan: ddouble else lnstd/num/ddouble/ln: (x : ddouble) -> ddouble((onestd/num/ddouble/one: ddouble +std/num/ddouble/(+): (x : ddouble, y : ddouble) -> ddouble xx: ddouble) /std/num/ddouble/(/): (x : ddouble, y : ddouble) -> ddouble (onestd/num/ddouble/one: ddouble -std/num/ddouble/(-): (x : ddouble, y : ddouble) -> ddouble xx: ddouble)).halfstd/num/ddouble/half: (x : ddouble) -> ddouble