|
| 1 | +unit scimath |
| 2 | + |
| 3 | +use { |
| 4 | + math.Math, |
| 5 | + scimath.SciMath |
| 6 | +} |
| 7 | + |
| 8 | +share Integration { |
| 9 | + share trapezoid(values: [], dx: int|float) :: int|float { |
| 10 | + if values.size < 2 { ~> 0 } |
| 11 | + sum: int|float = 0 |
| 12 | + for i of 0 to values.size - 2 { |
| 13 | + sum = sum + (values[i] + values[i + 1]) * 0.5 * dx |
| 14 | + } |
| 15 | + ~> sum |
| 16 | + } |
| 17 | +} |
| 18 | + |
| 19 | +share Interpolation { |
| 20 | + share linear(x0: int|float, y0: int|float, x1: int|float, y1: int|float, x: int|float) :: int|float { |
| 21 | + if x1 == x0 { ~> y0 } |
| 22 | + t := (x - x0) / (x1 - x0) |
| 23 | + ~> y0 + t * (y1 - y0) |
| 24 | + } |
| 25 | +} |
| 26 | + |
| 27 | +share Convolution { |
| 28 | + share discrete(a: [], b: []) :: [] { |
| 29 | + if a.size == 0 { ~> [] } |
| 30 | + if b.size == 0 { ~> [] } |
| 31 | + outSize := a.size + b.size - 1 |
| 32 | + outArr := [0 to outSize - 1] |
| 33 | + for i of 0 to outSize - 1 { |
| 34 | + outArr[i] = 0 |
| 35 | + } |
| 36 | + for i of 0 to a.size - 1 { |
| 37 | + for j of 0 to b.size - 1 { |
| 38 | + outArr[i + j] = outArr[i + j] + a[i] * b[j] |
| 39 | + } |
| 40 | + } |
| 41 | + ~> outArr |
| 42 | + } |
| 43 | +} |
| 44 | + |
| 45 | +share Spectral { |
| 46 | + share fft(values: []) :: [] { ~> SciMath.fft(values) } |
| 47 | + share ifft(values: []) :: [] { ~> SciMath.ifft(values) } |
| 48 | +} |
| 49 | + |
| 50 | +share Window { |
| 51 | + share hann(n: int) :: [] { |
| 52 | + if n <= 0 { ~> [] } |
| 53 | + if n == 1 { ~> [1] } |
| 54 | + outArr := [0 to n - 1] |
| 55 | + for i of 0 to n - 1 { |
| 56 | + outArr[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (n - 1)) |
| 57 | + } |
| 58 | + ~> outArr |
| 59 | + } |
| 60 | +} |
| 61 | + |
| 62 | +share ODE { |
| 63 | + share rk45(ode: [], y0: int|float, t0: int|float, t1: int|float, h: int|float) :: [] { |
| 64 | + ~> SciMath.rk45(ode, y0, t0, t1, h) |
| 65 | + } |
| 66 | + |
| 67 | + share dop853(ode: [], y0: int|float, t0: int|float, t1: int|float, h: int|float) :: [] { |
| 68 | + ~> SciMath.dop853(ode, y0, t0, t1, h) |
| 69 | + } |
| 70 | + |
| 71 | + share odeint(ode: [], y0: int|float, t0: int|float, t1: int|float, h: int|float) :: [] { |
| 72 | + ~> SciMath.odeint(ode, y0, t0, t1, h) |
| 73 | + } |
| 74 | +} |
| 75 | + |
| 76 | +share ErrorFunction { |
| 77 | + share erf(x: int|float) :: int|float { |
| 78 | + sign: int|float = 1 |
| 79 | + xx: int|float = x |
| 80 | + if x < 0 { |
| 81 | + sign = -1 |
| 82 | + xx = -x |
| 83 | + } |
| 84 | + a1: float = 0.254829592 |
| 85 | + a2: float = -0.284496736 |
| 86 | + a3: float = 1.421413741 |
| 87 | + a4: float = -1.453152027 |
| 88 | + a5: float = 1.061405429 |
| 89 | + p: float = 0.3275911 |
| 90 | + t := 1 / (1 + p * xx) |
| 91 | + y := 1 - (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t) * Math.exp(-xx * xx) |
| 92 | + ~> sign * y |
| 93 | + } |
| 94 | +} |
| 95 | + |
| 96 | +share Bessel { |
| 97 | + share j0(x: int|float) :: int|float { |
| 98 | + xx := x * x |
| 99 | + ~> 1 - xx / 4 + (xx * xx) / 64 |
| 100 | + } |
| 101 | +} |
| 102 | + |
| 103 | +share Airy { |
| 104 | + share ai(x: int|float) :: int|float { |
| 105 | + ~> Math.exp(-2 * Math.abs(x) / 3) |
| 106 | + } |
| 107 | +} |
| 108 | + |
| 109 | +share Elliptic { |
| 110 | + share k(m: int|float) :: int|float { |
| 111 | + if m < 0 { ~> 0 } |
| 112 | + if m >= 1 { ~> 0 } |
| 113 | + ~> (Math.PI / 2) * (1 + m / 4 + 9 * m * m / 64) |
| 114 | + } |
| 115 | +} |
| 116 | + |
| 117 | +share Bernoulli { |
| 118 | + share pmf(k: int, p: int|float) :: int|float { |
| 119 | + if k == 1 { ~> p } |
| 120 | + if k == 0 { ~> 1 - p } |
| 121 | + ~> 0 |
| 122 | + } |
| 123 | +} |
| 124 | + |
| 125 | +share Binomial { |
| 126 | + share pmf(k: int, n: int, p: int|float) :: int|float { |
| 127 | + if k < 0 { ~> 0 } |
| 128 | + if k > n { ~> 0 } |
| 129 | + c := Math.binomial(n, k) |
| 130 | + ~> c * Math.pow(p, k) * Math.pow(1 - p, n - k) |
| 131 | + } |
| 132 | +} |
| 133 | + |
| 134 | +share Geometric { |
| 135 | + share pmf(k: int, p: int|float) :: int|float { |
| 136 | + if k < 1 { ~> 0 } |
| 137 | + ~> Math.pow(1 - p, k - 1) * p |
| 138 | + } |
| 139 | +} |
| 140 | + |
| 141 | +share LogNormal { |
| 142 | + share pdf(x: int|float, mu: int|float, sigma: int|float) :: int|float { |
| 143 | + if x <= 0 { ~> 0 } |
| 144 | + if sigma <= 0 { ~> 0 } |
| 145 | + z := (Math.log(x) - mu) / sigma |
| 146 | + den := x * sigma * Math.sqrt(2 * Math.PI) |
| 147 | + ~> Math.exp(-0.5 * z * z) / den |
| 148 | + } |
| 149 | +} |
| 150 | + |
| 151 | +share Weibull { |
| 152 | + share pdf(x: int|float, k: int|float, lambda: int|float) :: int|float { |
| 153 | + if x < 0 { ~> 0 } |
| 154 | + if k <= 0 { ~> 0 } |
| 155 | + if lambda <= 0 { ~> 0 } |
| 156 | + xx := x / lambda |
| 157 | + if xx <= 0 { ~> 0 } |
| 158 | + ~> (k / lambda) * Math.exp((k - 1) * Math.log(xx)) * Math.exp(-Math.exp(k * Math.log(xx))) |
| 159 | + } |
| 160 | +} |
| 161 | + |
| 162 | +share Pareto { |
| 163 | + share pdf(x: int|float, xm: int|float, alpha: int|float) :: int|float { |
| 164 | + if x < xm { ~> 0 } |
| 165 | + if xm <= 0 { ~> 0 } |
| 166 | + if alpha <= 0 { ~> 0 } |
| 167 | + ~> alpha * Math.exp(alpha * Math.log(xm)) / Math.exp((alpha + 1) * Math.log(x)) |
| 168 | + } |
| 169 | +} |
0 commit comments