中期检查PR #4
|
|
@ -1,4 +0,0 @@
|
|||
bigint_go
|
||||
src
|
||||
*.wasm
|
||||
*.wat
|
||||
|
|
@ -1,37 +1,14 @@
|
|||
GitLink编程夏令营（GLCC），是在CCF中国计算机学会指导下，由CCF开源发展委员会（CCF ODC）举办的面向全国高校学生的暑期编程活动。活动将覆盖近千所高校，并联合各大开源基金会、开源企业、开源社区、开源专家，旨在鼓励青年学生通过参加真实的开源软件开发，提升自身技术能力，为开源社区输送优秀人才。为青年学生提供开放友好的交流平台，希望进一步推动国内开源社区的繁荣发展。
|
||||
# Pkg: waExample/main
|
||||
|
||||
凹语言开发组将作为指导组织参加本次的 [GitLink编程夏令营（GLCC 2023）](https://www.gitlink.org.cn/glcc/2023 "GitLink编程夏令营（GLCC 2023）")
|
||||
欢迎使用 Wa 语言!
|
||||
|
||||
## 目标
|
||||
使用国内开源社区创建的 [凹语言](https://wa-lang.org) 开发 天干地支码 编解码库。
|
||||
```
|
||||
+---+ +---+
|
||||
| o | | o |
|
||||
| +----+ |
|
||||
| |
|
||||
| \/\/ |
|
||||
| |
|
||||
+------------+
|
||||
```
|
||||
|
||||
## 难度
|
||||
中等
|
||||
|
||||
## 导师
|
||||
丁尔男
|
||||
电邮：<ending@wa-lang.org>
|
||||
|
||||
## 结果要求
|
||||
1. 60分线：提交**使用凹语言编写**的天干地支码编解码库；
|
||||
1. 80分线：除1外，额外提交测试用例、性能分析报告；
|
||||
1. 100分线：除2外，额外可编解码天干地支码的在线App（提交调用项目1成果）。
|
||||
开源创新，在很多时候体现为“使用有限的能力创造出更强的能力”，这也是我们对接题者的考察重点——既同学们需要使用凹语言的基本特性，从下至上的组合出本题所需的各种组件进而结题。
|
||||
|
||||
## 天干地支码简介
|
||||
天干地支纪年法，使用十天干（甲、乙、丙、丁、戊、己、庚、辛、壬、癸），十二地支（子、丑、寅、卯、辰、巳、午、未、申、酉、戌、亥）配合，每隔60年为一个轮回，可以将其类比为60进制（与二进制、十进制等对应）。天干地支码的本质既将输入的二进制数据串看作一个二进制的大整数B，将其转换为60进制的整数D，整数D的每一位对应一个干支码，最终所得的表达既为输入数据的天干地址码。每个干支码对应的10进制数如下表所示：
|
||||
|
||||
|||||||
|
||||
|---|---|---|---|---|---|
|
||||
|00 甲子|10 甲戌|20 甲申|30 甲午|40 甲辰|50 甲寅|
|
||||
|01 乙丑|11 乙亥|21 乙酉|31 乙未|41 乙巳|51 乙卯|
|
||||
|02 丙寅|12 丙子|22 丙戌|32 丙申|42 丙午|52 丙辰|
|
||||
|03 丁卯|13 丁丑|23 丁亥|33 丁酉|43 丁未|53 丁巳|
|
||||
|04 戊辰|14 戊寅|24 戊子|34 戊戌|44 戊申|54 戊午|
|
||||
|05 己巳|15 己卯|25 己丑|35 己亥|45 己酉|55 己未|
|
||||
|06 庚午|16 庚辰|26 庚寅|36 庚子|46 庚戌|56 庚申|
|
||||
|07 辛未|17 辛巳|27 辛卯|37 辛丑|47 辛亥|57 辛酉|
|
||||
|08 壬申|18 壬午|28 壬辰|38 壬寅|48 壬子|58 壬戌|
|
||||
|09 癸酉|19 癸未|29 癸巳|39 癸卯|49 癸丑|59 癸亥|
|
||||
|
||||
如：UTF8编码的“你好”，其对应的天干地支码为：“乙丑癸巳甲寅己亥丁卯甲申丁未甲午己巳”。
|
||||
|
|
|
|||
|
|
@ -0,0 +1,277 @@
|
|||
// Copyright 2009 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
// This file provides Go implementations of elementary multi-precision
|
||||
// arithmetic operations on word vectors. These have the suffix _g.
|
||||
// These are needed for platforms without assembly implementations of these routines.
|
||||
// This file also contains elementary operations that can be implemented
|
||||
// sufficiently efficiently in Go.
|
||||
|
||||
package big
|
||||
|
||||
import "math/bits"
|
||||
|
||||
// A Word represents a single digit of a multi-precision unsigned integer.
|
||||
type Word uint
|
||||
|
||||
const (
|
||||
_S = _W / 8 // word size in bytes
|
||||
|
||||
_W = bits.UintSize // word size in bits
|
||||
_B = 1 << _W // digit base
|
||||
_M = _B - 1 // digit mask
|
||||
)
|
||||
|
||||
// Many of the loops in this file are of the form
|
||||
// for i := 0; i < len(z) && i < len(x) && i < len(y); i++
|
||||
// i < len(z) is the real condition.
|
||||
// However, checking i < len(x) && i < len(y) as well is faster than
|
||||
// having the compiler do a bounds check in the body of the loop;
|
||||
// remarkably it is even faster than hoisting the bounds check
|
||||
// out of the loop, by doing something like
|
||||
// _, _ = x[len(z)-1], y[len(z)-1]
|
||||
// There are other ways to hoist the bounds check out of the loop,
|
||||
// but the compiler's BCE isn't powerful enough for them (yet?).
|
||||
// See the discussion in CL 164966.
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Elementary operations on words
|
||||
//
|
||||
// These operations are used by the vector operations below.
|
||||
|
||||
// z1<<_W + z0 = x*y
|
||||
func mulWW_g(x, y Word) (z1, z0 Word) {
|
||||
hi, lo := bits.Mul(uint(x), uint(y))
|
||||
return Word(hi), Word(lo)
|
||||
}
|
||||
|
||||
// z1<<_W + z0 = x*y + c
|
||||
func mulAddWWW_g(x, y, c Word) (z1, z0 Word) {
|
||||
hi, lo := bits.Mul(uint(x), uint(y))
|
||||
var cc uint
|
||||
lo, cc = bits.Add(lo, uint(c), 0)
|
||||
return Word(hi + cc), Word(lo)
|
||||
}
|
||||
|
||||
// nlz returns the number of leading zeros in x.
|
||||
// Wraps bits.LeadingZeros call for convenience.
|
||||
func nlz(x Word) uint {
|
||||
return uint(bits.LeadingZeros(uint(x)))
|
||||
}
|
||||
|
||||
// The resulting carry c is either 0 or 1.
|
||||
func addVV_g(z, x, y []Word) (c Word) {
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x) && i < len(y); i++ {
|
||||
zi, cc := bits.Add(uint(x[i]), uint(y[i]), uint(c))
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// The resulting carry c is either 0 or 1.
|
||||
func subVV_g(z, x, y []Word) (c Word) {
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x) && i < len(y); i++ {
|
||||
zi, cc := bits.Sub(uint(x[i]), uint(y[i]), uint(c))
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// The resulting carry c is either 0 or 1.
|
||||
func addVW_g(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
zi, cc := bits.Add(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// addVWlarge is addVW, but intended for large z.
|
||||
// The only difference is that we check on every iteration
|
||||
// whether we are done with carries,
|
||||
// and if so, switch to a much faster copy instead.
|
||||
// This is only a good idea for large z,
|
||||
// because the overhead of the check and the function call
|
||||
// outweigh the benefits when z is small.
|
||||
func addVWlarge(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
if c == 0 {
|
||||
copy(z[i:], x[i:])
|
||||
return
|
||||
}
|
||||
zi, cc := bits.Add(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
func subVW_g(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
zi, cc := bits.Sub(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// subVWlarge is to subVW as addVWlarge is to addVW.
|
||||
func subVWlarge(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
if c == 0 {
|
||||
copy(z[i:], x[i:])
|
||||
return
|
||||
}
|
||||
zi, cc := bits.Sub(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
func shlVU_g(z, x []Word, s uint) (c Word) {
|
||||
if s == 0 {
|
||||
copy(z, x)
|
||||
return
|
||||
}
|
||||
if len(z) == 0 {
|
||||
return
|
||||
}
|
||||
s &= _W - 1 // hint to the compiler that shifts by s don't need guard code
|
||||
ŝ := _W - s
|
||||
ŝ &= _W - 1 // ditto
|
||||
c = x[len(z)-1] >> ŝ
|
||||
for i := len(z) - 1; i > 0; i-- {
|
||||
z[i] = x[i]<<s | x[i-1]>>ŝ
|
||||
}
|
||||
z[0] = x[0] << s
|
||||
return
|
||||
}
|
||||
|
||||
func shrVU_g(z, x []Word, s uint) (c Word) {
|
||||
if s == 0 {
|
||||
copy(z, x)
|
||||
return
|
||||
}
|
||||
if len(z) == 0 {
|
||||
return
|
||||
}
|
||||
if len(x) != len(z) {
|
||||
// This is an invariant guaranteed by the caller.
|
||||
panic("len(x) != len(z)")
|
||||
}
|
||||
s &= _W - 1 // hint to the compiler that shifts by s don't need guard code
|
||||
ŝ := _W - s
|
||||
ŝ &= _W - 1 // ditto
|
||||
c = x[0] << ŝ
|
||||
for i := 1; i < len(z); i++ {
|
||||
z[i-1] = x[i-1]>>s | x[i]<<ŝ
|
||||
}
|
||||
z[len(z)-1] = x[len(z)-1] >> s
|
||||
return
|
||||
}
|
||||
|
||||
func mulAddVWW_g(z, x []Word, y, r Word) (c Word) {
|
||||
c = r
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
c, z[i] = mulAddWWW_g(x[i], y, c)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
func addMulVVW_g(z, x []Word, y Word) (c Word) {
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
z1, z0 := mulAddWWW_g(x[i], y, z[i])
|
||||
lo, cc := bits.Add(uint(z0), uint(c), 0)
|
||||
c, z[i] = Word(cc), Word(lo)
|
||||
c += z1
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// q = ( x1 << _W + x0 - r)/y. m = floor(( _B^2 - 1 ) / d - _B). Requiring x1<y.
|
||||
// An approximate reciprocal with a reference to "Improved Division by Invariant Integers
|
||||
// (IEEE Transactions on Computers, 11 Jun. 2010)"
|
||||
func divWW(x1, x0, y, m Word) (q, r Word) {
|
||||
s := nlz(y)
|
||||
if s != 0 {
|
||||
x1 = x1<<s | x0>>(_W-s)
|
||||
x0 <<= s
|
||||
y <<= s
|
||||
}
|
||||
d := uint(y)
|
||||
// We know that
|
||||
// m = ⎣(B^2-1)/d⎦-B
|
||||
// ⎣(B^2-1)/d⎦ = m+B
|
||||
// (B^2-1)/d = m+B+delta1 0 <= delta1 <= (d-1)/d
|
||||
// B^2/d = m+B+delta2 0 <= delta2 <= 1
|
||||
// The quotient we're trying to compute is
|
||||
// quotient = ⎣(x1*B+x0)/d⎦
|
||||
// = ⎣(x1*B*(B^2/d)+x0*(B^2/d))/B^2⎦
|
||||
// = ⎣(x1*B*(m+B+delta2)+x0*(m+B+delta2))/B^2⎦
|
||||
// = ⎣(x1*m+x1*B+x0)/B + x0*m/B^2 + delta2*(x1*B+x0)/B^2⎦
|
||||
// The latter two terms of this three-term sum are between 0 and 1.
|
||||
// So we can compute just the first term, and we will be low by at most 2.
|
||||
t1, t0 := bits.Mul(uint(m), uint(x1))
|
||||
_, c := bits.Add(t0, uint(x0), 0)
|
||||
t1, _ = bits.Add(t1, uint(x1), c)
|
||||
// The quotient is either t1, t1+1, or t1+2.
|
||||
// We'll try t1 and adjust if needed.
|
||||
qq := t1
|
||||
// compute remainder r=x-d*q.
|
||||
dq1, dq0 := bits.Mul(d, qq)
|
||||
r0, b := bits.Sub(uint(x0), dq0, 0)
|
||||
r1, _ := bits.Sub(uint(x1), dq1, b)
|
||||
// The remainder we just computed is bounded above by B+d:
|
||||
// r = x1*B + x0 - d*q.
|
||||
// = x1*B + x0 - d*⎣(x1*m+x1*B+x0)/B⎦
|
||||
// = x1*B + x0 - d*((x1*m+x1*B+x0)/B-alpha) 0 <= alpha < 1
|
||||
// = x1*B + x0 - x1*d/B*m - x1*d - x0*d/B + d*alpha
|
||||
// = x1*B + x0 - x1*d/B*⎣(B^2-1)/d-B⎦ - x1*d - x0*d/B + d*alpha
|
||||
// = x1*B + x0 - x1*d/B*⎣(B^2-1)/d-B⎦ - x1*d - x0*d/B + d*alpha
|
||||
// = x1*B + x0 - x1*d/B*((B^2-1)/d-B-beta) - x1*d - x0*d/B + d*alpha 0 <= beta < 1
|
||||
// = x1*B + x0 - x1*B + x1/B + x1*d + x1*d/B*beta - x1*d - x0*d/B + d*alpha
|
||||
// = x0 + x1/B + x1*d/B*beta - x0*d/B + d*alpha
|
||||
// = x0*(1-d/B) + x1*(1+d*beta)/B + d*alpha
|
||||
// < B*(1-d/B) + d*B/B + d because x0<B (and 1-d/B>0), x1<d, 1+d*beta<=B, alpha<1
|
||||
// = B - d + d + d
|
||||
// = B+d
|
||||
// So r1 can only be 0 or 1. If r1 is 1, then we know q was too small.
|
||||
// Add 1 to q and subtract d from r. That guarantees that r is <B, so
|
||||
// we no longer need to keep track of r1.
|
||||
if r1 != 0 {
|
||||
qq++
|
||||
r0 -= d
|
||||
}
|
||||
// If the remainder is still too large, increment q one more time.
|
||||
if r0 >= d {
|
||||
qq++
|
||||
r0 -= d
|
||||
}
|
||||
return Word(qq), Word(r0 >> s)
|
||||
}
|
||||
|
||||
// reciprocalWord return the reciprocal of the divisor. rec = floor(( _B^2 - 1 ) / u - _B). u = d1 << nlz(d1).
|
||||
func reciprocalWord(d1 Word) Word {
|
||||
u := uint(d1 << nlz(d1))
|
||||
x1 := ^u
|
||||
x0 := uint(_M)
|
||||
rec, _ := bits.Div(x1, x0, u) // (_B^2-1)/U-_B = (_B*(_M-C)+_M)/U
|
||||
return Word(rec)
|
||||
}
|
||||
|
|
@ -0,0 +1,54 @@
|
|||
// Copyright 2015 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
//go:build math_big_pure_go
|
||||
// +build math_big_pure_go
|
||||
|
||||
package big
|
||||
|
||||
func mulWW(x, y Word) (z1, z0 Word) {
|
||||
return mulWW_g(x, y)
|
||||
}
|
||||
|
||||
func addVV(z, x, y []Word) (c Word) {
|
||||
return addVV_g(z, x, y)
|
||||
}
|
||||
|
||||
func subVV(z, x, y []Word) (c Word) {
|
||||
return subVV_g(z, x, y)
|
||||
}
|
||||
|
||||
func addVW(z, x []Word, y Word) (c Word) {
|
||||
// TODO: remove indirect function call when golang.org/issue/30548 is fixed
|
||||
fn := addVW_g
|
||||
if len(z) > 32 {
|
||||
fn = addVWlarge
|
||||
}
|
||||
return fn(z, x, y)
|
||||
}
|
||||
|
||||
func subVW(z, x []Word, y Word) (c Word) {
|
||||
// TODO: remove indirect function call when golang.org/issue/30548 is fixed
|
||||
fn := subVW_g
|
||||
if len(z) > 32 {
|
||||
fn = subVWlarge
|
||||
}
|
||||
return fn(z, x, y)
|
||||
}
|
||||
|
||||
func shlVU(z, x []Word, s uint) (c Word) {
|
||||
return shlVU_g(z, x, s)
|
||||
}
|
||||
|
||||
func shrVU(z, x []Word, s uint) (c Word) {
|
||||
return shrVU_g(z, x, s)
|
||||
}
|
||||
|
||||
func mulAddVWW(z, x []Word, y, r Word) (c Word) {
|
||||
return mulAddVWW_g(z, x, y, r)
|
||||
}
|
||||
|
||||
func addMulVVW(z, x []Word, y Word) (c Word) {
|
||||
return addMulVVW_g(z, x, y)
|
||||
}
|
||||
|
|
@ -0,0 +1,697 @@
|
|||
// Copyright 2009 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
package big
|
||||
|
||||
import (
|
||||
"fmt"
|
||||
"internal/testenv"
|
||||
"math/bits"
|
||||
"math/rand"
|
||||
"strings"
|
||||
"testing"
|
||||
)
|
||||
|
||||
var isRaceBuilder = strings.HasSuffix(testenv.Builder(), "-race")
|
||||
|
||||
type funVV func(z, x, y []Word) (c Word)
|
||||
type argVV struct {
|
||||
z, x, y nat
|
||||
c Word
|
||||
}
|
||||
|
||||
var sumVV = []argVV{
|
||||
{},
|
||||
{nat{0}, nat{0}, nat{0}, 0},
|
||||
{nat{1}, nat{1}, nat{0}, 0},
|
||||
{nat{0}, nat{_M}, nat{1}, 1},
|
||||
{nat{80235}, nat{12345}, nat{67890}, 0},
|
||||
{nat{_M - 1}, nat{_M}, nat{_M}, 1},
|
||||
{nat{0, 0, 0, 0}, nat{_M, _M, _M, _M}, nat{1, 0, 0, 0}, 1},
|
||||
{nat{0, 0, 0, _M}, nat{_M, _M, _M, _M - 1}, nat{1, 0, 0, 0}, 0},
|
||||
{nat{0, 0, 0, 0}, nat{_M, 0, _M, 0}, nat{1, _M, 0, _M}, 1},
|
||||
}
|
||||
|
||||
func testFunVV(t *testing.T, msg string, f funVV, a argVV) {
|
||||
z := make(nat, len(a.z))
|
||||
c := f(z, a.x, a.y)
|
||||
for i, zi := range z {
|
||||
if zi != a.z[i] {
|
||||
t.Errorf("%s%+v\n\tgot z[%d] = %#x; want %#x", msg, a, i, zi, a.z[i])
|
||||
break
|
||||
}
|
||||
}
|
||||
if c != a.c {
|
||||
t.Errorf("%s%+v\n\tgot c = %#x; want %#x", msg, a, c, a.c)
|
||||
}
|
||||
}
|
||||
|
||||
func TestFunVV(t *testing.T) {
|
||||
for _, a := range sumVV {
|
||||
arg := a
|
||||
testFunVV(t, "addVV_g", addVV_g, arg)
|
||||
testFunVV(t, "addVV", addVV, arg)
|
||||
|
||||
arg = argVV{a.z, a.y, a.x, a.c}
|
||||
testFunVV(t, "addVV_g symmetric", addVV_g, arg)
|
||||
testFunVV(t, "addVV symmetric", addVV, arg)
|
||||
|
||||
arg = argVV{a.x, a.z, a.y, a.c}
|
||||
testFunVV(t, "subVV_g", subVV_g, arg)
|
||||
testFunVV(t, "subVV", subVV, arg)
|
||||
|
||||
arg = argVV{a.y, a.z, a.x, a.c}
|
||||
testFunVV(t, "subVV_g symmetric", subVV_g, arg)
|
||||
testFunVV(t, "subVV symmetric", subVV, arg)
|
||||
}
|
||||
}
|
||||
|
||||
// Always the same seed for reproducible results.
|
||||
var rnd = rand.New(rand.NewSource(0))
|
||||
|
||||
func rndW() Word {
|
||||
return Word(rnd.Int63()<<1 | rnd.Int63n(2))
|
||||
}
|
||||
|
||||
func rndV(n int) []Word {
|
||||
v := make([]Word, n)
|
||||
for i := range v {
|
||||
v[i] = rndW()
|
||||
}
|
||||
return v
|
||||
}
|
||||
|
||||
var benchSizes = []int{1, 2, 3, 4, 5, 1e1, 1e2, 1e3, 1e4, 1e5}
|
||||
|
||||
func BenchmarkAddVV(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
y := rndV(n)
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _W))
|
||||
for i := 0; i < b.N; i++ {
|
||||
addVV(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkSubVV(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
y := rndV(n)
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _W))
|
||||
for i := 0; i < b.N; i++ {
|
||||
subVV(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
type funVW func(z, x []Word, y Word) (c Word)
|
||||
type argVW struct {
|
||||
z, x nat
|
||||
y Word
|
||||
c Word
|
||||
}
|
||||
|
||||
var sumVW = []argVW{
|
||||
{},
|
||||
{nil, nil, 2, 2},
|
||||
{nat{0}, nat{0}, 0, 0},
|
||||
{nat{1}, nat{0}, 1, 0},
|
||||
{nat{1}, nat{1}, 0, 0},
|
||||
{nat{0}, nat{_M}, 1, 1},
|
||||
{nat{0, 0, 0, 0}, nat{_M, _M, _M, _M}, 1, 1},
|
||||
{nat{585}, nat{314}, 271, 0},
|
||||
}
|
||||
|
||||
var lshVW = []argVW{
|
||||
{},
|
||||
{nat{0}, nat{0}, 0, 0},
|
||||
{nat{0}, nat{0}, 1, 0},
|
||||
{nat{0}, nat{0}, 20, 0},
|
||||
|
||||
{nat{_M}, nat{_M}, 0, 0},
|
||||
{nat{_M << 1 & _M}, nat{_M}, 1, 1},
|
||||
{nat{_M << 20 & _M}, nat{_M}, 20, _M >> (_W - 20)},
|
||||
|
||||
{nat{_M, _M, _M}, nat{_M, _M, _M}, 0, 0},
|
||||
{nat{_M << 1 & _M, _M, _M}, nat{_M, _M, _M}, 1, 1},
|
||||
{nat{_M << 20 & _M, _M, _M}, nat{_M, _M, _M}, 20, _M >> (_W - 20)},
|
||||
}
|
||||
|
||||
var rshVW = []argVW{
|
||||
{},
|
||||
{nat{0}, nat{0}, 0, 0},
|
||||
{nat{0}, nat{0}, 1, 0},
|
||||
{nat{0}, nat{0}, 20, 0},
|
||||
|
||||
{nat{_M}, nat{_M}, 0, 0},
|
||||
{nat{_M >> 1}, nat{_M}, 1, _M << (_W - 1) & _M},
|
||||
{nat{_M >> 20}, nat{_M}, 20, _M << (_W - 20) & _M},
|
||||
|
||||
{nat{_M, _M, _M}, nat{_M, _M, _M}, 0, 0},
|
||||
{nat{_M, _M, _M >> 1}, nat{_M, _M, _M}, 1, _M << (_W - 1) & _M},
|
||||
{nat{_M, _M, _M >> 20}, nat{_M, _M, _M}, 20, _M << (_W - 20) & _M},
|
||||
}
|
||||
|
||||
func testFunVW(t *testing.T, msg string, f funVW, a argVW) {
|
||||
z := make(nat, len(a.z))
|
||||
c := f(z, a.x, a.y)
|
||||
for i, zi := range z {
|
||||
if zi != a.z[i] {
|
||||
t.Errorf("%s%+v\n\tgot z[%d] = %#x; want %#x", msg, a, i, zi, a.z[i])
|
||||
break
|
||||
}
|
||||
}
|
||||
if c != a.c {
|
||||
t.Errorf("%s%+v\n\tgot c = %#x; want %#x", msg, a, c, a.c)
|
||||
}
|
||||
}
|
||||
|
||||
func testFunVWext(t *testing.T, msg string, f funVW, f_g funVW, a argVW) {
|
||||
// using the result of addVW_g/subVW_g as golden
|
||||
z_g := make(nat, len(a.z))
|
||||
c_g := f_g(z_g, a.x, a.y)
|
||||
c := f(a.z, a.x, a.y)
|
||||
|
||||
for i, zi := range a.z {
|
||||
if zi != z_g[i] {
|
||||
t.Errorf("%s\n\tgot z[%d] = %#x; want %#x", msg, i, zi, z_g[i])
|
||||
break
|
||||
}
|
||||
}
|
||||
if c != c_g {
|
||||
t.Errorf("%s\n\tgot c = %#x; want %#x", msg, c, c_g)
|
||||
}
|
||||
}
|
||||
|
||||
func makeFunVW(f func(z, x []Word, s uint) (c Word)) funVW {
|
||||
return func(z, x []Word, s Word) (c Word) {
|
||||
return f(z, x, uint(s))
|
||||
}
|
||||
}
|
||||
|
||||
func TestFunVW(t *testing.T) {
|
||||
for _, a := range sumVW {
|
||||
arg := a
|
||||
testFunVW(t, "addVW_g", addVW_g, arg)
|
||||
testFunVW(t, "addVW", addVW, arg)
|
||||
|
||||
arg = argVW{a.x, a.z, a.y, a.c}
|
||||
testFunVW(t, "subVW_g", subVW_g, arg)
|
||||
testFunVW(t, "subVW", subVW, arg)
|
||||
}
|
||||
|
||||
shlVW_g := makeFunVW(shlVU_g)
|
||||
shlVW := makeFunVW(shlVU)
|
||||
for _, a := range lshVW {
|
||||
arg := a
|
||||
testFunVW(t, "shlVU_g", shlVW_g, arg)
|
||||
testFunVW(t, "shlVU", shlVW, arg)
|
||||
}
|
||||
|
||||
shrVW_g := makeFunVW(shrVU_g)
|
||||
shrVW := makeFunVW(shrVU)
|
||||
for _, a := range rshVW {
|
||||
arg := a
|
||||
testFunVW(t, "shrVU_g", shrVW_g, arg)
|
||||
testFunVW(t, "shrVU", shrVW, arg)
|
||||
}
|
||||
}
|
||||
|
||||
// Construct a vector comprising the same word, usually '0' or 'maximum uint'
|
||||
func makeWordVec(e Word, n int) []Word {
|
||||
v := make([]Word, n)
|
||||
for i := range v {
|
||||
v[i] = e
|
||||
}
|
||||
return v
|
||||
}
|
||||
|
||||
// Extended testing to addVW and subVW using various kinds of input data.
|
||||
// We utilize the results of addVW_g and subVW_g as golden reference to check
|
||||
// correctness.
|
||||
func TestFunVWExt(t *testing.T) {
|
||||
// 32 is the current threshold that triggers an optimized version of
|
||||
// calculation for large-sized vector, ensure we have sizes around it tested.
|
||||
var vwSizes = []int{0, 1, 3, 4, 5, 8, 9, 23, 31, 32, 33, 34, 35, 36, 50, 120}
|
||||
for _, n := range vwSizes {
|
||||
// vector of random numbers, using the result of addVW_g/subVW_g as golden
|
||||
x := rndV(n)
|
||||
y := rndW()
|
||||
z := make(nat, n)
|
||||
arg := argVW{z, x, y, 0}
|
||||
testFunVWext(t, "addVW, random inputs", addVW, addVW_g, arg)
|
||||
testFunVWext(t, "subVW, random inputs", subVW, subVW_g, arg)
|
||||
|
||||
// vector of random numbers, but make 'x' and 'z' share storage
|
||||
arg = argVW{x, x, y, 0}
|
||||
testFunVWext(t, "addVW, random inputs, sharing storage", addVW, addVW_g, arg)
|
||||
testFunVWext(t, "subVW, random inputs, sharing storage", subVW, subVW_g, arg)
|
||||
|
||||
// vector of maximum uint, to force carry flag set in each 'add'
|
||||
y = ^Word(0)
|
||||
x = makeWordVec(y, n)
|
||||
arg = argVW{z, x, y, 0}
|
||||
testFunVWext(t, "addVW, vector of max uint", addVW, addVW_g, arg)
|
||||
|
||||
// vector of '0', to force carry flag set in each 'sub'
|
||||
x = makeWordVec(0, n)
|
||||
arg = argVW{z, x, 1, 0}
|
||||
testFunVWext(t, "subVW, vector of zero", subVW, subVW_g, arg)
|
||||
}
|
||||
}
|
||||
|
||||
type argVU struct {
|
||||
d []Word // d is a Word slice, the input parameters x and z come from this array.
|
||||
l uint // l is the length of the input parameters x and z.
|
||||
xp uint // xp is the starting position of the input parameter x, x := d[xp:xp+l].
|
||||
zp uint // zp is the starting position of the input parameter z, z := d[zp:zp+l].
|
||||
s uint // s is the shift number.
|
||||
r []Word // r is the expected output result z.
|
||||
c Word // c is the expected return value.
|
||||
m string // message.
|
||||
}
|
||||
|
||||
var argshlVUIn = []Word{1, 2, 4, 8, 16, 32, 64, 0, 0, 0}
|
||||
var argshlVUr0 = []Word{1, 2, 4, 8, 16, 32, 64}
|
||||
var argshlVUr1 = []Word{2, 4, 8, 16, 32, 64, 128}
|
||||
var argshlVUrWm1 = []Word{1 << (_W - 1), 0, 1, 2, 4, 8, 16}
|
||||
|
||||
var argshlVU = []argVU{
|
||||
// test cases for shlVU
|
||||
{[]Word{1, _M, _M, _M, _M, _M, 3 << (_W - 2), 0}, 7, 0, 0, 1, []Word{2, _M - 1, _M, _M, _M, _M, 1<<(_W-1) + 1}, 1, "complete overlap of shlVU"},
|
||||
{[]Word{1, _M, _M, _M, _M, _M, 3 << (_W - 2), 0, 0, 0, 0}, 7, 0, 3, 1, []Word{2, _M - 1, _M, _M, _M, _M, 1<<(_W-1) + 1}, 1, "partial overlap by half of shlVU"},
|
||||
{[]Word{1, _M, _M, _M, _M, _M, 3 << (_W - 2), 0, 0, 0, 0, 0, 0, 0}, 7, 0, 6, 1, []Word{2, _M - 1, _M, _M, _M, _M, 1<<(_W-1) + 1}, 1, "partial overlap by 1 Word of shlVU"},
|
||||
{[]Word{1, _M, _M, _M, _M, _M, 3 << (_W - 2), 0, 0, 0, 0, 0, 0, 0, 0}, 7, 0, 7, 1, []Word{2, _M - 1, _M, _M, _M, _M, 1<<(_W-1) + 1}, 1, "no overlap of shlVU"},
|
||||
// additional test cases with shift values of 0, 1 and (_W-1)
|
||||
{argshlVUIn, 7, 0, 0, 0, argshlVUr0, 0, "complete overlap of shlVU and shift of 0"},
|
||||
{argshlVUIn, 7, 0, 0, 1, argshlVUr1, 0, "complete overlap of shlVU and shift of 1"},
|
||||
{argshlVUIn, 7, 0, 0, _W - 1, argshlVUrWm1, 32, "complete overlap of shlVU and shift of _W - 1"},
|
||||
{argshlVUIn, 7, 0, 1, 0, argshlVUr0, 0, "partial overlap by 6 Words of shlVU and shift of 0"},
|
||||
{argshlVUIn, 7, 0, 1, 1, argshlVUr1, 0, "partial overlap by 6 Words of shlVU and shift of 1"},
|
||||
{argshlVUIn, 7, 0, 1, _W - 1, argshlVUrWm1, 32, "partial overlap by 6 Words of shlVU and shift of _W - 1"},
|
||||
{argshlVUIn, 7, 0, 2, 0, argshlVUr0, 0, "partial overlap by 5 Words of shlVU and shift of 0"},
|
||||
{argshlVUIn, 7, 0, 2, 1, argshlVUr1, 0, "partial overlap by 5 Words of shlVU and shift of 1"},
|
||||
{argshlVUIn, 7, 0, 2, _W - 1, argshlVUrWm1, 32, "partial overlap by 5 Words of shlVU abd shift of _W - 1"},
|
||||
{argshlVUIn, 7, 0, 3, 0, argshlVUr0, 0, "partial overlap by 4 Words of shlVU and shift of 0"},
|
||||
{argshlVUIn, 7, 0, 3, 1, argshlVUr1, 0, "partial overlap by 4 Words of shlVU and shift of 1"},
|
||||
{argshlVUIn, 7, 0, 3, _W - 1, argshlVUrWm1, 32, "partial overlap by 4 Words of shlVU and shift of _W - 1"},
|
||||
}
|
||||
|
||||
var argshrVUIn = []Word{0, 0, 0, 1, 2, 4, 8, 16, 32, 64}
|
||||
var argshrVUr0 = []Word{1, 2, 4, 8, 16, 32, 64}
|
||||
var argshrVUr1 = []Word{0, 1, 2, 4, 8, 16, 32}
|
||||
var argshrVUrWm1 = []Word{4, 8, 16, 32, 64, 128, 0}
|
||||
|
||||
var argshrVU = []argVU{
|
||||
// test cases for shrVU
|
||||
{[]Word{0, 3, _M, _M, _M, _M, _M, 1 << (_W - 1)}, 7, 1, 1, 1, []Word{1<<(_W-1) + 1, _M, _M, _M, _M, _M >> 1, 1 << (_W - 2)}, 1 << (_W - 1), "complete overlap of shrVU"},
|
||||
{[]Word{0, 0, 0, 0, 3, _M, _M, _M, _M, _M, 1 << (_W - 1)}, 7, 4, 1, 1, []Word{1<<(_W-1) + 1, _M, _M, _M, _M, _M >> 1, 1 << (_W - 2)}, 1 << (_W - 1), "partial overlap by half of shrVU"},
|
||||
{[]Word{0, 0, 0, 0, 0, 0, 0, 3, _M, _M, _M, _M, _M, 1 << (_W - 1)}, 7, 7, 1, 1, []Word{1<<(_W-1) + 1, _M, _M, _M, _M, _M >> 1, 1 << (_W - 2)}, 1 << (_W - 1), "partial overlap by 1 Word of shrVU"},
|
||||
{[]Word{0, 0, 0, 0, 0, 0, 0, 0, 3, _M, _M, _M, _M, _M, 1 << (_W - 1)}, 7, 8, 1, 1, []Word{1<<(_W-1) + 1, _M, _M, _M, _M, _M >> 1, 1 << (_W - 2)}, 1 << (_W - 1), "no overlap of shrVU"},
|
||||
// additional test cases with shift values of 0, 1 and (_W-1)
|
||||
{argshrVUIn, 7, 3, 3, 0, argshrVUr0, 0, "complete overlap of shrVU and shift of 0"},
|
||||
{argshrVUIn, 7, 3, 3, 1, argshrVUr1, 1 << (_W - 1), "complete overlap of shrVU and shift of 1"},
|
||||
{argshrVUIn, 7, 3, 3, _W - 1, argshrVUrWm1, 2, "complete overlap of shrVU and shift of _W - 1"},
|
||||
{argshrVUIn, 7, 3, 2, 0, argshrVUr0, 0, "partial overlap by 6 Words of shrVU and shift of 0"},
|
||||
{argshrVUIn, 7, 3, 2, 1, argshrVUr1, 1 << (_W - 1), "partial overlap by 6 Words of shrVU and shift of 1"},
|
||||
{argshrVUIn, 7, 3, 2, _W - 1, argshrVUrWm1, 2, "partial overlap by 6 Words of shrVU and shift of _W - 1"},
|
||||
{argshrVUIn, 7, 3, 1, 0, argshrVUr0, 0, "partial overlap by 5 Words of shrVU and shift of 0"},
|
||||
{argshrVUIn, 7, 3, 1, 1, argshrVUr1, 1 << (_W - 1), "partial overlap by 5 Words of shrVU and shift of 1"},
|
||||
{argshrVUIn, 7, 3, 1, _W - 1, argshrVUrWm1, 2, "partial overlap by 5 Words of shrVU and shift of _W - 1"},
|
||||
{argshrVUIn, 7, 3, 0, 0, argshrVUr0, 0, "partial overlap by 4 Words of shrVU and shift of 0"},
|
||||
{argshrVUIn, 7, 3, 0, 1, argshrVUr1, 1 << (_W - 1), "partial overlap by 4 Words of shrVU and shift of 1"},
|
||||
{argshrVUIn, 7, 3, 0, _W - 1, argshrVUrWm1, 2, "partial overlap by 4 Words of shrVU and shift of _W - 1"},
|
||||
}
|
||||
|
||||
func testShiftFunc(t *testing.T, f func(z, x []Word, s uint) Word, a argVU) {
|
||||
// work on copy of a.d to preserve the original data.
|
||||
b := make([]Word, len(a.d))
|
||||
copy(b, a.d)
|
||||
z := b[a.zp : a.zp+a.l]
|
||||
x := b[a.xp : a.xp+a.l]
|
||||
c := f(z, x, a.s)
|
||||
for i, zi := range z {
|
||||
if zi != a.r[i] {
|
||||
t.Errorf("d := %v, %s(d[%d:%d], d[%d:%d], %d)\n\tgot z[%d] = %#x; want %#x", a.d, a.m, a.zp, a.zp+a.l, a.xp, a.xp+a.l, a.s, i, zi, a.r[i])
|
||||
break
|
||||
}
|
||||
}
|
||||
if c != a.c {
|
||||
t.Errorf("d := %v, %s(d[%d:%d], d[%d:%d], %d)\n\tgot c = %#x; want %#x", a.d, a.m, a.zp, a.zp+a.l, a.xp, a.xp+a.l, a.s, c, a.c)
|
||||
}
|
||||
}
|
||||
|
||||
func TestShiftOverlap(t *testing.T) {
|
||||
for _, a := range argshlVU {
|
||||
arg := a
|
||||
testShiftFunc(t, shlVU, arg)
|
||||
}
|
||||
|
||||
for _, a := range argshrVU {
|
||||
arg := a
|
||||
testShiftFunc(t, shrVU, arg)
|
||||
}
|
||||
}
|
||||
|
||||
func TestIssue31084(t *testing.T) {
|
||||
// compute 10^n via 5^n << n.
|
||||
const n = 165
|
||||
p := nat(nil).expNN(nat{5}, nat{n}, nil)
|
||||
p = p.shl(p, n)
|
||||
got := string(p.utoa(10))
|
||||
want := "1" + strings.Repeat("0", n)
|
||||
if got != want {
|
||||
t.Errorf("shl(%v, %v)\n\tgot %s\n\twant %s", p, n, got, want)
|
||||
}
|
||||
}
|
||||
|
||||
const issue42838Value = "159309191113245227702888039776771180559110455519261878607388585338616290151305816094308987472018268594098344692611135542392730712890625"
|
||||
|
||||
func TestIssue42838(t *testing.T) {
|
||||
const s = 192
|
||||
z, _, _, _ := nat(nil).scan(strings.NewReader(issue42838Value), 0, false)
|
||||
z = z.shl(z, s)
|
||||
got := string(z.utoa(10))
|
||||
want := "1" + strings.Repeat("0", s)
|
||||
if got != want {
|
||||
t.Errorf("shl(%v, %v)\n\tgot %s\n\twant %s", z, s, got, want)
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkAddVW(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
y := rndW()
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _S))
|
||||
for i := 0; i < b.N; i++ {
|
||||
addVW(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
// Benchmarking addVW using vector of maximum uint to force carry flag set
|
||||
func BenchmarkAddVWext(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
y := ^Word(0)
|
||||
x := makeWordVec(y, n)
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _S))
|
||||
for i := 0; i < b.N; i++ {
|
||||
addVW(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkSubVW(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
y := rndW()
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _S))
|
||||
for i := 0; i < b.N; i++ {
|
||||
subVW(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
// Benchmarking subVW using vector of zero to force carry flag set
|
||||
func BenchmarkSubVWext(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := makeWordVec(0, n)
|
||||
y := Word(1)
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _S))
|
||||
for i := 0; i < b.N; i++ {
|
||||
subVW(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
type funVWW func(z, x []Word, y, r Word) (c Word)
|
||||
type argVWW struct {
|
||||
z, x nat
|
||||
y, r Word
|
||||
c Word
|
||||
}
|
||||
|
||||
var prodVWW = []argVWW{
|
||||
{},
|
||||
{nat{0}, nat{0}, 0, 0, 0},
|
||||
{nat{991}, nat{0}, 0, 991, 0},
|
||||
{nat{0}, nat{_M}, 0, 0, 0},
|
||||
{nat{991}, nat{_M}, 0, 991, 0},
|
||||
{nat{0}, nat{0}, _M, 0, 0},
|
||||
{nat{991}, nat{0}, _M, 991, 0},
|
||||
{nat{1}, nat{1}, 1, 0, 0},
|
||||
{nat{992}, nat{1}, 1, 991, 0},
|
||||
{nat{22793}, nat{991}, 23, 0, 0},
|
||||
{nat{22800}, nat{991}, 23, 7, 0},
|
||||
{nat{0, 0, 0, 22793}, nat{0, 0, 0, 991}, 23, 0, 0},
|
||||
{nat{7, 0, 0, 22793}, nat{0, 0, 0, 991}, 23, 7, 0},
|
||||
{nat{0, 0, 0, 0}, nat{7893475, 7395495, 798547395, 68943}, 0, 0, 0},
|
||||
{nat{991, 0, 0, 0}, nat{7893475, 7395495, 798547395, 68943}, 0, 991, 0},
|
||||
{nat{0, 0, 0, 0}, nat{0, 0, 0, 0}, 894375984, 0, 0},
|
||||
{nat{991, 0, 0, 0}, nat{0, 0, 0, 0}, 894375984, 991, 0},
|
||||
{nat{_M << 1 & _M}, nat{_M}, 1 << 1, 0, _M >> (_W - 1)},
|
||||
{nat{_M<<1&_M + 1}, nat{_M}, 1 << 1, 1, _M >> (_W - 1)},
|
||||
{nat{_M << 7 & _M}, nat{_M}, 1 << 7, 0, _M >> (_W - 7)},
|
||||
{nat{_M<<7&_M + 1<<6}, nat{_M}, 1 << 7, 1 << 6, _M >> (_W - 7)},
|
||||
{nat{_M << 7 & _M, _M, _M, _M}, nat{_M, _M, _M, _M}, 1 << 7, 0, _M >> (_W - 7)},
|
||||
{nat{_M<<7&_M + 1<<6, _M, _M, _M}, nat{_M, _M, _M, _M}, 1 << 7, 1 << 6, _M >> (_W - 7)},
|
||||
}
|
||||
|
||||
func testFunVWW(t *testing.T, msg string, f funVWW, a argVWW) {
|
||||
z := make(nat, len(a.z))
|
||||
c := f(z, a.x, a.y, a.r)
|
||||
for i, zi := range z {
|
||||
if zi != a.z[i] {
|
||||
t.Errorf("%s%+v\n\tgot z[%d] = %#x; want %#x", msg, a, i, zi, a.z[i])
|
||||
break
|
||||
}
|
||||
}
|
||||
if c != a.c {
|
||||
t.Errorf("%s%+v\n\tgot c = %#x; want %#x", msg, a, c, a.c)
|
||||
}
|
||||
}
|
||||
|
||||
// TODO(gri) mulAddVWW and divWVW are symmetric operations but
|
||||
// their signature is not symmetric. Try to unify.
|
||||
|
||||
type funWVW func(z []Word, xn Word, x []Word, y Word) (r Word)
|
||||
type argWVW struct {
|
||||
z nat
|
||||
xn Word
|
||||
x nat
|
||||
y Word
|
||||
r Word
|
||||
}
|
||||
|
||||
func testFunWVW(t *testing.T, msg string, f funWVW, a argWVW) {
|
||||
z := make(nat, len(a.z))
|
||||
r := f(z, a.xn, a.x, a.y)
|
||||
for i, zi := range z {
|
||||
if zi != a.z[i] {
|
||||
t.Errorf("%s%+v\n\tgot z[%d] = %#x; want %#x", msg, a, i, zi, a.z[i])
|
||||
break
|
||||
}
|
||||
}
|
||||
if r != a.r {
|
||||
t.Errorf("%s%+v\n\tgot r = %#x; want %#x", msg, a, r, a.r)
|
||||
}
|
||||
}
|
||||
|
||||
func TestFunVWW(t *testing.T) {
|
||||
for _, a := range prodVWW {
|
||||
arg := a
|
||||
testFunVWW(t, "mulAddVWW_g", mulAddVWW_g, arg)
|
||||
testFunVWW(t, "mulAddVWW", mulAddVWW, arg)
|
||||
|
||||
if a.y != 0 && a.r < a.y {
|
||||
arg := argWVW{a.x, a.c, a.z, a.y, a.r}
|
||||
testFunWVW(t, "divWVW", divWVW, arg)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
var mulWWTests = []struct {
|
||||
x, y Word
|
||||
q, r Word
|
||||
}{
|
||||
{_M, _M, _M - 1, 1},
|
||||
// 32 bit only: {0xc47dfa8c, 50911, 0x98a4, 0x998587f4},
|
||||
}
|
||||
|
||||
func TestMulWW(t *testing.T) {
|
||||
for i, test := range mulWWTests {
|
||||
q, r := mulWW_g(test.x, test.y)
|
||||
if q != test.q || r != test.r {
|
||||
t.Errorf("#%d got (%x, %x) want (%x, %x)", i, q, r, test.q, test.r)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
var mulAddWWWTests = []struct {
|
||||
x, y, c Word
|
||||
q, r Word
|
||||
}{
|
||||
// TODO(agl): These will only work on 64-bit platforms.
|
||||
// {15064310297182388543, 0xe7df04d2d35d5d80, 13537600649892366549, 13644450054494335067, 10832252001440893781},
|
||||
// {15064310297182388543, 0xdab2f18048baa68d, 13644450054494335067, 12869334219691522700, 14233854684711418382},
|
||||
{_M, _M, 0, _M - 1, 1},
|
||||
{_M, _M, _M, _M, 0},
|
||||
}
|
||||
|
||||
func TestMulAddWWW(t *testing.T) {
|
||||
for i, test := range mulAddWWWTests {
|
||||
q, r := mulAddWWW_g(test.x, test.y, test.c)
|
||||
if q != test.q || r != test.r {
|
||||
t.Errorf("#%d got (%x, %x) want (%x, %x)", i, q, r, test.q, test.r)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
var divWWTests = []struct {
|
||||
x1, x0, y Word
|
||||
q, r Word
|
||||
}{
|
||||
{_M >> 1, 0, _M, _M >> 1, _M >> 1},
|
||||
{_M - (1 << (_W - 2)), _M, 3 << (_W - 2), _M, _M - (1 << (_W - 2))},
|
||||
}
|
||||
|
||||
const testsNumber = 1 << 16
|
||||
|
||||
func TestDivWW(t *testing.T) {
|
||||
i := 0
|
||||
for i, test := range divWWTests {
|
||||
rec := reciprocalWord(test.y)
|
||||
q, r := divWW(test.x1, test.x0, test.y, rec)
|
||||
if q != test.q || r != test.r {
|
||||
t.Errorf("#%d got (%x, %x) want (%x, %x)", i, q, r, test.q, test.r)
|
||||
}
|
||||
}
|
||||
//random tests
|
||||
for ; i < testsNumber; i++ {
|
||||
x1 := rndW()
|
||||
x0 := rndW()
|
||||
y := rndW()
|
||||
if x1 >= y {
|
||||
continue
|
||||
}
|
||||
rec := reciprocalWord(y)
|
||||
qGot, rGot := divWW(x1, x0, y, rec)
|
||||
qWant, rWant := bits.Div(uint(x1), uint(x0), uint(y))
|
||||
if uint(qGot) != qWant || uint(rGot) != rWant {
|
||||
t.Errorf("#%d got (%x, %x) want (%x, %x)", i, qGot, rGot, qWant, rWant)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkMulAddVWW(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
z := make([]Word, n+1)
|
||||
x := rndV(n)
|
||||
y := rndW()
|
||||
r := rndW()
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _W))
|
||||
for i := 0; i < b.N; i++ {
|
||||
mulAddVWW(z, x, y, r)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkAddMulVVW(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
y := rndW()
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _W))
|
||||
for i := 0; i < b.N; i++ {
|
||||
addMulVVW(z, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
func BenchmarkDivWVW(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
y := rndW()
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _W))
|
||||
for i := 0; i < b.N; i++ {
|
||||
divWVW(z, 0, x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkNonZeroShifts(b *testing.B) {
|
||||
for _, n := range benchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
x := rndV(n)
|
||||
s := uint(rand.Int63n(_W-2)) + 1 // avoid 0 and over-large shifts
|
||||
z := make([]Word, n)
|
||||
b.Run(fmt.Sprint(n), func(b *testing.B) {
|
||||
b.SetBytes(int64(n * _W))
|
||||
b.Run("shrVU", func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
_ = shrVU(z, x, s)
|
||||
}
|
||||
})
|
||||
b.Run("shlVU", func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
_ = shlVU(z, x, s)
|
||||
}
|
||||
})
|
||||
})
|
||||
}
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,816 @@
|
|||
// Copyright 2009 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
package big
|
||||
|
||||
import (
|
||||
"fmt"
|
||||
"runtime"
|
||||
"strings"
|
||||
"testing"
|
||||
)
|
||||
|
||||
var cmpTests = []struct {
|
||||
x, y nat
|
||||
r int
|
||||
}{
|
||||
{nil, nil, 0},
|
||||
{nil, nat(nil), 0},
|
||||
{nat(nil), nil, 0},
|
||||
{nat(nil), nat(nil), 0},
|
||||
{nat{0}, nat{0}, 0},
|
||||
{nat{0}, nat{1}, -1},
|
||||
{nat{1}, nat{0}, 1},
|
||||
{nat{1}, nat{1}, 0},
|
||||
{nat{0, _M}, nat{1}, 1},
|
||||
{nat{1}, nat{0, _M}, -1},
|
||||
{nat{1, _M}, nat{0, _M}, 1},
|
||||
{nat{0, _M}, nat{1, _M}, -1},
|
||||
{nat{16, 571956, 8794, 68}, nat{837, 9146, 1, 754489}, -1},
|
||||
{nat{34986, 41, 105, 1957}, nat{56, 7458, 104, 1957}, 1},
|
||||
}
|
||||
|
||||
func TestCmp(t *testing.T) {
|
||||
for i, a := range cmpTests {
|
||||
r := a.x.cmp(a.y)
|
||||
if r != a.r {
|
||||
t.Errorf("#%d got r = %v; want %v", i, r, a.r)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
type funNN func(z, x, y nat) nat
|
||||
type argNN struct {
|
||||
z, x, y nat
|
||||
}
|
||||
|
||||
var sumNN = []argNN{
|
||||
{},
|
||||
{nat{1}, nil, nat{1}},
|
||||
{nat{1111111110}, nat{123456789}, nat{987654321}},
|
||||
{nat{0, 0, 0, 1}, nil, nat{0, 0, 0, 1}},
|
||||
{nat{0, 0, 0, 1111111110}, nat{0, 0, 0, 123456789}, nat{0, 0, 0, 987654321}},
|
||||
{nat{0, 0, 0, 1}, nat{0, 0, _M}, nat{0, 0, 1}},
|
||||
}
|
||||
|
||||
var prodNN = []argNN{
|
||||
{},
|
||||
{nil, nil, nil},
|
||||
{nil, nat{991}, nil},
|
||||
{nat{991}, nat{991}, nat{1}},
|
||||
{nat{991 * 991}, nat{991}, nat{991}},
|
||||
{nat{0, 0, 991 * 991}, nat{0, 991}, nat{0, 991}},
|
||||
{nat{1 * 991, 2 * 991, 3 * 991, 4 * 991}, nat{1, 2, 3, 4}, nat{991}},
|
||||
{nat{4, 11, 20, 30, 20, 11, 4}, nat{1, 2, 3, 4}, nat{4, 3, 2, 1}},
|
||||
// 3^100 * 3^28 = 3^128
|
||||
{
|
||||
natFromString("11790184577738583171520872861412518665678211592275841109096961"),
|
||||
natFromString("515377520732011331036461129765621272702107522001"),
|
||||
natFromString("22876792454961"),
|
||||
},
|
||||
// z = 111....1 (70000 digits)
|
||||
// x = 10^(99*700) + ... + 10^1400 + 10^700 + 1
|
||||
// y = 111....1 (700 digits, larger than Karatsuba threshold on 32-bit and 64-bit)
|
||||
{
|
||||
natFromString(strings.Repeat("1", 70000)),
|
||||
natFromString("1" + strings.Repeat(strings.Repeat("0", 699)+"1", 99)),
|
||||
natFromString(strings.Repeat("1", 700)),
|
||||
},
|
||||
// z = 111....1 (20000 digits)
|
||||
// x = 10^10000 + 1
|
||||
// y = 111....1 (10000 digits)
|
||||
{
|
||||
natFromString(strings.Repeat("1", 20000)),
|
||||
natFromString("1" + strings.Repeat("0", 9999) + "1"),
|
||||
natFromString(strings.Repeat("1", 10000)),
|
||||
},
|
||||
}
|
||||
|
||||
func natFromString(s string) nat {
|
||||
x, _, _, err := nat(nil).scan(strings.NewReader(s), 0, false)
|
||||
if err != nil {
|
||||
panic(err)
|
||||
}
|
||||
return x
|
||||
}
|
||||
|
||||
func TestSet(t *testing.T) {
|
||||
for _, a := range sumNN {
|
||||
z := nat(nil).set(a.z)
|
||||
if z.cmp(a.z) != 0 {
|
||||
t.Errorf("got z = %v; want %v", z, a.z)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func testFunNN(t *testing.T, msg string, f funNN, a argNN) {
|
||||
z := f(nil, a.x, a.y)
|
||||
if z.cmp(a.z) != 0 {
|
||||
t.Errorf("%s%+v\n\tgot z = %v; want %v", msg, a, z, a.z)
|
||||
}
|
||||
}
|
||||
|
||||
func TestFunNN(t *testing.T) {
|
||||
for _, a := range sumNN {
|
||||
arg := a
|
||||
testFunNN(t, "add", nat.add, arg)
|
||||
|
||||
arg = argNN{a.z, a.y, a.x}
|
||||
testFunNN(t, "add symmetric", nat.add, arg)
|
||||
|
||||
arg = argNN{a.x, a.z, a.y}
|
||||
testFunNN(t, "sub", nat.sub, arg)
|
||||
|
||||
arg = argNN{a.y, a.z, a.x}
|
||||
testFunNN(t, "sub symmetric", nat.sub, arg)
|
||||
}
|
||||
|
||||
for _, a := range prodNN {
|
||||
arg := a
|
||||
testFunNN(t, "mul", nat.mul, arg)
|
||||
|
||||
arg = argNN{a.z, a.y, a.x}
|
||||
testFunNN(t, "mul symmetric", nat.mul, arg)
|
||||
}
|
||||
}
|
||||
|
||||
var mulRangesN = []struct {
|
||||
a, b uint64
|
||||
prod string
|
||||
}{
|
||||
{0, 0, "0"},
|
||||
{1, 1, "1"},
|
||||
{1, 2, "2"},
|
||||
{1, 3, "6"},
|
||||
{10, 10, "10"},
|
||||
{0, 100, "0"},
|
||||
{0, 1e9, "0"},
|
||||
{1, 0, "1"}, // empty range
|
||||
{100, 1, "1"}, // empty range
|
||||
{1, 10, "3628800"}, // 10!
|
||||
{1, 20, "2432902008176640000"}, // 20!
|
||||
{1, 100,
|
||||
"933262154439441526816992388562667004907159682643816214685929" +
|
||||
"638952175999932299156089414639761565182862536979208272237582" +
|
||||
"51185210916864000000000000000000000000", // 100!
|
||||
},
|
||||
}
|
||||
|
||||
func TestMulRangeN(t *testing.T) {
|
||||
for i, r := range mulRangesN {
|
||||
prod := string(nat(nil).mulRange(r.a, r.b).utoa(10))
|
||||
if prod != r.prod {
|
||||
t.Errorf("#%d: got %s; want %s", i, prod, r.prod)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// allocBytes returns the number of bytes allocated by invoking f.
|
||||
func allocBytes(f func()) uint64 {
|
||||
var stats runtime.MemStats
|
||||
runtime.ReadMemStats(&stats)
|
||||
t := stats.TotalAlloc
|
||||
f()
|
||||
runtime.ReadMemStats(&stats)
|
||||
return stats.TotalAlloc - t
|
||||
}
|
||||
|
||||
// TestMulUnbalanced tests that multiplying numbers of different lengths
|
||||
// does not cause deep recursion and in turn allocate too much memory.
|
||||
// Test case for issue 3807.
|
||||
func TestMulUnbalanced(t *testing.T) {
|
||||
defer runtime.GOMAXPROCS(runtime.GOMAXPROCS(1))
|
||||
x := rndNat(50000)
|
||||
y := rndNat(40)
|
||||
allocSize := allocBytes(func() {
|
||||
nat(nil).mul(x, y)
|
||||
})
|
||||
inputSize := uint64(len(x)+len(y)) * _S
|
||||
if ratio := allocSize / uint64(inputSize); ratio > 10 {
|
||||
t.Errorf("multiplication uses too much memory (%d > %d times the size of inputs)", allocSize, ratio)
|
||||
}
|
||||
}
|
||||
|
||||
// rndNat returns a random nat value >= 0 of (usually) n words in length.
|
||||
// In extremely unlikely cases it may be smaller than n words if the top-
|
||||
// most words are 0.
|
||||
func rndNat(n int) nat {
|
||||
return nat(rndV(n)).norm()
|
||||
}
|
||||
|
||||
// rndNat1 is like rndNat but the result is guaranteed to be > 0.
|
||||
func rndNat1(n int) nat {
|
||||
x := nat(rndV(n)).norm()
|
||||
if len(x) == 0 {
|
||||
x.setWord(1)
|
||||
}
|
||||
return x
|
||||
}
|
||||
|
||||
func BenchmarkMul(b *testing.B) {
|
||||
mulx := rndNat(1e4)
|
||||
muly := rndNat(1e4)
|
||||
b.ResetTimer()
|
||||
for i := 0; i < b.N; i++ {
|
||||
var z nat
|
||||
z.mul(mulx, muly)
|
||||
}
|
||||
}
|
||||
|
||||
func benchmarkNatMul(b *testing.B, nwords int) {
|
||||
x := rndNat(nwords)
|
||||
y := rndNat(nwords)
|
||||
var z nat
|
||||
b.ResetTimer()
|
||||
for i := 0; i < b.N; i++ {
|
||||
z.mul(x, y)
|
||||
}
|
||||
}
|
||||
|
||||
var mulBenchSizes = []int{10, 100, 1000, 10000, 100000}
|
||||
|
||||
func BenchmarkNatMul(b *testing.B) {
|
||||
for _, n := range mulBenchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
b.Run(fmt.Sprintf("%d", n), func(b *testing.B) {
|
||||
benchmarkNatMul(b, n)
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func TestNLZ(t *testing.T) {
|
||||
var x Word = _B >> 1
|
||||
for i := 0; i <= _W; i++ {
|
||||
if int(nlz(x)) != i {
|
||||
t.Errorf("failed at %x: got %d want %d", x, nlz(x), i)
|
||||
}
|
||||
x >>= 1
|
||||
}
|
||||
}
|
||||
|
||||
type shiftTest struct {
|
||||
in nat
|
||||
shift uint
|
||||
out nat
|
||||
}
|
||||
|
||||
var leftShiftTests = []shiftTest{
|
||||
{nil, 0, nil},
|
||||
{nil, 1, nil},
|
||||
{natOne, 0, natOne},
|
||||
{natOne, 1, natTwo},
|
||||
{nat{1 << (_W - 1)}, 1, nat{0}},
|
||||
{nat{1 << (_W - 1), 0}, 1, nat{0, 1}},
|
||||
}
|
||||
|
||||
func TestShiftLeft(t *testing.T) {
|
||||
for i, test := range leftShiftTests {
|
||||
var z nat
|
||||
z = z.shl(test.in, test.shift)
|
||||
for j, d := range test.out {
|
||||
if j >= len(z) || z[j] != d {
|
||||
t.Errorf("#%d: got: %v want: %v", i, z, test.out)
|
||||
break
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
var rightShiftTests = []shiftTest{
|
||||
{nil, 0, nil},
|
||||
{nil, 1, nil},
|
||||
{natOne, 0, natOne},
|
||||
{natOne, 1, nil},
|
||||
{natTwo, 1, natOne},
|
||||
{nat{0, 1}, 1, nat{1 << (_W - 1)}},
|
||||
{nat{2, 1, 1}, 1, nat{1<<(_W-1) + 1, 1 << (_W - 1)}},
|
||||
}
|
||||
|
||||
func TestShiftRight(t *testing.T) {
|
||||
for i, test := range rightShiftTests {
|
||||
var z nat
|
||||
z = z.shr(test.in, test.shift)
|
||||
for j, d := range test.out {
|
||||
if j >= len(z) || z[j] != d {
|
||||
t.Errorf("#%d: got: %v want: %v", i, z, test.out)
|
||||
break
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkZeroShifts(b *testing.B) {
|
||||
x := rndNat(800)
|
||||
|
||||
b.Run("Shl", func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
var z nat
|
||||
z.shl(x, 0)
|
||||
}
|
||||
})
|
||||
b.Run("ShlSame", func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
x.shl(x, 0)
|
||||
}
|
||||
})
|
||||
|
||||
b.Run("Shr", func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
var z nat
|
||||
z.shr(x, 0)
|
||||
}
|
||||
})
|
||||
b.Run("ShrSame", func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
x.shr(x, 0)
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
type modWTest struct {
|
||||
in string
|
||||
dividend string
|
||||
out string
|
||||
}
|
||||
|
||||
var modWTests32 = []modWTest{
|
||||
{"23492635982634928349238759823742", "252341", "220170"},
|
||||
}
|
||||
|
||||
var modWTests64 = []modWTest{
|
||||
{"6527895462947293856291561095690465243862946", "524326975699234", "375066989628668"},
|
||||
}
|
||||
|
||||
func runModWTests(t *testing.T, tests []modWTest) {
|
||||
for i, test := range tests {
|
||||
in, _ := new(Int).SetString(test.in, 10)
|
||||
d, _ := new(Int).SetString(test.dividend, 10)
|
||||
out, _ := new(Int).SetString(test.out, 10)
|
||||
|
||||
r := in.abs.modW(d.abs[0])
|
||||
if r != out.abs[0] {
|
||||
t.Errorf("#%d failed: got %d want %s", i, r, out)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestModW(t *testing.T) {
|
||||
if _W >= 32 {
|
||||
runModWTests(t, modWTests32)
|
||||
}
|
||||
if _W >= 64 {
|
||||
runModWTests(t, modWTests64)
|
||||
}
|
||||
}
|
||||
|
||||
var montgomeryTests = []struct {
|
||||
x, y, m string
|
||||
k0 uint64
|
||||
out32, out64 string
|
||||
}{
|
||||
{
|
||||
"0xffffffffffffffffffffffffffffffffffffffffffffffffe",
|
||||
"0xffffffffffffffffffffffffffffffffffffffffffffffffe",
|
||||
"0xfffffffffffffffffffffffffffffffffffffffffffffffff",
|
||||
1,
|
||||
"0x1000000000000000000000000000000000000000000",
|
||||
"0x10000000000000000000000000000000000",
|
||||
},
|
||||
{
|
||||
"0x000000000ffffff5",
|
||||
"0x000000000ffffff0",
|
||||
"0x0000000010000001",
|
||||
0xff0000000fffffff,
|
||||
"0x000000000bfffff4",
|
||||
"0x0000000003400001",
|
||||
},
|
||||
{
|
||||
"0x0000000080000000",
|
||||
"0x00000000ffffffff",
|
||||
"0x1000000000000001",
|
||||
0xfffffffffffffff,
|
||||
"0x0800000008000001",
|
||||
"0x0800000008000001",
|
||||
},
|
||||
{
|
||||
"0x0000000080000000",
|
||||
"0x0000000080000000",
|
||||
"0xffffffff00000001",
|
||||
0xfffffffeffffffff,
|
||||
"0xbfffffff40000001",
|
||||
"0xbfffffff40000001",
|
||||
},
|
||||
{
|
||||
"0x0000000080000000",
|
||||
"0x0000000080000000",
|
||||
"0x00ffffff00000001",
|
||||
0xfffffeffffffff,
|
||||
"0xbfffff40000001",
|
||||
"0xbfffff40000001",
|
||||
},
|
||||
{
|
||||
"0x0000000080000000",
|
||||
"0x0000000080000000",
|
||||
"0x0000ffff00000001",
|
||||
0xfffeffffffff,
|
||||
"0xbfff40000001",
|
||||
"0xbfff40000001",
|
||||
},
|
||||
{
|
||||
"0x3321ffffffffffffffffffffffffffff00000000000022222623333333332bbbb888c0",
|
||||
"0x3321ffffffffffffffffffffffffffff00000000000022222623333333332bbbb888c0",
|
||||
"0x33377fffffffffffffffffffffffffffffffffffffffffffff0000000000022222eee1",
|
||||
0xdecc8f1249812adf,
|
||||
"0x04eb0e11d72329dc0915f86784820fc403275bf2f6620a20e0dd344c5cd0875e50deb5",
|
||||
"0x0d7144739a7d8e11d72329dc0915f86784820fc403275bf2f61ed96f35dd34dbb3d6a0",
|
||||
},
|
||||
{
|
||||
"0x10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffffffffffffff00000000000022222223333333333444444444",
|
||||
"0x10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffffffffffffff999999999999999aaabbbbbbbbcccccccccccc",
|
||||
"0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff33377fffffffffffffffffffffffffffffffffffffffffffff0000000000022222eee1",
|
||||
0xdecc8f1249812adf,
|
||||
"0x5c0d52f451aec609b15da8e5e5626c4eaa88723bdeac9d25ca9b961269400410ca208a16af9c2fb07d7a11c7772cba02c22f9711078d51a3797eb18e691295293284d988e349fa6deba46b25a4ecd9f715",
|
||||
"0x92fcad4b5c0d52f451aec609b15da8e5e5626c4eaa88723bdeac9d25ca9b961269400410ca208a16af9c2fb07d799c32fe2f3cc5422f9711078d51a3797eb18e691295293284d8f5e69caf6decddfe1df6",
|
||||
},
|
||||
}
|
||||
|
||||
func TestMontgomery(t *testing.T) {
|
||||
one := NewInt(1)
|
||||
_B := new(Int).Lsh(one, _W)
|
||||
for i, test := range montgomeryTests {
|
||||
x := natFromString(test.x)
|
||||
y := natFromString(test.y)
|
||||
m := natFromString(test.m)
|
||||
for len(x) < len(m) {
|
||||
x = append(x, 0)
|
||||
}
|
||||
for len(y) < len(m) {
|
||||
y = append(y, 0)
|
||||
}
|
||||
|
||||
if x.cmp(m) > 0 {
|
||||
_, r := nat(nil).div(nil, x, m)
|
||||
t.Errorf("#%d: x > m (0x%s > 0x%s; use 0x%s)", i, x.utoa(16), m.utoa(16), r.utoa(16))
|
||||
}
|
||||
if y.cmp(m) > 0 {
|
||||
_, r := nat(nil).div(nil, x, m)
|
||||
t.Errorf("#%d: y > m (0x%s > 0x%s; use 0x%s)", i, y.utoa(16), m.utoa(16), r.utoa(16))
|
||||
}
|
||||
|
||||
var out nat
|
||||
if _W == 32 {
|
||||
out = natFromString(test.out32)
|
||||
} else {
|
||||
out = natFromString(test.out64)
|
||||
}
|
||||
|
||||
// t.Logf("#%d: len=%d\n", i, len(m))
|
||||
|
||||
// check output in table
|
||||
xi := &Int{abs: x}
|
||||
yi := &Int{abs: y}
|
||||
mi := &Int{abs: m}
|
||||
p := new(Int).Mod(new(Int).Mul(xi, new(Int).Mul(yi, new(Int).ModInverse(new(Int).Lsh(one, uint(len(m))*_W), mi))), mi)
|
||||
if out.cmp(p.abs.norm()) != 0 {
|
||||
t.Errorf("#%d: out in table=0x%s, computed=0x%s", i, out.utoa(16), p.abs.norm().utoa(16))
|
||||
}
|
||||
|
||||
// check k0 in table
|
||||
k := new(Int).Mod(&Int{abs: m}, _B)
|
||||
k = new(Int).Sub(_B, k)
|
||||
k = new(Int).Mod(k, _B)
|
||||
k0 := Word(new(Int).ModInverse(k, _B).Uint64())
|
||||
if k0 != Word(test.k0) {
|
||||
t.Errorf("#%d: k0 in table=%#x, computed=%#x\n", i, test.k0, k0)
|
||||
}
|
||||
|
||||
// check montgomery with correct k0 produces correct output
|
||||
z := nat(nil).montgomery(x, y, m, k0, len(m))
|
||||
z = z.norm()
|
||||
if z.cmp(out) != 0 {
|
||||
t.Errorf("#%d: got 0x%s want 0x%s", i, z.utoa(16), out.utoa(16))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
var expNNTests = []struct {
|
||||
x, y, m string
|
||||
out string
|
||||
}{
|
||||
{"0", "0", "0", "1"},
|
||||
{"0", "0", "1", "0"},
|
||||
{"1", "1", "1", "0"},
|
||||
{"2", "1", "1", "0"},
|
||||
{"2", "2", "1", "0"},
|
||||
{"10", "100000000000", "1", "0"},
|
||||
{"0x8000000000000000", "2", "", "0x40000000000000000000000000000000"},
|
||||
{"0x8000000000000000", "2", "6719", "4944"},
|
||||
{"0x8000000000000000", "3", "6719", "5447"},
|
||||
{"0x8000000000000000", "1000", "6719", "1603"},
|
||||
{"0x8000000000000000", "1000000", "6719", "3199"},
|
||||
{
|
||||
"2938462938472983472983659726349017249287491026512746239764525612965293865296239471239874193284792387498274256129746192347",
|
||||
"298472983472983471903246121093472394872319615612417471234712061",
|
||||
"29834729834729834729347290846729561262544958723956495615629569234729836259263598127342374289365912465901365498236492183464",
|
||||
"23537740700184054162508175125554701713153216681790245129157191391322321508055833908509185839069455749219131480588829346291",
|
||||
},
|
||||
{
|
||||
"11521922904531591643048817447554701904414021819823889996244743037378330903763518501116638828335352811871131385129455853417360623007349090150042001944696604737499160174391019030572483602867266711107136838523916077674888297896995042968746762200926853379",
|
||||
"426343618817810911523",
|
||||
"444747819283133684179",
|
||||
"42",
|
||||
},
|
||||
}
|
||||
|
||||
func TestExpNN(t *testing.T) {
|
||||
for i, test := range expNNTests {
|
||||
x := natFromString(test.x)
|
||||
y := natFromString(test.y)
|
||||
out := natFromString(test.out)
|
||||
|
||||
var m nat
|
||||
if len(test.m) > 0 {
|
||||
m = natFromString(test.m)
|
||||
}
|
||||
|
||||
z := nat(nil).expNN(x, y, m)
|
||||
if z.cmp(out) != 0 {
|
||||
t.Errorf("#%d got %s want %s", i, z.utoa(10), out.utoa(10))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkExp3Power(b *testing.B) {
|
||||
const x = 3
|
||||
for _, y := range []Word{
|
||||
0x10, 0x40, 0x100, 0x400, 0x1000, 0x4000, 0x10000, 0x40000, 0x100000, 0x400000,
|
||||
} {
|
||||
b.Run(fmt.Sprintf("%#x", y), func(b *testing.B) {
|
||||
var z nat
|
||||
for i := 0; i < b.N; i++ {
|
||||
z.expWW(x, y)
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func fibo(n int) nat {
|
||||
switch n {
|
||||
case 0:
|
||||
return nil
|
||||
case 1:
|
||||
return nat{1}
|
||||
}
|
||||
f0 := fibo(0)
|
||||
f1 := fibo(1)
|
||||
var f2 nat
|
||||
for i := 1; i < n; i++ {
|
||||
f2 = f2.add(f0, f1)
|
||||
f0, f1, f2 = f1, f2, f0
|
||||
}
|
||||
return f1
|
||||
}
|
||||
|
||||
var fiboNums = []string{
|
||||
"0",
|
||||
"55",
|
||||
"6765",
|
||||
"832040",
|
||||
"102334155",
|
||||
"12586269025",
|
||||
"1548008755920",
|
||||
"190392490709135",
|
||||
"23416728348467685",
|
||||
"2880067194370816120",
|
||||
"354224848179261915075",
|
||||
}
|
||||
|
||||
func TestFibo(t *testing.T) {
|
||||
for i, want := range fiboNums {
|
||||
n := i * 10
|
||||
got := string(fibo(n).utoa(10))
|
||||
if got != want {
|
||||
t.Errorf("fibo(%d) failed: got %s want %s", n, got, want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkFibo(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
fibo(1e0)
|
||||
fibo(1e1)
|
||||
fibo(1e2)
|
||||
fibo(1e3)
|
||||
fibo(1e4)
|
||||
fibo(1e5)
|
||||
}
|
||||
}
|
||||
|
||||
var bitTests = []struct {
|
||||
x string
|
||||
i uint
|
||||
want uint
|
||||
}{
|
||||
{"0", 0, 0},
|
||||
{"0", 1, 0},
|
||||
{"0", 1000, 0},
|
||||
|
||||
{"0x1", 0, 1},
|
||||
{"0x10", 0, 0},
|
||||
{"0x10", 3, 0},
|
||||
{"0x10", 4, 1},
|
||||
{"0x10", 5, 0},
|
||||
|
||||
{"0x8000000000000000", 62, 0},
|
||||
{"0x8000000000000000", 63, 1},
|
||||
{"0x8000000000000000", 64, 0},
|
||||
|
||||
{"0x3" + strings.Repeat("0", 32), 127, 0},
|
||||
{"0x3" + strings.Repeat("0", 32), 128, 1},
|
||||
{"0x3" + strings.Repeat("0", 32), 129, 1},
|
||||
{"0x3" + strings.Repeat("0", 32), 130, 0},
|
||||
}
|
||||
|
||||
func TestBit(t *testing.T) {
|
||||
for i, test := range bitTests {
|
||||
x := natFromString(test.x)
|
||||
if got := x.bit(test.i); got != test.want {
|
||||
t.Errorf("#%d: %s.bit(%d) = %v; want %v", i, test.x, test.i, got, test.want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
var stickyTests = []struct {
|
||||
x string
|
||||
i uint
|
||||
want uint
|
||||
}{
|
||||
{"0", 0, 0},
|
||||
{"0", 1, 0},
|
||||
{"0", 1000, 0},
|
||||
|
||||
{"0x1", 0, 0},
|
||||
{"0x1", 1, 1},
|
||||
|
||||
{"0x1350", 0, 0},
|
||||
{"0x1350", 4, 0},
|
||||
{"0x1350", 5, 1},
|
||||
|
||||
{"0x8000000000000000", 63, 0},
|
||||
{"0x8000000000000000", 64, 1},
|
||||
|
||||
{"0x1" + strings.Repeat("0", 100), 400, 0},
|
||||
{"0x1" + strings.Repeat("0", 100), 401, 1},
|
||||
}
|
||||
|
||||
func TestSticky(t *testing.T) {
|
||||
for i, test := range stickyTests {
|
||||
x := natFromString(test.x)
|
||||
if got := x.sticky(test.i); got != test.want {
|
||||
t.Errorf("#%d: %s.sticky(%d) = %v; want %v", i, test.x, test.i, got, test.want)
|
||||
}
|
||||
if test.want == 1 {
|
||||
// all subsequent i's should also return 1
|
||||
for d := uint(1); d <= 3; d++ {
|
||||
if got := x.sticky(test.i + d); got != 1 {
|
||||
t.Errorf("#%d: %s.sticky(%d) = %v; want %v", i, test.x, test.i+d, got, 1)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func testSqr(t *testing.T, x nat) {
|
||||
got := make(nat, 2*len(x))
|
||||
want := make(nat, 2*len(x))
|
||||
got = got.sqr(x)
|
||||
want = want.mul(x, x)
|
||||
if got.cmp(want) != 0 {
|
||||
t.Errorf("basicSqr(%v), got %v, want %v", x, got, want)
|
||||
}
|
||||
}
|
||||
|
||||
func TestSqr(t *testing.T) {
|
||||
for _, a := range prodNN {
|
||||
if a.x != nil {
|
||||
testSqr(t, a.x)
|
||||
}
|
||||
if a.y != nil {
|
||||
testSqr(t, a.y)
|
||||
}
|
||||
if a.z != nil {
|
||||
testSqr(t, a.z)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func benchmarkNatSqr(b *testing.B, nwords int) {
|
||||
x := rndNat(nwords)
|
||||
var z nat
|
||||
b.ResetTimer()
|
||||
for i := 0; i < b.N; i++ {
|
||||
z.sqr(x)
|
||||
}
|
||||
}
|
||||
|
||||
var sqrBenchSizes = []int{
|
||||
1, 2, 3, 5, 8, 10, 20, 30, 50, 80,
|
||||
100, 200, 300, 500, 800,
|
||||
1000, 10000, 100000,
|
||||
}
|
||||
|
||||
func BenchmarkNatSqr(b *testing.B) {
|
||||
for _, n := range sqrBenchSizes {
|
||||
if isRaceBuilder && n > 1e3 {
|
||||
continue
|
||||
}
|
||||
b.Run(fmt.Sprintf("%d", n), func(b *testing.B) {
|
||||
benchmarkNatSqr(b, n)
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkNatSetBytes(b *testing.B) {
|
||||
const maxLength = 128
|
||||
lengths := []int{
|
||||
// No remainder:
|
||||
8, 24, maxLength,
|
||||
// With remainder:
|
||||
7, 23, maxLength - 1,
|
||||
}
|
||||
n := make(nat, maxLength/_W) // ensure n doesn't need to grow during the test
|
||||
buf := make([]byte, maxLength)
|
||||
for _, l := range lengths {
|
||||
b.Run(fmt.Sprint(l), func(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
n.setBytes(buf[:l])
|
||||
}
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
func TestNatDiv(t *testing.T) {
|
||||
sizes := []int{
|
||||
1, 2, 5, 8, 15, 25, 40, 65, 100,
|
||||
200, 500, 800, 1500, 2500, 4000, 6500, 10000,
|
||||
}
|
||||
for _, i := range sizes {
|
||||
for _, j := range sizes {
|
||||
a := rndNat1(i)
|
||||
b := rndNat1(j)
|
||||
// the test requires b >= 2
|
||||
if len(b) == 1 && b[0] == 1 {
|
||||
b[0] = 2
|
||||
}
|
||||
// choose a remainder c < b
|
||||
c := rndNat1(len(b))
|
||||
if len(c) == len(b) && c[len(c)-1] >= b[len(b)-1] {
|
||||
c[len(c)-1] = 0
|
||||
c = c.norm()
|
||||
}
|
||||
// compute x = a*b+c
|
||||
x := nat(nil).mul(a, b)
|
||||
x = x.add(x, c)
|
||||
|
||||
var q, r nat
|
||||
q, r = q.div(r, x, b)
|
||||
if q.cmp(a) != 0 {
|
||||
t.Fatalf("wrong quotient: got %s; want %s for %s/%s", q.utoa(10), a.utoa(10), x.utoa(10), b.utoa(10))
|
||||
}
|
||||
if r.cmp(c) != 0 {
|
||||
t.Fatalf("wrong remainder: got %s; want %s for %s/%s", r.utoa(10), c.utoa(10), x.utoa(10), b.utoa(10))
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// TestIssue37499 triggers the edge case of divBasic where
|
||||
// the inaccurate estimate of the first word's quotient
|
||||
// happens at the very beginning of the loop.
|
||||
func TestIssue37499(t *testing.T) {
|
||||
// Choose u and v such that v is slightly larger than u >> N.
|
||||
// This tricks divBasic into choosing 1 as the first word
|
||||
// of the quotient. This works in both 32-bit and 64-bit settings.
|
||||
u := natFromString("0x2b6c385a05be027f5c22005b63c42a1165b79ff510e1706b39f8489c1d28e57bb5ba4ef9fd9387a3e344402c0a453381")
|
||||
v := natFromString("0x2b6c385a05be027f5c22005b63c42a1165b79ff510e1706c")
|
||||
|
||||
q := nat(nil).make(8)
|
||||
q.divBasic(u, v)
|
||||
q = q.norm()
|
||||
if s := string(q.utoa(16)); s != "fffffffffffffffffffffffffffffffffffffffffffffffb" {
|
||||
t.Fatalf("incorrect quotient: %s", s)
|
||||
}
|
||||
}
|
||||
|
||||
// TestIssue42552 triggers an edge case of recursive division
|
||||
// where the first division loop is never entered, and correcting
|
||||
// the remainder takes exactly two iterations in the final loop.
|
||||
func TestIssue42552(t *testing.T) {
|
||||
u := natFromString("0xc23b166884c3869092a520eceedeced2b00847bd256c9cf3b2c5e2227c15bd5e6ee7ef8a2f49236ad0eedf2c8a3b453cf6e0706f64285c526b372c4b1321245519d430540804a50b7ca8b6f1b34a2ec05cdbc24de7599af112d3e3c8db347e8799fe70f16e43c6566ba3aeb169463a3ecc486172deb2d9b80a3699c776e44fef20036bd946f1b4d054dd88a2c1aeb986199b0b2b7e58c42288824b74934d112fe1fc06e06b4d99fe1c5e725946b23210521e209cd507cce90b5f39a523f27e861f9e232aee50c3f585208b4573dcc0b897b6177f2ba20254fd5c50a033e849dee1b3a93bd2dc44ba8ca836cab2c2ae50e50b126284524fa0187af28628ff0face68d87709200329db1392852c8b8963fbe3d05fb1efe19f0ed5ca9fadc2f96f82187c24bb2512b2e85a66333a7e176605695211e1c8e0b9b9e82813e50654964945b1e1e66a90840396c7d10e23e47f364d2d3f660fa54598e18d1ca2ea4fe4f35a40a11f69f201c80b48eaee3e2e9b0eda63decf92bec08a70f731587d4ed0f218d5929285c8b2ccbc497e20db42de73885191fa453350335990184d8df805072f958d5354debda38f5421effaaafd6cb9b721ace74be0892d77679f62a4a126697cd35797f6858193da4ba1770c06aea2e5c59ec04b8ea26749e61b72ecdde403f3bc7e5e546cd799578cc939fa676dfd5e648576d4a06cbadb028adc2c0b461f145b2321f42e5e0f3b4fb898ecd461df07a6f5154067787bf74b5cc5c03704a1ce47494961931f0263b0aac32505102595957531a2de69dd71aac51f8a49902f81f21283dbe8e21e01e5d82517868826f86acf338d935aa6b4d5a25c8d540389b277dd9d64569d68baf0f71bd03dba45b92a7fc052601d1bd011a2fc6790a23f97c6fa5caeea040ab86841f268d39ce4f7caf01069df78bba098e04366492f0c2ac24f1bf16828752765fa523c9a4d42b71109d123e6be8c7b1ab3ccf8ea03404075fe1a9596f1bba1d267f9a7879ceece514818316c9c0583469d2367831fc42b517ea028a28df7c18d783d16ea2436cee2b15d52db68b5dfdee6b4d26f0905f9b030c911a04d078923a4136afea96eed6874462a482917353264cc9bee298f167ac65a6db4e4eda88044b39cc0b33183843eaa946564a00c3a0ab661f2c915e70bf0bb65bfbb6fa2eea20aed16bf2c1a1d00ec55fb4ff2f76b8e462ea70c19efa579c9ee78194b86708fdae66a9ce6e2cf3d366037798cfb50277ba6d2fd4866361022fd788ab7735b40b8b61d55e32243e06719e53992e9ac16c9c4b6e6933635c3c47c8f7e73e17dd54d0dd8aeba5d76de46894e7b3f9d3ec25ad78ee82297ba69905ea0fa094b8667faa2b8885e2187b3da80268aa1164761d7b0d6de206b676777348152b8ae1d4afed753bc63c739a5ca8ce7afb2b241a226bd9e502baba391b5b13f5054f070b65a9cf3a67063bfaa803ba390732cd03888f664023f888741d04d564e0b5674b0a183ace81452001b3fbb4214c77d42ca75376742c471e58f67307726d56a1032bd236610cbcbcd03d0d7a452900136897dc55bb3ce959d10d4e6a10fb635006bd8c41cd9ded2d3dfdd8f2e229590324a7370cb2124210b2330f4c56155caa09a2564932ceded8d92c79664dcdeb87faad7d3da006cc2ea267ee3df41e9677789cc5a8cc3b83add6491561b3047919e0648b1b2e97d7ad6f6c2aa80cab8e9ae10e1f75b1fdd0246151af709d259a6a0ed0b26bd711024965ecad7c41387de45443defce53f66612948694a6032279131c257119ed876a8e805dfb49576ef5c563574115ee87050d92d191bc761ef51d966918e2ef925639400069e3959d8fe19f36136e947ff430bf74e71da0aa5923b00000000")
|
||||
v := natFromString("0x838332321d443a3d30373d47301d47073847473a383d3030f25b3d3d3e00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002e00000000000000000041603038331c3d32f5303441e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e0e01c0a5459bfc7b9be9fcbb9d2383840464319434707303030f43a32f53034411c0a5459413820878787878787878787878787878787878787878787878787878787878787878787870630303a3a30334036605b923a6101f83638413943413960204337602043323801526040523241846038414143015238604060328452413841413638523c0240384141364036605b923a6101f83638413943413960204334602043323801526040523241846038414143015238604060328452413841413638523c02403841413638433030f25a8b83838383838383838383838383838383837d838383ffffffffffffffff838383838383838383000000000000000000030000007d26e27c7c8b83838383838383838383838383838383837d838383ffffffffffffffff83838383838383838383838383838383838383838383435960f535073030f3343200000000000000011881301938343030fa398383300000002300000000000000000000f11af4600c845252904141364138383c60406032414443095238010241414303364443434132305b595a15434160b042385341ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff47476043410536613603593a6005411c437405fcfcfcfcfcfcfc0000000000005a3b075815054359000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000")
|
||||
q := nat(nil).make(16)
|
||||
q.div(q, u, v)
|
||||
}
|
||||
|
|
@ -0,0 +1,884 @@
|
|||
// Copyright 2009 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
/*
|
||||
|
||||
Multi-precision division. Here be dragons.
|
||||
|
||||
Given u and v, where u is n+m digits, and v is n digits (with no leading zeros),
|
||||
the goal is to return quo, rem such that u = quo*v + rem, where 0 ≤ rem < v.
|
||||
That is, quo = ⌊u/v⌋ where ⌊x⌋ denotes the floor (truncation to integer) of x,
|
||||
and rem = u - quo·v.
|
||||
|
||||
|
||||
Long Division
|
||||
|
||||
Division in a computer proceeds the same as long division in elementary school,
|
||||
but computers are not as good as schoolchildren at following vague directions,
|
||||
so we have to be much more precise about the actual steps and what can happen.
|
||||
|
||||
We work from most to least significant digit of the quotient, doing:
|
||||
|
||||
• Guess a digit q, the number of v to subtract from the current
|
||||
section of u to zero out the topmost digit.
|
||||
• Check the guess by multiplying q·v and comparing it against
|
||||
the current section of u, adjusting the guess as needed.
|
||||
• Subtract q·v from the current section of u.
|
||||
• Add q to the corresponding section of the result quo.
|
||||
|
||||
When all digits have been processed, the final remainder is left in u
|
||||
and returned as rem.
|
||||
|
||||
For example, here is a sketch of dividing 5 digits by 3 digits (n=3, m=2).
|
||||
|
||||
q₂ q₁ q₀
|
||||
_________________
|
||||
v₂ v₁ v₀ ) u₄ u₃ u₂ u₁ u₀
|
||||
↓ ↓ ↓ | |
|
||||
[u₄ u₃ u₂]| |
|
||||
- [ q₂·v ]| |
|
||||
----------- ↓ |
|
||||
[ rem | u₁]|
|
||||
- [ q₁·v ]|
|
||||
----------- ↓
|
||||
[ rem | u₀]
|
||||
- [ q₀·v ]
|
||||
------------
|
||||
[ rem ]
|
||||
|
||||
Instead of creating new storage for the remainders and copying digits from u
|
||||
as indicated by the arrows, we use u's storage directly as both the source
|
||||
and destination of the subtractions, so that the remainders overwrite
|
||||
successive overlapping sections of u as the division proceeds, using a slice
|
||||
of u to identify the current section. This avoids all the copying as well as
|
||||
shifting of remainders.
|
||||
|
||||
Division of u with n+m digits by v with n digits (in base B) can in general
|
||||
produce at most m+1 digits, because:
|
||||
|
||||
• u < B^(n+m) [B^(n+m) has n+m+1 digits]
|
||||
• v ≥ B^(n-1) [B^(n-1) is the smallest n-digit number]
|
||||
• u/v < B^(n+m) / B^(n-1) [divide bounds for u, v]
|
||||
• u/v < B^(m+1) [simplify]
|
||||
|
||||
The first step is special: it takes the top n digits of u and divides them by
|
||||
the n digits of v, producing the first quotient digit and an n-digit remainder.
|
||||
In the example, q₂ = ⌊u₄u₃u₂ / v⌋.
|
||||
|
||||
The first step divides n digits by n digits to ensure that it produces only a
|
||||
single digit.
|
||||
|
||||
Each subsequent step appends the next digit from u to the remainder and divides
|
||||
those n+1 digits by the n digits of v, producing another quotient digit and a
|
||||
new n-digit remainder.
|
||||
|
||||
Subsequent steps divide n+1 digits by n digits, an operation that in general
|
||||
might produce two digits. However, as used in the algorithm, that division is
|
||||
guaranteed to produce only a single digit. The dividend is of the form
|
||||
rem·B + d, where rem is a remainder from the previous step and d is a single
|
||||
digit, so:
|
||||
|
||||
• rem ≤ v - 1 [rem is a remainder from dividing by v]
|
||||
• rem·B ≤ v·B - B [multiply by B]
|
||||
• d ≤ B - 1 [d is a single digit]
|
||||
• rem·B + d ≤ v·B - 1 [add]
|
||||
• rem·B + d < v·B [change ≤ to <]
|
||||
• (rem·B + d)/v < B [divide by v]
|
||||
|
||||
|
||||
Guess and Check
|
||||
|
||||
At each step we need to divide n+1 digits by n digits, but this is for the
|
||||
implementation of division by n digits, so we can't just invoke a division
|
||||
routine: we _are_ the division routine. Instead, we guess at the answer and
|
||||
then check it using multiplication. If the guess is wrong, we correct it.
|
||||
|
||||
How can this guessing possibly be efficient? It turns out that the following
|
||||
statement (let's call it the Good Guess Guarantee) is true.
|
||||
|
||||
If
|
||||
|
||||
• q = ⌊u/v⌋ where u is n+1 digits and v is n digits,
|
||||
• q < B, and
|
||||
• the topmost digit of v = vₙ₋₁ ≥ B/2,
|
||||
|
||||
then q̂ = ⌊uₙuₙ₋₁ / vₙ₋₁⌋ satisfies q ≤ q̂ ≤ q+2. (Proof below.)
|
||||
|
||||
That is, if we know the answer has only a single digit and we guess an answer
|
||||
by ignoring the bottom n-1 digits of u and v, using a 2-by-1-digit division,
|
||||
then that guess is at least as large as the correct answer. It is also not
|
||||
too much larger: it is off by at most two from the correct answer.
|
||||
|
||||
Note that in the first step of the overall division, which is an n-by-n-digit
|
||||
division, the 2-by-1 guess uses an implicit uₙ = 0.
|
||||
|
||||
Note that using a 2-by-1-digit division here does not mean calling ourselves
|
||||
recursively. Instead, we use an efficient direct hardware implementation of
|
||||
that operation.
|
||||
|
||||
Note that because q is u/v rounded down, q·v must not exceed u: u ≥ q·v.
|
||||
If a guess q̂ is too big, it will not satisfy this test. Viewed a different way,
|
||||
the remainder r̂ for a given q̂ is u - q̂·v, which must be positive. If it is
|
||||
negative, then the guess q̂ is too big.
|
||||
|
||||
This gives us a way to compute q. First compute q̂ with 2-by-1-digit division.
|
||||
Then, while u < q̂·v, decrement q̂; this loop executes at most twice, because
|
||||
q̂ ≤ q+2.
|
||||
|
||||
|
||||
Scaling Inputs
|
||||
|
||||
The Good Guess Guarantee requires that the top digit of v (vₙ₋₁) be at least B/2.
|
||||
For example in base 10, ⌊172/19⌋ = 9, but ⌊18/1⌋ = 18: the guess is wildly off
|
||||
because the first digit 1 is smaller than B/2 = 5.
|
||||
|
||||
We can ensure that v has a large top digit by multiplying both u and v by the
|
||||
right amount. Continuing the example, if we multiply both 172 and 19 by 3, we
|
||||
now have ⌊516/57⌋, the leading digit of v is now ≥ 5, and sure enough
|
||||
⌊51/5⌋ = 10 is much closer to the correct answer 9. It would be easier here
|
||||
to multiply by 4, because that can be done with a shift. Specifically, we can
|
||||
always count the number of leading zeros i in the first digit of v and then
|
||||
shift both u and v left by i bits.
|
||||
|
||||
Having scaled u and v, the value ⌊u/v⌋ is unchanged, but the remainder will
|
||||
be scaled: 172 mod 19 is 1, but 516 mod 57 is 3. We have to divide the remainder
|
||||
by the scaling factor (shifting right i bits) when we finish.
|
||||
|
||||
Note that these shifts happen before and after the entire division algorithm,
|
||||
not at each step in the per-digit iteration.
|
||||
|
||||
Note the effect of scaling inputs on the size of the possible quotient.
|
||||
In the scaled u/v, u can gain a digit from scaling; v never does, because we
|
||||
pick the scaling factor to make v's top digit larger but without overflowing.
|
||||
If u and v have n+m and n digits after scaling, then:
|
||||
|
||||
• u < B^(n+m) [B^(n+m) has n+m+1 digits]
|
||||
• v ≥ B^n / 2 [vₙ₋₁ ≥ B/2, so vₙ₋₁·B^(n-1) ≥ B^n/2]
|
||||
• u/v < B^(n+m) / (B^n / 2) [divide bounds for u, v]
|
||||
• u/v < 2 B^m [simplify]
|
||||
|
||||
The quotient can still have m+1 significant digits, but if so the top digit
|
||||
must be a 1. This provides a different way to handle the first digit of the
|
||||
result: compare the top n digits of u against v and fill in either a 0 or a 1.
|
||||
|
||||
|
||||
Refining Guesses
|
||||
|
||||
Before we check whether u < q̂·v, we can adjust our guess to change it from
|
||||
q̂ = ⌊uₙuₙ₋₁ / vₙ₋₁⌋ into the refined guess ⌊uₙuₙ₋₁uₙ₋₂ / vₙ₋₁vₙ₋₂⌋.
|
||||
Although not mentioned above, the Good Guess Guarantee also promises that this
|
||||
3-by-2-digit division guess is more precise and at most one away from the real
|
||||
answer q. The improvement from the 2-by-1 to the 3-by-2 guess can also be done
|
||||
without n-digit math.
|
||||
|
||||
If we have a guess q̂ = ⌊uₙuₙ₋₁ / vₙ₋₁⌋ and we want to see if it also equal to
|
||||
⌊uₙuₙ₋₁uₙ₋₂ / vₙ₋₁vₙ₋₂⌋, we can use the same check we would for the full division:
|
||||
if uₙuₙ₋₁uₙ₋₂ < q̂·vₙ₋₁vₙ₋₂, then the guess is too large and should be reduced.
|
||||
|
||||
Checking uₙuₙ₋₁uₙ₋₂ < q̂·vₙ₋₁vₙ₋₂ is the same as uₙuₙ₋₁uₙ₋₂ - q̂·vₙ₋₁vₙ₋₂ < 0,
|
||||
and
|
||||
|
||||
uₙuₙ₋₁uₙ₋₂ - q̂·vₙ₋₁vₙ₋₂ = (uₙuₙ₋₁·B + uₙ₋₂) - q̂·(vₙ₋₁·B + vₙ₋₂)
|
||||
[splitting off the bottom digit]
|
||||
= (uₙuₙ₋₁ - q̂·vₙ₋₁)·B + uₙ₋₂ - q̂·vₙ₋₂
|
||||
[regrouping]
|
||||
|
||||
The expression (uₙuₙ₋₁ - q̂·vₙ₋₁) is the remainder of uₙuₙ₋₁ / vₙ₋₁.
|
||||
If the initial guess returns both q̂ and its remainder r̂, then checking
|
||||
whether uₙuₙ₋₁uₙ₋₂ < q̂·vₙ₋₁vₙ₋₂ is the same as checking r̂·B + uₙ₋₂ < q̂·vₙ₋₂.
|
||||
|
||||
If we find that r̂·B + uₙ₋₂ < q̂·vₙ₋₂, then we can adjust the guess by
|
||||
decrementing q̂ and adding vₙ₋₁ to r̂. We repeat until r̂·B + uₙ₋₂ ≥ q̂·vₙ₋₂.
|
||||
(As before, this fixup is only needed at most twice.)
|
||||
|
||||
Now that q̂ = ⌊uₙuₙ₋₁uₙ₋₂ / vₙ₋₁vₙ₋₂⌋, as mentioned above it is at most one
|
||||
away from the correct q, and we've avoided doing any n-digit math.
|
||||
(If we need the new remainder, it can be computed as r̂·B + uₙ₋₂ - q̂·vₙ₋₂.)
|
||||
|
||||
The final check u < q̂·v and the possible fixup must be done at full precision.
|
||||
For random inputs, a fixup at this step is exceedingly rare: the 3-by-2 guess
|
||||
is not often wrong at all. But still we must do the check. Note that since the
|
||||
3-by-2 guess is off by at most 1, it can be convenient to perform the final
|
||||
u < q̂·v as part of the computation of the remainder r = u - q̂·v. If the
|
||||
subtraction underflows, decremeting q̂ and adding one v back to r is enough to
|
||||
arrive at the final q, r.
|
||||
|
||||
That's the entirety of long division: scale the inputs, and then loop over
|
||||
each output position, guessing, checking, and correcting the next output digit.
|
||||
|
||||
For a 2n-digit number divided by an n-digit number (the worst size-n case for
|
||||
division complexity), this algorithm uses n+1 iterations, each of which must do
|
||||
at least the 1-by-n-digit multiplication q̂·v. That's O(n) iterations of
|
||||
O(n) time each, so O(n²) time overall.
|
||||
|
||||
|
||||
Recursive Division
|
||||
|
||||
For very large inputs, it is possible to improve on the O(n²) algorithm.
|
||||
Let's call a group of n/2 real digits a (very) “wide digit”. We can run the
|
||||
standard long division algorithm explained above over the wide digits instead of
|
||||
the actual digits. This will result in many fewer steps, but the math involved in
|
||||
each step is more work.
|
||||
|
||||
Where basic long division uses a 2-by-1-digit division to guess the initial q̂,
|
||||
the new algorithm must use a 2-by-1-wide-digit division, which is of course
|
||||
really an n-by-n/2-digit division. That's OK: if we implement n-digit division
|
||||
in terms of n/2-digit division, the recursion will terminate when the divisor
|
||||
becomes small enough to handle with standard long division or even with the
|
||||
2-by-1 hardware instruction.
|
||||
|
||||
For example, here is a sketch of dividing 10 digits by 4, proceeding with
|
||||
wide digits corresponding to two regular digits. The first step, still special,
|
||||
must leave off a (regular) digit, dividing 5 by 4 and producing a 4-digit
|
||||
remainder less than v. The middle steps divide 6 digits by 4, guaranteed to
|
||||
produce two output digits each (one wide digit) with 4-digit remainders.
|
||||
The final step must use what it has: the 4-digit remainder plus one more,
|
||||
5 digits to divide by 4.
|
||||
|
||||
q₆ q₅ q₄ q₃ q₂ q₁ q₀
|
||||
_______________________________
|
||||
v₃ v₂ v₁ v₀ ) u₉ u₈ u₇ u₆ u₅ u₄ u₃ u₂ u₁ u₀
|
||||
↓ ↓ ↓ ↓ ↓ | | | | |
|
||||
[u₉ u₈ u₇ u₆ u₅]| | | | |
|
||||
- [ q₆q₅·v ]| | | | |
|
||||
----------------- ↓ ↓ | | |
|
||||
[ rem |u₄ u₃]| | |
|
||||
- [ q₄q₃·v ]| | |
|
||||
-------------------- ↓ ↓ |
|
||||
[ rem |u₂ u₁]|
|
||||
- [ q₂q₁·v ]|
|
||||
-------------------- ↓
|
||||
[ rem |u₀]
|
||||
- [ q₀·v ]
|
||||
------------------
|
||||
[ rem ]
|
||||
|
||||
An alternative would be to look ahead to how well n/2 divides into n+m and
|
||||
adjust the first step to use fewer digits as needed, making the first step
|
||||
more special to make the last step not special at all. For example, using the
|
||||
same input, we could choose to use only 4 digits in the first step, leaving
|
||||
a full wide digit for the last step:
|
||||
|
||||
q₆ q₅ q₄ q₃ q₂ q₁ q₀
|
||||
_______________________________
|
||||
v₃ v₂ v₁ v₀ ) u₉ u₈ u₇ u₆ u₅ u₄ u₃ u₂ u₁ u₀
|
||||
↓ ↓ ↓ ↓ | | | | | |
|
||||
[u₉ u₈ u₇ u₆]| | | | | |
|
||||
- [ q₆·v ]| | | | | |
|
||||
-------------- ↓ ↓ | | | |
|
||||
[ rem |u₅ u₄]| | | |
|
||||
- [ q₅q₄·v ]| | | |
|
||||
-------------------- ↓ ↓ | |
|
||||
[ rem |u₃ u₂]| |
|
||||
- [ q₃q₂·v ]| |
|
||||
-------------------- ↓ ↓
|
||||
[ rem |u₁ u₀]
|
||||
- [ q₁q₀·v ]
|
||||
---------------------
|
||||
[ rem ]
|
||||
|
||||
Today, the code in divRecursiveStep works like the first example. Perhaps in
|
||||
the future we will make it work like the alternative, to avoid a special case
|
||||
in the final iteration.
|
||||
|
||||
Either way, each step is a 3-by-2-wide-digit division approximated first by
|
||||
a 2-by-1-wide-digit division, just as we did for regular digits in long division.
|
||||
Because the actual answer we want is a 3-by-2-wide-digit division, instead of
|
||||
multiplying q̂·v directly during the fixup, we can use the quick refinement
|
||||
from long division (an n/2-by-n/2 multiply) to correct q to its actual value
|
||||
and also compute the remainder (as mentioned above), and then stop after that,
|
||||
never doing a full n-by-n multiply.
|
||||
|
||||
Instead of using an n-by-n/2-digit division to produce n/2 digits, we can add
|
||||
(not discard) one more real digit, doing an (n+1)-by-(n/2+1)-digit division that
|
||||
produces n/2+1 digits. That single extra digit tightens the Good Guess Guarantee
|
||||
to q ≤ q̂ ≤ q+1 and lets us drop long division's special treatment of the first
|
||||
digit. These benefits are discussed more after the Good Guess Guarantee proof
|
||||
below.
|
||||
|
||||
|
||||
How Fast is Recursive Division?
|
||||
|
||||
For a 2n-by-n-digit division, this algorithm runs a 4-by-2 long division over
|
||||
wide digits, producing two wide digits plus a possible leading regular digit 1,
|
||||
which can be handled without a recursive call. That is, the algorithm uses two
|
||||
full iterations, each using an n-by-n/2-digit division and an n/2-by-n/2-digit
|
||||
multiplication, along with a few n-digit additions and subtractions. The standard
|
||||
n-by-n-digit multiplication algorithm requires O(n²) time, making the overall
|
||||
algorithm require time T(n) where
|
||||
|
||||
T(n) = 2T(n/2) + O(n) + O(n²)
|
||||
|
||||
which, by the Bentley-Haken-Saxe theorem, ends up reducing to T(n) = O(n²).
|
||||
This is not an improvement over regular long division.
|
||||
|
||||
When the number of digits n becomes large enough, Karatsuba's algorithm for
|
||||
multiplication can be used instead, which takes O(n^log₂3) = O(n^1.6) time.
|
||||
(Karatsuba multiplication is implemented in func karatsuba in nat.go.)
|
||||
That makes the overall recursive division algorithm take O(n^1.6) time as well,
|
||||
which is an improvement, but again only for large enough numbers.
|
||||
|
||||
It is not critical to make sure that every recursion does only two recursive
|
||||
calls. While in general the number of recursive calls can change the time
|
||||
analysis, in this case doing three calls does not change the analysis:
|
||||
|
||||
T(n) = 3T(n/2) + O(n) + O(n^log₂3)
|
||||
|
||||
ends up being T(n) = O(n^log₂3). Because the Karatsuba multiplication taking
|
||||
time O(n^log₂3) is itself doing 3 half-sized recursions, doing three for the
|
||||
division does not hurt the asymptotic performance. Of course, it is likely
|
||||
still faster in practice to do two.
|
||||
|
||||
|
||||
Proof of the Good Guess Guarantee
|
||||
|
||||
Given numbers x, y, let us break them into the quotients and remainders when
|
||||
divided by some scaling factor S, with the added constraints that the quotient
|
||||
x/y and the high part of y are both less than some limit T, and that the high
|
||||
part of y is at least half as big as T.
|
||||
|
||||
x₁ = ⌊x/S⌋ y₁ = ⌊y/S⌋
|
||||
x₀ = x mod S y₀ = y mod S
|
||||
|
||||
x = x₁·S + x₀ 0 ≤ x₀ < S x/y < T
|
||||
y = y₁·S + y₀ 0 ≤ y₀ < S T/2 ≤ y₁ < T
|
||||
|
||||
And consider the two truncated quotients:
|
||||
|
||||
q = ⌊x/y⌋
|
||||
q̂ = ⌊x₁/y₁⌋
|
||||
|
||||
We will prove that q ≤ q̂ ≤ q+2.
|
||||
|
||||
The guarantee makes no real demands on the scaling factor S: it is simply the
|
||||
magnitude of the digits cut from both x and y to produce x₁ and y₁.
|
||||
The guarantee makes only limited demands on T: it must be large enough to hold
|
||||
the quotient x/y, and y₁ must have roughly the same size.
|
||||
|
||||
To apply to the earlier discussion of 2-by-1 guesses in long division,
|
||||
we would choose:
|
||||
|
||||
S = Bⁿ⁻¹
|
||||
T = B
|
||||
x = u
|
||||
x₁ = uₙuₙ₋₁
|
||||
x₀ = uₙ₋₂...u₀
|
||||
y = v
|
||||
y₁ = vₙ₋₁
|
||||
y₀ = vₙ₋₂...u₀
|
||||
|
||||
These simpler variables avoid repeating those longer expressions in the proof.
|
||||
|
||||
Note also that, by definition, truncating division ⌊x/y⌋ satisfies
|
||||
|
||||
x/y - 1 < ⌊x/y⌋ ≤ x/y.
|
||||
|
||||
This fact will be used a few times in the proofs.
|
||||
|
||||
Proof that q ≤ q̂:
|
||||
|
||||
q̂·y₁ = ⌊x₁/y₁⌋·y₁ [by definition, q̂ = ⌊x₁/y₁⌋]
|
||||
> (x₁/y₁ - 1)·y₁ [x₁/y₁ - 1 < ⌊x₁/y₁⌋]
|
||||
= x₁ - y₁ [distribute y₁]
|
||||
|
||||
So q̂·y₁ > x₁ - y₁.
|
||||
Since q̂·y₁ is an integer, q̂·y₁ ≥ x₁ - y₁ + 1.
|
||||
|
||||
q̂ - q = q̂ - ⌊x/y⌋ [by definition, q = ⌊x/y⌋]
|
||||
≥ q̂ - x/y [⌊x/y⌋ < x/y]
|
||||
= (1/y)·(q̂·y - x) [factor out 1/y]
|
||||
≥ (1/y)·(q̂·y₁·S - x) [y = y₁·S + y₀ ≥ y₁·S]
|
||||
≥ (1/y)·((x₁ - y₁ + 1)·S - x) [above: q̂·y₁ ≥ x₁ - y₁ + 1]
|
||||
= (1/y)·(x₁·S - y₁·S + S - x) [distribute S]
|
||||
= (1/y)·(S - x₀ - y₁·S) [-x = -x₁·S - x₀]
|
||||
> -y₁·S / y [x₀ < S, so S - x₀ < 0; drop it]
|
||||
≥ -1 [y₁·S ≤ y]
|
||||
|
||||
So q̂ - q > -1.
|
||||
Since q̂ - q is an integer, q̂ - q ≥ 0, or equivalently q ≤ q̂.
|
||||
|
||||
Proof that q̂ ≤ q+2:
|
||||
|
||||
x₁/y₁ - x/y = x₁·S/y₁·S - x/y [multiply left term by S/S]
|
||||
≤ x/y₁·S - x/y [x₁S ≤ x]
|
||||
= (x/y)·(y/y₁·S - 1) [factor out x/y]
|
||||
= (x/y)·((y - y₁·S)/y₁·S) [move -1 into y/y₁·S fraction]
|
||||
= (x/y)·(y₀/y₁·S) [y - y₁·S = y₀]
|
||||
= (x/y)·(1/y₁)·(y₀/S) [factor out 1/y₁]
|
||||
< (x/y)·(1/y₁) [y₀ < S, so y₀/S < 1]
|
||||
≤ (x/y)·(2/T) [y₁ ≥ T/2, so 1/y₁ ≤ 2/T]
|
||||
< T·(2/T) [x/y < T]
|
||||
= 2 [T·(2/T) = 2]
|
||||
|
||||
So x₁/y₁ - x/y < 2.
|
||||
|
||||
q̂ - q = ⌊x₁/y₁⌋ - q [by definition, q̂ = ⌊x₁/y₁⌋]
|
||||
= ⌊x₁/y₁⌋ - ⌊x/y⌋ [by definition, q = ⌊x/y⌋]
|
||||
≤ x₁/y₁ - ⌊x/y⌋ [⌊x₁/y₁⌋ ≤ x₁/y₁]
|
||||
< x₁/y₁ - (x/y - 1) [⌊x/y⌋ > x/y - 1]
|
||||
= (x₁/y₁ - x/y) + 1 [regrouping]
|
||||
< 2 + 1 [above: x₁/y₁ - x/y < 2]
|
||||
= 3
|
||||
|
||||
So q̂ - q < 3.
|
||||
Since q̂ - q is an integer, q̂ - q ≤ 2.
|
||||
|
||||
Note that when x/y < T/2, the bounds tighten to x₁/y₁ - x/y < 1 and therefore
|
||||
q̂ - q ≤ 1.
|
||||
|
||||
Note also that in the general case 2n-by-n division where we don't know that
|
||||
x/y < T, we do know that x/y < 2T, yielding the bound q̂ - q ≤ 4. So we could
|
||||
remove the special case first step of long division as long as we allow the
|
||||
first fixup loop to run up to four times. (Using a simple comparison to decide
|
||||
whether the first digit is 0 or 1 is still more efficient, though.)
|
||||
|
||||
Finally, note that when dividing three leading base-B digits by two (scaled),
|
||||
we have T = B² and x/y < B = T/B, a much tighter bound than x/y < T.
|
||||
This in turn yields the much tighter bound x₁/y₁ - x/y < 2/B. This means that
|
||||
⌊x₁/y₁⌋ and ⌊x/y⌋ can only differ when x/y is less than 2/B greater than an
|
||||
integer. For random x and y, the chance of this is 2/B, or, for large B,
|
||||
approximately zero. This means that after we produce the 3-by-2 guess in the
|
||||
long division algorithm, the fixup loop essentially never runs.
|
||||
|
||||
In the recursive algorithm, the extra digit in (2·⌊n/2⌋+1)-by-(⌊n/2⌋+1)-digit
|
||||
division has exactly the same effect: the probability of needing a fixup is the
|
||||
same 2/B. Even better, we can allow the general case x/y < 2T and the fixup
|
||||
probability only grows to 4/B, still essentially zero.
|
||||
|
||||
|
||||
References
|
||||
|
||||
There are no great references for implementing long division; thus this comment.
|
||||
Here are some notes about what to expect from the obvious references.
|
||||
|
||||
Knuth Volume 2 (Seminumerical Algorithms) section 4.3.1 is the usual canonical
|
||||
reference for long division, but that entire series is highly compressed, never
|
||||
repeating a necessary fact and leaving important insights to the exercises.
|
||||
For example, no rationale whatsoever is given for the calculation that extends
|
||||
q̂ from a 2-by-1 to a 3-by-2 guess, nor why it reduces the error bound.
|
||||
The proof that the calculation even has the desired effect is left to exercises.
|
||||
The solutions to those exercises provided at the back of the book are entirely
|
||||
calculations, still with no explanation as to what is going on or how you would
|
||||
arrive at the idea of doing those exact calculations. Nowhere is it mentioned
|
||||
that this test extends the 2-by-1 guess into a 3-by-2 guess. The proof of the
|
||||
Good Guess Guarantee is only for the 2-by-1 guess and argues by contradiction,
|
||||
making it difficult to understand how modifications like adding another digit
|
||||
or adjusting the quotient range affects the overall bound.
|
||||
|
||||
All that said, Knuth remains the canonical reference. It is dense but packed
|
||||
full of information and references, and the proofs are simpler than many other
|
||||
presentations. The proofs above are reworkings of Knuth's to remove the
|
||||
arguments by contradiction and add explanations or steps that Knuth omitted.
|
||||
But beware of errors in older printings. Take the published errata with you.
|
||||
|
||||
Brinch Hansen's “Multiple-length Division Revisited: a Tour of the Minefield”
|
||||
starts with a blunt critique of Knuth's presentation (among others) and then
|
||||
presents a more detailed and easier to follow treatment of long division,
|
||||
including an implementation in Pascal. But the algorithm and implementation
|
||||
work entirely in terms of 3-by-2 division, which is much less useful on modern
|
||||
hardware than an algorithm using 2-by-1 division. The proofs are a bit too
|
||||
focused on digit counting and seem needlessly complex, especially compared to
|
||||
the ones given above.
|
||||
|
||||
Burnikel and Ziegler's “Fast Recursive Division” introduced the key insight of
|
||||
implementing division by an n-digit divisor using recursive calls to division
|
||||
by an n/2-digit divisor, relying on Karatsuba multiplication to yield a
|
||||
sub-quadratic run time. However, the presentation decisions are made almost
|
||||
entirely for the purpose of simplifying the run-time analysis, rather than
|
||||
simplifying the presentation. Instead of a single algorithm that loops over
|
||||
quotient digits, the paper presents two mutually-recursive algorithms, for
|
||||
2n-by-n and 3n-by-2n. The paper also does not present any general (n+m)-by-n
|
||||
algorithm.
|
||||
|
||||
The proofs in the paper are remarkably complex, especially considering that
|
||||
the algorithm is at its core just long division on wide digits, so that the
|
||||
usual long division proofs apply essentially unaltered.
|
||||
*/
|
||||
|
||||
package big
|
||||
|
||||
import "math/bits"
|
||||
|
||||
// div returns q, r such that q = ⌊u/v⌋ and r = u%v = u - q·v.
|
||||
// It uses z and z2 as the storage for q and r.
|
||||
func (z nat) div(z2, u, v nat) (q, r nat) {
|
||||
if len(v) == 0 {
|
||||
panic("division by zero")
|
||||
}
|
||||
|
||||
if u.cmp(v) < 0 {
|
||||
q = z[:0]
|
||||
r = z2.set(u)
|
||||
return
|
||||
}
|
||||
|
||||
if len(v) == 1 {
|
||||
// Short division: long optimized for a single-word divisor.
|
||||
// In that case, the 2-by-1 guess is all we need at each step.
|
||||
var r2 Word
|
||||
q, r2 = z.divW(u, v[0])
|
||||
r = z2.setWord(r2)
|
||||
return
|
||||
}
|
||||
|
||||
q, r = z.divLarge(z2, u, v)
|
||||
return
|
||||
}
|
||||
|
||||
// divW returns q, r such that q = ⌊x/y⌋ and r = x%y = x - q·y.
|
||||
// It uses z as the storage for q.
|
||||
// Note that y is a single digit (Word), not a big number.
|
||||
func (z nat) divW(x nat, y Word) (q nat, r Word) {
|
||||
m := len(x)
|
||||
switch {
|
||||
case y == 0:
|
||||
panic("division by zero")
|
||||
case y == 1:
|
||||
q = z.set(x) // result is x
|
||||
return
|
||||
case m == 0:
|
||||
q = z[:0] // result is 0
|
||||
return
|
||||
}
|
||||
// m > 0
|
||||
z = z.make(m)
|
||||
r = divWVW(z, 0, x, y)
|
||||
q = z.norm()
|
||||
return
|
||||
}
|
||||
|
||||
// modW returns x % d.
|
||||
func (x nat) modW(d Word) (r Word) {
|
||||
// TODO(agl): we don't actually need to store the q value.
|
||||
var q nat
|
||||
q = q.make(len(x))
|
||||
return divWVW(q, 0, x, d)
|
||||
}
|
||||
|
||||
// divWVW overwrites z with ⌊x/y⌋, returning the remainder r.
|
||||
// The caller must ensure that len(z) = len(x).
|
||||
func divWVW(z []Word, xn Word, x []Word, y Word) (r Word) {
|
||||
r = xn
|
||||
if len(x) == 1 {
|
||||
qq, rr := bits.Div(uint(r), uint(x[0]), uint(y))
|
||||
z[0] = Word(qq)
|
||||
return Word(rr)
|
||||
}
|
||||
rec := reciprocalWord(y)
|
||||
for i := len(z) - 1; i >= 0; i-- {
|
||||
z[i], r = divWW(r, x[i], y, rec)
|
||||
}
|
||||
return r
|
||||
}
|
||||
|
||||
// div returns q, r such that q = ⌊uIn/vIn⌋ and r = uIn%vIn = uIn - q·vIn.
|
||||
// It uses z and u as the storage for q and r.
|
||||
// The caller must ensure that len(vIn) ≥ 2 (use divW otherwise)
|
||||
// and that len(uIn) ≥ len(vIn) (the answer is 0, uIn otherwise).
|
||||
func (z nat) divLarge(u, uIn, vIn nat) (q, r nat) {
|
||||
n := len(vIn)
|
||||
m := len(uIn) - n
|
||||
|
||||
// Scale the inputs so vIn's top bit is 1 (see “Scaling Inputs” above).
|
||||
// vIn is treated as a read-only input (it may be in use by another
|
||||
// goroutine), so we must make a copy.
|
||||
// uIn is copied to u.
|
||||
shift := nlz(vIn[n-1])
|
||||
vp := getNat(n)
|
||||
v := *vp
|
||||
shlVU(v, vIn, shift)
|
||||
u = u.make(len(uIn) + 1)
|
||||
u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift)
|
||||
|
||||
// The caller should not pass aliased z and u, since those are
|
||||
// the two different outputs, but correct just in case.
|
||||
if alias(z, u) {
|
||||
z = nil
|
||||
}
|
||||
q = z.make(m + 1)
|
||||
|
||||
// Use basic or recursive long division depending on size.
|
||||
if n < divRecursiveThreshold {
|
||||
q.divBasic(u, v)
|
||||
} else {
|
||||
q.divRecursive(u, v)
|
||||
}
|
||||
putNat(vp)
|
||||
|
||||
q = q.norm()
|
||||
|
||||
// Undo scaling of remainder.
|
||||
shrVU(u, u, shift)
|
||||
r = u.norm()
|
||||
|
||||
return q, r
|
||||
}
|
||||
|
||||
// divBasic implements long division as described above.
|
||||
// It overwrites q with ⌊u/v⌋ and overwrites u with the remainder r.
|
||||
// q must be large enough to hold ⌊u/v⌋.
|
||||
func (q nat) divBasic(u, v nat) {
|
||||
n := len(v)
|
||||
m := len(u) - n
|
||||
|
||||
qhatvp := getNat(n + 1)
|
||||
qhatv := *qhatvp
|
||||
|
||||
// Set up for divWW below, precomputing reciprocal argument.
|
||||
vn1 := v[n-1]
|
||||
rec := reciprocalWord(vn1)
|
||||
|
||||
// Compute each digit of quotient.
|
||||
for j := m; j >= 0; j-- {
|
||||
// Compute the 2-by-1 guess q̂.
|
||||
// The first iteration must invent a leading 0 for u.
|
||||
qhat := Word(_M)
|
||||
var ujn Word
|
||||
if j+n < len(u) {
|
||||
ujn = u[j+n]
|
||||
}
|
||||
|
||||
// ujn ≤ vn1, or else q̂ would be more than one digit.
|
||||
// For ujn == vn1, we set q̂ to the max digit M above.
|
||||
// Otherwise, we compute the 2-by-1 guess.
|
||||
if ujn != vn1 {
|
||||
var rhat Word
|
||||
qhat, rhat = divWW(ujn, u[j+n-1], vn1, rec)
|
||||
|
||||
// Refine q̂ to a 3-by-2 guess. See “Refining Guesses” above.
|
||||
vn2 := v[n-2]
|
||||
x1, x2 := mulWW(qhat, vn2)
|
||||
ujn2 := u[j+n-2]
|
||||
for greaterThan(x1, x2, rhat, ujn2) { // x1x2 > r̂ u[j+n-2]
|
||||
qhat--
|
||||
prevRhat := rhat
|
||||
rhat += vn1
|
||||
// If r̂ overflows, then
|
||||
// r̂ u[j+n-2]v[n-1] is now definitely > x1 x2.
|
||||
if rhat < prevRhat {
|
||||
break
|
||||
}
|
||||
// TODO(rsc): No need for a full mulWW.
|
||||
// x2 += vn2; if x2 overflows, x1++
|
||||
x1, x2 = mulWW(qhat, vn2)
|
||||
}
|
||||
}
|
||||
|
||||
// Compute q̂·v.
|
||||
qhatv[n] = mulAddVWW(qhatv[0:n], v, qhat, 0)
|
||||
qhl := len(qhatv)
|
||||
if j+qhl > len(u) && qhatv[n] == 0 {
|
||||
qhl--
|
||||
}
|
||||
|
||||
// Subtract q̂·v from the current section of u.
|
||||
// If it underflows, q̂·v > u, which we fix up
|
||||
// by decrementing q̂ and adding v back.
|
||||
c := subVV(u[j:j+qhl], u[j:], qhatv)
|
||||
if c != 0 {
|
||||
c := addVV(u[j:j+n], u[j:], v)
|
||||
// If n == qhl, the carry from subVV and the carry from addVV
|
||||
// cancel out and don't affect u[j+n].
|
||||
if n < qhl {
|
||||
u[j+n] += c
|
||||
}
|
||||
qhat--
|
||||
}
|
||||
|
||||
// Save quotient digit.
|
||||
// Caller may know the top digit is zero and not leave room for it.
|
||||
if j == m && m == len(q) && qhat == 0 {
|
||||
continue
|
||||
}
|
||||
q[j] = qhat
|
||||
}
|
||||
|
||||
putNat(qhatvp)
|
||||
}
|
||||
|
||||
// greaterThan reports whether the two digit numbers x1 x2 > y1 y2.
|
||||
// TODO(rsc): In contradiction to most of this file, x1 is the high
|
||||
// digit and x2 is the low digit. This should be fixed.
|
||||
func greaterThan(x1, x2, y1, y2 Word) bool {
|
||||
return x1 > y1 || x1 == y1 && x2 > y2
|
||||
}
|
||||
|
||||
// divRecursiveThreshold is the number of divisor digits
|
||||
// at which point divRecursive is faster than divBasic.
|
||||
const divRecursiveThreshold = 100
|
||||
|
||||
// divRecursive implements recursive division as described above.
|
||||
// It overwrites z with ⌊u/v⌋ and overwrites u with the remainder r.
|
||||
// z must be large enough to hold ⌊u/v⌋.
|
||||
// This function is just for allocating and freeing temporaries
|
||||
// around divRecursiveStep, the real implementation.
|
||||
func (z nat) divRecursive(u, v nat) {
|
||||
// Recursion depth is (much) less than 2 log₂(len(v)).
|
||||
// Allocate a slice of temporaries to be reused across recursion,
|
||||
// plus one extra temporary not live across the recursion.
|
||||
recDepth := 2 * bits.Len(uint(len(v)))
|
||||
tmp := getNat(3 * len(v))
|
||||
temps := make([]*nat, recDepth)
|
||||
|
||||
z.clear()
|
||||
z.divRecursiveStep(u, v, 0, tmp, temps)
|
||||
|
||||
// Free temporaries.
|
||||
for _, n := range temps {
|
||||
if n != nil {
|
||||
putNat(n)
|
||||
}
|
||||
}
|
||||
putNat(tmp)
|
||||
}
|
||||
|
||||
// divRecursiveStep is the actual implementation of recursive division.
|
||||
// It adds ⌊u/v⌋ to z and overwrites u with the remainder r.
|
||||
// z must be large enough to hold ⌊u/v⌋.
|
||||
// It uses temps[depth] (allocating if needed) as a temporary live across
|
||||
// the recursive call. It also uses tmp, but not live across the recursion.
|
||||
func (z nat) divRecursiveStep(u, v nat, depth int, tmp *nat, temps []*nat) {
|
||||
// u is a subsection of the original and may have leading zeros.
|
||||
// TODO(rsc): The v = v.norm() is useless and should be removed.
|
||||
// We know (and require) that v's top digit is ≥ B/2.
|
||||
u = u.norm()
|
||||
v = v.norm()
|
||||
if len(u) == 0 {
|
||||
z.clear()
|
||||
return
|
||||
}
|
||||
|
||||
// Fall back to basic division if the problem is now small enough.
|
||||
n := len(v)
|
||||
if n < divRecursiveThreshold {
|
||||
z.divBasic(u, v)
|
||||
return
|
||||
}
|
||||
|
||||
// Nothing to do if u is shorter than v (implies u < v).
|
||||
m := len(u) - n
|
||||
if m < 0 {
|
||||
return
|
||||
}
|
||||
|
||||
// We consider B digits in a row as a single wide digit.
|
||||
// (See “Recursive Division” above.)
|
||||
//
|
||||
// TODO(rsc): rename B to Wide, to avoid confusion with _B,
|
||||
// which is something entirely different.
|
||||
// TODO(rsc): Look into whether using ⌈n/2⌉ is better than ⌊n/2⌋.
|
||||
B := n / 2
|
||||
|
||||
// Allocate a nat for qhat below.
|
||||
if temps[depth] == nil {
|
||||
temps[depth] = getNat(n) // TODO(rsc): Can be just B+1.
|
||||
} else {
|
||||
*temps[depth] = temps[depth].make(B + 1)
|
||||
}
|
||||
|
||||
// Compute each wide digit of the quotient.
|
||||
//
|
||||
// TODO(rsc): Change the loop to be
|
||||
// for j := (m+B-1)/B*B; j > 0; j -= B {
|
||||
// which will make the final step a regular step, letting us
|
||||
// delete what amounts to an extra copy of the loop body below.
|
||||
j := m
|
||||
for j > B {
|
||||
// Divide u[j-B:j+n] (3 wide digits) by v (2 wide digits).
|
||||
// First make the 2-by-1-wide-digit guess using a recursive call.
|
||||
// Then extend the guess to the full 3-by-2 (see “Refining Guesses”).
|
||||
//
|
||||
// For the 2-by-1-wide-digit guess, instead of doing 2B-by-B-digit,
|
||||
// we use a (2B+1)-by-(B+1) digit, which handles the possibility that
|
||||
// the result has an extra leading 1 digit as well as guaranteeing
|
||||
// that the computed q̂ will be off by at most 1 instead of 2.
|
||||
|
||||
// s is the number of digits to drop from the 3B- and 2B-digit chunks.
|
||||
// We drop B-1 to be left with 2B+1 and B+1.
|
||||
s := (B - 1)
|
||||
|
||||
// uu is the up-to-3B-digit section of u we are working on.
|
||||
uu := u[j-B:]
|
||||
|
||||
// Compute the 2-by-1 guess q̂, leaving r̂ in uu[s:B+n].
|
||||
qhat := *temps[depth]
|
||||
qhat.clear()
|
||||
qhat.divRecursiveStep(uu[s:B+n], v[s:], depth+1, tmp, temps)
|
||||
qhat = qhat.norm()
|
||||
|
||||
// Extend to a 3-by-2 quotient and remainder.
|
||||
// Because divRecursiveStep overwrote the top part of uu with
|
||||
// the remainder r̂, the full uu already contains the equivalent
|
||||
// of r̂·B + uₙ₋₂ from the “Refining Guesses” discussion.
|
||||
// Subtracting q̂·vₙ₋₂ from it will compute the full-length remainder.
|
||||
// If that subtraction underflows, q̂·v > u, which we fix up
|
||||
// by decrementing q̂ and adding v back, same as in long division.
|
||||
|
||||
// TODO(rsc): Instead of subtract and fix-up, this code is computing
|
||||
// q̂·vₙ₋₂ and decrementing q̂ until that product is ≤ u.
|
||||
// But we can do the subtraction directly, as in the comment above
|
||||
// and in long division, because we know that q̂ is wrong by at most one.
|
||||
qhatv := tmp.make(3 * n)
|
||||
qhatv.clear()
|
||||
qhatv = qhatv.mul(qhat, v[:s])
|
||||
for i := 0; i < 2; i++ {
|
||||
e := qhatv.cmp(uu.norm())
|
||||
if e <= 0 {
|
||||
break
|
||||
}
|
||||
subVW(qhat, qhat, 1)
|
||||
c := subVV(qhatv[:s], qhatv[:s], v[:s])
|
||||
if len(qhatv) > s {
|
||||
subVW(qhatv[s:], qhatv[s:], c)
|
||||
}
|
||||
addAt(uu[s:], v[s:], 0)
|
||||
}
|
||||
if qhatv.cmp(uu.norm()) > 0 {
|
||||
panic("impossible")
|
||||
}
|
||||
c := subVV(uu[:len(qhatv)], uu[:len(qhatv)], qhatv)
|
||||
if c > 0 {
|
||||
subVW(uu[len(qhatv):], uu[len(qhatv):], c)
|
||||
}
|
||||
addAt(z, qhat, j-B)
|
||||
j -= B
|
||||
}
|
||||
|
||||
// TODO(rsc): Rewrite loop as described above and delete all this code.
|
||||
|
||||
// Now u < (v<<B), compute lower bits in the same way.
|
||||
// Choose shift = B-1 again.
|
||||
s := B - 1
|
||||
qhat := *temps[depth]
|
||||
qhat.clear()
|
||||
qhat.divRecursiveStep(u[s:].norm(), v[s:], depth+1, tmp, temps)
|
||||
qhat = qhat.norm()
|
||||
qhatv := tmp.make(3 * n)
|
||||
qhatv.clear()
|
||||
qhatv = qhatv.mul(qhat, v[:s])
|
||||
// Set the correct remainder as before.
|
||||
for i := 0; i < 2; i++ {
|
||||
if e := qhatv.cmp(u.norm()); e > 0 {
|
||||
subVW(qhat, qhat, 1)
|
||||
c := subVV(qhatv[:s], qhatv[:s], v[:s])
|
||||
if len(qhatv) > s {
|
||||
subVW(qhatv[s:], qhatv[s:], c)
|
||||
}
|
||||
addAt(u[s:], v[s:], 0)
|
||||
}
|
||||
}
|
||||
if qhatv.cmp(u.norm()) > 0 {
|
||||
panic("impossible")
|
||||
}
|
||||
c := subVV(u[0:len(qhatv)], u[0:len(qhatv)], qhatv)
|
||||
if c > 0 {
|
||||
c = subVW(u[len(qhatv):], u[len(qhatv):], c)
|
||||
}
|
||||
if c > 0 {
|
||||
panic("impossible")
|
||||
}
|
||||
|
||||
// Done!
|
||||
addAt(z, qhat.norm(), 0)
|
||||
}
|
||||
|
|
@ -0,0 +1,277 @@
|
|||
// Copyright 2009 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
// This file provides Go implementations of elementary multi-precision
|
||||
// arithmetic operations on word vectors. These have the suffix _g.
|
||||
// These are needed for platforms without assembly implementations of these routines.
|
||||
// This file also contains elementary operations that can be implemented
|
||||
// sufficiently efficiently in Go.
|
||||
|
||||
package mybig
|
||||
|
||||
import "math/bits"
|
||||
|
||||
// A Word represents a single digit of a multi-precision unsigned integer.
|
||||
type Word uint
|
||||
|
||||
const (
|
||||
_S = _W / 8 // word size in bytes
|
||||
|
||||
_W = bits.UintSize // word size in bits
|
||||
_B = 1 << _W // digit base
|
||||
_M = _B - 1 // digit mask
|
||||
)
|
||||
|
||||
// Many of the loops in this file are of the form
|
||||
// for i := 0; i < len(z) && i < len(x) && i < len(y); i++
|
||||
// i < len(z) is the real condition.
|
||||
// However, checking i < len(x) && i < len(y) as well is faster than
|
||||
// having the compiler do a bounds check in the body of the loop;
|
||||
// remarkably it is even faster than hoisting the bounds check
|
||||
// out of the loop, by doing something like
|
||||
// _, _ = x[len(z)-1], y[len(z)-1]
|
||||
// There are other ways to hoist the bounds check out of the loop,
|
||||
// but the compiler's BCE isn't powerful enough for them (yet?).
|
||||
// See the discussion in CL 164966.
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
// Elementary operations on words
|
||||
//
|
||||
// These operations are used by the vector operations below.
|
||||
|
||||
// z1<<_W + z0 = x*y
|
||||
func mulWW_g(x, y Word) (z1, z0 Word) {
|
||||
hi, lo := bits.Mul(uint(x), uint(y))
|
||||
return Word(hi), Word(lo)
|
||||
}
|
||||
|
||||
// z1<<_W + z0 = x*y + c
|
||||
func mulAddWWW_g(x, y, c Word) (z1, z0 Word) {
|
||||
hi, lo := bits.Mul(uint(x), uint(y))
|
||||
var cc uint
|
||||
lo, cc = bits.Add(lo, uint(c), 0)
|
||||
return Word(hi + cc), Word(lo)
|
||||
}
|
||||
|
||||
// nlz returns the number of leading zeros in x.
|
||||
// Wraps bits.LeadingZeros call for convenience.
|
||||
func nlz(x Word) uint {
|
||||
return uint(bits.LeadingZeros(uint(x)))
|
||||
}
|
||||
|
||||
// The resulting carry c is either 0 or 1.
|
||||
func addVV_g(z, x, y []Word) (c Word) {
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x) && i < len(y); i++ {
|
||||
zi, cc := bits.Add(uint(x[i]), uint(y[i]), uint(c))
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// The resulting carry c is either 0 or 1.
|
||||
func subVV_g(z, x, y []Word) (c Word) {
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x) && i < len(y); i++ {
|
||||
zi, cc := bits.Sub(uint(x[i]), uint(y[i]), uint(c))
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// The resulting carry c is either 0 or 1.
|
||||
func addVW_g(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
zi, cc := bits.Add(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// addVWlarge is addVW, but intended for large z.
|
||||
// The only difference is that we check on every iteration
|
||||
// whether we are done with carries,
|
||||
// and if so, switch to a much faster copy instead.
|
||||
// This is only a good idea for large z,
|
||||
// because the overhead of the check and the function call
|
||||
// outweigh the benefits when z is small.
|
||||
func addVWlarge(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
if c == 0 {
|
||||
copy(z[i:], x[i:])
|
||||
return
|
||||
}
|
||||
zi, cc := bits.Add(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
func subVW_g(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
zi, cc := bits.Sub(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// subVWlarge is to subVW as addVWlarge is to addVW.
|
||||
func subVWlarge(z, x []Word, y Word) (c Word) {
|
||||
c = y
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
if c == 0 {
|
||||
copy(z[i:], x[i:])
|
||||
return
|
||||
}
|
||||
zi, cc := bits.Sub(uint(x[i]), uint(c), 0)
|
||||
z[i] = Word(zi)
|
||||
c = Word(cc)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
func shlVU_g(z, x []Word, s uint) (c Word) {
|
||||
if s == 0 {
|
||||
copy(z, x)
|
||||
return
|
||||
}
|
||||
if len(z) == 0 {
|
||||
return
|
||||
}
|
||||
s &= _W - 1 // hint to the compiler that shifts by s don't need guard code
|
||||
ŝ := _W - s
|
||||
ŝ &= _W - 1 // ditto
|
||||
c = x[len(z)-1] >> ŝ
|
||||
for i := len(z) - 1; i > 0; i-- {
|
||||
z[i] = x[i]<<s | x[i-1]>>ŝ
|
||||
}
|
||||
z[0] = x[0] << s
|
||||
return
|
||||
}
|
||||
|
||||
func shrVU_g(z, x []Word, s uint) (c Word) {
|
||||
if s == 0 {
|
||||
copy(z, x)
|
||||
return
|
||||
}
|
||||
if len(z) == 0 {
|
||||
return
|
||||
}
|
||||
if len(x) != len(z) {
|
||||
// This is an invariant guaranteed by the caller.
|
||||
panic("len(x) != len(z)")
|
||||
}
|
||||
s &= _W - 1 // hint to the compiler that shifts by s don't need guard code
|
||||
ŝ := _W - s
|
||||
ŝ &= _W - 1 // ditto
|
||||
c = x[0] << ŝ
|
||||
for i := 1; i < len(z); i++ {
|
||||
z[i-1] = x[i-1]>>s | x[i]<<ŝ
|
||||
}
|
||||
z[len(z)-1] = x[len(z)-1] >> s
|
||||
return
|
||||
}
|
||||
|
||||
func mulAddVWW_g(z, x []Word, y, r Word) (c Word) {
|
||||
c = r
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
c, z[i] = mulAddWWW_g(x[i], y, c)
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
func addMulVVW_g(z, x []Word, y Word) (c Word) {
|
||||
// The comment near the top of this file discusses this for loop condition.
|
||||
for i := 0; i < len(z) && i < len(x); i++ {
|
||||
z1, z0 := mulAddWWW_g(x[i], y, z[i])
|
||||
lo, cc := bits.Add(uint(z0), uint(c), 0)
|
||||
c, z[i] = Word(cc), Word(lo)
|
||||
c += z1
|
||||
}
|
||||
return
|
||||
}
|
||||
|
||||
// q = ( x1 << _W + x0 - r)/y. m = floor(( _B^2 - 1 ) / d - _B). Requiring x1<y.
|
||||
// An approximate reciprocal with a reference to "Improved Division by Invariant Integers
|
||||
// (IEEE Transactions on Computers, 11 Jun. 2010)"
|
||||
func divWW(x1, x0, y, m Word) (q, r Word) {
|
||||
s := nlz(y)
|
||||
if s != 0 {
|
||||
x1 = x1<<s | x0>>(_W-s)
|
||||
x0 <<= s
|
||||
y <<= s
|
||||
}
|
||||
d := uint(y)
|
||||
// We know that
|
||||
// m = ⎣(B^2-1)/d⎦-B
|
||||
// ⎣(B^2-1)/d⎦ = m+B
|
||||
// (B^2-1)/d = m+B+delta1 0 <= delta1 <= (d-1)/d
|
||||
// B^2/d = m+B+delta2 0 <= delta2 <= 1
|
||||
// The quotient we're trying to compute is
|
||||
// quotient = ⎣(x1*B+x0)/d⎦
|
||||
// = ⎣(x1*B*(B^2/d)+x0*(B^2/d))/B^2⎦
|
||||
// = ⎣(x1*B*(m+B+delta2)+x0*(m+B+delta2))/B^2⎦
|
||||
// = ⎣(x1*m+x1*B+x0)/B + x0*m/B^2 + delta2*(x1*B+x0)/B^2⎦
|
||||
// The latter two terms of this three-term sum are between 0 and 1.
|
||||
// So we can compute just the first term, and we will be low by at most 2.
|
||||
t1, t0 := bits.Mul(uint(m), uint(x1))
|
||||
_, c := bits.Add(t0, uint(x0), 0)
|
||||
t1, _ = bits.Add(t1, uint(x1), c)
|
||||
// The quotient is either t1, t1+1, or t1+2.
|
||||
// We'll try t1 and adjust if needed.
|
||||
qq := t1
|
||||
// compute remainder r=x-d*q.
|
||||
dq1, dq0 := bits.Mul(d, qq)
|
||||
r0, b := bits.Sub(uint(x0), dq0, 0)
|
||||
r1, _ := bits.Sub(uint(x1), dq1, b)
|
||||
// The remainder we just computed is bounded above by B+d:
|
||||
// r = x1*B + x0 - d*q.
|
||||
// = x1*B + x0 - d*⎣(x1*m+x1*B+x0)/B⎦
|
||||
// = x1*B + x0 - d*((x1*m+x1*B+x0)/B-alpha) 0 <= alpha < 1
|
||||
// = x1*B + x0 - x1*d/B*m - x1*d - x0*d/B + d*alpha
|
||||
// = x1*B + x0 - x1*d/B*⎣(B^2-1)/d-B⎦ - x1*d - x0*d/B + d*alpha
|
||||
// = x1*B + x0 - x1*d/B*⎣(B^2-1)/d-B⎦ - x1*d - x0*d/B + d*alpha
|
||||
// = x1*B + x0 - x1*d/B*((B^2-1)/d-B-beta) - x1*d - x0*d/B + d*alpha 0 <= beta < 1
|
||||
// = x1*B + x0 - x1*B + x1/B + x1*d + x1*d/B*beta - x1*d - x0*d/B + d*alpha
|
||||
// = x0 + x1/B + x1*d/B*beta - x0*d/B + d*alpha
|
||||
// = x0*(1-d/B) + x1*(1+d*beta)/B + d*alpha
|
||||
// < B*(1-d/B) + d*B/B + d because x0<B (and 1-d/B>0), x1<d, 1+d*beta<=B, alpha<1
|
||||
// = B - d + d + d
|
||||
// = B+d
|
||||
// So r1 can only be 0 or 1. If r1 is 1, then we know q was too small.
|
||||
// Add 1 to q and subtract d from r. That guarantees that r is <B, so
|
||||
// we no longer need to keep track of r1.
|
||||
if r1 != 0 {
|
||||
qq++
|
||||
r0 -= d
|
||||
}
|
||||
// If the remainder is still too large, increment q one more time.
|
||||
if r0 >= d {
|
||||
qq++
|
||||
r0 -= d
|
||||
}
|
||||
return Word(qq), Word(r0 >> s)
|
||||
}
|
||||
|
||||
// reciprocalWord return the reciprocal of the divisor. rec = floor(( _B^2 - 1 ) / u - _B). u = d1 << nlz(d1).
|
||||
func reciprocalWord(d1 Word) Word {
|
||||
u := uint(d1 << nlz(d1))
|
||||
x1 := ^u
|
||||
x0 := uint(_M)
|
||||
rec, _ := bits.Div(x1, x0, u) // (_B^2-1)/U-_B = (_B*(_M-C)+_M)/U
|
||||
return Word(rec)
|
||||
}
|
||||
|
|
@ -0,0 +1,51 @@
|
|||
// Copyright 2015 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
package mybig
|
||||
|
||||
func mulWW(x, y Word) (z1, z0 Word) {
|
||||
return mulWW_g(x, y)
|
||||
}
|
||||
|
||||
func addVV(z, x, y []Word) (c Word) {
|
||||
return addVV_g(z, x, y)
|
||||
}
|
||||
|
||||
func subVV(z, x, y []Word) (c Word) {
|
||||
return subVV_g(z, x, y)
|
||||
}
|
||||
|
||||
func addVW(z, x []Word, y Word) (c Word) {
|
||||
// TODO: remove indirect function call when golang.org/issue/30548 is fixed
|
||||
fn := addVW_g
|
||||
if len(z) > 32 {
|
||||
fn = addVWlarge
|
||||
}
|
||||
return fn(z, x, y)
|
||||
}
|
||||
|
||||
func subVW(z, x []Word, y Word) (c Word) {
|
||||
// TODO: remove indirect function call when golang.org/issue/30548 is fixed
|
||||
fn := subVW_g
|
||||
if len(z) > 32 {
|
||||
fn = subVWlarge
|
||||
}
|
||||
return fn(z, x, y)
|
||||
}
|
||||
|
||||
func shlVU(z, x []Word, s uint) (c Word) {
|
||||
return shlVU_g(z, x, s)
|
||||
}
|
||||
|
||||
func shrVU(z, x []Word, s uint) (c Word) {
|
||||
return shrVU_g(z, x, s)
|
||||
}
|
||||
|
||||
func mulAddVWW(z, x []Word, y, r Word) (c Word) {
|
||||
return mulAddVWW_g(z, x, y, r)
|
||||
}
|
||||
|
||||
func addMulVVW(z, x []Word, y Word) (c Word) {
|
||||
return addMulVVW_g(z, x, y)
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,885 @@
|
|||
// Copyright 2009 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
/*
|
||||
|
||||
Multi-precision division. Here be dragons.
|
||||
|
||||
Given u and v, where u is n+m digits, and v is n digits (with no leading zeros),
|
||||
the goal is to return quo, rem such that u = quo*v + rem, where 0 ≤ rem < v.
|
||||
That is, quo = ⌊u/v⌋ where ⌊x⌋ denotes the floor (truncation to integer) of x,
|
||||
and rem = u - quo·v.
|
||||
|
||||
|
||||
Long Division
|
||||
|
||||
Division in a computer proceeds the same as long division in elementary school,
|
||||
but computers are not as good as schoolchildren at following vague directions,
|
||||
so we have to be much more precise about the actual steps and what can happen.
|
||||
|
||||
We work from most to least significant digit of the quotient, doing:
|
||||
|
||||
• Guess a digit q, the number of v to subtract from the current
|
||||
section of u to zero out the topmost digit.
|
||||
• Check the guess by multiplying q·v and comparing it against
|
||||
the current section of u, adjusting the guess as needed.
|
||||
• Subtract q·v from the current section of u.
|
||||
• Add q to the corresponding section of the result quo.
|
||||
|
||||
When all digits have been processed, the final remainder is left in u
|
||||
and returned as rem.
|
||||
|
||||
For example, here is a sketch of dividing 5 digits by 3 digits (n=3, m=2).
|
||||
|
||||
q₂ q₁ q₀
|
||||
_________________
|
||||
v₂ v₁ v₀ ) u₄ u₃ u₂ u₁ u₀
|
||||
↓ ↓ ↓ | |
|
||||
[u₄ u₃ u₂]| |
|
||||
- [ q₂·v ]| |
|
||||
----------- ↓ |
|
||||
[ rem | u₁]|
|
||||
- [ q₁·v ]|
|
||||
----------- ↓
|
||||
[ rem | u₀]
|
||||
- [ q₀·v ]
|
||||
------------
|
||||
[ rem ]
|
||||
|
||||
Instead of creating new storage for the remainders and copying digits from u
|
||||
as indicated by the arrows, we use u's storage directly as both the source
|
||||
and destination of the subtractions, so that the remainders overwrite
|
||||
successive overlapping sections of u as the division proceeds, using a slice
|
||||
of u to identify the current section. This avoids all the copying as well as
|
||||
shifting of remainders.
|
||||
|
||||
Division of u with n+m digits by v with n digits (in base B) can in general
|
||||
produce at most m+1 digits, because:
|
||||
|
||||
• u < B^(n+m) [B^(n+m) has n+m+1 digits]
|
||||
• v ≥ B^(n-1) [B^(n-1) is the smallest n-digit number]
|
||||
• u/v < B^(n+m) / B^(n-1) [divide bounds for u, v]
|
||||
• u/v < B^(m+1) [simplify]
|
||||
|
||||
The first step is special: it takes the top n digits of u and divides them by
|
||||
the n digits of v, producing the first quotient digit and an n-digit remainder.
|
||||
In the example, q₂ = ⌊u₄u₃u₂ / v⌋.
|
||||
|
||||
The first step divides n digits by n digits to ensure that it produces only a
|
||||
single digit.
|
||||
|
||||
Each subsequent step appends the next digit from u to the remainder and divides
|
||||
those n+1 digits by the n digits of v, producing another quotient digit and a
|
||||
new n-digit remainder.
|
||||
|
||||
Subsequent steps divide n+1 digits by n digits, an operation that in general
|
||||
might produce two digits. However, as used in the algorithm, that division is
|
||||
guaranteed to produce only a single digit. The dividend is of the form
|
||||
rem·B + d, where rem is a remainder from the previous step and d is a single
|
||||
digit, so:
|
||||
|
||||
• rem ≤ v - 1 [rem is a remainder from dividing by v]
|
||||
• rem·B ≤ v·B - B [multiply by B]
|
||||
• d ≤ B - 1 [d is a single digit]
|
||||
• rem·B + d ≤ v·B - 1 [add]
|
||||
• rem·B + d < v·B [change ≤ to <]
|
||||
• (rem·B + d)/v < B [divide by v]
|
||||
|
||||
|
||||
Guess and Check
|
||||
|
||||
At each step we need to divide n+1 digits by n digits, but this is for the
|
||||
implementation of division by n digits, so we can't just invoke a division
|
||||
routine: we _are_ the division routine. Instead, we guess at the answer and
|
||||
then check it using multiplication. If the guess is wrong, we correct it.
|
||||
|
||||
How can this guessing possibly be efficient? It turns out that the following
|
||||
statement (let's call it the Good Guess Guarantee) is true.
|
||||
|
||||
If
|
||||
|
||||
• q = ⌊u/v⌋ where u is n+1 digits and v is n digits,
|
||||
• q < B, and
|
||||
• the topmost digit of v = vₙ₋₁ ≥ B/2,
|
||||
|
||||
then q̂ = ⌊uₙuₙ₋₁ / vₙ₋₁⌋ satisfies q ≤ q̂ ≤ q+2. (Proof below.)
|
||||
|
||||
That is, if we know the answer has only a single digit and we guess an answer
|
||||
by ignoring the bottom n-1 digits of u and v, using a 2-by-1-digit division,
|
||||
then that guess is at least as large as the correct answer. It is also not
|
||||
too much larger: it is off by at most two from the correct answer.
|
||||
|
||||
Note that in the first step of the overall division, which is an n-by-n-digit
|
||||
division, the 2-by-1 guess uses an implicit uₙ = 0.
|
||||
|
||||
Note that using a 2-by-1-digit division here does not mean calling ourselves
|
||||
recursively. Instead, we use an efficient direct hardware implementation of
|
||||
that operation.
|
||||
|
||||
Note that because q is u/v rounded down, q·v must not exceed u: u ≥ q·v.
|
||||
If a guess q̂ is too big, it will not satisfy this test. Viewed a different way,
|
||||
the remainder r̂ for a given q̂ is u - q̂·v, which must be positive. If it is
|
||||
negative, then the guess q̂ is too big.
|
||||
|
||||
This gives us a way to compute q. First compute q̂ with 2-by-1-digit division.
|
||||
Then, while u < q̂·v, decrement q̂; this loop executes at most twice, because
|
||||
q̂ ≤ q+2.
|
||||
|
||||
|
||||
Scaling Inputs
|
||||
|
||||
The Good Guess Guarantee requires that the top digit of v (vₙ₋₁) be at least B/2.
|
||||
For example in base 10, ⌊172/19⌋ = 9, but ⌊18/1⌋ = 18: the guess is wildly off
|
||||
because the first digit 1 is smaller than B/2 = 5.
|
||||
|
||||
We can ensure that v has a large top digit by multiplying both u and v by the
|
||||
right amount. Continuing the example, if we multiply both 172 and 19 by 3, we
|
||||
now have ⌊516/57⌋, the leading digit of v is now ≥ 5, and sure enough
|
||||
⌊51/5⌋ = 10 is much closer to the correct answer 9. It would be easier here
|
||||
to multiply by 4, because that can be done with a shift. Specifically, we can
|
||||
always count the number of leading zeros i in the first digit of v and then
|
||||
shift both u and v left by i bits.
|
||||
|
||||
Having scaled u and v, the value ⌊u/v⌋ is unchanged, but the remainder will
|
||||
be scaled: 172 mod 19 is 1, but 516 mod 57 is 3. We have to divide the remainder
|
||||
by the scaling factor (shifting right i bits) when we finish.
|
||||
|
||||
Note that these shifts happen before and after the entire division algorithm,
|
||||
not at each step in the per-digit iteration.
|
||||
|
||||
Note the effect of scaling inputs on the size of the possible quotient.
|
||||
In the scaled u/v, u can gain a digit from scaling; v never does, because we
|
||||
pick the scaling factor to make v's top digit larger but without overflowing.
|
||||
If u and v have n+m and n digits after scaling, then:
|
||||
|
||||
• u < B^(n+m) [B^(n+m) has n+m+1 digits]
|
||||
• v ≥ B^n / 2 [vₙ₋₁ ≥ B/2, so vₙ₋₁·B^(n-1) ≥ B^n/2]
|
||||
• u/v < B^(n+m) / (B^n / 2) [divide bounds for u, v]
|
||||
• u/v < 2 B^m [simplify]
|
||||
|
||||
The quotient can still have m+1 significant digits, but if so the top digit
|
||||
must be a 1. This provides a different way to handle the first digit of the
|
||||
result: compare the top n digits of u against v and fill in either a 0 or a 1.
|
||||
|
||||
|
||||
Refining Guesses
|
||||
|
||||
Before we check whether u < q̂·v, we can adjust our guess to change it from
|
||||
q̂ = ⌊uₙuₙ₋₁ / vₙ₋₁⌋ into the refined guess ⌊uₙuₙ₋₁uₙ₋₂ / vₙ₋₁vₙ₋₂⌋.
|
||||
Although not mentioned above, the Good Guess Guarantee also promises that this
|
||||
3-by-2-digit division guess is more precise and at most one away from the real
|
||||
answer q. The improvement from the 2-by-1 to the 3-by-2 guess can also be done
|
||||
without n-digit math.
|
||||
|
||||
If we have a guess q̂ = ⌊uₙuₙ₋₁ / vₙ₋₁⌋ and we want to see if it also equal to
|
||||
⌊uₙuₙ₋₁uₙ₋₂ / vₙ₋₁vₙ₋₂⌋, we can use the same check we would for the full division:
|
||||
if uₙuₙ₋₁uₙ₋₂ < q̂·vₙ₋₁vₙ₋₂, then the guess is too large and should be reduced.
|
||||
|
||||
Checking uₙuₙ₋₁uₙ₋₂ < q̂·vₙ₋₁vₙ₋₂ is the same as uₙuₙ₋₁uₙ₋₂ - q̂·vₙ₋₁vₙ₋₂ < 0,
|
||||
and
|
||||
|
||||
uₙuₙ₋₁uₙ₋₂ - q̂·vₙ₋₁vₙ₋₂ = (uₙuₙ₋₁·B + uₙ₋₂) - q̂·(vₙ₋₁·B + vₙ₋₂)
|
||||
[splitting off the bottom digit]
|
||||
= (uₙuₙ₋₁ - q̂·vₙ₋₁)·B + uₙ₋₂ - q̂·vₙ₋₂
|
||||
[regrouping]
|
||||
|
||||
The expression (uₙuₙ₋₁ - q̂·vₙ₋₁) is the remainder of uₙuₙ₋₁ / vₙ₋₁.
|
||||
If the initial guess returns both q̂ and its remainder r̂, then checking
|
||||
whether uₙuₙ₋₁uₙ₋₂ < q̂·vₙ₋₁vₙ₋₂ is the same as checking r̂·B + uₙ₋₂ < q̂·vₙ₋₂.
|
||||
|
||||
If we find that r̂·B + uₙ₋₂ < q̂·vₙ₋₂, then we can adjust the guess by
|
||||
decrementing q̂ and adding vₙ₋₁ to r̂. We repeat until r̂·B + uₙ₋₂ ≥ q̂·vₙ₋₂.
|
||||
(As before, this fixup is only needed at most twice.)
|
||||
|
||||
Now that q̂ = ⌊uₙuₙ₋₁uₙ₋₂ / vₙ₋₁vₙ₋₂⌋, as mentioned above it is at most one
|
||||
away from the correct q, and we've avoided doing any n-digit math.
|
||||
(If we need the new remainder, it can be computed as r̂·B + uₙ₋₂ - q̂·vₙ₋₂.)
|
||||
|
||||
The final check u < q̂·v and the possible fixup must be done at full precision.
|
||||
For random inputs, a fixup at this step is exceedingly rare: the 3-by-2 guess
|
||||
is not often wrong at all. But still we must do the check. Note that since the
|
||||
3-by-2 guess is off by at most 1, it can be convenient to perform the final
|
||||
u < q̂·v as part of the computation of the remainder r = u - q̂·v. If the
|
||||
subtraction underflows, decremeting q̂ and adding one v back to r is enough to
|
||||
arrive at the final q, r.
|
||||
|
||||
That's the entirety of long division: scale the inputs, and then loop over
|
||||
each output position, guessing, checking, and correcting the next output digit.
|
||||
|
||||
For a 2n-digit number divided by an n-digit number (the worst size-n case for
|
||||
division complexity), this algorithm uses n+1 iterations, each of which must do
|
||||
at least the 1-by-n-digit multiplication q̂·v. That's O(n) iterations of
|
||||
O(n) time each, so O(n²) time overall.
|
||||
|
||||
|
||||
Recursive Division
|
||||
|
||||
For very large inputs, it is possible to improve on the O(n²) algorithm.
|
||||
Let's call a group of n/2 real digits a (very) “wide digit”. We can run the
|
||||
standard long division algorithm explained above over the wide digits instead of
|
||||
the actual digits. This will result in many fewer steps, but the math involved in
|
||||
each step is more work.
|
||||
|
||||
Where basic long division uses a 2-by-1-digit division to guess the initial q̂,
|
||||
the new algorithm must use a 2-by-1-wide-digit division, which is of course
|
||||
really an n-by-n/2-digit division. That's OK: if we implement n-digit division
|
||||
in terms of n/2-digit division, the recursion will terminate when the divisor
|
||||
becomes small enough to handle with standard long division or even with the
|
||||
2-by-1 hardware instruction.
|
||||
|
||||
For example, here is a sketch of dividing 10 digits by 4, proceeding with
|
||||
wide digits corresponding to two regular digits. The first step, still special,
|
||||
must leave off a (regular) digit, dividing 5 by 4 and producing a 4-digit
|
||||
remainder less than v. The middle steps divide 6 digits by 4, guaranteed to
|
||||
produce two output digits each (one wide digit) with 4-digit remainders.
|
||||
The final step must use what it has: the 4-digit remainder plus one more,
|
||||
5 digits to divide by 4.
|
||||
|
||||
q₆ q₅ q₄ q₃ q₂ q₁ q₀
|
||||
_______________________________
|
||||
v₃ v₂ v₁ v₀ ) u₉ u₈ u₇ u₆ u₅ u₄ u₃ u₂ u₁ u₀
|
||||
↓ ↓ ↓ ↓ ↓ | | | | |
|
||||
[u₉ u₈ u₇ u₆ u₅]| | | | |
|
||||
- [ q₆q₅·v ]| | | | |
|
||||
----------------- ↓ ↓ | | |
|
||||
[ rem |u₄ u₃]| | |
|
||||
- [ q₄q₃·v ]| | |
|
||||
-------------------- ↓ ↓ |
|
||||
[ rem |u₂ u₁]|
|
||||
- [ q₂q₁·v ]|
|
||||
-------------------- ↓
|
||||
[ rem |u₀]
|
||||
- [ q₀·v ]
|
||||
------------------
|
||||
[ rem ]
|
||||
|
||||
An alternative would be to look ahead to how well n/2 divides into n+m and
|
||||
adjust the first step to use fewer digits as needed, making the first step
|
||||
more special to make the last step not special at all. For example, using the
|
||||
same input, we could choose to use only 4 digits in the first step, leaving
|
||||
a full wide digit for the last step:
|
||||
|
||||
q₆ q₅ q₄ q₃ q₂ q₁ q₀
|
||||
_______________________________
|
||||
v₃ v₂ v₁ v₀ ) u₉ u₈ u₇ u₆ u₅ u₄ u₃ u₂ u₁ u₀
|
||||
↓ ↓ ↓ ↓ | | | | | |
|
||||
[u₉ u₈ u₇ u₆]| | | | | |
|
||||
- [ q₆·v ]| | | | | |
|
||||
-------------- ↓ ↓ | | | |
|
||||
[ rem |u₅ u₄]| | | |
|
||||
- [ q₅q₄·v ]| | | |
|
||||
-------------------- ↓ ↓ | |
|
||||
[ rem |u₃ u₂]| |
|
||||
- [ q₃q₂·v ]| |
|
||||
-------------------- ↓ ↓
|
||||
[ rem |u₁ u₀]
|
||||
- [ q₁q₀·v ]
|
||||
---------------------
|
||||
[ rem ]
|
||||
|
||||
Today, the code in divRecursiveStep works like the first example. Perhaps in
|
||||
the future we will make it work like the alternative, to avoid a special case
|
||||
in the final iteration.
|
||||
|
||||
Either way, each step is a 3-by-2-wide-digit division approximated first by
|
||||
a 2-by-1-wide-digit division, just as we did for regular digits in long division.
|
||||
Because the actual answer we want is a 3-by-2-wide-digit division, instead of
|
||||
multiplying q̂·v directly during the fixup, we can use the quick refinement
|
||||
from long division (an n/2-by-n/2 multiply) to correct q to its actual value
|
||||
and also compute the remainder (as mentioned above), and then stop after that,
|
||||
never doing a full n-by-n multiply.
|
||||
|
||||
Instead of using an n-by-n/2-digit division to produce n/2 digits, we can add
|
||||
(not discard) one more real digit, doing an (n+1)-by-(n/2+1)-digit division that
|
||||
produces n/2+1 digits. That single extra digit tightens the Good Guess Guarantee
|
||||
to q ≤ q̂ ≤ q+1 and lets us drop long division's special treatment of the first
|
||||
digit. These benefits are discussed more after the Good Guess Guarantee proof
|
||||
below.
|
||||
|
||||
|
||||
How Fast is Recursive Division?
|
||||
|
||||
For a 2n-by-n-digit division, this algorithm runs a 4-by-2 long division over
|
||||
wide digits, producing two wide digits plus a possible leading regular digit 1,
|
||||
which can be handled without a recursive call. That is, the algorithm uses two
|
||||
full iterations, each using an n-by-n/2-digit division and an n/2-by-n/2-digit
|
||||
multiplication, along with a few n-digit additions and subtractions. The standard
|
||||
n-by-n-digit multiplication algorithm requires O(n²) time, making the overall
|
||||
algorithm require time T(n) where
|
||||
|
||||
T(n) = 2T(n/2) + O(n) + O(n²)
|
||||
|
||||
which, by the Bentley-Haken-Saxe theorem, ends up reducing to T(n) = O(n²).
|
||||
This is not an improvement over regular long division.
|
||||
|
||||
When the number of digits n becomes large enough, Karatsuba's algorithm for
|
||||
multiplication can be used instead, which takes O(n^log₂3) = O(n^1.6) time.
|
||||
(Karatsuba multiplication is implemented in func karatsuba in nat.go.)
|
||||
That makes the overall recursive division algorithm take O(n^1.6) time as well,
|
||||
which is an improvement, but again only for large enough numbers.
|
||||
|
||||
It is not critical to make sure that every recursion does only two recursive
|
||||
calls. While in general the number of recursive calls can change the time
|
||||
analysis, in this case doing three calls does not change the analysis:
|
||||
|
||||
T(n) = 3T(n/2) + O(n) + O(n^log₂3)
|
||||
|
||||
ends up being T(n) = O(n^log₂3). Because the Karatsuba multiplication taking
|
||||
time O(n^log₂3) is itself doing 3 half-sized recursions, doing three for the
|
||||
division does not hurt the asymptotic performance. Of course, it is likely
|
||||
still faster in practice to do two.
|
||||
|
||||
|
||||
Proof of the Good Guess Guarantee
|
||||
|
||||
Given numbers x, y, let us break them into the quotients and remainders when
|
||||
divided by some scaling factor S, with the added constraints that the quotient
|
||||
x/y and the high part of y are both less than some limit T, and that the high
|
||||
part of y is at least half as big as T.
|
||||
|
||||
x₁ = ⌊x/S⌋ y₁ = ⌊y/S⌋
|
||||
x₀ = x mod S y₀ = y mod S
|
||||
|
||||
x = x₁·S + x₀ 0 ≤ x₀ < S x/y < T
|
||||
y = y₁·S + y₀ 0 ≤ y₀ < S T/2 ≤ y₁ < T
|
||||
|
||||
And consider the two truncated quotients:
|
||||
|
||||
q = ⌊x/y⌋
|
||||
q̂ = ⌊x₁/y₁⌋
|
||||
|
||||
We will prove that q ≤ q̂ ≤ q+2.
|
||||
|
||||
The guarantee makes no real demands on the scaling factor S: it is simply the
|
||||
magnitude of the digits cut from both x and y to produce x₁ and y₁.
|
||||
The guarantee makes only limited demands on T: it must be large enough to hold
|
||||
the quotient x/y, and y₁ must have roughly the same size.
|
||||
|
||||
To apply to the earlier discussion of 2-by-1 guesses in long division,
|
||||
we would choose:
|
||||
|
||||
S = Bⁿ⁻¹
|
||||
T = B
|
||||
x = u
|
||||
x₁ = uₙuₙ₋₁
|
||||
x₀ = uₙ₋₂...u₀
|
||||
y = v
|
||||
y₁ = vₙ₋₁
|
||||
y₀ = vₙ₋₂...u₀
|
||||
|
||||
These simpler variables avoid repeating those longer expressions in the proof.
|
||||
|
||||
Note also that, by definition, truncating division ⌊x/y⌋ satisfies
|
||||
|
||||
x/y - 1 < ⌊x/y⌋ ≤ x/y.
|
||||
|
||||
This fact will be used a few times in the proofs.
|
||||
|
||||
Proof that q ≤ q̂:
|
||||
|
||||
q̂·y₁ = ⌊x₁/y₁⌋·y₁ [by definition, q̂ = ⌊x₁/y₁⌋]
|
||||
> (x₁/y₁ - 1)·y₁ [x₁/y₁ - 1 < ⌊x₁/y₁⌋]
|
||||
= x₁ - y₁ [distribute y₁]
|
||||
|
||||
So q̂·y₁ > x₁ - y₁.
|
||||
Since q̂·y₁ is an integer, q̂·y₁ ≥ x₁ - y₁ + 1.
|
||||
|
||||
q̂ - q = q̂ - ⌊x/y⌋ [by definition, q = ⌊x/y⌋]
|
||||
≥ q̂ - x/y [⌊x/y⌋ < x/y]
|
||||
= (1/y)·(q̂·y - x) [factor out 1/y]
|
||||
≥ (1/y)·(q̂·y₁·S - x) [y = y₁·S + y₀ ≥ y₁·S]
|
||||
≥ (1/y)·((x₁ - y₁ + 1)·S - x) [above: q̂·y₁ ≥ x₁ - y₁ + 1]
|
||||
= (1/y)·(x₁·S - y₁·S + S - x) [distribute S]
|
||||
= (1/y)·(S - x₀ - y₁·S) [-x = -x₁·S - x₀]
|
||||
> -y₁·S / y [x₀ < S, so S - x₀ < 0; drop it]
|
||||
≥ -1 [y₁·S ≤ y]
|
||||
|
||||
So q̂ - q > -1.
|
||||
Since q̂ - q is an integer, q̂ - q ≥ 0, or equivalently q ≤ q̂.
|
||||
|
||||
Proof that q̂ ≤ q+2:
|
||||
|
||||
x₁/y₁ - x/y = x₁·S/y₁·S - x/y [multiply left term by S/S]
|
||||
≤ x/y₁·S - x/y [x₁S ≤ x]
|
||||
= (x/y)·(y/y₁·S - 1) [factor out x/y]
|
||||
= (x/y)·((y - y₁·S)/y₁·S) [move -1 into y/y₁·S fraction]
|
||||
= (x/y)·(y₀/y₁·S) [y - y₁·S = y₀]
|
||||
= (x/y)·(1/y₁)·(y₀/S) [factor out 1/y₁]
|
||||
< (x/y)·(1/y₁) [y₀ < S, so y₀/S < 1]
|
||||
≤ (x/y)·(2/T) [y₁ ≥ T/2, so 1/y₁ ≤ 2/T]
|
||||
< T·(2/T) [x/y < T]
|
||||
= 2 [T·(2/T) = 2]
|
||||
|
||||
So x₁/y₁ - x/y < 2.
|
||||
|
||||
q̂ - q = ⌊x₁/y₁⌋ - q [by definition, q̂ = ⌊x₁/y₁⌋]
|
||||
= ⌊x₁/y₁⌋ - ⌊x/y⌋ [by definition, q = ⌊x/y⌋]
|
||||
≤ x₁/y₁ - ⌊x/y⌋ [⌊x₁/y₁⌋ ≤ x₁/y₁]
|
||||
< x₁/y₁ - (x/y - 1) [⌊x/y⌋ > x/y - 1]
|
||||
= (x₁/y₁ - x/y) + 1 [regrouping]
|
||||
< 2 + 1 [above: x₁/y₁ - x/y < 2]
|
||||
= 3
|
||||
|
||||
So q̂ - q < 3.
|
||||
Since q̂ - q is an integer, q̂ - q ≤ 2.
|
||||
|
||||
Note that when x/y < T/2, the bounds tighten to x₁/y₁ - x/y < 1 and therefore
|
||||
q̂ - q ≤ 1.
|
||||
|
||||
Note also that in the general case 2n-by-n division where we don't know that
|
||||
x/y < T, we do know that x/y < 2T, yielding the bound q̂ - q ≤ 4. So we could
|
||||
remove the special case first step of long division as long as we allow the
|
||||
first fixup loop to run up to four times. (Using a simple comparison to decide
|
||||
whether the first digit is 0 or 1 is still more efficient, though.)
|
||||
|
||||
Finally, note that when dividing three leading base-B digits by two (scaled),
|
||||
we have T = B² and x/y < B = T/B, a much tighter bound than x/y < T.
|
||||
This in turn yields the much tighter bound x₁/y₁ - x/y < 2/B. This means that
|
||||
⌊x₁/y₁⌋ and ⌊x/y⌋ can only differ when x/y is less than 2/B greater than an
|
||||
integer. For random x and y, the chance of this is 2/B, or, for large B,
|
||||
approximately zero. This means that after we produce the 3-by-2 guess in the
|
||||
long division algorithm, the fixup loop essentially never runs.
|
||||
|
||||
In the recursive algorithm, the extra digit in (2·⌊n/2⌋+1)-by-(⌊n/2⌋+1)-digit
|
||||
division has exactly the same effect: the probability of needing a fixup is the
|
||||
same 2/B. Even better, we can allow the general case x/y < 2T and the fixup
|
||||
probability only grows to 4/B, still essentially zero.
|
||||
|
||||
|
||||
References
|
||||
|
||||
There are no great references for implementing long division; thus this comment.
|
||||
Here are some notes about what to expect from the obvious references.
|
||||
|
||||
Knuth Volume 2 (Seminumerical Algorithms) section 4.3.1 is the usual canonical
|
||||
reference for long division, but that entire series is highly compressed, never
|
||||
repeating a necessary fact and leaving important insights to the exercises.
|
||||
For example, no rationale whatsoever is given for the calculation that extends
|
||||
q̂ from a 2-by-1 to a 3-by-2 guess, nor why it reduces the error bound.
|
||||
The proof that the calculation even has the desired effect is left to exercises.
|
||||
The solutions to those exercises provided at the back of the book are entirely
|
||||
calculations, still with no explanation as to what is going on or how you would
|
||||
arrive at the idea of doing those exact calculations. Nowhere is it mentioned
|
||||
that this test extends the 2-by-1 guess into a 3-by-2 guess. The proof of the
|
||||
Good Guess Guarantee is only for the 2-by-1 guess and argues by contradiction,
|
||||
making it difficult to understand how modifications like adding another digit
|
||||
or adjusting the quotient range affects the overall bound.
|
||||
|
||||
All that said, Knuth remains the canonical reference. It is dense but packed
|
||||
full of information and references, and the proofs are simpler than many other
|
||||
presentations. The proofs above are reworkings of Knuth's to remove the
|
||||
arguments by contradiction and add explanations or steps that Knuth omitted.
|
||||
But beware of errors in older printings. Take the published errata with you.
|
||||
|
||||
Brinch Hansen's “Multiple-length Division Revisited: a Tour of the Minefield”
|
||||
starts with a blunt critique of Knuth's presentation (among others) and then
|
||||
presents a more detailed and easier to follow treatment of long division,
|
||||
including an implementation in Pascal. But the algorithm and implementation
|
||||
work entirely in terms of 3-by-2 division, which is much less useful on modern
|
||||
hardware than an algorithm using 2-by-1 division. The proofs are a bit too
|
||||
focused on digit counting and seem needlessly complex, especially compared to
|
||||
the ones given above.
|
||||
|
||||
Burnikel and Ziegler's “Fast Recursive Division” introduced the key insight of
|
||||
implementing division by an n-digit divisor using recursive calls to division
|
||||
by an n/2-digit divisor, relying on Karatsuba multiplication to yield a
|
||||
sub-quadratic run time. However, the presentation decisions are made almost
|
||||
entirely for the purpose of simplifying the run-time analysis, rather than
|
||||
simplifying the presentation. Instead of a single algorithm that loops over
|
||||
quotient digits, the paper presents two mutually-recursive algorithms, for
|
||||
2n-by-n and 3n-by-2n. The paper also does not present any general (n+m)-by-n
|
||||
algorithm.
|
||||
|
||||
The proofs in the paper are remarkably complex, especially considering that
|
||||
the algorithm is at its core just long division on wide digits, so that the
|
||||
usual long division proofs apply essentially unaltered.
|
||||
*/
|
||||
|
||||
package mybig
|
||||
|
||||
import "math/bits"
|
||||
|
||||
// div returns q, r such that q = ⌊u/v⌋ and r = u%v = u - q·v.
|
||||
// It uses z and z2 as the storage for q and r.
|
||||
func (z nat) div(z2, u, v nat) (q, r nat) {
|
||||
if len(v) == 0 {
|
||||
panic("division by zero")
|
||||
}
|
||||
|
||||
if u.cmp(v) < 0 {
|
||||
q = z[:0]
|
||||
r = z2.set(u)
|
||||
return
|
||||
}
|
||||
|
||||
if len(v) == 1 {
|
||||
// Short division: long optimized for a single-word divisor.
|
||||
// In that case, the 2-by-1 guess is all we need at each step.
|
||||
var r2 Word
|
||||
q, r2 = z.divW(u, v[0])
|
||||
r = z2.setWord(r2)
|
||||
return
|
||||
}
|
||||
|
||||
q, r = z.divLarge(z2, u, v)
|
||||
return
|
||||
}
|
||||
|
||||
// divW returns q, r such that q = ⌊x/y⌋ and r = x%y = x - q·y.
|
||||
// It uses z as the storage for q.
|
||||
// Note that y is a single digit (Word), not a big number.
|
||||
func (z nat) divW(x nat, y Word) (q nat, r Word) {
|
||||
m := len(x)
|
||||
switch {
|
||||
case y == 0:
|
||||
panic("division by zero")
|
||||
case y == 1:
|
||||
q = z.set(x) // result is x
|
||||
return
|
||||
case m == 0:
|
||||
q = z[:0] // result is 0
|
||||
return
|
||||
}
|
||||
// m > 0
|
||||
z = z.make(m)
|
||||
r = divWVW(z, 0, x, y)
|
||||
q = z.norm()
|
||||
return
|
||||
}
|
||||
|
||||
// modW returns x % d.
|
||||
/*
|
||||
func (x nat) modW(d Word) (r Word) {
|
||||
// TODO(agl): we don't actually need to store the q value.
|
||||
var q nat
|
||||
q = q.make(len(x))
|
||||
return divWVW(q, 0, x, d)
|
||||
} //*/
|
||||
|
||||
// divWVW overwrites z with ⌊x/y⌋, returning the remainder r.
|
||||
// The caller must ensure that len(z) = len(x).
|
||||
func divWVW(z []Word, xn Word, x []Word, y Word) (r Word) {
|
||||
r = xn
|
||||
if len(x) == 1 {
|
||||
qq, rr := bits.Div(uint(r), uint(x[0]), uint(y))
|
||||
z[0] = Word(qq)
|
||||
return Word(rr)
|
||||
}
|
||||
rec := reciprocalWord(y)
|
||||
for i := len(z) - 1; i >= 0; i-- {
|
||||
z[i], r = divWW(r, x[i], y, rec)
|
||||
}
|
||||
return r
|
||||
}
|
||||
|
||||
// div returns q, r such that q = ⌊uIn/vIn⌋ and r = uIn%vIn = uIn - q·vIn.
|
||||
// It uses z and u as the storage for q and r.
|
||||
// The caller must ensure that len(vIn) ≥ 2 (use divW otherwise)
|
||||
// and that len(uIn) ≥ len(vIn) (the answer is 0, uIn otherwise).
|
||||
func (z nat) divLarge(u, uIn, vIn nat) (q, r nat) {
|
||||
n := len(vIn)
|
||||
m := len(uIn) - n
|
||||
|
||||
// Scale the inputs so vIn's top bit is 1 (see “Scaling Inputs” above).
|
||||
// vIn is treated as a read-only input (it may be in use by another
|
||||
// goroutine), so we must make a copy.
|
||||
// uIn is copied to u.
|
||||
shift := nlz(vIn[n-1])
|
||||
vp := getNat(n)
|
||||
v := *vp
|
||||
shlVU(v, vIn, shift)
|
||||
u = u.make(len(uIn) + 1)
|
||||
u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift)
|
||||
|
||||
// The caller should not pass aliased z and u, since those are
|
||||
// the two different outputs, but correct just in case.
|
||||
if alias(z, u) {
|
||||
z = nil
|
||||
}
|
||||
q = z.make(m + 1)
|
||||
|
||||
// Use basic or recursive long division depending on size.
|
||||
if n < divRecursiveThreshold {
|
||||
q.divBasic(u, v)
|
||||
} else {
|
||||
q.divRecursive(u, v)
|
||||
}
|
||||
putNat(vp)
|
||||
|
||||
q = q.norm()
|
||||
|
||||
// Undo scaling of remainder.
|
||||
shrVU(u, u, shift)
|
||||
r = u.norm()
|
||||
|
||||
return q, r
|
||||
}
|
||||
|
||||
// divBasic implements long division as described above.
|
||||
// It overwrites q with ⌊u/v⌋ and overwrites u with the remainder r.
|
||||
// q must be large enough to hold ⌊u/v⌋.
|
||||
func (q nat) divBasic(u, v nat) {
|
||||
n := len(v)
|
||||
m := len(u) - n
|
||||
|
||||
qhatvp := getNat(n + 1)
|
||||
qhatv := *qhatvp
|
||||
|
||||
// Set up for divWW below, precomputing reciprocal argument.
|
||||
vn1 := v[n-1]
|
||||
rec := reciprocalWord(vn1)
|
||||
|
||||
// Compute each digit of quotient.
|
||||
for j := m; j >= 0; j-- {
|
||||
// Compute the 2-by-1 guess q̂.
|
||||
// The first iteration must invent a leading 0 for u.
|
||||
qhat := Word(_M)
|
||||
var ujn Word
|
||||
if j+n < len(u) {
|
||||
ujn = u[j+n]
|
||||
}
|
||||
|
||||
// ujn ≤ vn1, or else q̂ would be more than one digit.
|
||||
// For ujn == vn1, we set q̂ to the max digit M above.
|
||||
// Otherwise, we compute the 2-by-1 guess.
|
||||
if ujn != vn1 {
|
||||
var rhat Word
|
||||
qhat, rhat = divWW(ujn, u[j+n-1], vn1, rec)
|
||||
|
||||
// Refine q̂ to a 3-by-2 guess. See “Refining Guesses” above.
|
||||
vn2 := v[n-2]
|
||||
x1, x2 := mulWW(qhat, vn2)
|
||||
ujn2 := u[j+n-2]
|
||||
for greaterThan(x1, x2, rhat, ujn2) { // x1x2 > r̂ u[j+n-2]
|
||||
qhat--
|
||||
prevRhat := rhat
|
||||
rhat += vn1
|
||||
// If r̂ overflows, then
|
||||
// r̂ u[j+n-2]v[n-1] is now definitely > x1 x2.
|
||||
if rhat < prevRhat {
|
||||
break
|
||||
}
|
||||
// TODO(rsc): No need for a full mulWW.
|
||||
// x2 += vn2; if x2 overflows, x1++
|
||||
x1, x2 = mulWW(qhat, vn2)
|
||||
}
|
||||
}
|
||||
|
||||
// Compute q̂·v.
|
||||
qhatv[n] = mulAddVWW(qhatv[0:n], v, qhat, 0)
|
||||
qhl := len(qhatv)
|
||||
if j+qhl > len(u) && qhatv[n] == 0 {
|
||||
qhl--
|
||||
}
|
||||
|
||||
// Subtract q̂·v from the current section of u.
|
||||
// If it underflows, q̂·v > u, which we fix up
|
||||
// by decrementing q̂ and adding v back.
|
||||
c := subVV(u[j:j+qhl], u[j:], qhatv)
|
||||
if c != 0 {
|
||||
c := addVV(u[j:j+n], u[j:], v)
|
||||
// If n == qhl, the carry from subVV and the carry from addVV
|
||||
// cancel out and don't affect u[j+n].
|
||||
if n < qhl {
|
||||
u[j+n] += c
|
||||
}
|
||||
qhat--
|
||||
}
|
||||
|
||||
// Save quotient digit.
|
||||
// Caller may know the top digit is zero and not leave room for it.
|
||||
if j == m && m == len(q) && qhat == 0 {
|
||||
continue
|
||||
}
|
||||
q[j] = qhat
|
||||
}
|
||||
|
||||
putNat(qhatvp)
|
||||
}
|
||||
|
||||
// greaterThan reports whether the two digit numbers x1 x2 > y1 y2.
|
||||
// TODO(rsc): In contradiction to most of this file, x1 is the high
|
||||
// digit and x2 is the low digit. This should be fixed.
|
||||
func greaterThan(x1, x2, y1, y2 Word) bool {
|
||||
return x1 > y1 || x1 == y1 && x2 > y2
|
||||
}
|
||||
|
||||
// divRecursiveThreshold is the number of divisor digits
|
||||
// at which point divRecursive is faster than divBasic.
|
||||
const divRecursiveThreshold = 100
|
||||
|
||||
// divRecursive implements recursive division as described above.
|
||||
// It overwrites z with ⌊u/v⌋ and overwrites u with the remainder r.
|
||||
// z must be large enough to hold ⌊u/v⌋.
|
||||
// This function is just for allocating and freeing temporaries
|
||||
// around divRecursiveStep, the real implementation.
|
||||
func (z nat) divRecursive(u, v nat) {
|
||||
// Recursion depth is (much) less than 2 log₂(len(v)).
|
||||
// Allocate a slice of temporaries to be reused across recursion,
|
||||
// plus one extra temporary not live across the recursion.
|
||||
recDepth := 2 * bits.Len(uint(len(v)))
|
||||
tmp := getNat(3 * len(v))
|
||||
temps := make([]*nat, recDepth)
|
||||
|
||||
z.clear()
|
||||
z.divRecursiveStep(u, v, 0, tmp, temps)
|
||||
|
||||
// Free temporaries.
|
||||
for _, n := range temps {
|
||||
if n != nil {
|
||||
putNat(n)
|
||||
}
|
||||
}
|
||||
putNat(tmp)
|
||||
}
|
||||
|
||||
// divRecursiveStep is the actual implementation of recursive division.
|
||||
// It adds ⌊u/v⌋ to z and overwrites u with the remainder r.
|
||||
// z must be large enough to hold ⌊u/v⌋.
|
||||
// It uses temps[depth] (allocating if needed) as a temporary live across
|
||||
// the recursive call. It also uses tmp, but not live across the recursion.
|
||||
func (z nat) divRecursiveStep(u, v nat, depth int, tmp *nat, temps []*nat) {
|
||||
// u is a subsection of the original and may have leading zeros.
|
||||
// TODO(rsc): The v = v.norm() is useless and should be removed.
|
||||
// We know (and require) that v's top digit is ≥ B/2.
|
||||
u = u.norm()
|
||||
v = v.norm()
|
||||
if len(u) == 0 {
|
||||
z.clear()
|
||||
return
|
||||
}
|
||||
|
||||
// Fall back to basic division if the problem is now small enough.
|
||||
n := len(v)
|
||||
if n < divRecursiveThreshold {
|
||||
z.divBasic(u, v)
|
||||
return
|
||||
}
|
||||
|
||||
// Nothing to do if u is shorter than v (implies u < v).
|
||||
m := len(u) - n
|
||||
if m < 0 {
|
||||
return
|
||||
}
|
||||
|
||||
// We consider B digits in a row as a single wide digit.
|
||||
// (See “Recursive Division” above.)
|
||||
//
|
||||
// TODO(rsc): rename B to Wide, to avoid confusion with _B,
|
||||
// which is something entirely different.
|
||||
// TODO(rsc): Look into whether using ⌈n/2⌉ is better than ⌊n/2⌋.
|
||||
B := n / 2
|
||||
|
||||
// Allocate a nat for qhat below.
|
||||
if temps[depth] == nil {
|
||||
temps[depth] = getNat(n) // TODO(rsc): Can be just B+1.
|
||||
} else {
|
||||
*temps[depth] = temps[depth].make(B + 1)
|
||||
}
|
||||
|
||||
// Compute each wide digit of the quotient.
|
||||
//
|
||||
// TODO(rsc): Change the loop to be
|
||||
// for j := (m+B-1)/B*B; j > 0; j -= B {
|
||||
// which will make the final step a regular step, letting us
|
||||
// delete what amounts to an extra copy of the loop body below.
|
||||
j := m
|
||||
for j > B {
|
||||
// Divide u[j-B:j+n] (3 wide digits) by v (2 wide digits).
|
||||
// First make the 2-by-1-wide-digit guess using a recursive call.
|
||||
// Then extend the guess to the full 3-by-2 (see “Refining Guesses”).
|
||||
//
|
||||
// For the 2-by-1-wide-digit guess, instead of doing 2B-by-B-digit,
|
||||
// we use a (2B+1)-by-(B+1) digit, which handles the possibility that
|
||||
// the result has an extra leading 1 digit as well as guaranteeing
|
||||
// that the computed q̂ will be off by at most 1 instead of 2.
|
||||
|
||||
// s is the number of digits to drop from the 3B- and 2B-digit chunks.
|
||||
// We drop B-1 to be left with 2B+1 and B+1.
|
||||
s := (B - 1)
|
||||
|
||||
// uu is the up-to-3B-digit section of u we are working on.
|
||||
uu := u[j-B:]
|
||||
|
||||
// Compute the 2-by-1 guess q̂, leaving r̂ in uu[s:B+n].
|
||||
qhat := *temps[depth]
|
||||
qhat.clear()
|
||||
qhat.divRecursiveStep(uu[s:B+n], v[s:], depth+1, tmp, temps)
|
||||
qhat = qhat.norm()
|
||||
|
||||
// Extend to a 3-by-2 quotient and remainder.
|
||||
// Because divRecursiveStep overwrote the top part of uu with
|
||||
// the remainder r̂, the full uu already contains the equivalent
|
||||
// of r̂·B + uₙ₋₂ from the “Refining Guesses” discussion.
|
||||
// Subtracting q̂·vₙ₋₂ from it will compute the full-length remainder.
|
||||
// If that subtraction underflows, q̂·v > u, which we fix up
|
||||
// by decrementing q̂ and adding v back, same as in long division.
|
||||
|
||||
// TODO(rsc): Instead of subtract and fix-up, this code is computing
|
||||
// q̂·vₙ₋₂ and decrementing q̂ until that product is ≤ u.
|
||||
// But we can do the subtraction directly, as in the comment above
|
||||
// and in long division, because we know that q̂ is wrong by at most one.
|
||||
qhatv := tmp.make(3 * n)
|
||||
qhatv.clear()
|
||||
qhatv = qhatv.mul(qhat, v[:s])
|
||||
for i := 0; i < 2; i++ {
|
||||
e := qhatv.cmp(uu.norm())
|
||||
if e <= 0 {
|
||||
break
|
||||
}
|
||||
subVW(qhat, qhat, 1)
|
||||
c := subVV(qhatv[:s], qhatv[:s], v[:s])
|
||||
if len(qhatv) > s {
|
||||
subVW(qhatv[s:], qhatv[s:], c)
|
||||
}
|
||||
addAt(uu[s:], v[s:], 0)
|
||||
}
|
||||
if qhatv.cmp(uu.norm()) > 0 {
|
||||
panic("impossible")
|
||||
}
|
||||
c := subVV(uu[:len(qhatv)], uu[:len(qhatv)], qhatv)
|
||||
if c > 0 {
|
||||
subVW(uu[len(qhatv):], uu[len(qhatv):], c)
|
||||
}
|
||||
addAt(z, qhat, j-B)
|
||||
j -= B
|
||||
}
|
||||
|
||||
// TODO(rsc): Rewrite loop as described above and delete all this code.
|
||||
|
||||
// Now u < (v<<B), compute lower bits in the same way.
|
||||
// Choose shift = B-1 again.
|
||||
s := B - 1
|
||||
qhat := *temps[depth]
|
||||
qhat.clear()
|
||||
qhat.divRecursiveStep(u[s:].norm(), v[s:], depth+1, tmp, temps)
|
||||
qhat = qhat.norm()
|
||||
qhatv := tmp.make(3 * n)
|
||||
qhatv.clear()
|
||||
qhatv = qhatv.mul(qhat, v[:s])
|
||||
// Set the correct remainder as before.
|
||||
for i := 0; i < 2; i++ {
|
||||
if e := qhatv.cmp(u.norm()); e > 0 {
|
||||
subVW(qhat, qhat, 1)
|
||||
c := subVV(qhatv[:s], qhatv[:s], v[:s])
|
||||
if len(qhatv) > s {
|
||||
subVW(qhatv[s:], qhatv[s:], c)
|
||||
}
|
||||
addAt(u[s:], v[s:], 0)
|
||||
}
|
||||
}
|
||||
if qhatv.cmp(u.norm()) > 0 {
|
||||
panic("impossible")
|
||||
}
|
||||
c := subVV(u[0:len(qhatv)], u[0:len(qhatv)], qhatv)
|
||||
if c > 0 {
|
||||
c = subVW(u[len(qhatv):], u[len(qhatv):], c)
|
||||
}
|
||||
if c > 0 {
|
||||
panic("impossible")
|
||||
}
|
||||
|
||||
// Done!
|
||||
addAt(z, qhat.norm(), 0)
|
||||
}
|
||||
12
log/0710.md
12
log/0710.md
|
|
@ -1,12 +0,0 @@
|
|||
# 2023年7月10日会议纪要
|
||||
|
||||
## 与会人员
|
||||
柴树杉 丁尔男 温祖彤
|
||||
|
||||
## 会议时间
|
||||
21:30-22:10
|
||||
|
||||
## 会议内容
|
||||
- 提案代码回顾
|
||||
- 实现技术路线讨论：大整数运算模块参考 Go 语言 Big 包构建，争取达到标准库级别
|
||||
- 立项动员
|
||||
15
log/0712.md
15
log/0712.md
|
|
@ -1,15 +0,0 @@
|
|||
# 2023年7月12日会议纪要
|
||||
|
||||
## 与会人员
|
||||
丁尔男 温祖彤
|
||||
|
||||
## 会议时间
|
||||
19:30-21:00
|
||||
|
||||
## 会议内容
|
||||
- 凹语言中自定义类型方法的语法特点讲解：指针接收器、this
|
||||
- Big 包代码结构讲解
|
||||
- 实现难点：
|
||||
- nat 类型需要使用结构体包裹
|
||||
- 值接收器要更改为指针接收器，在此过程中需要仔细检查转换前后语义是否等价
|
||||
- 由于接收器变更为指针，能否继续采用链式调用存疑
|
||||
12
log/0714.md
12
log/0714.md
|
|
@ -1,12 +0,0 @@
|
|||
# 2023年7月14日会议纪要
|
||||
|
||||
## 与会人员
|
||||
丁尔男 温祖彤
|
||||
|
||||
## 会议时间
|
||||
20:00-20:10
|
||||
|
||||
## 会议内容
|
||||
- 确认了指针接收器方法仍可使用链式调用：方法返回 this
|
||||
- 基本确认 Big 包中 nat 类型的方法可直接转为指针接收器
|
||||
- 凹语言暂未支持 panic，对本课题的必要性待讨论
|
||||
|
|
@ -0,0 +1,241 @@
|
|||
// 大整数结构体
|
||||
type BigInt struct {
|
||||
sign :bool // 符号，true表示正数，false表示负数
|
||||
digits :[]int // 数字数组，按从低位到高位排列
|
||||
}
|
||||
|
||||
func stringToBigInt(s: string) => BigInt {
|
||||
digits := make([]int, len(s))
|
||||
for i := 0; i < len(s); i++ {
|
||||
digits[len(s)-1-i] = int(s[i] - '0')
|
||||
}
|
||||
return BigInt{true, digits}
|
||||
}
|
||||
|
||||
func intToBigInt(i: int) => BigInt {
|
||||
var digits: []int
|
||||
for i > 0 {
|
||||
digits = append(digits, i%10)
|
||||
i /= 10
|
||||
}
|
||||
return BigInt{true, digits}
|
||||
}
|
||||
|
||||
// 比较两个大整数的绝对值大小，返回-1表示a < b，0表示a = b，1表示a > b
|
||||
func compareAbs(a: BigInt, b: BigInt) => int {
|
||||
if len(a.digits) < len(b.digits) {
|
||||
return -1
|
||||
} else if len(a.digits) > len(b.digits) {
|
||||
return 1
|
||||
}
|
||||
|
||||
for i := len(a.digits) - 1; i >= 0; i-- {
|
||||
if a.digits[i] < b.digits[i] {
|
||||
return -1
|
||||
} else if a.digits[i] > b.digits[i] {
|
||||
return 1
|
||||
}
|
||||
}
|
||||
|
||||
return 0
|
||||
}
|
||||
|
||||
// 大整数加法
|
||||
func add(a: BigInt, b: BigInt) => BigInt {
|
||||
// 确保a是较大的数
|
||||
if len(a.digits) < len(b.digits) {
|
||||
a, b = b, a
|
||||
}
|
||||
|
||||
carry := 0
|
||||
result := make([]int, len(a.digits)+1)
|
||||
|
||||
for i := 0; i < len(b.digits); i++ {
|
||||
sum := a.digits[i] + b.digits[i] + carry
|
||||
result[i] = sum % 10
|
||||
carry = sum / 10
|
||||
}
|
||||
|
||||
for i := len(b.digits); i < len(a.digits); i++ {
|
||||
sum := a.digits[i] + carry
|
||||
result[i] = sum % 10
|
||||
carry = sum / 10
|
||||
}
|
||||
|
||||
if carry > 0 {
|
||||
result[len(a.digits)] = carry
|
||||
} else {
|
||||
result = result[:len(a.digits)]
|
||||
}
|
||||
|
||||
return BigInt{a.sign, result}
|
||||
}
|
||||
|
||||
// 大整数减法
|
||||
func subtract(a: BigInt, b: BigInt) => BigInt {
|
||||
// 确保a是较大的数
|
||||
if compareAbs(a, b) < 0 {
|
||||
a, b = b, a
|
||||
}
|
||||
|
||||
borrow := 0
|
||||
result := make([]int, len(a.digits))
|
||||
|
||||
for i := 0; i < len(b.digits); i++ {
|
||||
diff := a.digits[i] - b.digits[i] - borrow
|
||||
|
||||
if diff < 0 {
|
||||
diff += 10
|
||||
borrow = 1
|
||||
} else {
|
||||
borrow = 0
|
||||
}
|
||||
|
||||
result[i] = diff
|
||||
}
|
||||
|
||||
for i := len(b.digits); i < len(a.digits); i++ {
|
||||
diff := a.digits[i] - borrow
|
||||
|
||||
if diff < 0 {
|
||||
diff += 10
|
||||
borrow = 1
|
||||
} else {
|
||||
borrow = 0
|
||||
}
|
||||
|
||||
result[i] = diff
|
||||
}
|
||||
|
||||
// 去除结果中的前导零
|
||||
for len(result) > 0 && result[len(result)-1] == 0 {
|
||||
result = result[:len(result)-1]
|
||||
}
|
||||
|
||||
// 如果结果为零，则设置为正数
|
||||
if len(result) == 0 {
|
||||
return BigInt{true, result}
|
||||
}
|
||||
|
||||
return BigInt{a.sign, result}
|
||||
}
|
||||
|
||||
// 大整数乘法
|
||||
func multiply(a: BigInt, b: BigInt) => BigInt {
|
||||
result := BigInt{true, make([]int, len(a.digits)+len(b.digits))}
|
||||
|
||||
for i := 0; i < len(a.digits); i++ {
|
||||
carry := 0
|
||||
|
||||
for j := 0; j < len(b.digits); j++ {
|
||||
product := a.digits[i]*b.digits[j] + result.digits[i+j] + carry
|
||||
result.digits[i+j] = product % 10
|
||||
carry = product / 10
|
||||
}
|
||||
|
||||
result.digits[i+len(b.digits)] = carry
|
||||
}
|
||||
|
||||
// 去除结果中的前导零
|
||||
for len(result.digits) > 0 && result.digits[len(result.digits)-1] == 0 {
|
||||
result.digits = result.digits[:len(result.digits)-1]
|
||||
}
|
||||
|
||||
return result
|
||||
}
|
||||
|
||||
// 小整数除法
|
||||
func divmodint(a: BigInt, b: int) => (quotient: BigInt, remainder: int) {
|
||||
quotientDigits := make([]int, len(a.digits))
|
||||
remainder = 0
|
||||
|
||||
for i := len(a.digits) - 1; i >= 0; i-- {
|
||||
dividend := remainder*10 + a.digits[i]
|
||||
quotientDigits[i] = dividend / b
|
||||
remainder = dividend % b
|
||||
}
|
||||
|
||||
// 去除商中的前导零
|
||||
for len(quotientDigits) > 0 && quotientDigits[len(quotientDigits)-1] == 0 {
|
||||
quotientDigits = quotientDigits[:len(quotientDigits)-1]
|
||||
}
|
||||
|
||||
quotient = BigInt{a.sign, quotientDigits}
|
||||
return quotient, remainder
|
||||
}
|
||||
|
||||
// 大整数除法
|
||||
func divide(a: BigInt, b: BigInt) => (quotient: BigInt, remainder: BigInt) {
|
||||
// 处理特殊情况：除数为零
|
||||
if len(b.digits) == 0 {
|
||||
// panic("division by zero")
|
||||
println("division by zero")
|
||||
}
|
||||
|
||||
// 处理特殊情况：被除数为零
|
||||
if len(a.digits) == 0 {
|
||||
return BigInt{true, []int{}}, BigInt{true, []int{}}
|
||||
}
|
||||
|
||||
// 处理特殊情况：被除数小于除数
|
||||
if compareAbs(a, b) < 0 {
|
||||
return BigInt{true, []int{}}, a
|
||||
}
|
||||
|
||||
quotientDigits := make([]int, len(a.digits))
|
||||
remainderDigits := make([]int, len(a.digits))
|
||||
|
||||
divisor := b.digits[len(b.digits)-1]
|
||||
dividend := append(a.digits, 0)
|
||||
|
||||
for i := len(a.digits) - len(b.digits); i >= 0; i-- {
|
||||
quotientDigit := (dividend[i+len(b.digits)]*10 + dividend[i+len(b.digits)-1]) / divisor
|
||||
remainderDigit := (dividend[i+len(b.digits)]*10 + dividend[i+len(b.digits)-1]) % divisor
|
||||
|
||||
for quotientDigit > 0 && (quotientDigit*b.digits[len(b.digits)-2]) > (remainderDigit*10+dividend[i+len(b.digits)-2]) {
|
||||
quotientDigit--
|
||||
remainderDigit += divisor
|
||||
}
|
||||
|
||||
// 更新被除数
|
||||
for j := len(b.digits) - 1; j >= 0; j-- {
|
||||
dividend[i+j] -= quotientDigit * b.digits[j]
|
||||
}
|
||||
|
||||
if dividend[i+len(b.digits)-1] < 0 {
|
||||
// 借位
|
||||
dividend[i+len(b.digits)-1] += 10
|
||||
dividend[i+len(b.digits)-2]--
|
||||
}
|
||||
|
||||
quotientDigits[i] = quotientDigit
|
||||
remainderDigits[i] = remainderDigit
|
||||
}
|
||||
|
||||
// 去除商中的前导零
|
||||
for len(quotientDigits) > 0 && quotientDigits[len(quotientDigits)-1] == 0 {
|
||||
quotientDigits = quotientDigits[:len(quotientDigits)-1]
|
||||
}
|
||||
|
||||
// 去除余数中的前导零
|
||||
for len(remainderDigits) > 0 && remainderDigits[len(remainderDigits)-1] == 0 {
|
||||
remainderDigits = remainderDigits[:len(remainderDigits)-1]
|
||||
}
|
||||
|
||||
return BigInt{a.sign == b.sign, quotientDigits}, BigInt{a.sign, remainderDigits}
|
||||
}
|
||||
|
||||
// 将大整数转换为字符串
|
||||
func String(a: BigInt) => string {
|
||||
var digits: string
|
||||
|
||||
for i := 0; i < len(a.digits); i++ {
|
||||
digits = intToString(a.digits[i]) + digits
|
||||
}
|
||||
|
||||
if len(digits) == 0 {
|
||||
return "0"
|
||||
}
|
||||
|
||||
return digits
|
||||
}
|
||||
|
|
@ -0,0 +1,196 @@
|
|||
// 版权 @2023 waExample/main 作者。保留所有权利。
|
||||
|
||||
func main {
|
||||
x := EncodeBase60(stringToByteArray("你好"))
|
||||
println(x)
|
||||
c := DecodeBase60(x)
|
||||
printlnByteArray(c)
|
||||
// println(string(c))
|
||||
}
|
||||
|
||||
func getTianganDizhiByIndex(index: int) => string {
|
||||
var tiangandizhi = []string{
|
||||
"甲子", "乙丑", "丙寅", "丁卯", "戊辰", "己巳", "庚午", "辛未", "壬申", "癸酉",
|
||||
"甲戌", "乙亥", "丙子", "丁丑", "戊寅", "己卯", "庚辰", "辛巳", "壬午", "癸未",
|
||||
"甲申", "乙酉", "丙戌", "丁亥", "戊子", "己丑", "庚寅", "辛卯", "壬辰", "癸巳",
|
||||
"甲午", "乙未", "丙申", "丁酉", "戊戌", "己亥", "庚子", "辛丑", "壬寅", "癸卯",
|
||||
"甲辰", "乙巳", "丙午", "丁未", "戊申", "己酉", "庚戌", "辛亥", "壬子", "癸丑",
|
||||
"甲寅", "乙卯", "丙辰", "丁巳", "戊午", "己未", "庚申", "辛酉", "壬戌", "癸亥",
|
||||
}
|
||||
return tiangandizhi[index]
|
||||
}
|
||||
func getIndexByTianganDizhi(s: string) => int {
|
||||
var tiangandizhi = []string{
|
||||
"甲子", "乙丑", "丙寅", "丁卯", "戊辰", "己巳", "庚午", "辛未", "壬申", "癸酉",
|
||||
"甲戌", "乙亥", "丙子", "丁丑", "戊寅", "己卯", "庚辰", "辛巳", "壬午", "癸未",
|
||||
"甲申", "乙酉", "丙戌", "丁亥", "戊子", "己丑", "庚寅", "辛卯", "壬辰", "癸巳",
|
||||
"甲午", "乙未", "丙申", "丁酉", "戊戌", "己亥", "庚子", "辛丑", "壬寅", "癸卯",
|
||||
"甲辰", "乙巳", "丙午", "丁未", "戊申", "己酉", "庚戌", "辛亥", "壬子", "癸丑",
|
||||
"甲寅", "乙卯", "丙辰", "丁巳", "戊午", "己未", "庚申", "辛酉", "壬戌", "癸亥",
|
||||
}
|
||||
for i := 0; i < len(tiangandizhi); i++ {
|
||||
if s == tiangandizhi[i] {
|
||||
return i
|
||||
}
|
||||
}
|
||||
return -1
|
||||
}
|
||||
func printlnByteArray(byteArray: []byte) {
|
||||
print('[')
|
||||
for i := 0; i < len(byteArray); i++ {
|
||||
if (i != 0) {print(',')}
|
||||
print(int(byteArray[i]))
|
||||
|
||||
}
|
||||
print(']')
|
||||
println()
|
||||
}
|
||||
func stringToByteArray(str: string) => []byte {
|
||||
var result: []byte
|
||||
for i := 0; i < len(str); i++ {
|
||||
result = append(result, byte(str[i]))
|
||||
}
|
||||
return result
|
||||
}
|
||||
|
||||
func intToString(i: int) => string {
|
||||
var tr = []string{"0", "1", "2", "3", "4", "5", "6", "7", "8", "9"}
|
||||
var res: string
|
||||
for i > 0 {
|
||||
res = tr[i%10] + res
|
||||
i = i / 10
|
||||
}
|
||||
if len(res) == 0 {
|
||||
return "0"
|
||||
} else {
|
||||
return res
|
||||
}
|
||||
}
|
||||
|
||||
func byteArrayToString(byteArray: []byte) => string {
|
||||
var result: string
|
||||
var tr = []string{
|
||||
"0x00",
|
||||
"0x01", "0x02", "0x03", "0x04", "0x05", "0x06", "0x07", "0x08",
|
||||
"0x09", "0x0a", "0x0b", "0x0c", "0x0d", "0x0e", "0x0f", "0x10",
|
||||
"0x11", "0x12", "0x13", "0x14", "0x15", "0x16", "0x17", "0x18",
|
||||
"0x19", "0x1a", "0x1b", "0x1c", "0x1d", "0x1e", "0x1f", "0x20",
|
||||
"0x21", "0x22", "0x23", "0x24", "0x25", "0x26", "0x27", "0x28",
|
||||
"0x29", "0x2a", "0x2b", "0x2c", "0x2d", "0x2e", "0x2f", "0x30",
|
||||
"0x31", "0x32", "0x33", "0x34", "0x35", "0x36", "0x37", "0x38",
|
||||
"0x39", "0x3a", "0x3b", "0x3c", "0x3d", "0x3e", "0x3f", "0x40",
|
||||
"0x41", "0x42", "0x43", "0x44", "0x45", "0x46", "0x47", "0x48",
|
||||
"0x49", "0x4a", "0x4b", "0x4c", "0x4d", "0x4e", "0x4f", "0x50",
|
||||
"0x51", "0x52", "0x53", "0x54", "0x55", "0x56", "0x57", "0x58",
|
||||
"0x59", "0x5a", "0x5b", "0x5c", "0x5d", "0x5e", "0x5f", "0x60",
|
||||
"0x61", "0x62", "0x63", "0x64", "0x65", "0x66", "0x67", "0x68",
|
||||
"0x69", "0x6a", "0x6b", "0x6c", "0x6d", "0x6e", "0x6f", "0x70",
|
||||
"0x71", "0x72", "0x73", "0x74", "0x75", "0x76", "0x77", "0x78",
|
||||
"0x79", "0x7a", "0x7b", "0x7c", "0x7d", "0x7e", "0x7f", "0x80",
|
||||
"0x81", "0x82", "0x83", "0x84", "0x85", "0x86", "0x87", "0x88",
|
||||
"0x89", "0x8a", "0x8b", "0x8c", "0x8d", "0x8e", "0x8f", "0x90",
|
||||
"0x91", "0x92", "0x93", "0x94", "0x95", "0x96", "0x97", "0x98",
|
||||
"0x99", "0x9a", "0x9b", "0x9c", "0x9d", "0x9e", "0x9f", "0xa0",
|
||||
"0xa1", "0xa2", "0xa3", "0xa4", "0xa5", "0xa6", "0xa7", "0xa8",
|
||||
"0xa9", "0xaa", "0xab", "0xac", "0xad", "0xae", "0xaf", "0xb0",
|
||||
// 一次性生成后面所有行
|
||||
"0xb1", "0xb2", "0xb3", "0xb4", "0xb5", "0xb6", "0xb7", "0xb8",
|
||||
"0xb9", "0xba", "0xbb", "0xbc", "0xbd", "0xbe", "0xbf", "0xc0",
|
||||
"0xc1", "0xc2", "0xc3", "0xc4", "0xc5", "0xc6", "0xc7", "0xc8",
|
||||
"0xc9", "0xca", "0xcb", "0xcc", "0xcd", "0xce", "0xcf", "0xd0",
|
||||
"0xd1", "0xd2", "0xd3", "0xd4", "0xd5", "0xd6", "0xd7", "0xd8",
|
||||
"0xd9", "0xda", "0xdb", "0xdc", "0xdd", "0xde", "0xdf", "0xe0",
|
||||
"0xe1", "0xe2", "0xe3", "0xe4", "0xe5", "0xe6", "0xe7", "0xe8",
|
||||
"0xe9", "0xea", "0xeb", "0xec", "0xed", "0xee", "0xef", "0xf0",
|
||||
"0xf1", "0xf2", "0xf3", "0xf4", "0xf5", "0xf6", "0xf7", "0xf8",
|
||||
"0xf9", "0xfa", "0xfb", "0xfc", "0xfd", "0xfe", "0xff",
|
||||
}
|
||||
for i := 0; i < len(byteArray); i++ {
|
||||
result += tr[i]
|
||||
}
|
||||
return result
|
||||
}
|
||||
|
||||
func Encode(str: string) => string {
|
||||
var base64 = []string{
|
||||
"A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z", "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "+", "/",
|
||||
}
|
||||
|
||||
var result: string
|
||||
length := len(str)
|
||||
|
||||
// Convert string to byte array
|
||||
|
||||
// Encode
|
||||
for i := 0; i < length; i += 3 {
|
||||
var buffer: [4]int
|
||||
buffer[0] = int(str[i]) / 4
|
||||
buffer[1] = (int(str[i]) % 4) * 16
|
||||
if length-i > 1 {
|
||||
buffer[1] += int(str[i+1]) / 16
|
||||
buffer[2] += (int(str[i+1]) % 16) * 4
|
||||
if length-i > 2 {
|
||||
buffer[2] += int(str[i+2]) / 64
|
||||
buffer[3] += int(str[i+2]) % 64
|
||||
}
|
||||
}
|
||||
result += base64[buffer[0]]
|
||||
result += base64[buffer[1]]
|
||||
if length-i > 1 {
|
||||
result += base64[buffer[2]]
|
||||
}
|
||||
if length-i > 2 {
|
||||
result += base64[buffer[3]]
|
||||
}
|
||||
}
|
||||
switch length % 3 {
|
||||
case 1:
|
||||
result += "=="
|
||||
case 2:
|
||||
result += "="
|
||||
}
|
||||
return result
|
||||
}
|
||||
|
||||
func EncodeBase60(data: []byte) => string {
|
||||
printlnByteArray(data)
|
||||
x := intToBigInt(0)
|
||||
baseLength := intToBigInt(256)
|
||||
for i := 0; i < len(data); i++ {
|
||||
x = multiply(x, baseLength)
|
||||
x = add(x, intToBigInt(int(data[i])))
|
||||
}
|
||||
println(String(x))
|
||||
zero := intToBigInt(0)
|
||||
var result: string
|
||||
for compareAbs(x, zero) > 0 {
|
||||
quotient, remainder := divmodint(x, 60)
|
||||
result = getTianganDizhiByIndex(remainder) + result
|
||||
x = quotient
|
||||
}
|
||||
return result
|
||||
}
|
||||
|
||||
func DecodeBase60(s: string) => []byte {
|
||||
x := intToBigInt(0)
|
||||
baseLength := intToBigInt(60)
|
||||
for i := 0; i < len(s); i += 6 {
|
||||
k := getIndexByTianganDizhi(string(s[i : i+6]))
|
||||
if k == -1 {
|
||||
println("Decode Error!")
|
||||
return []byte{0}
|
||||
}
|
||||
x = multiply(x, baseLength)
|
||||
x = add(x, intToBigInt(k))
|
||||
// println(String(x))
|
||||
}
|
||||
println(String(x))
|
||||
var res: []byte
|
||||
zero := intToBigInt(0)
|
||||
for compareAbs(x, zero) > 0 {
|
||||
quotient, remainder := divmodint(x, 256)
|
||||
res = append([]byte{byte(remainder)}, res...)
|
||||
x = quotient
|
||||
}
|
||||
return res
|
||||
}
|
||||
|
|
@ -0,0 +1,16 @@
|
|||
// 版权 @2023 waExample/main 作者。保留所有权利。
|
||||
|
||||
// 打印素数
|
||||
func PrintPrime(max: int) {
|
||||
for n := 2; n <= max; n = n + 1 {
|
||||
isPrime := 1
|
||||
for i := 2; i*i <= n; i = i + 1 {
|
||||
if x := n % i; x == 0 {
|
||||
isPrime = 0
|
||||
}
|
||||
}
|
||||
if isPrime != 0 {
|
||||
println(n)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
// 版权 @2023 waExample/main 作者。保留所有权利。
|
||||
|
||||
// 打印整数
|
||||
func Println(x: int) {
|
||||
println(x)
|
||||
}
|
||||
|
|
@ -0,0 +1,5 @@
|
|||
// 版权 @2023 waExample/main 作者。保留所有权利。
|
||||
|
||||
func Println(x: int) {
|
||||
println(x)
|
||||
}
|
||||
|
|
@ -1,60 +0,0 @@
|
|||
type nat struct {
|
||||
data: []uint
|
||||
}
|
||||
|
||||
func nat.clear() {
|
||||
for i := range this.data {
|
||||
this.data[i] = 0
|
||||
}
|
||||
}
|
||||
|
||||
var (
|
||||
natOne = nat{data: []uint{1}}
|
||||
)
|
||||
|
||||
func nat.raw_print() {
|
||||
print("raw_nat: {")
|
||||
for i := range this.data {
|
||||
print(this.data[i])
|
||||
if i != len(this.data)-1 {
|
||||
print(",")
|
||||
}
|
||||
}
|
||||
print("}")
|
||||
}
|
||||
|
||||
func nat.norm() => nat {
|
||||
i := len(this.data)
|
||||
for i > 0 && this.data[i-1] == 0 {
|
||||
i--
|
||||
}
|
||||
return nat{data: this.data[0:i]}
|
||||
}
|
||||
|
||||
func nat.make(n: int) => nat {
|
||||
if n <= cap(this.data) {
|
||||
return nat{data: this.data[:n]} // reuse this.data
|
||||
}
|
||||
if n == 1 {
|
||||
// Most nats start small and stay that way; don't over-allocate.
|
||||
return nat{data: make([]uint, 1)}
|
||||
}
|
||||
// Choosing a good value for e has significant performance impact
|
||||
// because it increases the chance that a value can be reused.
|
||||
const e = 4 // extra capacity
|
||||
return nat{data: make([]uint, n, n+e)}
|
||||
}
|
||||
|
||||
func Test() {
|
||||
testNorm()
|
||||
}
|
||||
|
||||
func testNorm() {
|
||||
t := nat{data: []uint{1, 0, 0, 0}}
|
||||
t = t.norm()
|
||||
if len(t.data) == 1 && t.data[0] == 1 {
|
||||
println("NORM TEST PASS!")
|
||||
} else {
|
||||
println("NORM TEST FAILED!")
|
||||
}
|
||||
}
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
{
|
||||
"name": "wa_bigint_std",
|
||||
"name": "waExample/main",
|
||||
"pkgpath": "myapp",
|
||||
"version": "0.0.1",
|
||||
"authors": ["author", "author2"],
|
||||
|
|
|
|||
Loading…
Reference in New Issue