Sources/OrgCore/Compute/Calc.swift
600 lines · 24011 bytes
1import Foundation
2
3/// The part of Emacs Calc that table formulas use, with Calc's number semantics: exact
4/// integers, decimal floats rounded to `calc-internal-prec` digits after every operation, and
5/// Calc's display formats. Anything outside this domain throws `Calc.Unsupported`, and the
6/// caller hands the work to Emacs.
7public enum Calc {
8 public struct Unsupported: Error, Equatable, CustomStringConvertible {
9 public let reason: String
10 public var description: String { reason }
11 init(_ reason: String) { self.reason = reason }
12 }
13
14 public enum Value: Equatable, Sendable {
15 case int(Int)
16 /// Mantissa and exponent, normalized: at most `precision` digits, no trailing zeros.
17 case float(Int, Int)
18 case vector([Value])
19 }
20
21 public enum FloatFormat: Equatable, Sendable {
22 case float(Int), fix(Int), sci(Int), eng(Int)
23 }
24
25 public struct Modes: Equatable, Sendable {
26 public var precision = 12
27 public var format = FloatFormat.float(8)
28 public var degrees = true
29 public init() {}
30 }
31
32 /// `calc-eval` of `expression`, as the string Calc would return.
33 public static func evaluate(_ expression: String, modes: Modes = Modes()) throws -> String {
34 guard modes.precision == 12 else { throw Unsupported("calc precision") }
35 var parser = Parser(expression)
36 let tree = try parser.parse()
37 let value = try Evaluator(modes: modes).eval(tree)
38 return try format(value, modes.format)
39 }
40
41 // MARK: - Numbers
42
43 static let precision = 12
44
45 static func digits(_ x: Int128) -> Int {
46 var n = x.magnitude
47 var count = 1
48 while n >= 10 {
49 n /= 10
50 count += 1
51 }
52 return count
53 }
54
55 static func power10(_ n: Int) throws -> Int128 {
56 guard n >= 0, n <= 37 else { throw Unsupported("number size") }
57 var result: Int128 = 1
58 for _ in 0..<n { result *= 10 }
59 return result
60 }
61
62 /// `math-scale-rounding`: `x` times 10^n, rounding half away from zero when n < 0.
63 static func scaleRounding(_ x: Int128, _ n: Int) throws -> Int128 {
64 if n >= 0 {
65 let (result, overflow) = x.multipliedReportingOverflow(by: try power10(n))
66 if overflow { throw Unsupported("number size") }
67 return result
68 }
69 if x < 0 { return -(try scaleRounding(-x, n)) }
70 if -n > 38 { return 0 }
71 let truncated = -n - 1 > 37 ? 0 : x / (try power10(-n - 1))
72 return (truncated + 5) / 10
73 }
74
75 /// `math-make-float`.
76 static func makeFloat(_ mantissa: Int128, _ exponent: Int) throws -> Value {
77 var m = mantissa
78 var e = exponent
79 if m == 0 { return .float(0, 0) }
80 let excess = digits(m) - precision
81 if excess > 0 {
82 m = try scaleRounding(m, -excess)
83 e += excess
84 }
85 if m == 0 { return .float(0, 0) }
86 while m % 10 == 0 {
87 m /= 10
88 e += 1
89 }
90 return .float(Int(m), e)
91 }
92
93 static func parts(_ v: Value) throws -> (Int128, Int) {
94 switch v {
95 case .int(let i): return (Int128(i), 0)
96 case .float(let m, let e): return (Int128(m), e)
97 case .vector: throw Unsupported("vector arithmetic")
98 }
99 }
100
101 static func double(_ v: Value) throws -> Double {
102 switch v {
103 case .int(let i): return Double(i)
104 case .float(let m, let e): return Double(m) * pow(10, Double(e))
105 case .vector: throw Unsupported("vector arithmetic")
106 }
107 }
108
109 static func fromDouble(_ d: Double) throws -> Value {
110 guard d.isFinite else { throw Unsupported("infinite result") }
111 if d == 0 { return .float(0, 0) }
112 // Seventeen significant digits are exact for a Double; rounding to 12 follows.
113 let text = String(format: "%.16e", d)
114 let pieces = text.split(separator: "e")
115 let digitsText = pieces[0].replacingOccurrences(of: ".", with: "")
116 guard let mantissa = Int128(digitsText), let exponent = Int(pieces[1]) else { throw Unsupported("number") }
117 return try makeFloat(mantissa, exponent - 16)
118 }
119
120 static func add(_ a: Value, _ b: Value) throws -> Value {
121 if case .int(let x) = a, case .int(let y) = b {
122 let (sum, overflow) = x.addingReportingOverflow(y)
123 if overflow { throw Unsupported("integer size") }
124 return .int(sum)
125 }
126 let (m1, e1) = try parts(a)
127 let (m2, e2) = try parts(b)
128 if m1 == 0 { return try makeFloat(m2, e2) }
129 if m2 == 0 { return try makeFloat(m1, e1) }
130 let ediff = e1 - e2
131 if ediff >= 0 {
132 if ediff >= 2 * precision { return try makeFloat(m1, e1) }
133 return try makeFloat(try scaleRounding(m1, ediff) + m2, e2)
134 }
135 if -ediff >= 2 * precision { return try makeFloat(m2, e2) }
136 return try makeFloat(m1 + (try scaleRounding(m2, -ediff)), e1)
137 }
138
139 static func negate(_ a: Value) throws -> Value {
140 switch a {
141 case .int(let i):
142 guard i != .min else { throw Unsupported("integer size") }
143 return .int(-i)
144 case .float(let m, let e): return .float(-m, e)
145 case .vector: throw Unsupported("vector arithmetic")
146 }
147 }
148
149 static func multiply(_ a: Value, _ b: Value) throws -> Value {
150 if case .int(let x) = a, case .int(let y) = b {
151 let (product, overflow) = x.multipliedReportingOverflow(by: y)
152 if overflow { throw Unsupported("integer size") }
153 return .int(product)
154 }
155 let (m1, e1) = try parts(a)
156 let (m2, e2) = try parts(b)
157 let (product, overflow) = m1.multipliedReportingOverflow(by: m2)
158 if overflow { throw Unsupported("number size") }
159 return try makeFloat(product, e1 + e2)
160 }
161
162 static func divide(_ a: Value, _ b: Value) throws -> Value {
163 if case .int(let x) = a, case .int(let y) = b {
164 guard y != 0 else { throw Unsupported("division by zero") }
165 if x % y == 0 { return .int(x / y) }
166 }
167 let (m1, e1) = try parts(a)
168 let (m2, e2) = try parts(b)
169 guard m2 != 0 else { throw Unsupported("division by zero") }
170 // `math-div-float`: scale so the truncated quotient has one digit more than needed.
171 let ldiff = max(precision + 1 - (digits(m1) - digits(m2)), 0)
172 let quotient = (try scaleRounding(m1, ldiff)) / m2
173 return try makeFloat(quotient, e1 - ldiff - e2)
174 }
175
176 static func isInteger(_ v: Value) -> Bool {
177 if case .int = v { return true }
178 return false
179 }
180
181 /// `math-ipow`.
182 static func integerPower(_ a: Value, _ n: Int) throws -> Value {
183 if n < 0 { return try integerPower(try divide(.int(1), a), -n) }
184 if n == 0 { return .int(1) }
185 if n == 1 { return a }
186 let square = try multiply(a, a)
187 if n % 2 == 0 { return try integerPower(square, n / 2) }
188 return try multiply(a, try integerPower(square, n / 2))
189 }
190
191 static func power(_ a: Value, _ b: Value) throws -> Value {
192 if case .int(let n) = b { return try integerPower(a, n) }
193 let base = try double(a)
194 guard base >= 0 else { throw Unsupported("complex result") }
195 return try fromDouble(pow(base, try double(b)))
196 }
197
198 static func compare(_ a: Value, _ b: Value) throws -> Int {
199 let difference = try add(a, try negate(b))
200 let (m, _) = try parts(difference)
201 return m == 0 ? 0 : (m < 0 ? -1 : 1)
202 }
203
204 // MARK: - Display
205
206 /// `math-format-number`.
207 static func format(_ value: Value, _ floatFormat: FloatFormat) throws -> String {
208 switch value {
209 case .int(let i): return String(i)
210 case .vector: throw Unsupported("vector result")
211 case .float(let m, _) where m < 0:
212 guard case .float(let m, let e) = value else { fatalError() }
213 return "-" + (try format(.float(-m, e), floatFormat))
214 case .float(let m, var exp):
215 var mant = Int128(m)
216 var figs: Int
217 let kind: String
218 switch floatFormat {
219 case .float(let n): (figs, kind) = (n, "float")
220 case .fix(let n): (figs, kind) = (n, "fix")
221 case .sci(let n): (figs, kind) = (n, "sci")
222 case .eng(let n): (figs, kind) = (n, "eng")
223 }
224 if kind == "fix", figs < 0 || exp + digits(mant) > -figs {
225 if figs < 0 { figs = -figs }
226 mant = try scaleRounding(mant, exp + figs)
227 var str = String(mant)
228 if str.count <= figs { str = String(repeating: "0", count: figs + 1 - str.count) + str }
229 if figs > 0 {
230 return String(str.dropLast(figs)) + "." + String(str.suffix(figs))
231 }
232 return str + "."
233 }
234 if figs < 0 { figs += precision }
235 if figs > 0 {
236 let adj = figs - digits(mant)
237 if adj < 0 {
238 mant = try scaleRounding(mant, adj)
239 exp -= adj
240 }
241 }
242 var str = String(mant)
243 let len = str.count
244 let dpos = exp + len
245 // `calc-display-sci-high` 0, `calc-display-sci-low` -3.
246 if kind == "float", dpos <= precision, dpos >= -1 {
247 if dpos == 0 {
248 str = "0." + str
249 } else if exp <= 0, dpos > 0 {
250 str = String(str.prefix(dpos)) + "." + String(str.dropFirst(dpos))
251 } else if exp > 0 {
252 str += String(repeating: "0", count: exp) + "."
253 } else {
254 str = "0." + String(repeating: "0", count: -dpos) + str
255 }
256 return str
257 }
258 let eadj = exp + len
259 let scale = kind == "eng" ? 1 + ((eadj + 300002) % 3) : 1
260 if scale > str.count { str += String(repeating: "0", count: scale - str.count) }
261 if scale < str.count { str = String(str.prefix(scale)) + "." + String(str.dropFirst(scale)) }
262 return str + "e" + String(eadj - scale)
263 }
264 }
265
266 // MARK: - Parsing
267
268 indirect enum Node {
269 case number(Value)
270 case vector([Node])
271 case unary(Character, Node)
272 case binary(Character, Node, Node)
273 case call(String, [Node])
274 }
275
276 struct Parser {
277 let chars: [Character]
278 var i = 0
279
280 init(_ s: String) { chars = Array(s) }
281
282 mutating func skip() {
283 while i < chars.count, chars[i] == " " || chars[i] == "\t" { i += 1 }
284 }
285
286 mutating func peek() -> Character? {
287 skip()
288 return i < chars.count ? chars[i] : nil
289 }
290
291 mutating func parse() throws -> Node {
292 let node = try sum()
293 guard peek() == nil else { throw Unsupported("formula syntax") }
294 return node
295 }
296
297 // `+ -` < `/` < `*` < unary minus < `^`, as Calc's operator table orders them.
298 mutating func sum() throws -> Node {
299 var left = try quotient()
300 while let c = peek(), c == "+" || c == "-" {
301 i += 1
302 left = .binary(c, left, try quotient())
303 }
304 return left
305 }
306
307 mutating func quotient() throws -> Node {
308 var left = try product()
309 while peek() == "/" {
310 i += 1
311 left = .binary("/", left, try product())
312 }
313 return left
314 }
315
316 mutating func product() throws -> Node {
317 var left = try unary()
318 while peek() == "*" {
319 i += 1
320 left = .binary("*", left, try unary())
321 }
322 return left
323 }
324
325 mutating func unary() throws -> Node {
326 if let c = peek(), c == "-" || c == "+" {
327 i += 1
328 let operand = try unary()
329 return c == "-" ? .unary("-", operand) : operand
330 }
331 return try powerNode()
332 }
333
334 mutating func powerNode() throws -> Node {
335 let base = try primary()
336 if peek() == "^" {
337 i += 1
338 return .binary("^", base, try unary())
339 }
340 return base
341 }
342
343 mutating func primary() throws -> Node {
344 guard let c = peek() else { throw Unsupported("formula syntax") }
345 if c == "(" {
346 i += 1
347 let inner = try sum()
348 guard peek() == ")" else { throw Unsupported("formula syntax") }
349 i += 1
350 return inner
351 }
352 if c == "[" {
353 i += 1
354 var items: [Node] = []
355 if peek() != "]" {
356 items.append(try sum())
357 while peek() == "," {
358 i += 1
359 items.append(try sum())
360 }
361 }
362 guard peek() == "]" else { throw Unsupported("formula syntax") }
363 i += 1
364 return .vector(items)
365 }
366 if c.isASCII, c.isNumber || c == "." { return .number(try number()) }
367 if c.isLetter {
368 var name = ""
369 while i < chars.count, chars[i].isLetter || chars[i].isNumber || chars[i] == "_" {
370 name.append(chars[i])
371 i += 1
372 }
373 guard peek() == "(" else { throw Unsupported("variable \(name)") }
374 i += 1
375 var args: [Node] = []
376 if peek() != ")" {
377 args.append(try sum())
378 while peek() == "," {
379 i += 1
380 args.append(try sum())
381 }
382 }
383 guard peek() == ")" else { throw Unsupported("formula syntax") }
384 i += 1
385 return .call(name, args)
386 }
387 throw Unsupported("formula syntax")
388 }
389
390 /// Calc's number syntax: `7`, `007`, `1.5`, `.5`, `5.`, `1e3`, `1.5e-3`.
391 mutating func number() throws -> Value {
392 var whole = ""
393 while i < chars.count, chars[i].isASCII, chars[i].isNumber {
394 whole.append(chars[i])
395 i += 1
396 }
397 var fraction: String?
398 if i < chars.count, chars[i] == "." {
399 i += 1
400 var f = ""
401 while i < chars.count, chars[i].isASCII, chars[i].isNumber {
402 f.append(chars[i])
403 i += 1
404 }
405 fraction = f
406 }
407 var exponent: Int?
408 if i < chars.count, chars[i] == "e" || chars[i] == "E" {
409 var j = i + 1
410 var e = ""
411 if j < chars.count, chars[j] == "-" || chars[j] == "+" {
412 e.append(chars[j])
413 j += 1
414 }
415 var digitsSeen = false
416 while j < chars.count, chars[j].isASCII, chars[j].isNumber {
417 e.append(chars[j])
418 j += 1
419 digitsSeen = true
420 }
421 if digitsSeen {
422 exponent = Int(e)
423 i = j
424 }
425 }
426 guard !whole.isEmpty || !(fraction ?? "").isEmpty else { throw Unsupported("formula syntax") }
427 if fraction == nil, exponent == nil {
428 guard let value = Int(whole) else { throw Unsupported("integer size") }
429 return .int(value)
430 }
431 let digitsText = whole + (fraction ?? "")
432 guard let mantissa = Int128(digitsText.isEmpty ? "0" : digitsText) else { throw Unsupported("number size") }
433 return try makeFloat(mantissa, (exponent ?? 0) - (fraction?.count ?? 0))
434 }
435 }
436
437 // MARK: - Evaluation
438
439 struct Evaluator {
440 let modes: Modes
441
442 func eval(_ node: Node) throws -> Value {
443 switch node {
444 case .number(let v): return v
445 case .vector(let items): return .vector(try items.map(eval))
446 case .unary(_, let operand): return try negate(try eval(operand))
447 case .binary(let op, let l, let r):
448 let a = try eval(l)
449 let b = try eval(r)
450 switch op {
451 case "+": return try add(a, b)
452 case "-": return try add(a, try negate(b))
453 case "*": return try multiply(a, b)
454 case "/": return try divide(a, b)
455 default: return try power(a, b)
456 }
457 case .call(let name, let args):
458 return try call(name, try args.map(eval))
459 }
460 }
461
462 func elements(_ args: [Value]) throws -> [Value] {
463 guard args.count == 1 else { throw Unsupported("arguments") }
464 if case .vector(let items) = args[0] { return items }
465 return [args[0]]
466 }
467
468 func one(_ args: [Value]) throws -> Value {
469 guard args.count == 1 else { throw Unsupported("arguments") }
470 if case .vector = args[0] { throw Unsupported("vector argument") }
471 return args[0]
472 }
473
474 func radians(_ x: Double) -> Double { modes.degrees ? x * .pi / 180 : x }
475 func angle(_ x: Double) -> Double { modes.degrees ? x * 180 / .pi : x }
476
477 func call(_ name: String, _ args: [Value]) throws -> Value {
478 switch name {
479 case "vsum":
480 return try elements(args).reduce(Value.int(0)) { try add($0, $1) }
481 case "vprod":
482 return try elements(args).reduce(Value.int(1)) { try multiply($0, $1) }
483 case "vcount":
484 return .int(try elements(args).count)
485 case "vmean":
486 let items = try elements(args)
487 guard !items.isEmpty else { throw Unsupported("empty vector") }
488 return try divide(try items.reduce(Value.int(0)) { try add($0, $1) }, .int(items.count))
489 case "vmax", "vmin":
490 let items = try elements(args)
491 guard var best = items.first else { throw Unsupported("empty vector") }
492 for item in items.dropFirst() {
493 let order = try compare(item, best)
494 if name == "vmax" ? order > 0 : order < 0 { best = item }
495 }
496 return best
497 case "vmedian":
498 var items = try elements(args)
499 guard !items.isEmpty else { throw Unsupported("empty vector") }
500 try items.sort { try compare($0, $1) < 0 }
501 if items.count % 2 == 1 { return items[items.count / 2] }
502 return try divide(try add(items[items.count / 2 - 1], items[items.count / 2]), .int(2))
503 case "max", "min":
504 guard !args.isEmpty else { throw Unsupported("arguments") }
505 var best = args[0]
506 for item in args.dropFirst() {
507 let order = try compare(item, best)
508 if name == "max" ? order > 0 : order < 0 { best = item }
509 }
510 return best
511 case "abs":
512 let x = try one(args)
513 return try compare(x, .int(0)) < 0 ? try negate(x) : x
514 case "floor", "ceil", "trunc":
515 let x = try one(args)
516 if case .int = x { return x }
517 let d = try double(x)
518 let r = name == "floor" ? d.rounded(.down) : name == "ceil" ? d.rounded(.up) : d.rounded(.towardZero)
519 guard let i = Int(exactly: r) else { throw Unsupported("integer size") }
520 return .int(i)
521 case "round":
522 guard args.count == 1 || args.count == 2 else { throw Unsupported("arguments") }
523 let x = args[0]
524 if args.count == 2 {
525 guard case .int(let places) = args[1] else { throw Unsupported("arguments") }
526 guard case .float(let m, let e) = x else { return x }
527 return try makeFloat(try scaleRounding(Int128(m), e + places), -places)
528 }
529 guard case .float(let m, let e) = x else { return x }
530 guard let i = Int(exactly: try scaleRounding(Int128(m), e)) else { throw Unsupported("integer size") }
531 return .int(i)
532 case "sqrt":
533 let x = try one(args)
534 if case .int(let i) = x {
535 guard i >= 0 else { throw Unsupported("complex result") }
536 let root = Int(Double(i).squareRoot().rounded())
537 if root * root == i { return .int(root) }
538 }
539 let d = try double(x)
540 guard d >= 0 else { throw Unsupported("complex result") }
541 return try fromDouble(d.squareRoot())
542 case "exp":
543 let x = try one(args)
544 if x == .int(0) { return .int(1) }
545 return try fromDouble(Foundation.exp(try double(x)))
546 case "ln":
547 let x = try one(args)
548 if x == .int(1) { return .int(0) }
549 let d = try double(x)
550 guard d > 0 else { throw Unsupported("logarithm domain") }
551 return try fromDouble(Foundation.log(d))
552 case "log10":
553 let x = try one(args)
554 if case .int(let i) = x, i > 0 {
555 var n = i
556 var k = 0
557 while n % 10 == 0 {
558 n /= 10
559 k += 1
560 }
561 if n == 1 { return .int(k) }
562 }
563 let d = try double(x)
564 guard d > 0 else { throw Unsupported("logarithm domain") }
565 return try fromDouble(Foundation.log10(d))
566 case "sin", "cos", "tan":
567 let x = try one(args)
568 if modes.degrees, case .int(let i) = x, i % 90 == 0 {
569 let quarter = ((i / 90) % 4 + 4) % 4
570 switch name {
571 case "sin": return .int([0, 1, 0, -1][quarter])
572 case "cos": return .int([1, 0, -1, 0][quarter])
573 default:
574 guard quarter % 2 == 0 else { throw Unsupported("infinite result") }
575 return .int(0)
576 }
577 }
578 if !modes.degrees, x == .int(0) { return .int(name == "cos" ? 1 : 0) }
579 let r = radians(try double(x))
580 return try fromDouble(name == "sin" ? Foundation.sin(r) : name == "cos" ? Foundation.cos(r) : Foundation.tan(r))
581 case "arcsin", "arccos", "arctan":
582 let x = try one(args)
583 if modes.degrees, case .int(let i) = x {
584 switch (name, i) {
585 case ("arcsin", -1), ("arcsin", 0), ("arcsin", 1): return .int(90 * i)
586 case ("arccos", -1), ("arccos", 0), ("arccos", 1): return .int(90 - 90 * i)
587 case ("arctan", -1), ("arctan", 0), ("arctan", 1): return .int(45 * i)
588 default: break
589 }
590 }
591 let d = try double(x)
592 if name != "arctan", abs(d) > 1 { throw Unsupported("complex result") }
593 let r = name == "arcsin" ? Foundation.asin(d) : name == "arccos" ? Foundation.acos(d) : Foundation.atan(d)
594 return try fromDouble(angle(r))
595 default:
596 throw Unsupported("function \(name)")
597 }
598 }
599 }
600}