Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

No max operation in formulas #797

Merged
merged 15 commits into from
May 14, 2024
Merged

No max operation in formulas #797

Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
0349dee
no max
AYadrov Apr 23, 2024
3265969
missed max
AYadrov Apr 23, 2024
d71ee37
set lower boundary of precisions for each operation
AYadrov Apr 23, 2024
69d3dcc
oops 20->5 ampl bits
AYadrov Apr 23, 2024
717d730
few attempts to find bugs
AYadrov Apr 23, 2024
b198eae
typo
AYadrov Apr 23, 2024
fa5edae
Merge branch 'main' into artem-no-max-formulas
AYadrov Apr 24, 2024
12718d1
merged main
AYadrov Apr 29, 2024
47c5656
atan2 fomrula change, max is the right operation
AYadrov May 1, 2024
31bbe25
pow fix + super smart slack
AYadrov May 3, 2024
712b2f9
trash
AYadrov May 3, 2024
6d9a683
trying to fix overflows at ival-pow
AYadrov May 6, 2024
ed2f238
rolling back to the previous commit
AYadrov May 7, 2024
4c1f625
comments typos
AYadrov May 7, 2024
4292dd7
final version before merging
AYadrov May 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
89 changes: 52 additions & 37 deletions src/correct-round.rkt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@
[prec-old (in-vector (if (equal? 1 current-iter) vstart-precs vprecs))]
[prec-new (in-vector vprecs-new)]
[n (in-naturals)])
(vector-set! vrepeats n (and (equal? prec-new prec-old)
(vector-set! vrepeats n (and (<= prec-new prec-old)
(andmap (lambda (x) (or (< x varc) (vector-ref vrepeats (- x varc))))
(rest instr)))))

Expand All @@ -190,17 +190,17 @@
(define exps-from-above (vector-ref vprecs-new (- n varc))) ; vprecs-new is shifted by varc elements from vregs
(define new-exponents (get-exponent op output srcs))

(define final-parent-precision (min (*max-mpfr-prec*)
(+ exps-from-above
(vector-ref vstart-precs (- n varc)))))
(define final-parent-precision (max (+ exps-from-above
(vector-ref vstart-precs (- n varc)))
(*base-tuning-precision*)))

; This case is weird. if we have a cancellation in fma -> ival-mult in fma should be in higher precision
(when (equal? op ival-fma)
(set! final-parent-precision (+ final-parent-precision (first new-exponents))))

(when (equal? final-parent-precision (*max-mpfr-prec*)) ; Early stopping
(when (>= final-parent-precision (*max-mpfr-prec*)) ; Early stopping
(*sampling-iteration* (*max-sampling-iterations*)))
(vector-set! vprecs-new (- n varc) final-parent-precision)
(vector-set! vprecs-new (- n varc) (min final-parent-precision (*max-mpfr-prec*)))

(for ([x (in-list tail-registers)]
[new-exp (in-list new-exponents)]
Expand All @@ -222,7 +222,7 @@
; log[Г*]'x = log[Г*]'y = log[Г/]'x = log[Г/]'y = 1
; log[Гsqrt] = 0.5
; log[Гcbrt] = 0.3
(make-list (length srcs) 1)]
(make-list (length srcs) 0)] ; assume that *ampl-bits* already introduces this 1 bit

[(ival-add ival-sub)
; log[Г+]'x = log[x] - log[x + y] + 1 (1 for mantissa approximation when dividing)
Expand All @@ -244,31 +244,46 @@
0
(get-slack)))

(list (+ (max 0 (+ (- x-exp out-exp) slack)) 1) ; exponent per x
(+ (max 0 (+ (- y-exp out-exp) slack)) 1))] ; exponent per y
(list (+ (- x-exp out-exp) slack 1) ; exponent per x
(+ (- y-exp out-exp) slack 1))] ; exponent per y

