StatProfilerHTML.jl report
Generated on tor 10 okt 2019 11:38:33
File source code
Line Exclusive Inclusive Code
1 # This file is a part of Julia, but is derived from
2 # https://github.com/google/double-conversion which has the following license
3 #
4 # Copyright 2006-2014, the V8 project authors. All rights reserved.
5 # Redistribution and use in source and binary forms, with or without
6 # modification, are permitted provided that the following conditions are
7 # met:
8 #
9 # * Redistributions of source code must retain the above copyright
10 # notice, this list of conditions and the following disclaimer.
11 # * Redistributions in binary form must reproduce the above
12 # copyright notice, this list of conditions and the following
13 # disclaimer in the documentation and/or other materials provided
14 # with the distribution.
15 # * Neither the name of Google Inc. nor the names of its
16 # contributors may be used to endorse or promote products derived
17 # from this software without specific prior written permission.
18 #
19 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20 # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21 # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22 # A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
23 # OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
24 # SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
25 # LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
26 # DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
27 # THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
28 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
29 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30
31 module Bignums
32
33 import Base: ==, <
34
35 export Bignum
36
37 const kMaxSignificantBits = 3584
38
39 const Chunk = UInt32
40 const DoubleChunk = UInt64
41
42 const kChunkSize = sizeof(Chunk) * 8
43 const kDoubleChunkSize = sizeof(DoubleChunk) * 8
44 # With bigit size of 28 we loose some bits, but a double still fits easily
45 # into two chunks, and more importantly we can use the Comba multiplication.
46 const kBigitSize = 28
47 const kBigitMask = Chunk((1 << kBigitSize) - 1)
48 # Every instance allocates kBigitLength chunks on the stack. Bignums cannot
49 # grow. There are no checks if the stack-allocated space is sufficient.
50 const kBigitCapacity = div(kMaxSignificantBits, kBigitSize)
51
52 mutable struct Bignum
53 bigits::Vector{UInt32}
54 used_digits::Int32
55 exponent::Int32
56 function Bignum()
57 bigits = Vector{UInt32}(undef, kBigitCapacity)
58 @inbounds for i = 1:kBigitCapacity
59 bigits[i] = 0
60 end
61 new(bigits,0,0)
62 end
63 end
64
65 ==(a::Bignum,b::Bignum) = compare(a,b) == 0
66 <(a::Bignum,b::Bignum) = compare(a,b) < 0
67
68 times10!(x::Bignum) = multiplybyuint32!(x,UInt32(10))
69
70 plusequal(a,b,c) = pluscompare(a,b,c) == 0
71 pluslessequal(a,b,c) = pluscompare(a,b,c) <= 0
72 plusless(a,b,c) = pluscompare(a,b,c) < 0
73 lessequal(a::Bignum,b::Bignum) = compare(a,b) <= 0
74 less(a::Bignum,b::Bignum) = compare(a,b) < 0
75
76 bigitlength(x::Bignum) = x.used_digits + x.exponent
77
78 bitsize(value) = 8 * sizeof(value)
79
80 function zero!(x::Bignum)
81 for i = 1:x.used_digits
82 @inbounds x.bigits[i] = 0
83 end
84 x.used_digits = 0
85 x.exponent = 0
86 return
87 end
88
89 function clamp!(x::Bignum)
90 @inbounds while (x.used_digits > 0 && x.bigits[x.used_digits] == 0)
91 x.used_digits -= 1
92 end
93 x.used_digits == 0 && (x.exponent = 0)
94 return
95 end
96
97 isclamped(x::Bignum) = x.used_digits == 0 || x.bigits[x.used_digits] != 0
98
99 function align!(x::Bignum,other::Bignum)
100 @inbounds if x.exponent > other.exponent
101 zero_digits = x.exponent - other.exponent
102 for i = x.used_digits:-1:1
103 x.bigits[i + zero_digits] = x.bigits[i]
104 end
105 for i = 1:zero_digits
106 x.bigits[i] = 0
107 end
108 x.used_digits += zero_digits
109 x.exponent -= zero_digits
110 end
111 return
112 end
113
114 function bigitshiftleft!(x::Bignum,shift_amount)
115 carry::UInt32 = 0
116 @inbounds begin
117 for i = 1:x.used_digits
118 new_carry::Chunk = x.bigits[i] >> (kBigitSize - shift_amount)
119 x.bigits[i] = ((x.bigits[i] << shift_amount) + carry) & kBigitMask
120 carry = new_carry
121 end
122 if carry != 0
123 x.bigits[x.used_digits+1] = carry
124 x.used_digits += 1
125 end
126 end
127 return
128 end
129
130 function subtracttimes!(x::Bignum,other::Bignum,factor)
131 if factor < 3
132 for i = 1:factor
133 subtractbignum!(x,other)
134 end
135 return
136 end
137 borrow::Chunk = 0
138 exponent_diff = other.exponent - x.exponent
139 @inbounds begin
140 for i = 1:other.used_digits
141 product::DoubleChunk = DoubleChunk(factor) * other.bigits[i]
142 remove::DoubleChunk = borrow + product
143 difference::Chunk = (x.bigits[i+exponent_diff] - (remove & kBigitMask)) % Chunk
144 x.bigits[i+exponent_diff] = difference & kBigitMask
145 borrow = ((difference >> (kChunkSize - 1)) + (remove >> kBigitSize)) % Chunk
146 end
147 for i = (other.used_digits + exponent_diff + 1):x.used_digits
148 borrow == 0 && return
149 difference::Chunk = x.bigits[i] - borrow
150 x.bigits[i] = difference & kBigitMask
151 borrow = difference >> (kChunkSize - 1)
152 end
153 end
154 clamp!(x)
155 end
156
157 function assignuint16!(x::Bignum,value::UInt16)
158 zero!(x)
159 value == 0 && return
160 x.bigits[1] = value
161 x.used_digits = 1
162 return
163 end
164
165 const kUInt64Size = 64
166 function assignuint64!(x::Bignum,value::UInt64)
167 zero!(x)
168 value == 0 && return
169 needed_bigits = div(kUInt64Size,kBigitSize) + 1
170 @inbounds for i = 1:needed_bigits
171 x.bigits[i] = value & kBigitMask
172 value >>= kBigitSize
173 end
174 x.used_digits = needed_bigits
175 clamp!(x)
176 end
177
178 function assignbignum!(x::Bignum,other::Bignum)
179 x.exponent = other.exponent
180 @inbounds begin
181 for i = 1:other.used_digits
182 x.bigits[i] = other.bigits[i]
183 end
184 for i = (other.used_digits+1):x.used_digits
185 x.bigits[i] = 0
186 end
187 end
188 x.used_digits = other.used_digits
189 return
190 end
191
192 function adduint64!(x::Bignum,operand::UInt64)
193 operand == 0 && return
194 other = Bignum()
195 assignuint64!(other,operand)
196 addbignum!(x,other)
197 end
198
199 function addbignum!(x::Bignum,other::Bignum)
200 align!(x,other)
201 carry::Chunk = 0
202 bigit_pos = other.exponent - x.exponent
203 @inbounds for i = 1:other.used_digits
204 sum::Chunk = x.bigits[bigit_pos+1] + other.bigits[i] + carry
205 x.bigits[bigit_pos+1] = sum & kBigitMask
206 carry = sum >> kBigitSize
207 bigit_pos += 1
208 end
209 @inbounds while carry != 0
210 sum = x.bigits[bigit_pos+1] + carry
211 x.bigits[bigit_pos+1] = sum & kBigitMask
212 carry = sum >> kBigitSize
213 bigit_pos += 1
214 end
215 x.used_digits = max(bigit_pos,x.used_digits)
216 return
217 end
218
219 function subtractbignum!(x::Bignum,other::Bignum)
220 align!(x,other)
221 offset = other.exponent - x.exponent
222 borrow = Chunk(0)
223 @inbounds begin
224 for i = 1:other.used_digits
225 difference = x.bigits[i+offset] - other.bigits[i] - borrow
226 x.bigits[i+offset] = difference & kBigitMask
227 borrow = difference >> (kChunkSize - 1)
228 end
229 i = other.used_digits+1
230 while borrow != 0
231 difference = x.bigits[i+offset] - borrow
232 x.bigits[i+offset] = difference & kBigitMask
233 borrow = difference >> (kChunkSize - 1)
234 i += 1
235 end
236 end
237 clamp!(x)
238 end
239
240 function shiftleft!(x::Bignum,shift_amount)
241 x.used_digits == 0 && return
242 x.exponent += div(shift_amount,kBigitSize)
243 local_shift = shift_amount % kBigitSize
244 bigitshiftleft!(x,local_shift)
245 end
246
247 function multiplybyuint32!(x::Bignum,factor::UInt32)
248 factor == 1 && return
249 if factor == 0
250 zero!(x)
251 return
252 end
253 x.used_digits == 0 && return
254 carry::DoubleChunk = 0
255 @inbounds begin
256 for i = 1:x.used_digits
257 product::DoubleChunk = (factor % DoubleChunk) * x.bigits[i] + carry
258 x.bigits[i] = (product & kBigitMask) % Chunk
259 carry = product >> kBigitSize
260 end
261 while carry != 0
262 x.bigits[x.used_digits+1] = carry & kBigitMask
263 x.used_digits += 1
264 carry >>= kBigitSize
265 end
266 end
267 return
268 end
269
270 function multiplybyuint64!(x::Bignum,factor::UInt64)
271 factor == 1 && return
272 if factor == 0
273 zero!(x)
274 return
275 end
276 carry::UInt64 = 0
277 low::UInt64 = factor & 0xFFFFFFFF
278 high::UInt64 = factor >> 32
279 @inbounds begin
280 for i = 1:x.used_digits
281 product_low::UInt64 = low * x.bigits[i]
282 product_high::UInt64 = high * x.bigits[i]
283 tmp::UInt64 = (carry & kBigitMask) + product_low
284 x.bigits[i] = tmp & kBigitMask
285 carry = (carry >> kBigitSize) + (tmp >> kBigitSize) +
286 (product_high << (32 - kBigitSize))
287 end
288 while carry != 0
289 x.bigits[x.used_digits+1] = carry & kBigitMask
290 x.used_digits += 1
291 carry >>= kBigitSize
292 end
293 end
294 return
295 end
296
297 const kFive27 = UInt64(0x6765c793fa10079d)
298 const kFive1 = UInt16(5)
299 const kFive2 = UInt16(kFive1 * 5)
300 const kFive3 = UInt16(kFive2 * 5)
301 const kFive4 = UInt16(kFive3 * 5)
302 const kFive5 = UInt16(kFive4 * 5)
303 const kFive6 = UInt16(kFive5 * 5)
304 const kFive7 = UInt32(kFive6 * 5)
305 const kFive8 = UInt32(kFive7 * 5)
306 const kFive9 = UInt32(kFive8 * 5)
307 const kFive10 = UInt32(kFive9 * 5)
308 const kFive11 = UInt32(kFive10 * 5)
309 const kFive12 = UInt32(kFive11 * 5)
310 const kFive13 = UInt32(kFive12 * 5)
311 const kFive1_to_12 = UInt32[kFive1, kFive2, kFive3, kFive4, kFive5, kFive6,
312 kFive7, kFive8, kFive9, kFive10, kFive11, kFive12]
313 function multiplybypoweroften!(x::Bignum,exponent)
314 exponent == 0 && return
315 x.used_digits == 0 && return
316 remaining_exponent = exponent
317 while remaining_exponent >= 27
318 multiplybyuint64!(x,kFive27)
319 remaining_exponent -= 27
320 end
321 while remaining_exponent >= 13
322 multiplybyuint32!(x,kFive13)
323 remaining_exponent -= 13
324 end
325 remaining_exponent > 0 && multiplybyuint32!(x,
326 kFive1_to_12[remaining_exponent])
327 shiftleft!(x,exponent)
328 end
329
330 function square!(x::Bignum)
331 product_length = 2 * x.used_digits
332 (1 << (2 * (kChunkSize - kBigitSize))) <= x.used_digits && error("unimplemented")
333 accumulator::DoubleChunk = 0
334 copy_offset = x.used_digits
335 @inbounds begin
336 for i = 1:x.used_digits
337 x.bigits[copy_offset + i] = x.bigits[i]
338 end
339 for i = 1:x.used_digits
340 bigit_index1 = i-1
341 bigit_index2 = 0
342 while bigit_index1 >= 0
343 chunk1::Chunk = x.bigits[copy_offset + bigit_index1 + 1]
344 chunk2::Chunk = x.bigits[copy_offset + bigit_index2 + 1]
345 accumulator += (chunk1 % DoubleChunk) * chunk2
346 bigit_index1 -= 1
347 bigit_index2 += 1
348 end
349 x.bigits[i] = (accumulator % Chunk) & kBigitMask
350 accumulator >>= kBigitSize
351 end
352 for i = x.used_digits+1:product_length
353 bigit_index1 = x.used_digits - 1
354 bigit_index2 = i - bigit_index1 - 1
355 while bigit_index2 < x.used_digits
356 chunk1::Chunk = x.bigits[copy_offset + bigit_index1 + 1]
357 chunk2::Chunk = x.bigits[copy_offset + bigit_index2 + 1]
358 accumulator += (chunk1 % DoubleChunk) * chunk2
359 bigit_index1 -= 1
360 bigit_index2 += 1
361 end
362 x.bigits[i] = (accumulator % Chunk) & kBigitMask
363 accumulator >>= kBigitSize
364 end
365 end
366 x.used_digits = product_length
367 x.exponent *= 2
368 clamp!(x)
369 end
370
371 function assignpoweruint16!(x::Bignum,base::UInt16,power_exponent::Int)
372 if power_exponent == 0
373 assignuint16!(x,UInt16(1))
374 return
375 end
376 zero!(x)
377 shifts::Int = 0
378 while base & UInt16(1) == UInt16(0)
379 base >>= UInt16(1)
380 shifts += 1
381 end
382 bit_size::Int = 0
383 tmp_base::Int= base
384 while tmp_base != 0
385 tmp_base >>= 1
386 bit_size += 1
387 end
388 final_size = bit_size * power_exponent
389 mask::Int = 1
390 while power_exponent >= mask
391 mask <<= 1
392 end
393 mask >>= 2
394 this_value::UInt64 = base
395 delayed_multiplication = false
396 max_32bits::UInt64 = 0xFFFFFFFF
397 while mask != 0 && this_value <= max_32bits
398 this_value *= this_value
399 if (power_exponent & mask) != 0
400 base_bits_mask::UInt64 = ~(UInt64(1) << (64 - bit_size) - 1)
401 high_bits_zero = (this_value & base_bits_mask) == 0
402 if high_bits_zero
403 this_value *= base
404 else
405 delayed_multiplication = true
406 end
407 end
408 mask >>= 1
409 end
410 assignuint64!(x,this_value)
411 delayed_multiplication && multiplybyuint32!(x,UInt32(base))
412 while mask != 0
413 square!(x)
414 (power_exponent & mask) != 0 && multiplybyuint32!(x,UInt32(base))
415 mask >>= 1
416 end
417 shiftleft!(x,shifts * power_exponent)
418 end
419
420
1 (2.38%) samples spent in dividemodulointbignum!
function dividemodulointbignum!(x::Bignum,other::Bignum)
421 1 (2.38%) 1 (2.38%) bigitlength(x) < bigitlength(other) && return UInt16(0)
422 align!(x,other)
423 result::UInt16 = 0
424 @inbounds begin
425 while bigitlength(x) > bigitlength(other)
426 result += x.bigits[x.used_digits] % UInt16
427 subtracttimes!(x,other,x.bigits[x.used_digits])
428 end
429 this_bigit::Chunk = x.bigits[x.used_digits]
430 other_bigit::Chunk = other.bigits[other.used_digits]
431 if other.used_digits == 1
432 quotient = reinterpret(Int32,div(this_bigit,other_bigit))
433 x.bigits[x.used_digits] = this_bigit - other_bigit * reinterpret(UInt32,quotient)
434 result += quotient % UInt16
435 clamp!(x)
436 return result
437 end
438 end
439 division_estimate = reinterpret(Int32,div(this_bigit,other_bigit+Chunk(1)))
440 result += division_estimate % UInt16
441 subtracttimes!(x,other,division_estimate)
442 other_bigit * (division_estimate+1) > this_bigit && return result
443 while lessequal(other, x)
444 subtractbignum!(x,other)
445 result += UInt16(1)
446 end
447 return result
448 end
449
450 function pluscompare(a::Bignum,b::Bignum,c::Bignum)
451 bigitlength(a) < bigitlength(b) && return pluscompare(b,a,c)
452 bigitlength(a) + 1 < bigitlength(c) && return -1
453 bigitlength(a) > bigitlength(c) && return 1
454 a.exponent >= bigitlength(b) && bigitlength(a) < bigitlength(c) && return -1
455 borrow::Chunk = 0
456 min_exponent = min(a.exponent,b.exponent,c.exponent)
457 for i = (bigitlength(c)-1):-1:min_exponent
458 chunk_a::Chunk = bigitat(a,i)
459 chunk_b::Chunk = bigitat(b,i)
460 chunk_c::Chunk = bigitat(c,i)
461 sum::Chunk = chunk_a + chunk_b
462 if sum > chunk_c + borrow
463 return 1
464 else
465 borrow = chunk_c + borrow - sum
466 borrow > 1 && return -1
467 borrow <<= kBigitSize
468 end
469 end
470 borrow == 0 && return 0
471 return -1
472 end
473
474
1 (2.38%) samples spent in compare
function compare(a::Bignum,b::Bignum)
475 1 (2.38%) 1 (2.38%) bigit_length_a = bigitlength(a)
476 bigit_length_b = bigitlength(b)
477 bigit_length_a < bigit_length_b && return -1
478 bigit_length_a > bigit_length_b && return 1
479 for i = (bigit_length_a-1):-1:min(a.exponent,b.exponent)
480 bigit_a::Chunk = bigitat(a,i)
481 bigit_b::Chunk = bigitat(b,i)
482 bigit_a < bigit_b && return -1
483 bigit_a > bigit_b && return 1
484 end
485 return 0
486 end
487
488 function bigitat(x::Bignum,index)
489 index >= bigitlength(x) && return Chunk(0)
490 index < x.exponent && return Chunk(0)
491 @inbounds ret = x.bigits[index - x.exponent+1]::Chunk
492 return ret
493 end
494
495 end # module