(simple?) problem with multiplication

[Paul]

Hello,

I need a fast way to calculate (a * b) % B, with the following
constraints:

- a, b, B are int
- B = 10^N , N>0 arbitrary (B always fits in an int)
- the result must be valid even if the multiplication overflows
- portable code, no assumptions on the sizeof(int) and no types larger
than int can be used.

Any help will be much appreciated.

If you convert the integer into an array such that 543 was {5,4,3}, or in
a vector you can do it.
Here is an example of how to overcome such a problem using vectors, it may
not be perfect but it should give you the idea:
#include <iostream>
#include <vector>

typedef unsigned int uint;
typedef std::vector<uint> v_uint;

v_uint foo(v_uint v1, v_uint v2, uint N){
uint pow = N;
v_uint retv;
uint carry=0;
uint temp_prod =0;

//v1 must contain the largest.
for(int i=0; i<v1.size() && pow>=10; i++, pow/=10){
for (int j=i; j>=0; j-- ){
if(i<v2.size()){
temp_prod += v1[i-j]* v2[j];
}
if( i>v2.size() ){
if( (i-j)<v2.size() )
temp_prod+= v1[j]* v2[i-j];
}
}
temp_prod +=carry;
retv.push_back( temp_prod%10 );
carry = temp_prod/10;
temp_prod=0;
}
return retv;
}

int main(){
uint N = 100000;
v_uint a(5, 3);
v_uint b(4,5);
v_uint v = foo(a,b, N);

for(int i=0; i<v.size() ; i++){
std::cout<< v << std::endl;
}

std::cout<< (33333*5555)%100000;

}

Previous code was a bit buggy I will include a revised version which
includes a couple of helper functions.
Obviously there is room for improvement but it seems to work ok.

#include <iostream>
#include <vector>

typedef unsigned int uint;
typedef std::vector<uint> v_uint;

v_uint int_to_vect(uint arg, uint base){
v_uint retv;
while(arg){
retv.push_back(arg%base);
arg/=base;
}
return retv;
}

uint vect_to_int(v_uint v, uint B){
uint temp=0;
uint pow=1;

for (int i=0; i<v.size(); i++, pow*=B){
temp+=v*pow;
}

return temp;
}

uint power_mag(uint B, uint pow){
uint temp=(pow)?B:1;
for (int i=1; i<pow; i++){
if(temp*B>temp)
temp*=B;
}
return temp;
}


uint foo(uint a, uint b, uint B, uint p){
v_uint v1= int_to_vect((a>=b)?a:b, B);
v_uint v2= int_to_vect((a>=b)?b:a, B);
int v1s= v1.size();
int v2s= v2.size();
uint pow = power_mag(B,p);
v_uint retv;
uint carry=0;
uint temp_prod =0;

for(int i=0; i<(v1s+v2s-1) && pow>=B; i++, pow/=B){
for (int j=i; j>=0; j-- ){
if(i<v2s){
temp_prod += v1[i-j]* v2[j];
}
if( i>=v2s && i<v1s && (i-j)<v2s ){
temp_prod += v1[j]* v2[i-j];
}
if(i>=v1s && j<(v2s-1)-(i-v1s) ){
temp_prod+= v1[v1s-1-j] * v2[j+(i-(v1s-1))];
}
}
temp_prod +=carry;
retv.push_back( temp_prod%B );
carry = temp_prod/B;
temp_prod=0;
}
while(pow>=B && carry){
retv.push_back(carry%B);
carry/=B;
pow/=B;
}
return vect_to_int(retv,B);
}


int main(){
uint B = 8; //Base for divisor
uint pow = 6; //Power for divisor

uint x = foo(96,9677,B, pow);

//Example (9677*96)%8^6
}



Please let me know if you find a better way of doing it, It is alot simpler
if B==2 but unfortunately your requirement was B==10.

 

[gwowen]

Hello,

I need a fast way to calculate (a * b) % B, with the following
constraints:

- a, b, B are int
- B = 10^N , N>0 arbitrary (B always fits in an int)
- the result must be valid even if the multiplication overflows
- portable code, no assumptions on the sizeof(int) and no types larger
than int can be used.

Any help will be much appreciated.

Here's a sketch solution. It's not well tested, but I think the ideas
are sound... It uses unsigneds rather than ints, and takes N as a
template parameter, both for convenience, but nothing bigger.

[implementation - modint.hpp]
#ifndef MODINT_HPP
#define MODINT_HPP

#include<climits>