[(ival-pow)
; log[Гpow]'x = log[y]
; log[Гpow]'y = log[y] + log[log(x)]
; ^^^^^^^^^^^
; always less than 30
; log[Гpow]'y = log[y] + log[log(x)] <= log[y] + |log[x]|

(define x (first srcs))
(define y (second srcs))
(define outlo (ival-lo output))
(define outhi (ival-hi output))

(define x-exp (ival-max-log2-approx x))
(define y-exp (ival-max-log2-approx y))
(define slack (if (>= x-exp 3) ; if x > 2 (actually 2.718),
30 ; then at least 1 additional bit is needed
0))

(list (max 0 y-exp) ; exponent per x
(max 0 (+ y-exp slack)))] ; exponent per y

; when output crosses zero and x is negative - means that y was fractional and not fixed
; solution - add more slack for y to converge
(define y-slack (if (and (not (equal? (mpfr-sign outlo) (mpfr-sign outhi)))
(bfnegative? (ival-lo x))
(bfnegative? (ival-hi x)))
(get-slack)
0))

(define xlo-exp (mpfr-exp (ival-lo x)))
(define xhi-exp (mpfr-exp (ival-hi x)))
(define x-slack ; it is likely to be a handling case from IEEE, x is not close enough
(if (or (equal? xlo-exp -9223372036854775807)
(equal? xhi-exp -9223372036854775807) ; if interval contains 0
(and (> 2 xlo-exp) (<= 2 xhi-exp))) ; if interval crosses 1
(get-slack)
0))

(list (+ y-exp x-slack) ; exponent per x
(+ y-exp (abs x-exp) y-slack))] ; exponent per y

[(ival-exp)
; log[Гexp] = log[x]
(define x (car srcs))
(define x-exp (ival-max-log2-approx x))
(list (max 0 x-exp))]
(list x-exp)]

[(ival-tan)
; log[Гtan] = log[x] - log[sin(x)*cos(x)] <= log[x] + |log[tan(x)]| + 1
Expand All @@ -288,7 +303,7 @@
(get-slack) ; tan is (-inf, +inf) or around zero (but x != 0)
0))

(list (max 0 (+ x-exp out-exp-abs 1 slack)))]
(list (+ x-exp out-exp-abs 1 slack))]

[(ival-sin ival-cos ival-sinh ival-cosh)
; log[Гcos] = log[x] + log[sin(x)] - log[cos(x)] <= log[x] - log[cos(x)] + 1
Expand All @@ -308,7 +323,7 @@
(get-slack) ; Condition of uncertainty
0))

(list (+ (max 0 (+ (- x-exp out-exp) slack)) 1))]
(list (+ (- x-exp out-exp) slack 1))]

[(ival-log ival-log2 ival-log10)
; log[Гlog] = log[1/logx] = -log[log(x)]
Expand All @@ -324,7 +339,7 @@
0
(get-slack))) ; output crosses 0 - uncertainty

(list (max 0 (+ out-exp 1 slack)))]
(list (+ out-exp 1 slack))]

[(ival-asin ival-acos)
; log[Гasin] = log[x] - log[1-x^2]/2 - log[asin(x)] + 1
Expand All @@ -339,17 +354,17 @@
(get-slack) ; assumes that log[1-x^2]/2 is equal to slack
0))

(list (max 0 (+ (- x-exp out-exp) 1 slack)))]
(list (+ (- x-exp out-exp) 1 slack))]

[(ival-atan)
; log[Гatan] = log[x] - log[x^2+1] - log[atan(x)] <= -|log[x]| - log[atan(x)] <= 0
; log[Гatan] = log[x] - log[x^2+1] - log[atan(x)] <= -|log[x]| - log[atan(x)] + 1 <= 0
(define x (first srcs))

(define x-exp-abs (max (abs (log2-approx (ival-hi x)))
(abs (log2-approx (ival-lo x)))))
(define out-exp (ival-min-log2-approx output))

(list (max 0 (- (- x-exp-abs) out-exp)))]
(list (+ (- (- x-exp-abs) out-exp) 1))]

[(ival-fmod ival-remainder)
; x mod y = x - y*q, where q is rnd_down(x/y)
Expand All @@ -375,8 +390,8 @@
x-slack
(+ x-slack (get-slack)))) ; y crosses zero

