referrerpolicy=no-referrer-when-downgrade
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

// This file is part of Substrate.

// Copyright (C) Parity Technologies (UK) Ltd.
// SPDX-License-Identifier: Apache-2.0

// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// 	http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

/// Simple u32 power of 2 function - simply uses a bit shift
macro_rules! pow2 {
	($n:expr) => {
		1_u32 << $n
	};
}

/// Returns the k_th per_million taylor term for a log2 function
fn taylor_term(k: u32, y_num: u128, y_den: u128) -> u32 {
	let _2_div_ln_2: u128 = 2_885_390u128;

	if k == 0 {
		(_2_div_ln_2 * (y_num).pow(1) / (y_den).pow(1)).try_into().unwrap()
	} else {
		let mut res = _2_div_ln_2 * (y_num).pow(3) / (y_den).pow(3);
		for _ in 1..k {
			res = res * (y_num).pow(2) / (y_den).pow(2);
		}
		res /= 2 * k as u128 + 1;

		res.try_into().unwrap()
	}
}

/// Performs a log2 operation using a rational fraction
///
/// result = log2(p/q) where p/q is bound to [1, 1_000_000]
/// Where:
/// * q represents the numerator of the rational fraction input
/// * p represents the denominator of the rational fraction input
/// * result represents a per-million output of log2
pub fn log2(p: u32, q: u32) -> u32 {
	assert!(p >= q); // keep p/q bound to [1, inf)
	assert!(p <= u32::MAX / 2);

	// This restriction should not be mandatory. But function is only tested and used for this.
	assert!(p <= 1_000_000);
	assert!(q <= 1_000_000);

	// log2(1) = 0
	if p == q {
		return 0
	}

	// find the power of 2 where q * 2^n <= p < q * 2^(n+1)
	let mut n = 0u32;
	while (p < pow2!(n) * q) || (p >= pow2!(n + 1) * q) {
		n += 1;
		assert!(n < 32); // cannot represent 2^32 in u32
	}
	assert!(p < pow2!(n + 1) * q);

	let y_num: u32 = p - pow2!(n) * q;
	let y_den: u32 = p + pow2!(n) * q;

	// Loop through each Taylor series coefficient until it reaches 10^-6
	let mut res = n * 1_000_000u32;
	let mut k = 0;
	loop {
		let term = taylor_term(k, y_num.into(), y_den.into());
		if term == 0 {
			break
		}

		res += term;
		k += 1;
	}

	res
}

#[test]
fn test_log() {
	let div = 1_000;
	for p in 0..=div {
		for q in 1..=p {
			let p: u32 = (1_000_000 as u64 * p as u64 / div as u64).try_into().unwrap();
			let q: u32 = (1_000_000 as u64 * q as u64 / div as u64).try_into().unwrap();

			let res = -(log2(p, q) as i64);
			let expected = ((q as f64 / p as f64).log(2.0) * 1_000_000 as f64).round() as i64;
			assert!((res - expected).abs() <= 6);
		}
	}
}

#[test]
#[should_panic]
fn test_log_p_must_be_greater_than_q() {
	let p: u32 = 1_000;
	let q: u32 = 1_001;
	let _ = log2(p, q);
}

#[test]
#[should_panic]
fn test_log_p_upper_bound() {
	let p: u32 = 1_000_001;
	let q: u32 = 1_000_000;
	let _ = log2(p, q);
}

#[test]
#[should_panic]
fn test_log_q_limit() {
	let p: u32 = 1_000_000;
	let q: u32 = 0;
	let _ = log2(p, q);
}

#[test]
fn test_log_of_one_boundary() {
	let p: u32 = 1_000_000;
	let q: u32 = 1_000_000;
	assert_eq!(log2(p, q), 0);
}

#[test]
fn test_log_of_largest_input() {
	let p: u32 = 1_000_000;
	let q: u32 = 1;
	let expected = 19_931_568;
	let tolerance = 100;
	assert!((log2(p, q) as i32 - expected as i32).abs() < tolerance);
}