1 # SPDX-License-Identifier: LGPL-3-or-later
2 # Copyright 2022 Jacob Lifshay programmerjake@gmail.com
4 # Funded by NLnet Assure Programme 2021-02-052, https://nlnet.nl/assure part
5 # of Horizon 2020 EU Programme 957073.
7 from collections
import defaultdict
8 from dataclasses
import dataclass
, field
, fields
, replace
12 from fractions
import Fraction
13 from types
import FunctionType
14 from functools
import lru_cache
15 from nmigen
.hdl
.ast
import Signal
, unsigned
, Mux
, signed
16 from nmigen
.hdl
.dsl
import Module
, Elaboratable
17 from nmigen
.hdl
.mem
import Memory
18 from nmutil
.clz
import CLZ
21 from functools
import cached_property
23 from cached_property
import cached_property
25 # fix broken IDE type detection for cached_property
26 from typing
import TYPE_CHECKING
, Any
28 from functools
import cached_property
34 def cache_on_self(func
):
35 """like `functools.cached_property`, except for methods. unlike
36 `lru_cache` the cache is per-class instance rather than a global cache
39 assert isinstance(func
, FunctionType
), \
40 "non-plain methods are not supported"
42 cache_name
= func
.__name
__ + "__cache"
44 def wrapper(self
, *args
, **kwargs
):
45 # specifically access through `__dict__` to bypass frozen=True
46 cache
= self
.__dict
__.get(cache_name
, _NOT_FOUND
)
47 if cache
is _NOT_FOUND
:
48 self
.__dict
__[cache_name
] = cache
= {}
49 key
= (args
, *kwargs
.items())
50 retval
= cache
.get(key
, _NOT_FOUND
)
51 if retval
is _NOT_FOUND
:
52 retval
= func(self
, *args
, **kwargs
)
56 wrapper
.__doc
__ = func
.__doc
__
61 class RoundDir(enum
.Enum
):
64 NEAREST_TIES_UP
= enum
.auto()
65 ERROR_IF_INEXACT
= enum
.auto()
68 @dataclass(frozen
=True)
73 def __post_init__(self
):
74 # called by the autogenerated __init__
75 assert isinstance(self
.bits
, int)
76 assert isinstance(self
.frac_wid
, int) and self
.frac_wid
>= 0
80 """convert `value` to a fixed-point number with enough fractional
81 bits to preserve its value."""
82 if isinstance(value
, FixedPoint
):
84 if isinstance(value
, int):
85 return FixedPoint(value
, 0)
86 if isinstance(value
, str):
88 neg
= value
.startswith("-")
89 if neg
or value
.startswith("+"):
91 if value
.startswith(("0x", "0X")) and "." in value
:
101 raise ValueError("too many `.` in string")
106 if not digit
.isalnum():
107 raise ValueError("invalid hexadecimal digit")
109 bits |
= int("0x" + digit
, base
=16)
111 bits
= int(value
, base
=0)
115 return FixedPoint(bits
, frac_wid
)
117 if isinstance(value
, float):
118 n
, d
= value
.as_integer_ratio()
119 log2_d
= d
.bit_length() - 1
120 assert d
== 1 << log2_d
, ("d isn't a power of 2 -- won't ever "
121 "fail with float being IEEE 754")
122 return FixedPoint(n
, log2_d
)
123 raise TypeError("can't convert type to FixedPoint")
126 def with_frac_wid(value
, frac_wid
, round_dir
=RoundDir
.ERROR_IF_INEXACT
):
127 """convert `value` to the nearest fixed-point number with `frac_wid`
128 fractional bits, rounding according to `round_dir`."""
129 assert isinstance(frac_wid
, int) and frac_wid
>= 0
130 assert isinstance(round_dir
, RoundDir
)
131 if isinstance(value
, Fraction
):
132 numerator
= value
.numerator
133 denominator
= value
.denominator
135 value
= FixedPoint
.cast(value
)
136 numerator
= value
.bits
137 denominator
= 1 << value
.frac_wid
139 numerator
= -numerator
140 denominator
= -denominator
141 bits
, remainder
= divmod(numerator
<< frac_wid
, denominator
)
142 if round_dir
== RoundDir
.DOWN
:
144 elif round_dir
== RoundDir
.UP
:
147 elif round_dir
== RoundDir
.NEAREST_TIES_UP
:
148 if remainder
* 2 >= denominator
:
150 elif round_dir
== RoundDir
.ERROR_IF_INEXACT
:
152 raise ValueError("inexact conversion")
154 assert False, "unimplemented round_dir"
155 return FixedPoint(bits
, frac_wid
)
157 def to_frac_wid(self
, frac_wid
, round_dir
=RoundDir
.ERROR_IF_INEXACT
):
158 """convert to the nearest fixed-point number with `frac_wid`
159 fractional bits, rounding according to `round_dir`."""
160 return FixedPoint
.with_frac_wid(self
, frac_wid
, round_dir
)
163 # use truediv to get correct result even when bits
164 # and frac_wid are huge
165 return float(self
.bits
/ (1 << self
.frac_wid
))
167 def as_fraction(self
):
168 return Fraction(self
.bits
, 1 << self
.frac_wid
)
171 """compare self with rhs, returning a positive integer if self is
172 greater than rhs, zero if self is equal to rhs, and a negative integer
173 if self is less than rhs."""
174 rhs
= FixedPoint
.cast(rhs
)
175 common_frac_wid
= max(self
.frac_wid
, rhs
.frac_wid
)
176 lhs
= self
.to_frac_wid(common_frac_wid
)
177 rhs
= rhs
.to_frac_wid(common_frac_wid
)
178 return lhs
.bits
- rhs
.bits
180 def __eq__(self
, rhs
):
181 return self
.cmp(rhs
) == 0
183 def __ne__(self
, rhs
):
184 return self
.cmp(rhs
) != 0
186 def __gt__(self
, rhs
):
187 return self
.cmp(rhs
) > 0
189 def __lt__(self
, rhs
):
190 return self
.cmp(rhs
) < 0
192 def __ge__(self
, rhs
):
193 return self
.cmp(rhs
) >= 0
195 def __le__(self
, rhs
):
196 return self
.cmp(rhs
) <= 0
199 """return the fractional part of `self`.
200 that is `self - math.floor(self)`.
202 fract_mask
= (1 << self
.frac_wid
) - 1
203 return FixedPoint(self
.bits
& fract_mask
, self
.frac_wid
)
207 return "-" + str(-self
)
209 frac_digit_count
= (self
.frac_wid
+ digit_bits
- 1) // digit_bits
210 fract
= self
.fract().to_frac_wid(frac_digit_count
* digit_bits
)
211 frac_str
= hex(fract
.bits
)[2:].zfill(frac_digit_count
)
212 return hex(math
.floor(self
)) + "." + frac_str
215 return f
"FixedPoint.with_frac_wid({str(self)!r}, {self.frac_wid})"
217 def __add__(self
, rhs
):
218 rhs
= FixedPoint
.cast(rhs
)
219 common_frac_wid
= max(self
.frac_wid
, rhs
.frac_wid
)
220 lhs
= self
.to_frac_wid(common_frac_wid
)
221 rhs
= rhs
.to_frac_wid(common_frac_wid
)
222 return FixedPoint(lhs
.bits
+ rhs
.bits
, common_frac_wid
)
224 def __radd__(self
, lhs
):
226 return self
.__add
__(lhs
)
229 return FixedPoint(-self
.bits
, self
.frac_wid
)
231 def __sub__(self
, rhs
):
232 rhs
= FixedPoint
.cast(rhs
)
233 common_frac_wid
= max(self
.frac_wid
, rhs
.frac_wid
)
234 lhs
= self
.to_frac_wid(common_frac_wid
)
235 rhs
= rhs
.to_frac_wid(common_frac_wid
)
236 return FixedPoint(lhs
.bits
- rhs
.bits
, common_frac_wid
)
238 def __rsub__(self
, lhs
):
240 return -self
.__sub
__(lhs
)
242 def __mul__(self
, rhs
):
243 rhs
= FixedPoint
.cast(rhs
)
244 return FixedPoint(self
.bits
* rhs
.bits
, self
.frac_wid
+ rhs
.frac_wid
)
246 def __rmul__(self
, lhs
):
248 return self
.__mul
__(lhs
)
251 return self
.bits
>> self
.frac_wid
253 def div(self
, rhs
, frac_wid
, round_dir
=RoundDir
.ERROR_IF_INEXACT
):
254 assert isinstance(frac_wid
, int) and frac_wid
>= 0
255 assert isinstance(round_dir
, RoundDir
)
256 rhs
= FixedPoint
.cast(rhs
)
257 return FixedPoint
.with_frac_wid(self
.as_fraction()
261 def sqrt(self
, round_dir
=RoundDir
.ERROR_IF_INEXACT
):
262 assert isinstance(round_dir
, RoundDir
)
264 raise ValueError("can't compute sqrt of negative number")
267 retval
= FixedPoint(0, self
.frac_wid
)
268 int_part_wid
= self
.bits
.bit_length() - self
.frac_wid
269 first_bit_index
= -(-int_part_wid
// 2) # division rounds up
270 last_bit_index
= -self
.frac_wid
271 for bit_index
in range(first_bit_index
, last_bit_index
- 1, -1):
272 trial
= retval
+ FixedPoint(1 << (bit_index
+ self
.frac_wid
),
274 if trial
* trial
<= self
:
276 if round_dir
== RoundDir
.DOWN
:
278 elif round_dir
== RoundDir
.UP
:
279 if retval
* retval
< self
:
280 retval
+= FixedPoint(1, self
.frac_wid
)
281 elif round_dir
== RoundDir
.NEAREST_TIES_UP
:
282 half_way
= retval
+ FixedPoint(1, self
.frac_wid
+ 1)
283 if half_way
* half_way
<= self
:
284 retval
+= FixedPoint(1, self
.frac_wid
)
285 elif round_dir
== RoundDir
.ERROR_IF_INEXACT
:
286 if retval
* retval
!= self
:
287 raise ValueError("inexact sqrt")
289 assert False, "unimplemented round_dir"
292 def rsqrt(self
, round_dir
=RoundDir
.ERROR_IF_INEXACT
):
293 """compute the reciprocal-sqrt of `self`"""
294 assert isinstance(round_dir
, RoundDir
)
296 raise ValueError("can't compute rsqrt of negative number")
298 raise ZeroDivisionError("can't compute rsqrt of zero")
299 retval
= FixedPoint(0, self
.frac_wid
)
300 first_bit_index
= -(-self
.frac_wid
// 2) # division rounds up
301 last_bit_index
= -self
.frac_wid
302 for bit_index
in range(first_bit_index
, last_bit_index
- 1, -1):
303 trial
= retval
+ FixedPoint(1 << (bit_index
+ self
.frac_wid
),
305 if trial
* trial
* self
<= 1:
307 if round_dir
== RoundDir
.DOWN
:
309 elif round_dir
== RoundDir
.UP
:
310 if retval
* retval
* self
< 1:
311 retval
+= FixedPoint(1, self
.frac_wid
)
312 elif round_dir
== RoundDir
.NEAREST_TIES_UP
:
313 half_way
= retval
+ FixedPoint(1, self
.frac_wid
+ 1)
314 if half_way
* half_way
* self
<= 1:
315 retval
+= FixedPoint(1, self
.frac_wid
)
316 elif round_dir
== RoundDir
.ERROR_IF_INEXACT
:
317 if retval
* retval
* self
!= 1:
318 raise ValueError("inexact rsqrt")
320 assert False, "unimplemented round_dir"
324 class ParamsNotAccurateEnough(Exception):
325 """raised when the parameters aren't accurate enough to have goldschmidt
329 def _assert_accuracy(condition
, msg
="not accurate enough"):
332 raise ParamsNotAccurateEnough(msg
)
335 @dataclass(frozen
=True, unsafe_hash
=True)
336 class GoldschmidtDivParamsBase
:
337 """parameters for a Goldschmidt division algorithm, excluding derived
342 """bit-width of the input divisor and the result.
343 the input numerator is `2 * io_width`-bits wide.
347 """number of bits of additional precision used inside the algorithm."""
350 """the number of address bits used in the lookup-table."""
353 """the number of data bits used in the lookup-table."""
356 """the total number of iterations of the division algorithm's loop"""
359 @dataclass(frozen
=True, unsafe_hash
=True)
360 class GoldschmidtDivParams(GoldschmidtDivParamsBase
):
361 """parameters for a Goldschmidt division algorithm.
362 Use `GoldschmidtDivParams.get` to find a efficient set of parameters.
365 # tuple to be immutable, repr=False so repr() works for debugging even when
366 # __post_init__ hasn't finished running yet
367 table
: "tuple[FixedPoint, ...]" = field(init
=False, repr=False)
368 """the lookup-table"""
370 ops
: "tuple[GoldschmidtDivOp, ...]" = field(init
=False, repr=False)
371 """the operations needed to perform the goldschmidt division algorithm."""
373 def _shrink_bound(self
, bound
, round_dir
):
374 """prevent fractions from having huge numerators/denominators by
375 rounding to a `FixedPoint` and converting back to a `Fraction`.
377 This is intended only for values used to compute bounds, and not for
378 values that end up in the hardware.
380 assert isinstance(bound
, (Fraction
, int))
381 assert round_dir
is RoundDir
.DOWN
or round_dir
is RoundDir
.UP
, \
382 "you shouldn't use that round_dir on bounds"
383 frac_wid
= self
.io_width
* 4 + 100 # should be enough precision
384 fixed
= FixedPoint
.with_frac_wid(bound
, frac_wid
, round_dir
)
385 return fixed
.as_fraction()
387 def _shrink_min(self
, min_bound
):
388 """prevent fractions used as minimum bounds from having huge
389 numerators/denominators by rounding down to a `FixedPoint` and
390 converting back to a `Fraction`.
392 This is intended only for values used to compute bounds, and not for
393 values that end up in the hardware.
395 return self
._shrink
_bound
(min_bound
, RoundDir
.DOWN
)
397 def _shrink_max(self
, max_bound
):
398 """prevent fractions used as maximum bounds from having huge
399 numerators/denominators by rounding up to a `FixedPoint` and
400 converting back to a `Fraction`.
402 This is intended only for values used to compute bounds, and not for
403 values that end up in the hardware.
405 return self
._shrink
_bound
(max_bound
, RoundDir
.UP
)
408 def table_addr_count(self
):
409 """number of distinct addresses in the lookup-table."""
410 # used while computing self.table, so can't just do len(self.table)
411 return 1 << self
.table_addr_bits
413 def table_input_exact_range(self
, addr
):
414 """return the range of inputs as `Fraction`s used for the table entry
415 with address `addr`."""
416 assert isinstance(addr
, int)
417 assert 0 <= addr
< self
.table_addr_count
418 _assert_accuracy(self
.io_width
>= self
.table_addr_bits
)
419 addr_shift
= self
.io_width
- self
.table_addr_bits
420 min_numerator
= (1 << self
.io_width
) + (addr
<< addr_shift
)
421 denominator
= 1 << self
.io_width
422 values_per_table_entry
= 1 << addr_shift
423 max_numerator
= min_numerator
+ values_per_table_entry
- 1
424 min_input
= Fraction(min_numerator
, denominator
)
425 max_input
= Fraction(max_numerator
, denominator
)
426 min_input
= self
._shrink
_min
(min_input
)
427 max_input
= self
._shrink
_max
(max_input
)
428 assert 1 <= min_input
<= max_input
< 2
429 return min_input
, max_input
431 def table_value_exact_range(self
, addr
):
432 """return the range of values as `Fraction`s used for the table entry
433 with address `addr`."""
434 min_input
, max_input
= self
.table_input_exact_range(addr
)
435 # division swaps min/max
436 min_value
= 1 / max_input
437 max_value
= 1 / min_input
438 min_value
= self
._shrink
_min
(min_value
)
439 max_value
= self
._shrink
_max
(max_value
)
440 assert 0.5 < min_value
<= max_value
<= 1
441 return min_value
, max_value
443 def table_exact_value(self
, index
):
444 min_value
, max_value
= self
.table_value_exact_range(index
)
448 def __post_init__(self
):
449 # called by the autogenerated __init__
450 _assert_accuracy(self
.io_width
>= 1, "io_width out of range")
451 _assert_accuracy(self
.extra_precision
>= 0,
452 "extra_precision out of range")
453 _assert_accuracy(self
.table_addr_bits
>= 1,
454 "table_addr_bits out of range")
455 _assert_accuracy(self
.table_data_bits
>= 1,
456 "table_data_bits out of range")
457 _assert_accuracy(self
.iter_count
>= 1, "iter_count out of range")
459 for addr
in range(1 << self
.table_addr_bits
):
460 table
.append(FixedPoint
.with_frac_wid(self
.table_exact_value(addr
),
461 self
.table_data_bits
,
463 # we have to use object.__setattr__ since frozen=True
464 object.__setattr
__(self
, "table", tuple(table
))
465 object.__setattr
__(self
, "ops", tuple(self
.__make
_ops
()))
468 def expanded_width(self
):
469 """the total number of bits of precision used inside the algorithm."""
470 return self
.io_width
+ self
.extra_precision
473 def n_d_f_int_wid(self
):
474 """the number of bits in the integer part of `state.n`, `state.d`, and
475 `state.f` during the main iteration loop.
480 def n_d_f_total_wid(self
):
481 """the total number of bits (both integer and fraction bits) in
482 `state.n`, `state.d`, and `state.f` during the main iteration loop.
484 return self
.n_d_f_int_wid
+ self
.expanded_width
487 def max_neps(self
, i
):
488 """maximum value of `neps[i]`.
489 `neps[i]` is defined to be `n[i] * N_prime[i - 1] * F_prime[i - 1]`.
491 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
492 return Fraction(1, 1 << self
.expanded_width
)
495 def max_deps(self
, i
):
496 """maximum value of `deps[i]`.
497 `deps[i]` is defined to be `d[i] * D_prime[i - 1] * F_prime[i - 1]`.
499 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
500 return Fraction(1, 1 << self
.expanded_width
)
503 def max_feps(self
, i
):
504 """maximum value of `feps[i]`.
505 `feps[i]` is defined to be `f[i] * (2 - D_prime[i - 1])`.
507 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
508 # zero, because the computation of `F_prime[i]` in
509 # `GoldschmidtDivOp.MulDByF.run(...)` is exact.
514 """minimum and maximum values of `e[0]`
515 (the relative error in `F_prime[-1]`)
519 for addr
in range(self
.table_addr_count
):
520 # `F_prime[-1] = (1 - e[0]) / B`
521 # => `e[0] = 1 - B * F_prime[-1]`
522 min_b
, max_b
= self
.table_input_exact_range(addr
)
523 f_prime_m1
= self
.table
[addr
].as_fraction()
524 assert min_b
>= 0 and f_prime_m1
>= 0, \
525 "only positive quadrant of interval multiplication implemented"
526 min_product
= min_b
* f_prime_m1
527 max_product
= max_b
* f_prime_m1
528 # negation swaps min/max
529 cur_min_e0
= 1 - max_product
530 cur_max_e0
= 1 - min_product
531 min_e0
= min(min_e0
, cur_min_e0
)
532 max_e0
= max(max_e0
, cur_max_e0
)
533 min_e0
= self
._shrink
_min
(min_e0
)
534 max_e0
= self
._shrink
_max
(max_e0
)
535 return min_e0
, max_e0
539 """minimum value of `e[0]` (the relative error in `F_prime[-1]`)
541 min_e0
, max_e0
= self
.e0_range
546 """maximum value of `e[0]` (the relative error in `F_prime[-1]`)
548 min_e0
, max_e0
= self
.e0_range
552 def max_abs_e0(self
):
553 """maximum value of `abs(e[0])`."""
554 return max(abs(self
.min_e0
), abs(self
.max_e0
))
557 def min_abs_e0(self
):
558 """minimum value of `abs(e[0])`."""
563 """maximum value of `n[i]` (the relative error in `N_prime[i]`
564 relative to the previous iteration)
566 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
569 # `n[0] = neps[0] / ((1 - e[0]) * (A / B))`
570 # `n[0] <= 2 * neps[0] / (1 - e[0])`
572 assert self
.max_e0
< 1 and self
.max_neps(0) >= 0, \
573 "only one quadrant of interval division implemented"
574 retval
= 2 * self
.max_neps(0) / (1 - self
.max_e0
)
577 # `n[1] <= neps[1] / ((1 - f[0]) * (1 - pi[0] - delta[0]))`
578 min_mpd
= 1 - self
.max_pi(0) - self
.max_delta(0)
579 assert self
.max_f(0) <= 1 and min_mpd
>= 0, \
580 "only one quadrant of interval multiplication implemented"
581 prod
= (1 - self
.max_f(0)) * min_mpd
582 assert self
.max_neps(1) >= 0 and prod
> 0, \
583 "only one quadrant of interval division implemented"
584 retval
= self
.max_neps(1) / prod
587 # `0 <= n[i] <= 2 * max_neps[i] / (1 - pi[i - 1] - delta[i - 1])`
588 min_mpd
= 1 - self
.max_pi(i
- 1) - self
.max_delta(i
- 1)
589 assert self
.max_neps(i
) >= 0 and min_mpd
> 0, \
590 "only one quadrant of interval division implemented"
591 retval
= self
.max_neps(i
) / min_mpd
593 return self
._shrink
_max
(retval
)
597 """maximum value of `d[i]` (the relative error in `D_prime[i]`
598 relative to the previous iteration)
600 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
603 # `d[0] = deps[0] / (1 - e[0])`
605 assert self
.max_e0
< 1 and self
.max_deps(0) >= 0, \
606 "only one quadrant of interval division implemented"
607 retval
= self
.max_deps(0) / (1 - self
.max_e0
)
610 # `d[1] <= deps[1] / ((1 - f[0]) * (1 - delta[0] ** 2))`
611 assert self
.max_f(0) <= 1 and self
.max_delta(0) <= 1, \
612 "only one quadrant of interval multiplication implemented"
613 divisor
= (1 - self
.max_f(0)) * (1 - self
.max_delta(0) ** 2)
614 assert self
.max_deps(1) >= 0 and divisor
> 0, \
615 "only one quadrant of interval division implemented"
616 retval
= self
.max_deps(1) / divisor
619 # `0 <= d[i] <= max_deps[i] / (1 - delta[i - 1])`
620 assert self
.max_deps(i
) >= 0 and self
.max_delta(i
- 1) < 1, \
621 "only one quadrant of interval division implemented"
622 retval
= self
.max_deps(i
) / (1 - self
.max_delta(i
- 1))
624 return self
._shrink
_max
(retval
)
628 """maximum value of `f[i]` (the relative error in `F_prime[i]`
629 relative to the previous iteration)
631 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
634 # `f[0] = feps[0] / (1 - delta[0])`
636 assert self
.max_delta(0) < 1 and self
.max_feps(0) >= 0, \
637 "only one quadrant of interval division implemented"
638 retval
= self
.max_feps(0) / (1 - self
.max_delta(0))
642 retval
= self
.max_feps(1)
645 # `f[i] <= max_feps[i]`
646 retval
= self
.max_feps(i
)
648 return self
._shrink
_max
(retval
)
651 def max_delta(self
, i
):
652 """ maximum value of `delta[i]`.
653 `delta[i]` is defined in Definition 4 of paper.
655 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
657 # `delta[0] = abs(e[0]) + 3 * d[0] / 2`
658 retval
= self
.max_abs_e0
+ Fraction(3, 2) * self
.max_d(0)
660 # `delta[i] = delta[i - 1] ** 2 + f[i - 1]`
661 prev_max_delta
= self
.max_delta(i
- 1)
662 assert prev_max_delta
>= 0
663 retval
= prev_max_delta
** 2 + self
.max_f(i
- 1)
665 # `delta[i]` has to be smaller than one otherwise errors would go off
667 _assert_accuracy(retval
< 1)
669 return self
._shrink
_max
(retval
)
673 """ maximum value of `pi[i]`.
674 `pi[i]` is defined right below Theorem 5 of paper.
676 assert isinstance(i
, int) and 0 <= i
< self
.iter_count
677 # `pi[i] = 1 - (1 - n[i]) * prod`
678 # where `prod` is the product of,
679 # for `j` in `0 <= j < i`, `(1 - n[j]) / (1 + d[j])`
680 min_prod
= Fraction(1)
682 max_n_j
= self
.max_n(j
)
683 max_d_j
= self
.max_d(j
)
684 assert max_n_j
<= 1 and max_d_j
> -1, \
685 "only one quadrant of interval division implemented"
686 min_prod
*= (1 - max_n_j
) / (1 + max_d_j
)
687 max_n_i
= self
.max_n(i
)
688 assert max_n_i
<= 1 and min_prod
>= 0, \
689 "only one quadrant of interval multiplication implemented"
690 retval
= 1 - (1 - max_n_i
) * min_prod
691 return self
._shrink
_max
(retval
)
694 def max_n_shift(self
):
695 """ maximum value of `state.n_shift`.
697 # input numerator is `2*io_width`-bits
698 max_n
= (1 << (self
.io_width
* 2)) - 1
700 # normalize so 1 <= n < 2
708 """ maximum value of, for all `i`, `max_n(i)` and `max_d(i)`
711 for i
in range(self
.iter_count
):
712 n_hat
= max(n_hat
, self
.max_n(i
), self
.max_d(i
))
713 return self
._shrink
_max
(n_hat
)
715 def __make_ops(self
):
716 """ Goldschmidt division algorithm.
719 Even, G., Seidel, P. M., & Ferguson, W. E. (2003).
720 A Parametric Error Analysis of Goldschmidt's Division Algorithm.
721 https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.90.1238&rep=rep1&type=pdf
723 yields: GoldschmidtDivOp
724 the operations needed to perform the division.
726 # establish assumptions of the paper's error analysis (section 3.1):
728 # 1. normalize so A (numerator) and B (denominator) are in [1, 2)
729 yield GoldschmidtDivOp
.Normalize
731 # 2. ensure all relative errors from directed rounding are <= 1 / 4.
732 # the assumption is met by multipliers with > 4-bits precision
733 _assert_accuracy(self
.expanded_width
> 4)
735 # 3. require `abs(e[0]) + 3 * d[0] / 2 + f[0] < 1 / 2`.
736 _assert_accuracy(self
.max_abs_e0
+ 3 * self
.max_d(0) / 2
737 + self
.max_f(0) < Fraction(1, 2))
739 # 4. the initial approximation F'[-1] of 1/B is in [1/2, 1].
740 # (B is the denominator)
742 for addr
in range(self
.table_addr_count
):
743 f_prime_m1
= self
.table
[addr
]
744 _assert_accuracy(0.5 <= f_prime_m1
<= 1)
746 yield GoldschmidtDivOp
.FEqTableLookup
748 # we use Setting I (section 4.1 of the paper):
749 # Require `n[i] <= n_hat` and `d[i] <= n_hat` and `f[i] = 0`:
750 # the conditions on n_hat are satisfied by construction.
751 for i
in range(self
.iter_count
):
752 _assert_accuracy(self
.max_f(i
) == 0)
753 yield GoldschmidtDivOp
.MulNByF
754 if i
!= self
.iter_count
- 1:
755 yield GoldschmidtDivOp
.MulDByF
756 yield GoldschmidtDivOp
.FEq2MinusD
758 # relative approximation error `p(N_prime[i])`:
759 # `p(N_prime[i]) = (A / B - N_prime[i]) / (A / B)`
760 # `0 <= p(N_prime[i])`
761 # `p(N_prime[i]) <= (2 * i) * n_hat \`
762 # ` + (abs(e[0]) + 3 * n_hat / 2) ** (2 ** i)`
763 i
= self
.iter_count
- 1 # last used `i`
764 # compute power manually to prevent huge intermediate values
765 power
= self
._shrink
_max
(self
.max_abs_e0
+ 3 * self
.n_hat
/ 2)
767 power
= self
._shrink
_max
(power
* power
)
769 max_rel_error
= (2 * i
) * self
.n_hat
+ power
771 min_a_over_b
= Fraction(1, 2)
772 max_a_over_b
= Fraction(2)
773 max_allowed_abs_error
= max_a_over_b
/ (1 << self
.max_n_shift
)
774 max_allowed_rel_error
= max_allowed_abs_error
/ min_a_over_b
776 _assert_accuracy(max_rel_error
< max_allowed_rel_error
,
777 f
"not accurate enough: max_rel_error={max_rel_error}"
778 f
" max_allowed_rel_error={max_allowed_rel_error}")
780 yield GoldschmidtDivOp
.CalcResult
783 def default_cost_fn(self
):
784 """ calculate the estimated cost on an arbitrary scale of implementing
785 goldschmidt division with the specified parameters. larger cost
786 values mean worse parameters.
788 This is the default cost function for `GoldschmidtDivParams.get`.
792 rom_cells
= self
.table_data_bits
<< self
.table_addr_bits
793 cost
= float(rom_cells
)
795 if op
== GoldschmidtDivOp
.MulNByF \
796 or op
== GoldschmidtDivOp
.MulDByF
:
797 mul_cost
= self
.expanded_width
** 2
798 mul_cost
*= self
.expanded_width
.bit_length()
800 cost
+= 5e7
* self
.iter_count
804 @lru_cache(maxsize
=1 << 16)
805 def __cached_new(base_params
):
806 assert isinstance(base_params
, GoldschmidtDivParamsBase
)
807 # can't use dataclasses.asdict, since it's recursive and will also give
808 # child class fields too, which we don't want.
810 for field
in fields(GoldschmidtDivParamsBase
):
811 kwargs
[field
.name
] = getattr(base_params
, field
.name
)
813 return GoldschmidtDivParams(**kwargs
), None
814 except ParamsNotAccurateEnough
as e
:
818 def __raise(e
): # type: (ParamsNotAccurateEnough) -> Any
822 def cached_new(base_params
, handle_error
=__raise
):
823 assert isinstance(base_params
, GoldschmidtDivParamsBase
)
824 params
, error
= GoldschmidtDivParams
.__cached
_new
(base_params
)
828 return handle_error(error
)
831 def get(io_width
, cost_fn
=default_cost_fn
, max_table_addr_bits
=12):
832 """ find efficient parameters for a goldschmidt division algorithm
833 with `params.io_width == io_width`.
837 bit-width of the input divisor and the result.
838 the input numerator is `2 * io_width`-bits wide.
839 cost_fn: Callable[[GoldschmidtDivParams], float]
840 return the estimated cost on an arbitrary scale of implementing
841 goldschmidt division with the specified parameters. larger cost
842 values mean worse parameters.
843 max_table_addr_bits: int
844 maximum allowable value of `table_addr_bits`
846 assert isinstance(io_width
, int) and io_width
>= 1
847 assert callable(cost_fn
)
850 last_error_params
= None
852 def cached_new(base_params
):
854 nonlocal last_error
, last_error_params
856 last_error_params
= base_params
859 retval
= GoldschmidtDivParams
.cached_new(base_params
, handle_error
)
861 logging
.debug(f
"GoldschmidtDivParams.get: err: {base_params}")
863 logging
.debug(f
"GoldschmidtDivParams.get: ok: {base_params}")
866 @lru_cache(maxsize
=None)
867 def get_cost(base_params
):
868 params
= cached_new(base_params
)
871 retval
= cost_fn(params
)
872 logging
.debug(f
"GoldschmidtDivParams.get: cost={retval}: {params}")
875 # start with parameters big enough to always work.
876 initial_extra_precision
= io_width
* 2 + 4
877 initial_params
= GoldschmidtDivParamsBase(
879 extra_precision
=initial_extra_precision
,
880 table_addr_bits
=min(max_table_addr_bits
, io_width
),
881 table_data_bits
=io_width
+ initial_extra_precision
,
882 iter_count
=1 + io_width
.bit_length())
884 if cached_new(initial_params
) is None:
885 raise ValueError(f
"initial goldschmidt division algorithm "
886 f
"parameters are invalid: {initial_params}"
889 # find good initial `iter_count`
890 params
= initial_params
891 for iter_count
in range(1, initial_params
.iter_count
):
892 trial_params
= replace(params
, iter_count
=iter_count
)
893 if cached_new(trial_params
) is not None:
894 params
= trial_params
897 # now find `table_addr_bits`
898 cost
= get_cost(params
)
899 for table_addr_bits
in range(1, max_table_addr_bits
):
900 trial_params
= replace(params
, table_addr_bits
=table_addr_bits
)
901 trial_cost
= get_cost(trial_params
)
902 if trial_cost
< cost
:
903 params
= trial_params
907 # check one higher `iter_count` to see if it has lower cost
908 for table_addr_bits
in range(1, max_table_addr_bits
+ 1):
909 trial_params
= replace(params
,
910 table_addr_bits
=table_addr_bits
,
911 iter_count
=params
.iter_count
+ 1)
912 trial_cost
= get_cost(trial_params
)
913 if trial_cost
< cost
:
914 params
= trial_params
918 # now shrink `table_data_bits`
920 trial_params
= replace(params
,
921 table_data_bits
=params
.table_data_bits
- 1)
922 trial_cost
= get_cost(trial_params
)
923 if trial_cost
< cost
:
924 params
= trial_params
929 # and shrink `extra_precision`
931 trial_params
= replace(params
,
932 extra_precision
=params
.extra_precision
- 1)
933 trial_cost
= get_cost(trial_params
)
934 if trial_cost
< cost
:
935 params
= trial_params
940 return cached_new(params
)
944 """count leading zeros -- handy for debugging."""
945 assert isinstance(wid
, int)
946 assert isinstance(v
, int) and 0 <= v
< (1 << wid
)
947 return (1 << wid
).bit_length() - v
.bit_length()
951 class GoldschmidtDivOp(enum
.Enum
):
952 Normalize
= "n, d, n_shift = normalize(n, d)"
953 FEqTableLookup
= "f = table_lookup(d)"
956 FEq2MinusD
= "f = 2 - d"
957 CalcResult
= "result = unnormalize_and_round(n)"
959 def run(self
, params
, state
):
960 assert isinstance(params
, GoldschmidtDivParams
)
961 assert isinstance(state
, GoldschmidtDivState
)
962 expanded_width
= params
.expanded_width
963 table_addr_bits
= params
.table_addr_bits
964 if self
== GoldschmidtDivOp
.Normalize
:
965 # normalize so 1 <= d < 2
966 # can easily be done with count-leading-zeros and left shift
968 state
.n
= (state
.n
* 2).to_frac_wid(expanded_width
)
969 state
.d
= (state
.d
* 2).to_frac_wid(expanded_width
)
972 # normalize so 1 <= n < 2
974 state
.n
= (state
.n
* 0.5).to_frac_wid(expanded_width
)
976 elif self
== GoldschmidtDivOp
.FEqTableLookup
:
977 # compute initial f by table lookup
979 d_m_1
= d_m_1
.to_frac_wid(table_addr_bits
, RoundDir
.DOWN
)
980 assert 0 <= d_m_1
.bits
< (1 << params
.table_addr_bits
)
981 state
.f
= params
.table
[d_m_1
.bits
]
982 elif self
== GoldschmidtDivOp
.MulNByF
:
983 assert state
.f
is not None
984 n
= state
.n
* state
.f
985 state
.n
= n
.to_frac_wid(expanded_width
, round_dir
=RoundDir
.DOWN
)
986 elif self
== GoldschmidtDivOp
.MulDByF
:
987 assert state
.f
is not None
988 d
= state
.d
* state
.f
989 state
.d
= d
.to_frac_wid(expanded_width
, round_dir
=RoundDir
.UP
)
990 elif self
== GoldschmidtDivOp
.FEq2MinusD
:
991 state
.f
= (2 - state
.d
).to_frac_wid(expanded_width
)
992 elif self
== GoldschmidtDivOp
.CalcResult
:
993 assert state
.n_shift
is not None
994 # scale to correct value
995 n
= state
.n
* (1 << state
.n_shift
)
997 state
.quotient
= math
.floor(n
)
998 state
.remainder
= state
.orig_n
- state
.quotient
* state
.orig_d
999 if state
.remainder
>= state
.orig_d
:
1001 state
.remainder
-= state
.orig_d
1003 assert False, f
"unimplemented GoldschmidtDivOp: {self}"
1005 def gen_hdl(self
, params
, state
, sync_rom
):
1006 # FIXME: finish getting hdl/simulation to work
1007 """generate the hdl for this operation.
1010 params: GoldschmidtDivParams
1011 the goldschmidt division parameters.
1012 state: GoldschmidtDivHDLState
1013 the input/output state
1015 true if the rom should be read synchronously rather than
1016 combinatorially, incurring an extra clock cycle of latency.
1018 assert isinstance(params
, GoldschmidtDivParams
)
1019 assert isinstance(state
, GoldschmidtDivHDLState
)
1021 expanded_width
= params
.expanded_width
1022 table_addr_bits
= params
.table_addr_bits
1023 if self
== GoldschmidtDivOp
.Normalize
:
1024 # normalize so 1 <= d < 2
1025 assert state
.d
.width
== params
.io_width
1026 assert state
.n
.width
== 2 * params
.io_width
1027 d_leading_zeros
= CLZ(params
.io_width
)
1028 m
.submodules
.d_leading_zeros
= d_leading_zeros
1029 m
.d
.comb
+= d_leading_zeros
.sig_in
.eq(state
.d
)
1030 d_shift_out
= Signal
.like(state
.d
)
1031 m
.d
.comb
+= d_shift_out
.eq(state
.d
<< d_leading_zeros
.lz
)
1032 state
.d
= Signal(params
.n_d_f_total_wid
)
1033 m
.d
.comb
+= state
.d
.eq(d_shift_out
<< (params
.extra_precision
1034 + params
.n_d_f_int_wid
))
1036 # normalize so 1 <= n < 2
1037 n_leading_zeros
= CLZ(2 * params
.io_width
)
1038 m
.submodules
.n_leading_zeros
= n_leading_zeros
1039 m
.d
.comb
+= n_leading_zeros
.sig_in
.eq(state
.n
)
1040 n_shift_s_v
= (params
.io_width
+ d_leading_zeros
.lz
1041 - n_leading_zeros
.lz
)
1042 n_shift_s
= Signal
.like(n_shift_s_v
)
1043 state
.n_shift
= Signal(d_leading_zeros
.lz
.width
)
1045 n_shift_s
.eq(n_shift_s_v
),
1046 state
.n_shift
.eq(Mux(n_shift_s
< 0, 0, n_shift_s
)),
1048 n
= Signal(params
.n_d_f_total_wid
)
1049 shifted_n
= state
.n
<< state
.n_shift
1050 fixed_shift
= params
.expanded_width
- state
.n
.width
1051 m
.d
.comb
+= n
.eq(shifted_n
<< fixed_shift
)
1053 elif self
== GoldschmidtDivOp
.FEqTableLookup
:
1054 assert state
.d
.width
== params
.n_d_f_total_wid
, "invalid d width"
1055 # compute initial f by table lookup
1057 # extra bit for table entries == 1.0
1058 table_width
= 1 + params
.table_data_bits
1059 table
= Memory(width
=table_width
, depth
=len(params
.table
),
1060 init
=[i
.bits
for i
in params
.table
])
1061 addr
= state
.d
[:-params
.n_d_f_int_wid
][-table_addr_bits
:]
1063 table_read
= table
.read_port()
1064 m
.d
.comb
+= table_read
.addr
.eq(addr
)
1065 state
.insert_pipeline_register()
1067 table_read
= table
.read_port(domain
="comb")
1068 m
.d
.comb
+= table_read
.addr
.eq(addr
)
1069 m
.submodules
.table_read
= table_read
1070 state
.f
= Signal(params
.n_d_f_int_wid
+ params
.expanded_width
)
1071 data_shift
= (table_width
- params
.table_data_bits
1072 + params
.expanded_width
)
1073 m
.d
.comb
+= state
.f
.eq(table_read
.data
<< data_shift
)
1074 elif self
== GoldschmidtDivOp
.MulNByF
:
1075 assert state
.n
.width
== params
.n_d_f_total_wid
, "invalid n width"
1076 assert state
.f
is not None
1077 assert state
.f
.width
== params
.n_d_f_total_wid
, "invalid f width"
1078 n
= Signal
.like(state
.n
)
1079 m
.d
.comb
+= n
.eq((state
.n
* state
.f
) >> params
.expanded_width
)
1081 elif self
== GoldschmidtDivOp
.MulDByF
:
1082 assert state
.d
.width
== params
.n_d_f_total_wid
, "invalid d width"
1083 assert state
.f
is not None
1084 assert state
.f
.width
== params
.n_d_f_total_wid
, "invalid f width"
1085 d
= Signal
.like(state
.d
)
1086 m
.d
.comb
+= d
.eq((state
.d
* state
.f
) >> params
.expanded_width
)
1088 elif self
== GoldschmidtDivOp
.FEq2MinusD
:
1089 assert state
.d
.width
== params
.n_d_f_total_wid
, "invalid d width"
1090 f
= Signal
.like(state
.d
)
1091 m
.d
.comb
+= f
.eq((2 << params
.expanded_width
) - state
.d
)
1093 elif self
== GoldschmidtDivOp
.CalcResult
:
1094 assert state
.n
.width
== params
.n_d_f_total_wid
, "invalid n width"
1095 assert state
.n_shift
is not None
1096 # scale to correct value
1097 n
= state
.n
* (1 << state
.n_shift
)
1098 q_approx
= Signal(params
.io_width
)
1099 # extra bit for if it's bigger than orig_d
1100 r_approx
= Signal(params
.io_width
+ 1)
1101 adjusted_r
= Signal(signed(1 + params
.io_width
))
1103 q_approx
.eq((state
.n
<< state
.n_shift
)
1104 >> params
.expanded_width
),
1105 r_approx
.eq(state
.orig_n
- q_approx
* state
.orig_d
),
1106 adjusted_r
.eq(r_approx
- state
.orig_d
),
1108 state
.quotient
= Signal(params
.io_width
)
1109 state
.remainder
= Signal(params
.io_width
)
1111 with m
.If(adjusted_r
>= 0):
1113 state
.quotient
.eq(q_approx
+ 1),
1114 state
.remainder
.eq(adjusted_r
),
1118 state
.quotient
.eq(q_approx
),
1119 state
.remainder
.eq(r_approx
),
1122 assert False, f
"unimplemented GoldschmidtDivOp: {self}"
1126 class GoldschmidtDivState
:
1128 """original numerator"""
1131 """original denominator"""
1134 """numerator -- N_prime[i] in the paper's algorithm 2"""
1137 """denominator -- D_prime[i] in the paper's algorithm 2"""
1139 f
: "FixedPoint | None" = None
1140 """current factor -- F_prime[i] in the paper's algorithm 2"""
1142 quotient
: "int | None" = None
1143 """final quotient"""
1145 remainder
: "int | None" = None
1146 """final remainder"""
1148 n_shift
: "int | None" = None
1149 """amount the numerator needs to be left-shifted at the end of the
1154 def goldschmidt_div(n
, d
, params
):
1155 """ Goldschmidt division algorithm.
1158 Even, G., Seidel, P. M., & Ferguson, W. E. (2003).
1159 A Parametric Error Analysis of Goldschmidt's Division Algorithm.
1160 https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.90.1238&rep=rep1&type=pdf
1164 numerator. a `2*width`-bit unsigned integer.
1165 must be less than `d << width`, otherwise the quotient wouldn't
1166 fit in `width` bits.
1168 denominator. a `width`-bit unsigned integer. must not be zero.
1170 the bit-width of the inputs/outputs. must be a positive integer.
1172 returns: tuple[int, int]
1173 the quotient and remainder. a tuple of two `width`-bit unsigned
1176 assert isinstance(params
, GoldschmidtDivParams
)
1177 assert isinstance(d
, int) and 0 < d
< (1 << params
.io_width
)
1178 assert isinstance(n
, int) and 0 <= n
< (d
<< params
.io_width
)
1180 # this whole algorithm is done with fixed-point arithmetic where values
1181 # have `width` fractional bits
1183 state
= GoldschmidtDivState(
1186 n
=FixedPoint(n
, params
.io_width
),
1187 d
=FixedPoint(d
, params
.io_width
),
1190 for op
in params
.ops
:
1191 op
.run(params
, state
)
1193 assert state
.quotient
is not None
1194 assert state
.remainder
is not None
1196 return state
.quotient
, state
.remainder
1199 @dataclass(eq
=False)
1200 class GoldschmidtDivHDLState
:
1202 """The HDL Module"""
1205 """original numerator"""
1208 """original denominator"""
1211 """numerator -- N_prime[i] in the paper's algorithm 2"""
1214 """denominator -- D_prime[i] in the paper's algorithm 2"""
1216 f
: "Signal | None" = None
1217 """current factor -- F_prime[i] in the paper's algorithm 2"""
1219 quotient
: "Signal | None" = None
1220 """final quotient"""
1222 remainder
: "Signal | None" = None
1223 """final remainder"""
1225 n_shift
: "Signal | None" = None
1226 """amount the numerator needs to be left-shifted at the end of the
1230 old_signals
: "defaultdict[str, list[Signal]]" = field(repr=False,
1233 __signal_name_prefix
: "str" = field(default
="state_", repr=False,
1236 def __post_init__(self
):
1237 # called by the autogenerated __init__
1238 self
.old_signals
= defaultdict(list)
1240 def __setattr__(self
, name
, value
):
1241 assert isinstance(name
, str)
1242 if name
.startswith("_"):
1243 return super().__setattr
__(name
, value
)
1245 old_signals
= self
.old_signals
[name
]
1246 except AttributeError:
1247 # haven't yet finished __post_init__
1248 return super().__setattr
__(name
, value
)
1249 assert name
!= "m" and name
!= "old_signals", f
"can't write to {name}"
1250 assert isinstance(value
, Signal
)
1251 value
.name
= f
"{self.__signal_name_prefix}{name}_{len(old_signals)}"
1252 old_signal
= getattr(self
, name
, None)
1253 if old_signal
is not None:
1254 assert isinstance(old_signal
, Signal
)
1255 old_signals
.append(old_signal
)
1256 return super().__setattr
__(name
, value
)
1258 def insert_pipeline_register(self
):
1259 old_prefix
= self
.__signal
_name
_prefix
1261 for field
in fields(GoldschmidtDivHDLState
):
1262 if field
.name
.startswith("_") or field
.name
== "m":
1264 old_sig
= getattr(self
, field
.name
, None)
1267 assert isinstance(old_sig
, Signal
)
1268 new_sig
= Signal
.like(old_sig
)
1269 setattr(self
, field
.name
, new_sig
)
1270 self
.m
.d
.sync
+= new_sig
.eq(old_sig
)
1272 self
.__signal
_name
_prefix
= old_prefix
1275 class GoldschmidtDivHDL(Elaboratable
):
1276 # FIXME: finish getting hdl/simulation to work
1277 """ Goldschmidt division algorithm.
1280 Even, G., Seidel, P. M., & Ferguson, W. E. (2003).
1281 A Parametric Error Analysis of Goldschmidt's Division Algorithm.
1282 https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.90.1238&rep=rep1&type=pdf
1285 params: GoldschmidtDivParams
1286 the goldschmidt division algorithm parameters.
1287 pipe_reg_indexes: list[int]
1288 the operation indexes where pipeline registers should be inserted.
1289 duplicate values mean multiple registers should be inserted for
1290 that operation index -- this is useful to allow yosys to spread a
1291 multiplication across those multiple pipeline stages.
1293 true if the rom should be read synchronously rather than
1294 combinatorially, incurring an extra clock cycle of latency.
1295 n: Signal(unsigned(2 * params.io_width))
1296 input numerator. a `2 * params.io_width`-bit unsigned integer.
1297 must be less than `d << params.io_width`, otherwise the quotient
1298 wouldn't fit in `params.io_width` bits.
1299 d: Signal(unsigned(params.io_width))
1300 input denominator. a `params.io_width`-bit unsigned integer.
1302 q: Signal(unsigned(params.io_width))
1303 output quotient. only valid when `n < (d << params.io_width)`.
1304 r: Signal(unsigned(params.io_width))
1305 output remainder. only valid when `n < (d << params.io_width)`.
1309 def total_pipeline_registers(self
):
1310 """the total number of pipeline registers"""
1311 return len(self
.pipe_reg_indexes
) + self
.sync_rom
1313 def __init__(self
, params
, pipe_reg_indexes
=(), sync_rom
=False):
1314 assert isinstance(params
, GoldschmidtDivParams
)
1315 assert isinstance(sync_rom
, bool)
1316 self
.params
= params
1317 self
.pipe_reg_indexes
= sorted(int(i
) for i
in pipe_reg_indexes
)
1318 self
.sync_rom
= sync_rom
1319 self
.n
= Signal(unsigned(2 * params
.io_width
))
1320 self
.d
= Signal(unsigned(params
.io_width
))
1321 self
.q
= Signal(unsigned(params
.io_width
))
1322 self
.r
= Signal(unsigned(params
.io_width
))
1324 def elaborate(self
, platform
):
1326 state
= GoldschmidtDivHDLState(
1334 pipe_reg_indexes
= list(reversed(self
.pipe_reg_indexes
))
1336 for op_index
, op
in enumerate(self
.params
.ops
):
1337 while len(pipe_reg_indexes
) > 0 \
1338 and pipe_reg_indexes
[-1] <= op_index
:
1339 pipe_reg_indexes
.pop()
1340 state
.insert_pipeline_register()
1341 op
.gen_hdl(self
.params
, state
, self
.sync_rom
)
1343 while len(pipe_reg_indexes
) > 0:
1344 pipe_reg_indexes
.pop()
1345 state
.insert_pipeline_register()
1347 m
.d
.comb
+= self
.q
.eq(state
.quotient
)
1348 m
.d
.comb
+= self
.r
.eq(state
.remainder
)
1352 GOLDSCHMIDT_SQRT_RSQRT_TABLE_ADDR_INT_WID
= 2
1356 def goldschmidt_sqrt_rsqrt_table(table_addr_bits
, table_data_bits
):
1357 """Generate the look-up table needed for Goldschmidt's square-root and
1358 reciprocal-square-root algorithm.
1361 table_addr_bits: int
1362 the number of address bits for the look-up table.
1363 table_data_bits: int
1364 the number of data bits for the look-up table.
1366 assert isinstance(table_addr_bits
, int) and \
1367 table_addr_bits
>= GOLDSCHMIDT_SQRT_RSQRT_TABLE_ADDR_INT_WID
1368 assert isinstance(table_data_bits
, int) and table_data_bits
>= 1
1370 table_len
= 1 << table_addr_bits
1371 for addr
in range(table_len
):
1373 value
= FixedPoint(0, table_data_bits
)
1374 elif (addr
<< 2) < table_len
:
1375 value
= None # table entries should be unused
1377 table_addr_frac_wid
= table_addr_bits
1378 table_addr_frac_wid
-= GOLDSCHMIDT_SQRT_RSQRT_TABLE_ADDR_INT_WID
1379 max_input_value
= FixedPoint(addr
+ 1, table_addr_bits
- 2)
1380 max_frac_wid
= max(max_input_value
.frac_wid
, table_data_bits
)
1381 value
= max_input_value
.to_frac_wid(max_frac_wid
)
1382 value
= value
.rsqrt(RoundDir
.DOWN
)
1383 value
= value
.to_frac_wid(table_data_bits
, RoundDir
.DOWN
)
1386 # tuple for immutability
1389 # FIXME: add code to calculate error bounds and check that the algorithm will
1390 # actually work (like in the goldschmidt division algorithm).
1391 # FIXME: add code to calculate a good set of parameters based on the error
1395 def goldschmidt_sqrt_rsqrt(radicand
, io_width
, frac_wid
, extra_precision
,
1396 table_addr_bits
, table_data_bits
, iter_count
):
1397 """Goldschmidt's square-root and reciprocal-square-root algorithm.
1399 uses algorithm based on second method at:
1400 https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Goldschmidt%E2%80%99s_algorithm
1403 radicand: FixedPoint(frac_wid=frac_wid)
1404 the input value to take the square-root and reciprocal-square-root of.
1406 the number of bits in the input (`radicand`) and output values.
1408 the number of fraction bits in the input (`radicand`) and output
1410 extra_precision: int
1411 the number of bits of internal extra precision.
1412 table_addr_bits: int
1413 the number of address bits for the look-up table.
1414 table_data_bits: int
1415 the number of data bits for the look-up table.
1417 returns: tuple[FixedPoint, FixedPoint]
1418 the square-root and reciprocal-square-root, rounded down to the
1419 nearest representable value. If `radicand == 0`, then the
1420 reciprocal-square-root value returned is zero.
1422 assert (isinstance(radicand
, FixedPoint
)
1423 and radicand
.frac_wid
== frac_wid
1424 and 0 <= radicand
.bits
< (1 << io_width
))
1425 assert isinstance(io_width
, int) and io_width
>= 1
1426 assert isinstance(frac_wid
, int) and 0 <= frac_wid
< io_width
1427 assert isinstance(extra_precision
, int) and extra_precision
>= io_width
1428 assert isinstance(table_addr_bits
, int) and table_addr_bits
>= 1
1429 assert isinstance(table_data_bits
, int) and table_data_bits
>= 1
1430 assert isinstance(iter_count
, int) and iter_count
>= 0
1431 expanded_frac_wid
= frac_wid
+ extra_precision
1432 s
= radicand
.to_frac_wid(expanded_frac_wid
)
1433 sqrt_rshift
= extra_precision
1434 rsqrt_rshift
= extra_precision
1435 while s
!= 0 and s
< 1:
1436 s
= (s
* 4).to_frac_wid(expanded_frac_wid
)
1440 s
= s
.div(4, expanded_frac_wid
)
1443 table
= goldschmidt_sqrt_rsqrt_table(table_addr_bits
=table_addr_bits
,
1444 table_data_bits
=table_data_bits
)
1445 # core goldschmidt sqrt/rsqrt algorithm:
1447 table_addr_frac_wid
= table_addr_bits
1448 table_addr_frac_wid
-= GOLDSCHMIDT_SQRT_RSQRT_TABLE_ADDR_INT_WID
1449 addr
= s
.to_frac_wid(table_addr_frac_wid
, RoundDir
.DOWN
)
1450 assert 0 <= addr
.bits
< (1 << table_addr_bits
), "table addr out of range"
1451 f
= table
[addr
.bits
]
1452 assert f
is not None, "accessed invalid table entry"
1453 # use with_frac_wid to fix IDE type deduction
1454 f
= FixedPoint
.with_frac_wid(f
, expanded_frac_wid
, RoundDir
.DOWN
)
1455 x
= (s
* f
).to_frac_wid(expanded_frac_wid
, RoundDir
.DOWN
)
1456 h
= (f
* 0.5).to_frac_wid(expanded_frac_wid
, RoundDir
.DOWN
)
1457 for _
in range(iter_count
):
1459 f
= (1.5 - x
* h
).to_frac_wid(expanded_frac_wid
, RoundDir
.DOWN
)
1460 x
= (x
* f
).to_frac_wid(expanded_frac_wid
, RoundDir
.DOWN
)
1461 h
= (h
* f
).to_frac_wid(expanded_frac_wid
, RoundDir
.DOWN
)
1463 # now `x` is approximately `sqrt(s)` and `r` is approximately `rsqrt(s)`
1465 sqrt
= FixedPoint(x
.bits
>> sqrt_rshift
, frac_wid
)
1466 rsqrt
= FixedPoint(r
.bits
>> rsqrt_rshift
, frac_wid
)
1468 next_sqrt
= FixedPoint(sqrt
.bits
+ 1, frac_wid
)
1469 if next_sqrt
* next_sqrt
<= radicand
:
1472 next_rsqrt
= FixedPoint(rsqrt
.bits
+ 1, frac_wid
)
1473 if next_rsqrt
* next_rsqrt
* radicand
<= 1 and radicand
!= 0: