Karatsuba multiplication #95

dlesnoff · 2022-01-29T20:56:54Z

This is still a draft.
I get a strange error in one test of the tliterals.nim file. I have no clue why this test is wrong when I switch the karatsuba multiplication.

dlesnoff · 2022-01-31T12:01:36Z

Lack tests for multiplication.

konsumlamm

Note: I haven't looked at karatsubaMultiplication itself yet.

src/bigints.nim

konsumlamm · 2022-02-11T13:24:48Z

src/bigints.nim

+    if bl <= karatsubaTreshold:
+      karatsubaMultiplication(a, c, b)
+    else:
+      unsignedMultiplication(a, c, b)


Suggested change

if bl <= karatsubaTreshold:

karatsubaMultiplication(a, c, b)

else:

unsignedMultiplication(a, c, b)

if bl > karatsubaTreshold:

karatsubaMultiplication(a, c, b)

else:

unsignedMultiplication(a, c, b)

konsumlamm · 2022-02-11T13:25:07Z

src/bigints.nim

+    if cl <= karatsubaTreshold:
+      karatsubaMultiplication(a, b, c)
+    else:
+      unsignedMultiplication(a, b, c)


Suggested change

if cl <= karatsubaTreshold:

karatsubaMultiplication(a, b, c)

else:

unsignedMultiplication(a, b, c)

if cl > karatsubaTreshold:

karatsubaMultiplication(a, b, c)

else:

unsignedMultiplication(a, b, c)

konsumlamm · 2022-02-11T14:14:11Z

src/bigints.nim

+  if bl == 1:
+    # base case : multiply the only limb with each limb of second term
+    scalarMultiplication(a, c, b.limbs[0])
+    return 
+  if cl == 1:
+    scalarMultiplication(a, b, c.limbs[0])
+    return


Suggested change

if bl == 1:

# base case : multiply the only limb with each limb of second term

scalarMultiplication(a, c, b.limbs[0])

return

if cl == 1:

scalarMultiplication(a, b, c.limbs[0])

return

This is already done in unsignedMultiplication afaict.

The idea was to reduce the number of operations if we know there is only one limb.
We do not have to do a whole for loop.
But you are right, we can do a scalarMultiplication with unsignedMultiplication.

src/bigints.nim

konsumlamm · 2022-02-11T18:15:07Z

src/bigints.nim

+  if bl < karatsubaTreshold:
+    if cl <= bl:
+      unsignedMultiplication(a, b, c)
+    else:
+      unsignedMultiplication(a, c, b)
+    return
+  if cl < karatsubaTreshold:
+    if bl <= cl:
+      unsignedMultiplication(a, c, b)
+    else:
+      unsignedMultiplication(a, b, c)
+    return


I think it makes more sense to move this logic into unsignedMultiplication. Then we avoid checking this twice when calling karatsubaMultiplication from multiplication and we can call unsignedMultiplication inside karatsubaMultiplication (instead of karatsubaMultiplication itself).

I have fixed many bugs in the version you just reviewed, I am sorry I should have warned you, there will be many changes to this version.
I made karatsubaMultiplication and multiplication completely independent to test them independently and compare obtained results.

It's fine, I didn't expect this to be the final version anyway. As I said, I haven't looked at karatsubaMultiplication in depth yet, since I want to better understand the algorithm first.

konsumlamm · 2022-02-11T18:18:34Z

src/bigints.nim

+  low_b.limbs = b.limbs[0 .. k-1]
+  high_b.limbs = b.limbs[k .. ^1]
+  low_c.limbs = c.limbs[0 .. k-1]
+  high_c.limbs = c.limbs[k .. ^1]


These all create new seqs, which isn't very efficient. Perhaps we should use openArray[uint32] instead (then this can use toOpenArray for O(1) slicing) or make our own "sliceable seq" to avoid copies.

Ok, I did not thougt about this. I used seq because I expect to join the results into a seq after.

@konsumlamm openArray are only used for procs arguments, when we want either a seq or an array as argument.
Can't we use some pointers here ? Do we have to create a new structure ?
Internal structure is a seq, if we use another structure, we would have to make a copy anyway or change the internal structure.
We can not use arrays neither, since we do not know the size at compile time. (It depends on the variable k, which value is known at runtime).

An openArray is basically a pointer and a length, it doesn't create a copy. Anything that doesn't create a copy but just modifies indices/pointers should be good.

I do not see how we can call multiplication then on those pointers, without making a multiplication directly on arrays.
We need to initialize BigInts

These all create new seqs, which isn't very efficient. Perhaps we should use openArray[uint32] instead (then this can use toOpenArray for O(1) slicing) or make our own "sliceable seq" to avoid copies.

The algorithm wants the value of the polynomial corresponding to each parts of the slice.
The best way so far that I conceive is to modify the BigInt's limbs field from:

type BigInt* = object limbs: seq[uint32] isNegative: bool