template <unsigned N>
class modint
{
private:
unsigned val_;
static modint overflowmod_;

public:
void operator+=(const modint& rhs) {
unsigned increment = (rhs.val()%N);
if(val_ > UINT_MAX - increment){ // overflow
increment -= N-val_; // add enough to make val_ == N
~= 0
val_ = 0;
}
val_ = (val_ + increment) % N;
}

void rightshift() { val_ /= 2;}
void leftshift() {
const bool overflow = (val_ > UINT_MAX/2);
val_ *= 2;
val_ %= N;
if(overflow) *this += modint(overflowmod_);
}

unsigned val() const { return val_;}
modint(unsigned x) : val_(x % N) {}
};

template <unsigned N>
modint<N> modint<N>:verflowmod_ = (1 + (UINT_MAX % N));

// peasants arithmetic, mod N
template <unsigned N>
modint<N> operator*(modint<N> r1,modint<N> r2)
{
modint<N> tot = 0;
while(r2.val() != 0){
if(r2.val() % 2 == 1) {
tot += r1;
}
r2.rightshift();
r1.leftshift();
}
return tot;
}
#endif

[test program]

#include <iostream>
#include "modint.hpp"

int main(void)
{
const unsigned BIGVAL = 0xFFFFFFFF;
modint<BIGVAL> a = BIGVAL-3, b = BIGVAL-9;

std::cout << a.val() << " " << b.val() << std::endl;
std::cout << (a * b).val() << std::endl;
}

 

[Paul]

Hello,

I need a fast way to calculate (a * b) % B, with the following
constraints:

- a, b, B are int
- B = 10^N , N>0 arbitrary (B always fits in an int)
- the result must be valid even if the multiplication overflows
- portable code, no assumptions on the sizeof(int) and no types larger
than int can be used.

Any help will be much appreciated.

Here's a sketch solution. It's not well tested, but I think the ideas
are sound... It uses unsigneds rather than ints, and takes N as a
template parameter, both for convenience, but nothing bigger.

[implementation - modint.hpp]
#ifndef MODINT_HPP
#define MODINT_HPP

#include<climits>

template <unsigned N>
class modint
{
private:
unsigned val_;
static modint overflowmod_;

public:
void operator+=(const modint& rhs) {
unsigned increment = (rhs.val()%N);
if(val_ > UINT_MAX - increment){ // overflow
increment -= N-val_; // add enough to make val_ == N
~= 0
val_ = 0;
}
val_ = (val_ + increment) % N;
}

void rightshift() { val_ /= 2;}
void leftshift() {
const bool overflow = (val_ > UINT_MAX/2);
val_ *= 2;
val_ %= N;
if(overflow) *this += modint(overflowmod_);
}

unsigned val() const { return val_;}
modint(unsigned x) : val_(x % N) {}
};

template <unsigned N>
modint<N> modint<N>:verflowmod_ = (1 + (UINT_MAX % N));

// peasants arithmetic, mod N
template <unsigned N>
modint<N> operator*(modint<N> r1,modint<N> r2)
{
modint<N> tot = 0;
while(r2.val() != 0){
if(r2.val() % 2 == 1) {
tot += r1;
}
r2.rightshift();
r1.leftshift();
}
return tot;
}
#endif

[test program]

#include <iostream>
#include "modint.hpp"

int main(void)
{
const unsigned BIGVAL = 0xFFFFFFFF;
modint<BIGVAL> a = BIGVAL-3, b = BIGVAL-9;

std::cout << a.val() << " " << b.val() << std::endl;
std::cout << (a * b).val() << std::endl;
}

-----------------------------------------------------------

How is this used?

With the expression (a*b)% B^N

Is BIGVAL N or 10^N?

 

[SG]

In case this hasn't been mentioned already (I just skimmed this
thread):

(a*b)%B = ((a%B)*(b%B))%B

The last expression can be computed with a modified double-and-add
algorithm for multiplication that includes the modulo reduction.

Cheers!
SG

 

[Kai-Uwe Bux]

Yes, but it doesn't solve the problem. :-( Suppose a*b overflows and
that B is greater than both a and b. Then a%B is a, and b%B is b; then
(a%B)*(b%B) is a*b, and it overflows.

Another approach is brute force: write double-precision multiplication
and division routines; not too hard, but involves tedious details;

That seems to be vicious snipping on your part. SG said:

(a*b)%B = ((a%B)*(b%B))%B

The last expression can be computed with a modified double-and-add
algorithm for multiplication that includes the modulo reduction.

So, he did not rely on (a%B)*(b%B) not overflowing.

That said, a version of double-and-add for modulo arithmetic has been
proposed elsethread; and that version at least did assume that 2*(a%B) does
not overflow.