(list (max 0 (+ (- x-exp out-exp) 1 x-slack)) ; exponent per x
(max 0 (+ (- x-exp out-exp) 1 y-slack)))] ; exponent per y
(list (+ (- x-exp out-exp) 1 x-slack) ; exponent per x
(+ (- x-exp out-exp) 1 y-slack))] ; exponent per y

[(ival-fma)
; log[Гfma] = log[ max(x*y, -z) / fma(x,y,z)] ~ max(log[x] + log[y], log[z]) - log[fma(x,y,z)] + 1
Expand All @@ -399,7 +414,7 @@

(define lhs-exp (max (+ x-exp y-exp) ; max(log[x] + log[y], log[z])
z-exp))
(make-list 3 (max 0 (+ (- lhs-exp out-exp) 1 slack)))]
(make-list 3 (+ (- lhs-exp out-exp) 1 slack))] ; exponents per arguments

[(ival-hypot)
; hypot = sqrt(x^2+y^2)
Expand All @@ -413,7 +428,7 @@
(define y-exp (ival-max-log2-approx y))
(define out-exp (ival-min-log2-approx output))

(make-list 2 (max 0 (+ (* 2 (- (+ 1 (max x-exp y-exp)) out-exp)) 1)))]
(make-list 2 (+ (* 2 (- (+ 1 (max x-exp y-exp)) out-exp)) 1))]

; Currently log1p has a very poor approximation
[(ival-log1p)
Expand All @@ -432,7 +447,7 @@
(get-slack) ; if x in negative
0))

(list (max 0 (+ (- x-exp out-exp) 1 slack)))]
(list (+ (- x-exp out-exp) 1 slack))]

; Currently expm1 has a very poor solution for negative values
[(ival-expm1)
Expand All @@ -442,18 +457,18 @@
(define x-exp (ival-max-log2-approx x))
(define out-exp (ival-min-log2-approx output))

(list (max 0 (+ 1 x-exp) (+ 2 (- x-exp out-exp))))]
(list (max (+ 1 x-exp) (+ 2 (- x-exp out-exp))))]

[(ival-atan2)
; log[Гatan2]'x = log[Гatan2]'y = log[xy / ((x^2+y^2)*atan2)] <= log[x] + log[y] - 2*max[logx, logy] - log[atan2]
; log[Гatan2]'x = log[Гatan2]'y = log[xy / ((x^2+y^2)*atan2)] <= log[x] + log[y] - 2*max[logx, logy] - log[atan2] + 1
(define x (first srcs))
(define y (second srcs))

(define x-exp (ival-max-log2-approx x))
(define y-exp (ival-max-log2-approx y))
(define out-exp (ival-min-log2-approx output))

(make-list 2 (max 0 (- (+ x-exp y-exp) (* 2 (max x-exp y-exp)) out-exp)))]
(make-list 2 (+ (- (+ x-exp y-exp) (* 2 (max x-exp y-exp)) out-exp) 1))]

; Currently has a poor implementation
[(ival-tanh)
Expand All @@ -462,7 +477,7 @@
(define x-exp (ival-min-log2-approx x))
(define out-exp (ival-max-log2-approx output))

(list (max 0 (+ (- x-exp) out-exp)))]
(list (+ (- x-exp) out-exp))]

[(ival-atanh)
; log[Гarctanh] = log[x / ((1-x^2) * atanh)] = 1 if x < 0.5, otherwise slack
Expand All @@ -478,15 +493,15 @@
[(ival-acosh)
; log[Гacosh] = log[x / (sqrt(x-1) * sqrt(x+1) * acosh)] <= -log[acosh] + slack
(define out-exp (ival-min-log2-approx output))
(define slack (if (< out-exp 2) ; when acosh(x) < 1
(define slack (if (< out-exp 2) ; when acosh(x) < 1
(get-slack)
0))

(list (max 0 (+ (- out-exp) slack)))]
(list (+ (- out-exp) slack))]
; TODO
[(ival-erfc ival-erf ival-lgamma ival-tgamma)
(list (get-slack))]
[else (make-list (length srcs) 0)]))
[else (make-list (length srcs) 0)])) ; exponents for argumetns

(define (log2-approx x)
(define exp (mpfr-exp x))
Expand Down