to

type BigInt* = object limbs: ref seq[uint32] isNegative: bool

Otherwise, we will have to reimplement addition, subtraction and base case multiplication for another container.

seqs are openarrays, if you implement for openarrays, it's implemented for seq and arrays.

We have to change the whole code.

Addition, subtraction and multiplication take bigints with a sequence parameter as input.

For the Karatsuba algorithm as well as some others, we need to manipulate these sequences through a pointer and get a pointer for parts of the sequence.

We also need to operate on those slices of the sequence, i.e. get the value associated with each part of the slice, and add, subtract, and multiply these slices.

That's why we need a seq-like container with the possibility of getting a reference to each value of the seq.

None of the openarray types enables this.

As @konsumlamm said, I think openArray and toOpenArray is good enough to implement karatsuba multiplication.
You can write recursive procedure like this using openArray and toOpenArray:

proc sum(x: openArray[int]): int = case x.len: of 0: 0 of 1: x[0] of 2: x[0] + x[1] else: let mid = x.len div 2 sum(toOpenArray(x, 0, mid - 1)) + sum(toOpenArray(x, mid, x.high)) echo sum([1, 2, 3, 4, 5])

This is a part of karatsubaMultiplication in your current PR:

var low_b, high_b, low_c, high_c: BigInt # Decompose `b` and `c` in two parts of (almost) equal length low_b.limbs = b.limbs[0 .. k-1] high_b.limbs = b.limbs[k .. ^1] low_c.limbs = c.limbs[0 .. k-1] high_c.limbs = c.limbs[k .. ^1] # subtractive version of Karatsuba's algorithm to limit carry handling var lowProduct, highProduct, add3, add4, add5, middleTerm: BigInt = zero multiplication(lowProduct, low_b, low_c) multiplication(highProduct, high_b, high_c) add3 = low_b - high_b add4 = high_c - low_c

Above code would be written using toOpenArray like:

# Decompose `b` and `c` in two parts of (almost) equal length template low_b = toOpenArray(b.limbs, 0, k - 1) template high_b = toOpenArray(b.limbs, k, b.limbs.high) template low_c = toOpenArray(c.limbs, 0, k - 1) template high_c = toOpenArray(c.limbs, k, c.limbs.high) # subtractive version of Karatsuba's algorithm to limit carry handling var lowProduct, highProduct, add3, add4, add5, middleTerm: BigInt = zero multiplication(lowProduct, low_b, low_c) multiplication(highProduct, high_b, high_c) # This code requires subtraction proc that takes 2 `openArray` add3 = low_b - high_b add4 = high_c - low_c

Thank you for the time taken and the detailed changes.

# This code requires subtraction proc that takes 2 `openArray`

I just want to point out that scalarMultiplication and unsignedMultiplication will also have to take openArrays:

if bl == 1: scalarMultiplication(a, c, b.limbs[0]) # b and c are openArrays a.isNegative = b.isNegative xor c.isNegative return ... if bl < karatsubaThreshold: if cl <= bl: unsignedMultiplication(a, b, c) # b and c are openArrays here else: unsignedMultiplication(a, c, b) # same a.isNegative = b.isNegative xor c.isNegative return ...

I can make a multiplication(a: var Bigint, b, c: openArray[uint32]) = proc, and avoid to rewrite addition and shr for openArrays.

I read the unsignedMultiplication and scalarMultiplication proc algorithms again and they effectively don't need parameters to be BigInts.
The substract proc will be quite delicate to convert for OpenArray parameters though.

I will look into it.

konsumlamm · 2022-02-11T18:28:03Z

src/bigints.nim

+  # limit carry handling in opposition to the additive version
+  var
+    lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero
+  karatsubaMultiplication(lowProduct, low_b, low_c)


low_b and low_c aren't necessarily normalized afaict, so that may cause problems.

This might be the problem I encounter at compilation time. I have not done any test when the operands does not have the same limbs sequence’s length.

src/bigints.nim

Co-authored-by: konsumlamm <[email protected]>

Add recent changes in the library.

Karatsuba only works in running time. Strange errors when the tests are executed as static. I commented for the present. Tests should be moved from main to tbigints.nim

demotomohiro · 2022-09-16T23:57:27Z

src/bigints.nim

+  # Decompose `b` and `c` in two parts of (almost) equal length
+  low_b.limbs = b.limbs[0 .. k-1]
+  high_b.limbs = b.limbs[k .. ^1]
+  low_c.limbs = c.limbs[0 .. k-1]


When c.limbs.len is smaller than half of b.limbs.len, k is larger than c.limbs.len and low_c.limbs = c.limbs[0 .. k-1] cause IndexDefect.
It would be better to add tests that multiply a large value by a small value.

I can simply change n to the min of bl and cl in this case.
I want to add these tests, but I am waiting for my other PRs to be merged. In these, I have added random generation and benchmarks. I am especially waiting for #112 . With initRandomBigInt, I will be able to generate bigints of a specific size for tests.

dlesnoff · 2022-10-11T08:15:41Z

Thanks. I am not sure I can continue working on this PR, I am terribly busy. I will try to move away the tests under the when isMainModule instruction block. Maybe I will test your suggestion, but I do not think I will have time to do much debugging. Furthermore, I am really concern by the choice of uint32 vs uint64 (as suggested by @mratsim). This seems to be a big change for the library, that we need to tackle beforehand.

konsumlamm · 2022-10-16T11:52:08Z

Furthermore, I am really concern by the choice of uint32 vs uint64 (as suggested by @mratsim). This seems to be a big change for the library, that we need to tackle beforehand.

Changing that requires proper support for addition and multiplication of uint64 without overflow (and I don't think writing that in assembly is a good solution). I don't see what that has to do with Karatsuba multiplication though, it doesn't really change the algorithm.

dlesnoff · 2022-11-07T08:28:51Z

src/bigints.nim

@@ -420,6 +419,30 @@ func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} =
      inc pos
  normalize(a)

+func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} =
+  # always called with bl >= cl


This comment needs to be removed, since cl == 1 in this case.

dlesnoff added 5 commits January 14, 2022 21:30

First implementation of untested Karatsuba

bdabfe8

Call Karatsuba with a treshold - does not compile

244a6a8

Made karatsuba treshold a global const

f2a48e2

Merge branch 'master' into karatsuba

9e4efd6

Fixed some expressions

f3bcc9b

dlesnoff marked this pull request as draft January 29, 2022 20:58

dlesnoff force-pushed the karatsuba branch from 07433ff to 95067aa Compare January 31, 2022 10:32

Karatsuba multiplication now works

bc118cb

dlesnoff force-pushed the karatsuba branch from 95067aa to bc118cb Compare January 31, 2022 10:35

dlesnoff marked this pull request as ready for review January 31, 2022 10:38

add tests

438a26e

dlesnoff marked this pull request as draft January 31, 2022 12:01

Add randomized tests and fixed call to karatsuba

5dd2512

konsumlamm reviewed Feb 11, 2022

View reviewed changes

dlesnoff and others added 8 commits February 11, 2022 21:09

Remove overflow error in scalar multiplication

dc69b3c

Co-authored-by: konsumlamm <[email protected]>

treshold -> threshold

fd4d321

Co-authored-by: konsumlamm <[email protected]>

Many changes I forgot to commit

bb05ee5

Merge review commits

1327726

Merge branch 'master' into karatsuba

d7bfc9f

Add recent changes in the library.

Add tests and last recommandations of the review

290a1d1

Karatsuba only works in running time. Strange errors when the tests are executed as static. I commented for the present. Tests should be moved from main to tbigints.nim

Remove echo's, convert proc into func again

334537a

Merge branch 'nim-lang:master' into karatsuba

463ece4

demotomohiro suggested changes Sep 17, 2022

View reviewed changes

dlesnoff commented Nov 7, 2022

View reviewed changes

dlesnoff mentioned this pull request Feb 7, 2024

for loops for bigint is very slow #143

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Karatsuba multiplication #95

Karatsuba multiplication #95

dlesnoff commented Jan 29, 2022

dlesnoff commented Jan 31, 2022

konsumlamm left a comment

konsumlamm Feb 11, 2022

konsumlamm Feb 11, 2022

konsumlamm Feb 11, 2022

dlesnoff Feb 12, 2022

konsumlamm Feb 11, 2022

dlesnoff Feb 11, 2022 •

edited

Loading

konsumlamm Feb 12, 2022

konsumlamm Feb 11, 2022

dlesnoff Feb 11, 2022

dlesnoff Mar 27, 2022

konsumlamm Mar 30, 2022

dlesnoff May 3, 2022

dlesnoff Nov 3, 2022

mratsim Nov 3, 2022

dlesnoff Nov 3, 2022

demotomohiro Nov 7, 2022

dlesnoff Nov 7, 2022

konsumlamm Feb 11, 2022

dlesnoff Feb 12, 2022

demotomohiro Sep 16, 2022

dlesnoff Sep 19, 2022

dlesnoff commented Oct 11, 2022 via email •

edited

Loading

konsumlamm commented Oct 16, 2022

dlesnoff Nov 7, 2022

Karatsuba multiplication #95

Are you sure you want to change the base?

Karatsuba multiplication #95

Conversation

dlesnoff commented Jan 29, 2022

dlesnoff commented Jan 31, 2022

konsumlamm left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlesnoff Feb 11, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlesnoff commented Oct 11, 2022 via email • edited Loading

konsumlamm commented Oct 16, 2022

Choose a reason for hiding this comment

dlesnoff Feb 11, 2022 •

edited

Loading

dlesnoff commented Oct 11, 2022 via email •

edited

Loading