Best,

Kai-Uwe Bux

 

[Kai-Uwe Bux]

"Vicious"? Really?
Yes.


Um, if the expression ((a%B)*(b%B)) overflows the behavior is
undefined. "The last expression" is that result %B. If that's not your
interpretation, please allow for the possibility that others may have
read it differently without being "vicious".

A fair interpretation has to take into account the full material of the
text. How does your interpretation deal with the "can be compute with a
modified double-and-add algorithm for multiplication" part? That part makes
it perfectly clear that SG did not intend (a%B)*(b%B) to be evaluated as a C
expression. Your interpretation (forced upon readers of your posting by
leaving out the part it ignored) makes it sound as though he did. That is
why I think your quoting was vicious: it left out a part that is crucial to
the interpretation of the part that was quoted.


Best,

Kai-Uwe Bux

 

[Paul]

Okay. The signal to bullshit ratio in this newsgroup is fast approaching
zero. It's not worth the time spent reading it any more. Have a good life.

Firstly I'm sorry to see that Pete has spat the dummy out.

Secondly I have found another solution to the given problem, for unsigned
integers only.
Multiplication of two matrices where:
m1 == size of smallest int (number of decimal digits).
m2 == size of largest int + (size of smallest int -1) x size of smallest
int.

For example with 43296x21 we would multiply:
6 0
9 6
2 9 [6x2 array.]
3 2
4 3
0 4
multiplied by....
1
2 [2x1 array]

The resulting array would be a 6x1 array containing:
[6 21 20 7 10 8]

Which we would then simply loop through taking modulo and carrying result.
I don't feel too well tonight so early to bed for me, I may do an algorithm
for this tomorrow if I feel better.

Note: m1*m2 != m2*m1.
HTH.

 

[Alf P. Steinbach /Usenet]

* Kai-Uwe Bux, on 24.06.2011 20:28 -> Pete Becker:

... I think your quoting was vicious: it left out a part that is crucial to
the interpretation of the part that was quoted.

I think the word you're looking may be "inadvertently", not "vicious".


Cheers & hth.,

- Alf

 

[Kai-Uwe Bux]

My double-and-add version did have an overflow check. It's a simple
problem - its only the details are a pain.

Oops: I didn't see yours, whence I was not specific enough. The above should
have specifically referred to this proposal by Anonymous:

replace c with B in the below code, and long long with int, then test it

/* this function calculates (a*b)%c taking into account that a*b might
overflow */
long long mulmod(long long a,long long b,long long c){
long long x = 0,y=a%c;
while(b > 0){
if(b%2 == 1){
x = (x+y)%c;
}
y = (y*2)%c;
b /= 2;
}
return x%c;
}


Best,

Kai-Uwe Bux

 

[Paul]

My double-and-add version did have an overflow check. It's a simple
problem - its only the details are a pain.

Its so simple that you cannot provide a solution? lol very good.

 

[Paul]

A fair interpretation has to take into account the full material of the
text. How does your interpretation deal with the "can be compute with a
modified double-and-add algorithm for multiplication" part?

My interpretation of that is that it's unproven bullshit.

Lets see this magical "modified double and add algorithm" that, for you,
solves the problem

 

[SG]

This was (the uninteresting) part of the solution. You snipped the
rest. Obviously the whole line is not even well-formed C nor C++ code
since the left side is an rvalue. I simply used C's %-syntax for a
_mathematical_equation_. You seem to have ignored my reference to
"modified double-and-add algorithm" (the interesting part of the
solution).

"Vicious"? Really?

Yes, it was. You've overlooked something.

Um, if the expression ((a%B)*(b%B)) overflows the behavior is

If you copy&paste it into a source code file, of course. But at no
place I suggested to compute

((a%B)*(b%B))

I suggested to compute

((a%B)*(b%B))%B <-- as a mathematical expression, not a C
expression

VIA a _modified_double-and-add_algorithm_. And by that I don't mean
copy&paste this thing into source code. I expected people to use their
brain and read the whole post. It's not my job to implement this
algorithm and show it here for what looks like a homework question.

As for the reduced signal-to-bullshit ratio, I'm afraid you're not
completely innocent in this case.

Cheers!
SG

 

[SG]

[...]
My interpretation of that is that it's unproven bullshit.

Lets see this magical "modified double and add algorithm" that, for you,
solves the problem

Hmmm. It seems like I expected too much from you guys in terms of
reading comprehension, google skills, and education.

"double and add" for multiplication is the little brother of
"square and multiply" for computing powers.

Here it is:
http://en.wikipedia.org/wiki/Multiplication_algorithm#Shift_and_add

Also, I explained "modified" with a "that includes the mod-reduction"
in my initial post. In case you have no idea what the purpose of this
modification is: It's to avoid _overflows_.

*sigh*

Cheers!
SG

 

[Paul]

[...]
My interpretation of that is that it's unproven bullshit.

Lets see this magical "modified double and add algorithm" that, for you,
solves the problem

Hmmm. It seems like I expected too much from you guys in terms of
reading comprehension, google skills, and education.

"double and add" for multiplication is the little brother of
"square and multiply" for computing powers.

Here it is:
http://en.wikipedia.org/wiki/Multiplication_algorithm#Shift_and_add

Also, I explained "modified" with a "that includes the mod-reduction"
in my initial post. In case you have no idea what the purpose of this
modification is: It's to avoid _overflows_.

*sigh*

I know what shift and add means, but I don't know what your modified shift
and add means.
I want to see your modified shift and add algorithm that solves this
problem..

 

[Paul]

[...]
My interpretation of that is that it's unproven bullshit.

Lets see this magical "modified double and add algorithm" that, for you,
solves the problem

Hmmm. It seems like I expected too much from you guys in terms of
reading comprehension, google skills, and education.

"double and add" for multiplication is the little brother of
"square and multiply" for computing powers.

Here it is:
http://en.wikipedia.org/wiki/Multiplication_algorithm#Shift_and_add

Also, I explained "modified" with a "that includes the mod-reduction"
in my initial post. In case you have no idea what the purpose of this
modification is: It's to avoid _overflows_.

*sigh*

I know what shift and add means, but I don't know what your modified shift
and add means.
I want to see your modified shift and add algorithm that solves this
problem..

 

[Paul]

Okay. The signal to bullshit ratio in this newsgroup is fast approaching
zero. It's not worth the time spent reading it any more. Have a good
life.

Firstly I'm sorry to see that Pete has spat the dummy out.

Secondly I have found another solution to the given problem, for unsigned
integers only.
Multiplication of two matrices where:
m1 == size of smallest int (number of decimal digits).
m2 == size of largest int + (size of smallest int -1) x size of smallest
int.

For example with 43296x21 we would multiply:
6 0
9 6
2 9 [6x2 array.]
3 2
4 3
0 4
multiplied by....
1
2 [2x1 array]

The resulting array would be a 6x1 array containing:
[6 21 20 7 10 8]

Which we would then simply loop through taking modulo and carrying result.
I don't feel too well tonight so early to bed for me, I may do an
algorithm for this tomorrow if I feel better.

Note: m1*m2 != m2*m1.
HTH.

Ok this seems to do the trick:

#include <iostream>
#include <vector>

typedef unsigned int uint;
typedef std::vector<uint> v_uint;

v_uint int_to_vect(uint arg, uint B=10){
v_uint retv;
while(arg){
retv.push_back(arg%B);
arg/=B;
}
return retv;
}


uint mod_ab(uint a, uint b, uint exp=1, uint B=10){
v_uint v1= int_to_vect((a>=b)?a:b, B);
v_uint v2= int_to_vect((a>=b)?b:a, B);
int v1s= v1.size();
int v2s= (v2.size())-1;
uint retv=0;
uint mult=1;/*Annoying variable*/

for(int i=v2s; i>0; --i){
v1.push_back(0);
v1.insert(v1.begin(), 0);
}
v_uint vresult;
uint temp=0;
for(int i=0; i<(v1s+v2s) && i<exp; ++i, mult*=B){
for(int j=v2s; j>=0; --j){
temp+=v1[i+v2s-j]*v2[j];
}
retv += (temp%B)*mult;
temp = temp/B;
}
return retv;
}

int main(){
uint a=-1, b=-1;
std::cout<< mod_ab(a,b,3);
}

I haven't fully checked it as I just did it.

One thing that was annoying me was the need for the varaible mult.
All I need it for is to calculate B^i, which I was having problems
simplifying when i is zero.

I think it works for all bases up to max_value of uint but haven't checked
all that , it would probably be a good idea to make B unsigned and limit
that to max_base 255.

 

[Paul]

*sigh*

Put Paul in a killfile, and sigh no more.[0]

[0] 87.9% of respondents said their quality of life improved.

gwowen has tried and failed to solve this problem, Its his loss if he
killfiled me because obviously I am the superior programmer because I solved
the problem and he didn't <shrug>.

You seem to think you can produce some modified double and add algorithm,
which I think is doubtfull. If you can lets see it, if you can't just
*sigh*. Killfiling me won't help you either way