File Exchange

image thumbnail

Variable Precision Integer Arithmetic

version 1.45 (2.86 MB) by

Arithmetic with integers of fully arbitrary size. Arrays and vectors of vpi numbers are supported.

4.84091
56 Ratings

80 Downloads

Updated

View License

Editor's Note: This file was selected as MATLAB Central Pick of the Week

Every once in a while, I've wanted to do arithmetic with large integers with magnitude exceeding that which can fit into MATLAB's standard data types. Since I don't have the symbolic toolbox, the simple solution was to write it in MATLAB. I did that in these tools which are entirely written in MATLAB, so there is no need for compiled code.
Arithmetic is simple with the vpi tools.
A = vpi(17)^17
ans =
827240261886336764177

17 + A^17
ans =
39786732894291535047752038041559739510060813980024082
30012867731573722066105737100731556603857745946047229
53759676529121155309750944582301597489457676380805029
59227566911971103003303064782118652210655457390045806
99039190393572334521701109889855832341416056005878848
49943142324389193616484809157960034059531548585473213
36465170635561696613297503569949729314
 
There are also a few nice add ons, for example a tool to compute exact binomial coefficients for large arguments, or large factorials, or convert binary numbers with thousands of digits to decimal (vpi) form.

For example, the existing nchoosek function in matlab gets upset for even reasonably small binomial coefficients.

nchoosek(100,50)
Warning: Result may not be exact. Coefficient is greater than 1.000000e+15
 and is only accurate to 15 digits.
> In nchoosek at 66
ans =
   1.0089e+29

However, nchoosek has no such issues on vpi numbers.

nchoosek(vpi(100),50)
ans =
100891344545564193334812497256

Similarly, the computation of factorial(171) will cause an overflow. While I'll admit that there are many good ways to avoid this problem, the factvpi function has no problems at all.

factorial(171)
ans =
   Inf

factorial(vpi(171))
ans =
12410180702176678234248405241031039926166055775016931
85388951803611996075221691752992751978120487585576464
95950167038705280988985869071076733124203221848436431
04735778899685482782907545415619648521534683180442932
39598173696899657235903947616152278558180061176365108
428800000000000000000000000000000000000000000

There are now GCD and LCM tools, both of which can accept more than two input arguments.

lcm(vpi(123452356),12344332,65364467)
ans =
3557547184310976844988

I've also put in some tools that can test for primality. For example, the Mersenne prime:

p = vpi(2)^127 - 1
ans =
170141183460469231731687303715884105727

isprime(p)
ans =
     1

Factoring of vpi numbers is now implemented.

factor(vpi('1234567890123456789'))
ans =
       3 3 101 3541 3607 3803 27961

Vectors or arrays of vpi numbers work very nicely now.

A = vpi(eye(3))*3 + 1
A =
   4 1 1
   1 4 1
   1 1 4

A^17
ans =
   5642305908354 5642176768191 5642176768191
   5642176768191 5642305908354 5642176768191
   5642176768191 5642176768191 5642305908354

Dozens of other tools are also included. I've even included a tool just for fun that will convert a number into a readable text version of it as a large integer.

vpi2english(vpi('12000000110022987'))
ans =
twelve quadrillion, one hundred ten million, twenty two thousand, nine hundred eighty seven

For those Project Euler solvers around, vpi makes many of the problems easy to solve.

Addenda - Ben Petschel has graciously given me code for unique and sortrows operations on vpi arrays. Many thanks to Ben!

Comments and Ratings (163)

Mohsin Shah

I am writing a code for Paillier Cryptosystem for image encryption. The encryption of Paillier Crypto goes like this: C = (vpi(g)^pixelvalues)*vpi(r)^N mod N^2, where N is large number. The resultant C image is displayed in the workspace of my MATLAB as 256x256 vpi (256x256 being image dimensions). When I use imshow to see the encrypted image C, MATLAB gives this error:
Error using imageDisplayValidateParams
Expected input number 1, I, to be one of these types:
numeric, logical
Instead its type was vpi.

How to display image that is processed with vpi?

Mohsin Shah

Paul

Paul (view profile)

Jan Sieber

I have been using the tool for teaching purposes for a while. Thanks a lot for the excellent tool and a demo for OO programming in matlab. A small problem, though (not sure where else to contact):
>> nextprime(vpi(2)^78)
Error using nextprime (line 90)
The maximum value of N (for numeric input) allowed is 2^46.

isnumeric returns true for vpi's but nextprime gives an error for large numeric inputs. Am I making a silly mistake? It's easy to "fix", but I would like students to download vpi directly from the file exchange.

David Blake

Aaaah sorry, I wasn't constructing vpi variables! Silly noob mistake.

This works:

INT = vpi('920314723904709')
INT2 = vpi(123498374987439873'')
INT3 = INT + INT2

I'd remove my earlier comment if I could (how embarassing!) but I can't see how to do that!

David Blake

Hi, I'm a new MatLab user. I can't see how to use this code. All I want to do is add 2 large integers. The integers are in string format.

So I have:

INT = '920314723904709'

INT2 = '123498374987439873'

INT3 = INT + INT2

And I get an error message saying:

Error using +
Matrix dimensions must agree.

What am I missing please?

Arafat Asghar

Great work John. Keep it up. Is there any way I can convert a vpi number to a character string?

aiman saad

you are amazing , thank you very much, how did you do that????

John D'Errico

John D'Errico (view profile)

I could not have said it better than did Erich. This is perhaps the most common mistake I see people make, that of not recognizing when an intermediate computation is done as a double, versus VPI (or HPF, for my other high precision tool.)

Erich

Erich (view profile)

John,

The error message should give you a clue: Your 17*17*...*17 is all performed in double precision, so that the resulting "integer" cannot be accurately represented in double precision. Try this instead:

B=vpi(17)^17
>> B/(vpi(17)*17*17*17*17*17*17*17*17*17*17*17*17*17*17)
ans =
    289

Casting the first 17 to a vpi causes all further multiplication to be converted to vpi, so the resulting integer is exact.

Here is another hint:
>> vpi(17*17*17*17*17*17*17*17*17*17*17*17*17*17*17)
Error using vpi (line 160)
If N is a double, it may be no larger than 2^53 - 1

John BG

John BG (view profile)

hi Mr D'Ericco

I get an error when attempting to reach 17 from 17^17:

>> B=vpi(17)^17
B =
    827240261886336764177
 
>> B/(17*17*17)
ans =
    168377826559400929
>> B/(17*17*17*17*17*17)
ans =
    34271896307633
>> B/(17*17*17*17*17*17*17*17*17)
ans =
    6975757441
>> B/(17*17*17*17*17*17*17*17*17*17*17*17)
ans =
    1419857
>> B/(17*17*17*17*17*17*17*17*17*17*17*17*17*17*17)
Error using vpi (line 160)
If N is a double, it may be no larger than 2^53 - 1
Error in vpi/quotient (line 135)
  denominator = vpi(denominator);
Error in / (line 41)
  Q = quotient(numerator,denominator);

shouldn't it stand straight without returning error until 17 times 17?

thanks for time and attention
awaiting answer

John
jgb2012@sky.com

See Long Wong

John D'Errico

John D'Errico (view profile)

I've never used Simulink, nor do do I have Simulink to test it. I can only presume that if other MATLAB functionality is available in Simulink, then VPI would work similarly.

Can I use it in simulink?

Afefa Asiri

Yu Ang Tan

Well done! Thanks for sharing this awesome package!

Guangxi Liu

Ander Biguri

Ander Biguri (view profile)

Ren

Ren (view profile)

good work!

John D'Errico

John D'Errico (view profile)

If you have a case where the algorithm I used fails, feel free to send it to me.

p = vpi('124233534465235455634245356345246456233758457665875654647634645745657567464774575875658845645746364745747465766463352435473');

Now, square it.

p2 = p*p
p2 =
    15433971085724845807529053262922462001895202948051311590074369833049
509154900758731761178232321235094091820663015565827704498284337773042059
357095682140561450981630534851931113170834845007762750960372854097063386
898516121202278514094760628733729

And then form the square root.

q= sqrt(p2)
q =
    12423353446523545563424535634524645623375845766587565464763464574565
7567464774575875658845645746364745747465766463352435473

p - q
ans =
    0

As you can see, it found the exact square root, using a reasonable and efficient method.

John Klempir

Is there a way to allow it to calculate the exact square root? I have a perfect square, but the integer is big that the code doesn't even attempt an exact square root?

Hi,
Is it possible to use logical functions (XOR, AND,OR) on vpi numbers ? I need to do such operations in cryptography where operands are 192-bit size. thanks.

John D'Errico

John D'Errico (view profile)

Of course it is free. The file exchange provides tools that are free for your use. VPI does not come with MATLAB directly because it was written by a third party. In this case, I am that third party.

Just download it, install it as directed, then use it and enjoy it.

Is this free to download? R2014a which I am using does not have vpi apparently...

John D'Errico

John D'Errico (view profile)

Hi Nicolas,

VPI uses a Pollard rho algorithm for factoring. While it is significantly better than MATLAB's builtin factor tool, which can only handle numbers as large as about 10 digits or so, the scheme I used is often ineffectual on numbers that are the product of two seriously large primes. Of course, that is a difficult problem in general, else many of the encryption schemes in common use today would fail miserably.

One of my goals has been to re-write factor to use a better scheme, but until then I should put in a warning message when factor has failed to resolve a composite number fully into its prime factors.

So you did nothing wrong here, except to give factor a number too large for it to handle.

Hi.

First, great tools!

Second, what am I doing wrong?

>> p1 = vpi('40094690950920881030683735292761468389214899724061')
p1 =
40094690950920881030683735292761468389214899724061
>> p2 = vpi('37975227936943673922808872755445627854565536638199')
p2 = 37975227936943673922808872755445627854565536638199

>> factor(p1*p2)
ans =
    1522605027922533360535618378132637429718068114961380688657908494580122963258952897654000350692006139

It should give p1 and p2...

Thanks!

Thanks sir.
if someone is interested, I answer my question about using vpi in simulink (see below comment). I converted the vpi number into a matrix of bits (actually a bit is stored in double).

here are the steps :
seed = vpi('123456789')
N0 = randint(vpi(seed),[1,1])
N1 = vpi2bin(N0) % Convert vpi from decimal to binary form
N2 = [N1 - '0'] % Convert number to string (vector)
N3 = vec2mat(N2,8) % Convert vector into matrix of 8 columns.

John D'Errico

John D'Errico (view profile)

N = vpi('123456789012345678901234567890');

vpi2bin(N)
ans =
1100011101110100100001111111101101100001101110011111000001110111001001110001111110000101011010010

I suppose in hindsight, I should have used dec2bin as the name for this.

Hi,
How can I convert a vpi number to binary ?
Thanks.

Giuseppe Cardillo

Dear John,
I saw that you use the Miller-Rabin algorithm to test primality. considering that are known the strong pseudoprimes to the first 9 prime bases (http://arxiv.org/pdf/1207.0063.pdf), I suggest a little implementation to speed up your Miller-Rabin algorithm:

if N<2047
    witnesses=2;
elseif N<1373653
    witnesses=[2 3];
elseif N<25326001
    witnesses=[2 3 5];
elseif N<3215031751
    witnesses=[2 3 5 7];
elseif N<2152302898747
    witnesses=[2 3 5 7 11];
elseif N<3474749660383
    witnesses=[2 3 5 7 11 13];
elseif N<341550071728321
    witnesses=[2 3 5 7 11 13 17 19];
elseif N<3825123056546413051
    witnesses=[2 3 5 7 11 13 17 19 23];
else
    witnesses=[2 3 5 7 11 13 17 19 23 29 31 37];
end

Regards

John D'Errico

John D'Errico (view profile)

Hi Suma,

I recall seeing exactly this same question a few days ago on some forum. Maybe it was on Answers, but I'm not going to search for it. The responses there told you that the issue was floating point arithmetic.

Regardless, you are trying to solve a FLOATING point issue with an INTEGER arithmetic tool. Perhaps you need to consider that there is a problem with that whole idea. I'm not at all sure what answer you expect to see from the computation you did.

Worse, the problem is not even with eig, but still with your understanding of floating point arithmetic. (Please don't try throwing this same problem at my HPF tool, as it won't solve the problem either.)

SO. The fact is, it looks like you have a singular matrix. Floating point trash can creep in on problems like this to yield slightly negative numbers, on the order of eps. Throwing random pieces of code at the results of eig will not fix the problem by trying to increase the precision of MATLAB.

HOWEVER, if you know the eigenvalues can never be negative, then why not just use max? (By the way, please don't ask again about this question here, as I do NOT do consulting in the comments of my code, especially when the problem has no relevance to my code.)

Suma

Suma (view profile)

hi,

I have a 2x2 matrix A, its eigen value returned me negative zeros,eg-
 -0.000000000000000 and -0.000000000000000
  1.000000000000000 0.999999999999999

So I tried to use the vpi function like-
lamda(:,k) = vpi(eig(A(:,:,k)))

but this returned me error message-

Warning: Converted a real float to a vpi integer form
> In vpi.vpi at 172
  In vpi.vpi at 138

obviously my vpi implementation is wrong. can vpi be used in such cases?

Jean-Luc

Edgar,

I suspect the Simulink toolbox is calling isnumeric on a VPI object. Try making a file isnumeric.m in the @vpi folder with something like:

function val = isnumeric(INT)
val = true;

Jean-Luc

Tétef

Tétef (view profile)

John D'Errico

John D'Errico (view profile)

Edgar - I'm sorry, but I know absolutely nothing about Simulink, so I can offer no help here. Perhaps someone else can do so, or you can contact MathWorks and ask them.

Hi,
thanks Mr D'Errico for sharing,
If I may ask, a question to anyone : within Simulink, when I use an Embedded MATLAB function and I call from it another function that outputs a vpi number (Z), I always get the following error : "MATLAB Function Interface Error: MATLAB expression 'Z' is not numeric." what should I do to avoid it the error ? I still need a vpi number from that function.
thanks.

Cavin Dsouza

Cavin Dsouza (view profile)

Remarkable Sir

Jean-Luc

Hey John: thanks for implementing min and max in the latest version!

Nice submission!
Suggest the implementation of generating random numbers in Z_p^*, where p is a vpi number. Thanks!

John D'Errico

John D'Errico (view profile)

You install it the same way you would install any other submission from the FEX. The same for Windoze or any operating system. Unpack the zip file, then put the TOP level directory (VariablePrecisionIntegers) on your search path.

how do i install it on windows ? i have matlab 2008b 64bit, thx.

John D'Errico

John D'Errico (view profile)

Oops. Fixed totient(1) & uploaded.

the cyclist

Great submission. I use it rarely, but it is a life-saver when I need it.

I did find one foible, I believe. Your totient function outputs totient(1) as 0, but both the Wikipedia page you cite and the Sloane Online Encyclopedia of Integer Sequences (http://oeis.org/A000010) say that totient(1) = 1.

I used it for prime number factorization. Works great!

Seth Kosowsky

Thanks John. I did implement w/ repmat and as expected the loss of speed is huge. Still, I have no other way of doing this problem (finding n-digit numbers in base-n with unique digits and special divisibility properties) above base 13. While the repmat is slow, the vpi is doing its job.

Also, I might suggest for future work an equivalent 'isvpi' function.

thanks again.
-Seth

John D'Errico

John D'Errico (view profile)

As far as I know, bsxfun will not work with these tools. That tool is out of my control, as compiled code. I suppose I could overload bsxfun for vpi, but there are a variety of bsxfun replacements on the FEX. One of them might work, as long as they are written for arrays of a general class that supports addition and multiplication.

You can always just use repmat to solve the problem, although it will not offer any gain in speed, since it will allocate the additional memory that bsxfun avoids.

And, finally, the possibility of a new, faster version of vpi is on the horizon.

Seth Kosowsky

Is it possible to get bsxfun to work with the vpi tool box? when I blindly try it, matlab reports bsxfun not defined for vpi types.

Thanks.

kevin

kevin (view profile)

Juliano

Great quality, very easy to integrate. I needed to check for some bugs today that could be related to lack of precision of the standard numeric types, and VPI did the job.

kevin

kevin (view profile)

thank you John!

John D'Errico

John D'Errico (view profile)

I'm sorry, but I don't do consulting in the comments of a submission, as otherwise these things go on forever. More to the point, this very much appears to be a homework problem of some sort, and I don't do homework. I will only repeat my last statement, that analysis always trumps brute force. At the very least, it will allow you to reduce your search space. And I already gave you one hint. Spend some time with pencil and paper, BEFORE you throw brute force search at the problem.

kevin

kevin (view profile)

I ever tried to use "round(sqrt(temps))==sqrt(temps)" as criterion for perfect square, but finally obtained some pseudo solutions for m=17 very soon without any warning message .

kevin

kevin (view profile)

% code I am now using,but it does not work well for m=11
close all
clear
clc
N=1e7;
result=vpi(zeros(1,2));
counts=0;
test256=[0;1;4;9;16;17;25;33;36;41;49;57;64;65;68;73;81;89;97;100;105;113;121;129;132;137;144;145;153;161;164;169;177;185;193;196;201;209;217;225;228;233;241;249];
test10=[0;1;4;5;6;9];
test05=[0;1;4];
test07=[0;1;2;4];
test09=[0;1;4;7];
test13=[0;1;3;4;9;10;12];
test17=[0;1;2;4;8;9;13;15;16];

tic

for ii=2:N
    if counts>0
        break
    end
    for jj=1:ii-1
     if counts>0
     break
     end
        temps=(ii*jj)*(ii+jj)*(ii-jj)/11; % here the m.

       if round(temps)==temps
           if~isempty(find(test256==mod(temps,256)))
               if mod(temps,3)<2 &&(mod(temps,5)<2||mod(temps,5)==4)
       if ~isempty(find(test10==mod(temps,10))) || (mod(temps,4)<2)
           if ~isempty(find(test17==mod(temps,17)))
             if ~isempty(find(test13==mod(temps,13)))
                  if ~isempty(find(test09==mod(temps,9)))
         if temps>1e10
           [root,excess]=sqrt(vpi(temps));
         else
             xxx=sqrt(temps);
             excess=xxx-round(xxx);
         end
      if excess==0 % isSqr(temps)%mod(temps,temp)==0
           if counts>0
             break
            end
          numerator=(ii^2+jj^2)*(ii^2+jj^2);
          if temps>2e15
          temps=vpi(temps);
          numerator=vpi(ii^2+jj^2)*vpi(ii^2+jj^2);
          end
          denominator=temps*4;
           % if gcd(numerator,denominator)==1
               counts=counts+1;
               result(counts,:)=[numerator,denominator]
               toc
               break
       % end % disable
          end
           end
           end
           end
       end
               end
           end
    end
    end
end
toc

John D'Errico

John D'Errico (view profile)

Hi Chen - I'm sorry, but I don't have any reason or need or even the spare time at this point to "test" vpi. I know that it works. And, to be honest, this feels like some sort of homework problem, or maybe a Project Euler (like) problem, found on some web site. Problems like this are not a test of software, they are a test of your analytical skill to reduce it to something more easily solved in a reasonable time. So I'm not going to solve it for you either. I'll give you at most a hint.

Consider the moduli of the set of perfect squares, mod 194. Note that you can compute that modulus even without vpi if you are careful for quite large numbers. What does it tell you if two such elements have the same modulus? Does it help you to find a solution? There is still some search required, but this may reduce the extent of your search.

kevin

kevin (view profile)

hi John

The problem mainly lies the limitation of VPI toolbox's calculation capability.

e.g.

when determining whether a number is a square number, I used "round(sqrt(m)) == sqrt(m)" which does not work for larger number m.

Would you please try to find a rational number n which is subjected to:
√n, √n+97, √n-97 are all rational numbers ? I believe this would be a good verification method for the toolbox's calculating capability.

Michaela

Wow ... I design floating point hardware, and this is just what the doctor ordered. Big Thanks xxx!

Tristan

Very rad ! Worked right out of the box for what i needed, thanks!

John D'Errico

John D'Errico (view profile)

Maybe we are looking at this the wrong way.

Suppose your goal is to generate random numbers, WITHOUT replacement in the interval [0,2^200 - 1].

Yes, the odds of a collision are infinitesimal, but the point is as easily made for numbers in a smaller interval. Here I'll use 20 "digit" numbers, with each digit in base 1024.

n = 500000;
spares = 100;
B = floor(rand(n + spares,20)*1024);

Eliminate any collisions, although there essentially won't ever be any for numbers this large. Unique sorts them of course, so we now extract the number we originally wanted.

B = unique(B,'rows');
ind = randperm(n+spares);
B = B(ind(1:n),:);

Having done this, you can easily convert them into VPI numbers if you want them in a big decimal form for some reason, or simply leave then as they are. Or, suppose you wanted to use a decimal form, with bounds of [0,1e1000-1]. Here we could have used numbers with 100 base 1e10 digits. The same approach applies.

Generating a few spares is incredibly cheap here, so there is no real problem.

In fact, as long as you can factor the upper limit (+1) into a list of small primes, all of which are less than 2^53, one can always use a scheme as I have described above. The only problem is if you wanted to generate samples that lie in an interval like this:

 [0, 12323242124232456456345346]

Of course, I've chosen the number 12323242124232456456345347 carefully, to be a large prime number. This is why my randint code is slow, as it allows any general upper limit. (I can certainly speed that up with some thought though, since I chose an algorithm that may not have been optimal for speed. I'll need to leave that for the future re-write.)

Michael

For this particular experiment, I'm clearly on the cusp of breaking into larger problems with a 2^60 size. Ideally my problem size gets up to 100-200 which, until now, didn't seem possible. It may not be. My problem is not in storing big integers (tried using uint64), but in manipulating them with existing functions. Much support has been added for these numbers, but many functions, like binary/decimal conversion, are limited by the 2^53 size. This is where your tools come in perfectly--if you'll recall you found me trying to convert binary to decimal with some inferior software.

I know a lot of people use this page to gripe about your software not being perfect for what they want -- I was trying to express my gratitude and contribute something that might be useful to you or someone. We're not all out to get you, John. Peace, love, and thanks for the software.

Peacefully Unforced,

Mike

John D'Errico

John D'Errico (view profile)

While I'm happy to see that some find VPI useful, nothing forces you to use a tool that is designed to operate on numbers with hundreds or even many many thousands of decimal digits. If your goal really is nothing larger than that will fit into a 64 bit integer, then use uint64 for your operations. (Only 64 bits is amazingly small in context of the problems VPI is designed to handle.) VPI obviously cannot compete in speed with arithmetic on those numbers, as this is an issue of apples versus oranges.

One day I expect to get around to re-writing VPI, but even then I will not reasonably expect more than about 10-1 speedup in any tests I've done, and even that will force some potential compromises on my part.

Use the right tool for your problem. In this case, uint64 seems like a better choice.

Michael

Here are some stats on the big numbers. I think I will just write my code to ignore duplicates since it's not time effective to remove them. Sorting takes a little while, but I think I can work out of order just fine. So I'm good for now I think, just thought you'd be interested in some run times, and that you'd tell me if something looks off.

>> N=vpi(uint64(2^60))
>> tic,vsamples = randint(N,[500000,1]);toc
Elapsed time is 1138.346514 seconds.
>> tic, vsamplesU=unique(vsamples);toc
Elapsed time is 1811.413225 seconds.
>> tic,vsamplesS=sort(vsamples);toc
Elapsed time is 1833.296390 seconds.

Michael

I see also that unique() was contributed to the toolbox, so with only ~55 collisions per 500,000 samples, I could easily remove the duplicates using unique. I'm guessing that will be more efficient than the re-sampling approach.

John D'Errico

John D'Errico (view profile)

Michael - you do have the benefit of large numbers here. I saw in your question online that you were looking at really large numbers. So suppose we choose to sample 500,000 numbers from the set 1:2^50, but allowing replacements? What are the odds of a collision? This is just a large scale birthday paradox problem. We can compute the probability of a collision in that case:

>> 1-10^sum(log10(1-(1:500000)*2^-50))
ans =
   0.00011102

So the probability of a collision happening is fairly small here.

Had we allowed the random sampling with replacement to go as high as 2^60, which I recall was the basis of your question, the probability of any collision becomes quite a bit smaller. It is not impossible, but certainly a low probability event. You may choose to simply ignore the possibility completely. If you insist on a fully zero probability of any collision, this too is possible. Resampling is easiest. For example, suppose we choose to sample from the set of integers 1:2^50? As long as we stay under 2^53, we can use doubles. Thus

n = 500000;
m = 2^50;
A = ceil(rand(1,n)*m);
A = unique(A);
k = n - numel(A);
while k > 0
  A = [A,ceil(rand(1,k)*m)];
  A = unique(A);
  k = n - numel(A);
end
A = A(randperm(n));

Of course, when the limit exceeds 2^53-1, you could do a similar resampling scheme on vpi numbers using the randint function to generate the samples. unique is defined for vpi.

Michael

I've just been diving into the documentation (thank you for all of that so much, by the way). I see we can generate vectors of random numbers with replacement, but is it possible (or would it be difficult) to generate random numbers from a population without replacement?

Michael

Oh my God I love you.

John D'Errico

John D'Errico (view profile)

Well, I have no idea what is the value of l, since you don't say. However, I tried l = 7. To no surprise, it worked fine, giving these results:

    4
   16
   256
   65536
   4294967296
   18446744073709551616
   340282366920938463463374607431768211456

My guess is, that you are doing something again without thinking about what you did.

LOOK at the output!!!!!!! See that the value of t that is being displayed is a double precision number, NOT a VPI number. So regardless of what you tell me you did, this output is not from the code you posted.

Think about what you are doing BEFORE you ask a question. PLEASE. I am pleading with you.

Nanthini

hi john,
        i want to call a function recursively.. in that vpi numbers are involved...while executing that function code i got following error please suggest me some idea to overcome this problem....thanks in advance....

My code is
n=vpi('59623408128589294906667477489466511394227058875719821124229053681229189216330915080633066080885115044018179752252075074762255857650359062592367314714255870395157878516147169752676083317321479715959914788386131373552660119411485042306204724685302709882300932537368996737625591854420932024858921774454130558350942691706925108621964426007163860852250193402906989838459886518533364378045747452657637521597746398695785988555098978535368265198866036728222641091714452140558062069296423468874123503682189797544270604765655771');
m=vpi(2);
r=vpi(1);
c(1)=mod(power(m,r),n);
t(1)=c(1);
for i=1:l
    c(i+1)=vpi(c(i)).*vpi(c(i));
    if lt(c(i+1),n)==1
        t(i+1)=mod(c(i+1),n);
    else
        t(i+1)=c(i+1);
    end
        disp(t(i+1));
end

The output and error is:
   2

     4

    16

   256

       65536

  4.2950e+009

??? Error using ==> vpi.vpi at 141
If N is a double, it may be no larger than 2^53 - 1

Error in ==> vpi.lt at 46
  INT1 = vpi(INT1);

Error in ==> mdxp at 19
    if lt(c(i+1),n)==1
 

Nanthini

Thanks john...ya you are right here after i will go through the helps as much as i can...

John D'Errico

John D'Errico (view profile)

Sigh. I see I was half asleep when I responded there.

This will do the correct thing, but I was playing around with different examples, and editing my response, so it shows the wrong result in that response.

>> 1+order(vpi(123456))
ans =
     6

And of course, if you actually extract the digits themselves, then digits will succeed.

>> numel(digits(vpi(123456)))
ans =
     6

disp would also work, but you must subtract off the first few blanks. There are always four leading blanks.

>> numel(deblank(disp(vpi(123456)))) - 4
ans =
     6

John D'Errico

John D'Errico (view profile)

Scan through all of the functions in VPI. There really are not that many. Read through the helps. You might find the tools to do this, and you might find other tools that you never even knew existed in there.

Simplest is to use the order function, which returns the power of 10 associated with the highest order digit. Adding 1 to that will tell you what you want.

>> 1+order(vpi(123456))
ans =
     5

Or, since disp actually returns the digits of the number as a string variable...

>> numel(disp(vpi(123456)))
ans =
    10

Or, you could count the digits themselves, extracted into a numeric vector...

>> numel(digits(vpi(123456)))
ans =
     6

Finally, you could compute the log to the base 10, then take the ceiling of that, but since log10 returns a floating point number, this may actually be inaccurate, and there would be further issues for non-positive numbers.

Nanthini

Hi john,
Is there any possibility to find the length of vpi number???
For example n=vpi('59623408128589294906667477489466511394227058875719821124229053681229189216330915080633066080885115044018179752252075074762255857650359062592367314714255870395157878516147169752676083317321479715959914788386131373552660119411485042306204724685302709882300932537368996737625591854420932024858921774454130558350942691706925108621964426007163860852250193402906989838459886518533364378045747452657637521597746398695785988555098978535368265198866036728222641091714452140558062069296423468874123503682189797544270604765655771')
if suppose,i give length(n) it shows ans as 1...please suggest me some idea...

John D'Errico

John D'Errico (view profile)

powermod does not work on character input, only on vpi input. Only the vpi function itself is able to convert a character string of decimal digits into a vpi number. If you cannot bother to convert the number to vpi form, it won't correct your mistake, guessing that you intended the input to be treated as a number. I won't/can't write prescient code that corrects all possible errors. Computers do what you tell them to do. They follow your instructions precisely, not knowing what you wanted to do in your mind when you wrote down incorrect or ambiguous instructions.

Nanthini

I am using powermod function to get the result,to compute mod(M^E,N)where M,E,N should be vpi number right???but for M,i am converting char into ASCII value for example 'hai' is my input after converting to ASCII and also append padding bits it become like this 104097105000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
If i pass this number into powermod function where N=vpi('59623408128589294906667477489466511394227058875719821124229053681229189216330915080633066080885115044018179752252075074762255857650359062592367314714255870395157878516147169752676083317321479715959914788386131373552660119411485042306204724685302709882300932537368996737625591854420932024858921774454130558350942691706925108621964426007163860852250193402906989838459886518533364378045747452657637521597746398695785988555098978535368265198866036728222641091714452140558062069296423468874123503682189797544270604765655771');
E=vpi('877598789786798798477798374982794987349');
then i got following error
??? Error using ==> mtimes
Inner matrix dimensions must agree.

Error in ==> vpi.powermod at 56
    a2 = mod(a*a,n);

Error in ==> encryption at 13
    c=powermod(m,e,n);

I hope the problem here is i have not converted ASCII value into vpi number
How to resolve this error suggest me some ideas.thanks john

 

This is something truly great. Thank you for all your effort you've put into this great "product".

John D'Errico

John D'Errico (view profile)

You are overloading the comments by asking too many questions here, most of which you have not really thought about, or read the relevant help first. Please stop. Send me e-mail instead when you have a problem, or simply read the help.

Regardless, I have no idea what you are asking or what you might have done wrong, since it does work exactly as it should.

>> a=vpi('367366666786786668787987897979797879879');

>> nextprime(a)
ans =
   367366666786786668787987897979797880011

Nanthini

thanks john...if i generate vpi number assigned to some variable like a=vpi('367366666786786668787987897979797879879')
if want to use this number as input to other function such nextprime then if i will give variable 'a' as input its right?? if i give like that it doesn't work??
give me your suggestions for these questions
thanks in advance....

John D'Errico

John D'Errico (view profile)

Read the help for input.

>> inp = input('Please input a vpi number: ','s');
Please input a vpi number: 53634523453456456465766566879675668997876768797865767867567685453433233465463446
>> vpi(inp)
ans =
    53634523453456456465766566879675668997876768797865767867567685453433233465463446

Or, use inputdlg to get input from a dialog box.

Nanthini

hi john,If there is anyway to get interactive input from user for vpi
(i.e)input('enter the vpi no') how to assign that given input to vpi function???
and also one more question,if i generate vpi number assigned to some variable like a=vpi('367366666786786668787987897979797879879')
if want to use this number as input to other function such nextprime then if i will give variable 'a' as input its right?? if i give like that it doesn't work??
give me your suggestions for these questions
thanks in advance....

John D'Errico

John D'Errico (view profile)

To be honest, I never thought of min and max. Yes, it will be easy enough to write those utilities. In fact, I'd not written a min/max facility in the high precision floating point tool I've written either. Just a blind spot on my part.

Jean-Luc

Hi John,

Very nice package, thanks!

A small question: is there a reason that vpi doesn't overload min and max? It seems easy to do but that's why I'm surprised they're missing.

Thanks,

Jean-Luc

Jean-Luc

Nanthini

hi,i used that num2strexact function through the link sent by you,i can't compile it because it requires c/c++ compiler..i have dev-c++ compiler...while running that function it asking to mention the compiler name used. i gave it as dev-c++ compiler. but for that it shows the following error message:
mex('C:\Users\Nanthini\Documents\MATLAB\num2strexact.c')
 
 
Select a compiler:
 
[0] None
 
Compiler: Dev-C++ IDE
Please select from 0-0
Compiler: 0
 
  mex: No compiler selected. No action taken.
 
**************************************************************************
  Warning: The MATLAB C and Fortran API has changed to support MATLAB
           variables with more than 2^32-1 elements. In the near future
           you will be required to update your code to utilize the new
           API. You can find more information about this at:
           http://www.mathworks.com/support/solutions/en/data/1-5C27B9/?solution=1-5C27B9
           Building with the -largeArrayDims option enables the new API.
**************************************************************************
 
 
  C:\PROGRA~1\MATLAB\R2010A\BIN\MEX.PL: Error: No compiler options file could be found to compile source code. Please run "mex -setup" to rectify.
 
 
Well, *that* didn't work ... now trying it with mwSize defined ...
 
 
mex('-DDEFINEMWSIZE','C:\Users\Nanthini\Documents\MATLAB\num2strexact.c')
C:\PROGRA~1\MATLAB\R2010A\BIN\MEX.PL: Error: No compiler options file could be found to compile source code. Please run "mex -setup" to rectify.
 
Hmmm ... That didn't work either.
 
The mex command failed. This may be because you have already run
mex -setup and selected a non-C compiler, such as Fortran. If this
is the case, then rerun mex -setup and select a C/C++ compiler.
 
??? Error using ==> num2strexact at 155
Unable to compile num2strexact.c

Error in ==> check at 5
d=num2strexact(v);
 

why i am getting this error what i have to do to overcome this error..please give your suggestions..
thanks in advance....

  

John D'Errico

John D'Errico (view profile)

That number is not really an integer at all in matlab, so converting to an integer format seems a bit overkill. Or, perhaps you can think of it as MANY integers, all of which are impossible to distinguish from each other. Computing exact integer results for an integer that you cannot resolve seems a waste of time.

>> p = 2.9600e+024;
>> eps(p)
ans =
   536870912

You can think of this as roughly the magnitude of the smallest number that can be added to p and get a different result from that addition. Or think of it as the granularity of p, how much slop or uncertainty there is in determining that value. When eps(p) is greater than 1, then you are unable to effectively resolve p as an integer. Thus

>> eps(2^53 - 1)
ans =
     1

>> eps(2^53)
ans =
     2

As far as how many bits are used for p, log2 tells us that we would need over 80 bits of precision to define p as an integer.

>> log2(p)
ans =
       81.292

But ALL double precision numbers in MATLAB use an internal IEEE format that has 52 bits devoted to the mantissa, a sign bit, and effectively 11 bits in the exponent, to fit into 64 bits for the overall number.

If you insist on converting this number to vpi, you could probably use the function num2strexact, also found on the file exchange.

http://www.mathworks.com/matlabcentral/fileexchange/22239-num2strexact-exact-version-of-num2str

Once p is available as a string, trim off anything that is provided beyond the decimal point, and then vpi can interpret it as an integer.

Nanthini

2960000000000000100663296, can you tell me the number of bits occupied for the above mentioned number...
Thanks in advance...

John D'Errico

John D'Errico (view profile)

For ANY number represented in double precision that is larger than 2^53 - 1, you cannot convert it to a vpi number directly, because you have not provided ALL of the integer digits of that number. The double precision number you have in p is a FLOATING POINT NUMBER. In fact, the actual number stored in MATLAB when you specify p = 2.9600e+024 is 2960000000000000100663296. (Computed using another tool of mine.) This is because MATLAB stores numbers using only a limited number of binary bits. That number is not stored as the decimal integer you think it is. For this reason, I prevent you from creating a vpi number from anything when the number exceeds 2^53 - 1.

>> 2^53 - 1
ans =
   9.0072e+15

Here is a simple way to test that p as you have entered it is not exactly represented as an integer in matlab:

>> p = 2.9600e+024;
>> p == (p + 1)
ans =
     1

You can test for equality too. But note that that test is actually a fuzzy one.

>> isequal(2960000000000000100663296,p)
ans =
     1

>> isequal(2960000000000000100034423,p)
ans =
     1

>> isequal(2960000000000000000000002,p)
ans =
     1

>> isequal(2959999999999999999999375,p)
ans =
     1

See that when you add 1 to p, the number has not changed! In fact, I could change it by quite a lot and still not change it at all as far as MATLAB knows. This is not a bug in MATLAB as many novices think, but a fact of life when you use double precision. It is also one of the reasons I wrote vpi, because the precision of a number in matlab is limited in double precision.

In order to create that corresponding exact vpi number, you need to supply all of the digits of that number. This is why I allow you to provide a string input for vpi as an integer, so you can create numbers larger than 2^53-1.

>> p = vpi('296000000000000000000000');
>> gcd(67,p)
ans =
    1

Now however, try the add 1 test on the true vpi number.

>> p = vpi('296000000000000000000000')
p =
    296000000000000000000000

>> p + 1
ans =
    296000000000000000000001

>> p == (p + 1)
ans =
     0

Of course this all works exactly as arithmetic should, because vpi is handling the details.

Nanthini

i tried to find gcd of following two numbers but i got the following error will please explain why i got this error
e=67
p= 2.9600e+024

??? Error using ==> vpi.vpi at 141
If N is a double, it may be no larger than 2^53 - 1

Error in ==> test at 23
z=gcd(e,vpi(p));

thanks in advance..

Nanthini

thank you so much.... i got the output...:)

Nanthini

sorry i am using MATLAB R2011a mistakenly
i typed as R2021a...

John D'Errico

John D'Errico (view profile)

I tried your test case.

>> int1=vpi('4536788992900999999999999999999999999423456789072434546775885959');

>> int2=vpi('2345678998765432112345678987654323456789755654433245567787889899');

>> t=gcd(int1,int2)
t =
    1

As I expected, it works perfectly fine. Now, admittedly, I don't have version R2021a, because this is only 2011, and that version will not be released for another 10 years or so. I do have release R2011b. :)

My guess is that you have done something inappropriate when you downloaded the files. Most likely, you took some or all of the m-files from the @vpi directory. I'm at a complete random guess, based on the class error message, that probably cannot happen unless you moved some files around. Perhaps you put the @vpi directory on your search path, another no-no, which the readme file explicitly says not to do. If I might repeat what is said in the file ReadAboutVPI.rtf:

"Do NOT put the @vpi directory on your search path!!! Do NOT take any functions from that directory, as this will break the operation of the vpi code!!!"

Best seems to completely delete the VPI files you have downloaded, and re-install, this time following the instructions provided.

Nanthini

hello i want find gcd of two larger numbers for that i used vpi but i got some error...will please explain why i got this kind of error
my input
int1=vpi('4536788992900999999999999999999999999423456789072434546775885959');
int2=vpi('2345678998765432112345678987654323456789755654433245567787889899');
t=gcd(int1,int2);
disp(t);

but i got error like this
??? Error using ==> class
The CLASS function must be called from a class constructor.
Error in ==> vpi at 95
INT = class(INT,'vpi');
Error in ==> check at 1
int1=vpi('4536788992900999999999999999999999999423456789072434546775885959');

My MATLAB version is R2021a....

John D'Errico

John D'Errico (view profile)

M = vpi(2)^100;
M =
    1267650600228229401496703205376

Modular operations work as you would expect.

mod((M+1)^263-1,M)
ans =
    0

Modular multiplicative inverses work, as long as the inverse exists, minv should compute an inverse.

P = minv(3,M)
P =
    845100400152152934331135470251

mod(3*P,M)
ans =
    1

Lea Jade

Hello, I am new to matlab and so far I have had an amazing time using this toolbox and now I would like to use it to perform arithmetic within a galois field. The maximum power MATLAB supports is 2^16, how do I use VPI in order to have values like 2^100 and still be able to perform the arithmetic with the GF field. Thanks alot

Sarla

Sarla (view profile)

:) that disp was a later addition,since i was trying to forcibly stuff an enormous number into a vpi variable. i think the error is in the variable out. i am still running a few tests.

actually i rather like the randint function. it made my life slightly simpler.

thank you!!

John D'Errico

John D'Errico (view profile)

No. You don't need to convert an existing vpi number into vpi form using disp.

If k and PointOrder are already vpi numbers, then there is no need to use disp on them. However, the error message that is returned from vpi in your original question indicates that prod (or one of the arguments in that line of code) was NOT a vpi number.

What happens when two numbers are added (where one is a vpi number) is this:

I attempt to convert both numbers into vpi form, then I do the add. Try an example.

>> A = vpi(123456);
>> A + 12
ans =
    123468

vpi looks at this and realizes that A is already a vpi number, but that 12 is not. So it internally converts 12 into vpi form, then adds the digits, doing any necessary carries. All of this works nicely and completely transparently to you. But what happens when we try this on a really large number added to A?

>> A + 1e37
??? Error using ==> vpi.vpi at 160
If N is a double, it may be no larger than 2^53 - 1

Error in ==> vpi.plus at 61
      INT2 = vpi(INT2);

See that when vpi tries to convert 1e37 into a vpi, it fails. This is because 1e37 is too large to represent in MATLAB as an exact floating point integer, a FLINT. I don't (can't) know ALL of the digits of 1e37, since only roughly 16 of them are available. The limit for this is in fact 2^53-1 as a double.

>> 2^53-1
ans =
   9.0072e+15

Anything larger than that, and vpi gets a headache, since MATLAB does not store sufficient information to recover all of the digits to the left of the decimal point.

The error message that you showed says that one of the arguments to plus was indeed not a vpi number, and that it was too large to convert exactly to a vpi number. When this happens, vpi must return an error. What else can it do?

As far as randint goes, it should be reasonably efficient as these things go, but I had to be quite careful to be positive that the digits were truly uniformly distributed. I'll look again to see if it can be improved for a bit more speed though. If I try too hard to gain speed there, I may end up damaging the pseudo-randomness of the result. For example, suppose I wanted to generate a uniformly distributed random integer as large as 1e100 as a vpi?

I cannot simply generate a single random number using rand, and then scale it up to 1e100. The digits would not be uniformly distributed. All numbers in that range would not be possible and equally likely. So vpi/randint actually generates those digits independently, in a way that will result in a uniformly distributed number, and where all such numbers are equally likely in theory. But if you want a number with 100 random digits, vpi must generate roughly 100 random numbers along the way. This makes it slower than you like, but it has the desired properties.

Sarla

Sarla (view profile)

ok i should have been clear. my mistake. both invk and prod (i'll change the name) are vpi variables. those values were arrived at after some calculations(i copied those values from the command terminal). therefore i do not have the option of providing those values to vpi as strings.
my actual code is like this:
  invk=minv(k,PointOrder)% k and PointOrder are user defined vpi variables
    
  prod1=vpi(da*r)% da and r are user defined vpi variables

    s=vpi(mod(vpi(disp(invk))*vpi((disp(out)+disp(prod1))),vpi(PointOrder)))

in a case like the above one what can i do?

also regarding the randint function i am using it to generate private keys when dealing with vpi variables! but my slow processor takes too long to perform certain calculations. so now i use a mix of randint and randomnum in my program depending upon the precision i need.

John D'Errico

John D'Errico (view profile)

To answer the question about other random number generators, randint generates uniformly distributed random VPI integers. As such I have found it of particular value for factoring and primality testing. That was in fact why I wrote randint.

To generate numbers with other distributions, one needs to specify the distribution itself. And if you are working in vpi form, the numbers must be integer. Since there are not too many probability distributions in common use where the variates are integer (ok, Poisson is one) I'm not sure what I might provide that enough people would find of use.

John D'Errico

John D'Errico (view profile)

When you write things like this:

invk = 5431182971659226082032497739807769284499946472968309173822

prod =
66634193303502359521887744541301070677802959088568215854843250

Matlab creates DOUBLE precision numbers called prod and invk. (BY the way, the name prod is a poor choice, because it overloads the prod function, preventing you from later using the prod function.) Next, you try to do the computation:

s=vpi(mod((vpi(invk)*vpi(out))+prod)),vpi(PointOrder)))

Here matlab tries to add prod to a vpi number, converted from invk. In order to do this, it must convert invk and prod into vpi form. But those numbers are too large to convert exactly into a vpi number. Oops. This is your error. See that

invk =
   5.4312e+57

log2(invk)
ans =
       191.79

So invk is a double precision number, with far too many digits to represent as an integer in standard matlab. The limit in Matlab for a floating point integer (sometimes called a FLINT) is 2^53-1. But vpi can handle the numbers, if you call pass them in properly.

You can define them as vpi numbers directly, calling vpi with a character string of digits.

invk = vpi(' 5431182971659226082032497739807769284499946472968309173822')

prod = vpi('
66634193303502359521887744541301070677802959088568215854843250')

By providing the numbers to vpi as character strings, it can convert them directly into VPI form with NO loss in precision.

Sarla

Sarla (view profile)

hi
i am having an issue with performing certain calculations on vpi numbers.
i have a vpi variable
invk = 5431182971659226082032497739807769284499946472968309173822

i have another variable
prod =
66634193303502359521887744541301070677802959088568215854843250

now i have this expression
  s=vpi(mod((vpi(invk)*vpi(out))
+prod)),vpi(PointOrder)))

now the program works fine until the value of prod. but it comes to a complete halt when it comes to s.
and i get an error of the form :
??? Error using ==> vpi.vpi at 141
If N is a double, it may be no larger than 2^53 - 1

Error in ==> vpi.plus at 51
      INT1 = vpi(INT1);

Error in ==> dsaforlargenumbers at 98
    s=vpi(mod(vpi(invk)*vpi((out+prod)),PointOrder))

could you help me out?

Sarla

Sarla (view profile)

yup you are right. i did not really explore isprime for vpi.

is there any function other than randint to generate pseudo random numbers for vpi?

John D'Errico

John D'Errico (view profile)

You could count the number of factors of a number to determine primality. This would be efficient IF the number was prime, since the virtually first thing factor does is a call to isprime. When it is not prime however, factor will take longer, since factoring numbers is relatively time consuming.

So in the end, just call isprime if all you wish to do is determine primality of the number.

John D'Errico

John D'Errico (view profile)

Gabriel - The exponent in powermod can be rather large, a vpi number itself. For example...

tic,powermod(vpi('124142565235674676762435233523'),vpi('1234562425532452245524323146735634323443256767474343446454564657645466435556453557454646578567453543658567536457567346878067852'),vpi('137457798574858')),toc

ans =
    101407725757733

Of course, since powermod is the workhorse in many of the tools of vpi, such as isprime, modroot and factor, it needs to be efficient.

Gabriel

@Sarla

Use isprime(vpi(p))?
It's probably faster than factoring the whole thing. It return 1 when the number is prime.

Sarla

Sarla (view profile)

:) so even the stupidest of questions have their uses!
i realized that i could use the factor in vpi to test whether or not a number is prime. like this:

fac=factor(vpi(p));
if length(fac)==1
    thetruth=1;
else
    thetruth=0;
end
i.e. if the length were 1 the number would be prime! this way i don't really have to bother with primality tests!

Gabriel

Gabriel

-_-

Lesson learned... Thank you =)
Is there a limit to the exponent in powermod?

John D'Errico

John D'Errico (view profile)

Gabriel. Yes, you can use the code as written. But you might also be less than happy having done so. I would point out that there already exists a powermod tool in vpi!!!!! How does yours compare?

I saved the code you wrote, calling it powermod2, then ran it on a simple comparison.

>> tic,R = mod(vpi(2).^12345,vpi(9937)),toc
R =
    8315
Elapsed time is 0.457709 seconds.

>> tic,R = powermod(vpi(2),12345,vpi(9937)),toc
R =
    8315
Elapsed time is 0.024054 seconds.

>> tic,R = powermod2(vpi(2),12345,vpi(9937)),toc
R =
    7516
Elapsed time is 0.068513 seconds.

See that my own powermod code produces the correct, same result.

Your code is slower. And it gets the answer wrong. Looking to see if a code exists to solve your problem is always a good idea. Testing your code for correctness is ALWAYS just as important.

Gabriel

Thank you for this great package. I am definitely not a pro programmer but I wrote a little function because I needed something like powermod from the Matlab symbolic toolbox. I might have made mistakes, but it seemed a lot faster to use this to compute the modulo of a power.
function result = powermod(x,power,modulus)
intermediateResult = x;
flag = 1;
for i=1:32
   intermediateResult = mod(intermediateResult^2,modulus);
   if(bitget(power,i))
      if(flag)
         moduloResult = intermediateResult;
         flag = 0;
      else
         moduloResult = moduloResult*intermediateResult;
      end
   end
end
result = mod(moduloResult,modulus);
end

John D'Errico

John D'Errico (view profile)

You got me thinking that I should include a num2str function for exactly that purpose though, which would have been more obvious to a user.

Sarla

Sarla (view profile)

:)that worked perfectly. sorry. should have worked on it a bit more before posting up a silly question like that. thank you!

John D'Errico

John D'Errico (view profile)

Sarla - you are wrong that there is no conversion. Try this:

>> N = vpi(12)^37
N =
    8505622499821102144576131684114829934592

>> nchar = disp(N)
nchar =
    8505622499821102144576131684114829934592

>> whos nchar
  Name Size Bytes Class Attributes
  nchar 1x44 88 char

See that disp does the conversion, and if you provide an output argument, it gives you exactly what you want.

Sarla

Sarla (view profile)

hi
i was actually trying to display a vpi value in my matlab gui's text box. however i was unable to do that. i think there is no conversion from vpi back to string. so i keep getting an error message. is there any way i can work around this?

John D'Errico

John D'Errico (view profile)

Finding the factors of a significantly large number, in this case, one with over 600 decimal digits, can be extremely difficult. If these problems were easy, then the schemes used in modern methods of encryption would be easily broken.

The factoring scheme in vpi/factor is able to handle numbers with many digits, but there are limits, and I have not coded up the best possible method, but only a reasonably good one. (It is still far better than that offered by the default factor method in MATLAB. One day I hope to offer some of the better methods as an alternative.) If a very large factor is returned, then you can always test to see if it is truly prime using isprime, or if factor simply gave up. For example,

F = factor(vpi(2)^113 - 11)
F =
    3 7 7 101 699440541326170624170606362123

isprime(F)
ans =
     1 1 1 1 1

Zakir

Zakir (view profile)

I want to factorize very large no
Example:
a = vpi(2);
b = a^2067;
c = b-7;
How to sactorize it?
When i factorize by factor function then matlab display only those factors which are less then 2^32;
and the remaning value as a factor?
What is the Problem?
and How to handle it Please?

Zakir

Zakir (view profile)

Very Very Nice Job
Thanks for This.

Sarla

Sarla (view profile)

okay. thank you anyway!

John D'Errico

John D'Errico (view profile)

Hmm. I'd not implemented inf (or for that matter, nan) into vpi. For example, there seemed no reason to define an inf, as a vpi number has essentially no dynamic range restrictions on it anyway. I'll take a look, as admittedly I have occasionally thought that these extensions would be appropriate.

Sarla

Sarla (view profile)

hi
i seem to be facing a problem assigning inf to a vpi variable.

i have a variable P declared as vpi, and another condition wherein if this variable ever equals inf then a certain value is returned.
the problem is whenever i try to write something of the form
if P==inf
i get an error of the form

Error using ==> vpi.vpi at 139
N must be a finite number

Error in ==> vpi.eq at 44
  INT2 = vpi(INT2);

any suggestions as to how i can get around this?
p.s. the inf is sort of essential for the program.

Sarla

Sarla (view profile)

this submission is wonderful for working around with huge numbers, and after struggling for days with implementing cryptographic algorithms and being thwarted at every point, i finally stumbled upon this. thank you. bless you!

Evgeny Pr

Evgeny Pr (view profile)

John D'Errico

John D'Errico (view profile)

When you download this toolbox, or ANY toolbox that creates a new class, NEVER move the included files around. Do NOT put the @vpi directory on your search path. Do NOT place the downloaded directory in the Mathworks provided hierarchy of toolboxes. Use some other directory that you provide to store all downloaded files.

All that you should do is unzip the download zip file, and then put the VariablePrecisionIntegers directory on your search path.

This file would solve a lot of problems for me if I only could get it to work. No matter what I do, I keep getting this error message:

??? Error using ==> class
The CLASS function must be called from a class constructor.

Error in ==> vpi at 181
INT = class(INT,'vpi');

What am I doing wrong. I've tested just copy pasting some of your own examples and all I get is that message. What am I doing wrong? Do I need to have more files than the vpi.m in the same directory as the file I'm working with?

/Gustaf Lindberg

Joseph Kirk

Joseph Kirk (view profile)

This is quite possibly the best submission I have downloaded from the File Exchange. Great work John!

Harriet Fell

Rob Campbell

Rob Campbell (view profile)

John D'Errico

John D'Errico (view profile)

Thanks Ian. I've sent in an update to fix the legendresymbol bug. John

Ian

Ian (view profile)

Great work!

I got this as a reference and guide as I attempt to write my own Quadratic Sieve.

I noticed an issue on your legendresymbol. I think you have the 'iseven' check under p==2 backwards.

for instance: legendresymbol(87463,2) = 0. I expect '1'.

Thanks for such a great toolbox.

Nathan Greco

Guillaume

Just a small word to say that I really like this package, very handy !

However, a typo in the docs says that randint generates between 1 and N, which is wrong. It generates numbers between 0 and N.

Moreover, I would love a fast function similar to biterr (that tells the bit error rate of 2 integers). At the moment, I compute it with some sum and vpi2bin function, but it is quite slow.

Thanks hor this work !

John D'Errico

John D'Errico (view profile)

One of the things I never overloaded was colon. I briefly considered doing something like that, but logically, there is no need to do so. Consider...

Suppose the arguments to colon are small numbers? In that case, you can simply use double to convert the numbers to normal integers, then use colon on that. Even then, don't expect a loop like

for i = 1:2^53
   stuff
end

to terminate any time soon in matlab.

The second case is if the numbers are true vpi numbers, that way for a good reason. For example, in the case you pose, this loop will not terminate before the universe undergoes heat death, even on the fastest super-computer in the world, even on a cloud of the 1000 fastest super-computers in the world. You must appreciate that 2^1024 is a number with over 300 zeros in it.

log10(vpi(2)^1024)
ans =
       308.25

In fact, vpi(2)^1024 slightly exceeds the capability of vpi2english!

vpi2english(vpi(2)^1000)
ans =
ten novemnonagintillion, seven hundred fifteen octononagintillion, eighty six septennonagintillion, seventy one sexnonagintillion, eight hundred sixty two quinnonagintillion, six hundred seventy three quattuornonagintillion, two hundred nine trenonagintillion, four hundred eighty four duononagintillion, two hundred fifty unnonagintillion, four hundred ninety nonagintillion, six hundred novemoctogintillion, eighteen octooctogintillion, one hundred five septoctogintillion, six hundred fourteen sexoctogintillion, forty eight quinoctogintillion, one hundred seventeen quattuoroctogintillion, fifty five treoctogintillion, three hundred thirty six duooctogintillion, seventy four unoctogintillion, four hundred thirty seven octogintillion, five hundred three novemseptuagintillion, eight hundred eighty three octoseptuagintillion, seven hundred three septenseptuagintillion, five hundred ten sexseptuagintillion, five hundred eleven quinseptuagintillion, ...

It continues that way for quite a while too. This is a seriously BIG number. It is larger than the number of elementary particles in the entire universe, by MANY orders of magnitude.

So suppose we made it a wee bit simpler. How about

for i = vpi(2)^1023:vpi(2)^1024
    stuff
end

In theory, this loop would terminate in 1/2 the time as would your loop. So what? It will still take an essentially infinite amount of time.

So consider a much simpler problem, at least in theory.

for i = 1:vpi(1e21)
   stuff
end

That number only has 21 zeros in it. Even on a supercomputer capable of 1 petaflop performance, this will take on the order of some days in CPU time to execute. On your laptop, expect it to take years of time.

Ok, so maybe that one was too much to do. How about this simple loop?

N = vpi(10)^21;
for i = N:(N+100)
   stuff
end

While I might consider allowing you to do this, nothing stops you from doing something almost as simple.

N = vpi(10)^21;
for j = 0:100
   i = N + j;
   stuff
end

This last case is the only one I have shown that is possible on a normal computer in a finite amount of time, and it is trivially possible already as I show.

So this is why I never provided the facility to do as you wanted. There was never any reason.

Regards,
John

Hi Jhon,

I'm new using your vpi toolbox, and it caught my attention. Well, making tests I realized that using vpi type in a "for" of the form i=1:vpi(2)^1024, the result sends me the message: "Undefined function or method 'colon' for input arguments of type 'vpi'." .

 And my question is, the current version of vpi toolbox don't handle this task ? Or there is another way to solve this ?

In advance thanks for your attention.

John D'Errico

John D'Errico (view profile)

Sorry about that. Only the binomfactors code uses consolidator, and only nchoosek calls binomfactors. I never caught it in testing, since consolidator is on my search path.

Anyway, and would be more efficient to use a call to unique there, then accumarray anyway, since no tolerance is employed.

I've submitted a new release. It will be there on Monday morning.

John,

The demo_vpi.m needs your consolidator function:

http://www.mathworks.com/matlabcentral/fileexchange/8354-consolidator

Great package.

Derek O'Connor

Khaled Hamed

More fun with vpi2english if you use, for example,

speak(vpi2english(vpi(2^52)))

The function "speak" is summission File ID: #4769 named "Speak".

There are a also a couple more "speak" submissions around.

Khaled Hamed

Another good reason to use higher preciesion: the famous subtractive cancellation error

>> x=77617;y=33096;z=33375*y^6+100*x^2*(11*x^2*y^2-y^6-121*y^4-2)+550*y^8

z = -1.5112e+023

>> x=sym(77617);y=sym(33096);z=33375*y^6+100*x^2*(11*x^2*y^2-y^6-121*y^4-2)+550*y^8

z = -200

>>x=vpa(77617);y=vpa(33096);z=100*x^2*(11*x^2*y^2-y^6-121*y^4-2)+33375*y^6+550*y^8
 
z = -200.0

>> x=vpi(77617);y=vpi(33096);z=33375*y^6+100*x^2*(11*x^2*y^2-y^6-121*y^4-2)+550*y^8

z = -200

vpi is as accurate as symbolic and vpa, roughly four times faster than sym and eight times faster than vpa

John D'Errico

John D'Errico (view profile)

Hi Gwendolyn,

I was going to submit a new version today again. The simplest thing to do is to save that file in a form that is compatible with older releases that go back to version 7.

John

Hi John,
I like your submission, and it works fine for my puposes even in R14. Using your old factor routines, as the new factor command uses a mat file (_primeslist_), which the old matlab version can not read. Is it possible to provide a compatible mat file for MATLAB R14?

Thanks

John D'Errico

John D'Errico (view profile)

Ben - I've submitted an update, that now handles the pathological cases for plus and times.

John D'Errico

John D'Errico (view profile)

John - I think the problem is shown by your path. Are you using an older release of MATLAB? Perhaps release 11? If so, I'll not be at all surprised at seeing a problem. I would not expect vpi to run on older versions, and I have not tested it there.

If you actually have a more recent release, then please contact me via mail and I will help to resolve the problem.

John Edwards

Hi John,

I am a beginner with Matlab. However, your creation caught my eye and I downloaded the package. I used very simple examples but MATLAB kept returning the same error:
---------------------------------
A = vpi(17)^17
??? lseif ~ischar(N) &&
                    |
Missing variable or function.

Syntax error in ==> C:\MATLABR11\toolbox\VariablePrecisionIntegers\@vpi\vpi.m
On line 51 ==> elseif ~ischar(N) && ~isnumeric(N)
----------------------------

I am not sure how to resolve this error. Any ideas?

Thanks

John Edwards

Ben Petschel

John, you've done a great service to the mathematical community (esp. number theorists and Project Euler enthusiasts!)

In the vpi/add carry operation you could improve the performance of pathological cases such as 1+vpi(repmat('9',1,1e6)) by changing the while loop to a for loop.

Neha P

Neha P (view profile)

Dear Errico,

The tool box you have created in order to store large integer numbers has been very helpful for our project. Currently we have encountered a problem to store large floating numbers.
I would request you to give me any information regarding the above problem.

Lissa

Lissa (view profile)

John D'Errico

John D'Errico (view profile)

I just uploaded a new release, fixing the bug in vpi2english.

Great submission, but there seems to be a bug with vpi2english for numbers from 1 to 99:

>> vpi2english(vpi('1'))
??? Index exceeds matrix dimensions.

Error in ==> vpi.vpi2base at 114
    B(i) = base2dec(fliplr(cdigits(ind)),10);

Error in ==> vpi.vpi2english at 53
base1000 = vpi2base(N,1000);

John D'Errico

John D'Errico (view profile)

Sorry. I'll fix the demo.

Giampiero Campa

Giampiero Campa (view profile)

very cool,
however i think there's something wrong with line 129 of the demo:
>> iseven(2^127)
it should not be just a double, and yet is too large to be a vpi ...
it gives an error in r2009a anyway.

perhaps you meant iseven(vpi(2)^127) ...
Giampy

John D'Errico

John D'Errico (view profile)

I've done a more careful investigation into the question. Can I in fact use an integer format for the internal storage?

As it turns out, I cannot do so, at least not without seriously compromising the capabilities of this toolbox. First, the question is, can I use INT8 for storage? The answer here is no. While Bill checked to see if it would work, he apparently did not do a multiply between a pair of vpi numbers. vpi uses a nice trick to multiply a pair of vpi numbers. It uses conv, whereby it gains a massive increase in its multiplication speed. And since multiplication is necessary for these tools, the loss of conv would be rather bad for performance.

Sadly, conv does not support int8 inputs. Even worse, int8 would cause a dramatic restriction on the size of the numbers that vpi can handle. try this simple test:

N = [9 9];
conv(N,N)
ans =
    81 162 81

So were I to store the number 99 as a vpi number that uses int8 as a storage mechanism, then I could not actually multiply 99*99 using conv, even if int8 was supported as an input to conv. The conv result would overflow int8 storage.

This suggests an alternative. Perhaps I could use single as the internal storage for the digits. In fact, this offers some minor gain. Use of single would cut the memory requirements for vpi numbers in half. Even better, in my tests the single operations on integers are actually a bit faster than they are for doubles.

So for a few moments, I thought I could at least use single. The problem is, it too will cause a problem. This would limit vpi number to about a maximum of roughly 100,000 digits. Again, the problem is numeric overflow in multiplication using conv. While numbers that large might seem large, I'd still prefer to have essentially no limit on the number of digits.

One other option is to use a smaller form for storage, but then to automatically convert from int8 to doubles and then back (internally) for all operations. I could do this, but it seems risky, and likely to cause bugs.

So for the present time, I'm forced to stick with the internal storage as it is.

John D'Errico

John D'Errico (view profile)

Thanks Bill. I'll fix the examples, and look into the question of changing the internal storage format. It is something I've considered before. My initial reason for leaving them all as doubles internally was to get the maximum speed. I want to minimize the number of events where I need to force a type conversion between int and double.

Bill McKeeman

I am impressed the scope and presentation of this toolbox. Very nice work. I have a couple of nits to comment upon.
The last example in vpi.m seems to have lost its text.
I recommend that the examples be such that, if cut and pasted to the command line, run. I did this with my World Time Zones FX entry in response to a criticism and was pleased at the result.
I experimented with replacing
 INT.digits = double(fliplr(N))-48;
with
  INT.digits = int8(fliplr(N))-int8('0');
The package still seems to run and takes 8x less storage for the vpi numbers.

Dear Alain,

Thanks for your reply and correct answer mean(x) = 4.661273948537807e-005, (or
sum(x) = 7680*4.661273948537807e-005 = 0.3579858392477036 (16 digits) )
which is the mean of the data at https://hpcrd.lbl.gov/SCG/ocean/NRS/etaana.dat.

I meant to ask for the sum(x) not mean(x) because this sum(x) was discussed in
the paper by Yun He and Chris H.Q. Ding, ``Using Accurate Arithmetics
to Improve Numerical Reproducibility and Stability in Parallel Applications'',
Journal of Supercomputing, Vol 18, No 3, Mar 2001.

 http://crd.lbl.gov/~yunhe/papers/HPA_JSC.pdf

In what follows the 7680 d.p.f.p. numbers in etaana.dat are stored in ssh.mat

They give exact sum(ssh) = 0.35798583924770355224609375 (26 digits).

I think you should publish your code and also how the condition of any sum(x)
can be calculated. I have met many people who are surprised that a
simple summation can give very bad results.

Matlab gives sum(ssh) = 3.441476821899414e+001, which has a relative error
of 9.513444009773050e+001, nearly 2 orders of magnitude wrong.

Using the 1-norm my condition number for sum(ssh) is
k_1 = n*sumexact(abs(ssh))/abs(sumexact(ssh)) ~ 10^21.

If I use the Matlab sum (ssh) I get k_1 ~ 10^16,
which is not correct but still gives a warning that
this problem is very ill-conditioned.

Here is my Matlab implementation of D. Priest's method,
which gives the correct result, rounded to about 16 digits :

===========================================================
function snew = SumP(xorig);

% Using Douglas M. Priest's thesis recommendation.
% D.M. Priest: "On properties of floating point arithmetics:
% Numerical stability and the cost of accurate computations",
% Mathematics Department,
% University of California at Berkeley, CA, USA (1992).

% This function has complexity O(n log n) because of the sort.
% This is the price we pay to get the correct sum.

% Written by Derek O'Connor March 2006.

% --- Sort xorig in decreasing abs(x) order

    [absx absord] = sort(abs(xorig),'descend');
    x = xorig(absord); % |x1| >= |x2| >= ... >= |xn|

% --------- Compensated Summation --------
    sold = x(1);
    c = 0;
    for k = 2:length(x)
       y = x(k) + c;
       snew = sold + y;
       u = snew - sold;
       c = y - u;
       sold = snew;
    end;
% --------- end of SumP -----------------

============================================================

Here is the result using Priest's method:

>> SumP(ssh)
ans =
    3.579858392477036e-001

Here is the result using Rump's Intlab 5.5 set to 30-digit precision:

>> p = 30; display(SumRump(ssh,p),p)
long ans =
  3.579858392477035522460937500000 * 10^-0001

It seems that this Sea Surface Heights problem needs at least 26-digit
precision to obtain the exact sum.

Intlab is here : http://www.ti3.tu-harburg.de/rump/intlab/

Regards,

Derek O'Connor

Alain Barraud

vpi is an excellent toolbox, thanks a lot.
I shoud add a comment to derek O Coonor question. I have found 4.661273948537807e-005 as the mean of the data in "etaana" without any toolbox just using floating poitn computation within Matlab. I have used a function I have written some years ago which computes full precision sums of any condition number (here 1.5e17) in a very simple and efficient way.

MATLABician

Dear Kevin,

You "don't understand the utility of this function"? First of all, this submission is not a function, it's a toolbox. I suggest you read John's reply to get a glimpse of its utility, but even before that I suggest, as he did, that you learn some mathematics. That might help you to understand the utility of this toolbox, his other submission MIN2, MAX2, whose utility, too, you have had trouble understanding, and many other of his submissions for that matter - before you go around commenting on their utility. And after that, I suggest you go through the code. I'm sure you can learn one or two things about MATLAB.

Omid

PS. As always, great submission John. Thanks for sharing this with the community.

Kenneth Eaton

Kenneth Eaton (view profile)

What a marvelous tool! Thank you John!

Dear Kevin,

Please calculate the average of the following set of 7680 floating point numbers to 3 "significant" digits.

https://hpcrd.lbl.gov/SCG/ocean/NRS/etaana.dat

Please post your result and tell us how you did it.

Yours sincerely,

Derek O'Connor

John D'Errico

John D'Errico (view profile)

Dear Kevin,

(Sigh.) Just because you personally don't ever need more than 3 digits in any computation does not mean that no one else will ever have such a need. And there is NO reason to trash the useful work of another just because your own point of view is so specialized.

This toolbox does not have a specific purpose. As a toolbox, it should NOT be so specialized that others cannot use it for their purposes. And as an author of dozens of toolboxes over the years, I am continually surprised at how others have used the utilities I've offered, in a variety of ways. But if you need a reason for something to exist, I'll offer a few such reasons.

Surely the most well publicized reason to work with large integers is public key encryption. Here one must find large prime numbers (with very many digits) and work with them.

http://en.wikipedia.org/wiki/RSA

On the other end of encryption is the task of breaking such a code. Here tools are needed to factor large integers. While the specific tool I've provided to factor integers is only capable of factoring integers with a few dozen digits with any regularity, it is far better than the very limited tools provided in MATLAB itself. I've also begun to work on other factoring tools, with the potential to solve larger factoring problems. Admittedly, the best such tools would use multiple CPUs in such a task, parcelling out the work. But, this is something that is easily done in the factoring of large integers. Read about the quadratic sieve methods for factoring.

http://en.wikipedia.org/wiki/Quadratic_sieve

These are very useful methods that can succeed in factoring huge numbers. Again, they require the ability to use multiple CPUs, all working in parallel, but this is something that is now becoming quite feasible. In fact, I've even written a simple working quadratic sieve factoring code, although it is far below the standards I put on my own work to release.

This toolbox is not (yet) capable of working with the speed necessary to factor numbers with several hundred digits. But it can serve as a useful testbed for algorithmic development in that field, and as an easy, inexpensive way to learn some techniques for such work.

Are there needs to work on the breaking of codes? Perhaps World War II might have lasted a bit longer had those mathematicians working with Alan Turing not done their jobs. Read about codes and cyphers.

http://www.codesandciphers.org.uk/

The National Security Agency is one such place where many mathematicians work.

http://www.nsa.gov/CAREERS/

If this set of tools existed for no reason other than primality testing and factoring, I'd suggest it has some merit. But there are other tools I've provided inside that may offer value. Solutions of linear Diophantine equations are frequently requested on the newsgroup, so I added a simple tool to do so. It is one that I will look to improve (given time) with future releases. As well I've begun work on Pell equation solvers. These last tools are of interest to number theorists, so perhaps they are boring to you. All of them tend to use continued fraction approximations to numbers, approximations to fractions, approximations to square roots, etc. Number theory is itself a topic that offers careers to mathematicians worldwide.

http://en.wikipedia.org/wiki/Number_theory

As well, number theory is perhaps one of the prettiest branches of mathematics. It offers a way for a student to become interested in mathematics, with easily formulated problems with correspondingly beautiful solutions.

In my own case, I spent my first several years as a student mostly interested in number theory, before gravitating to more prosaic fields like numerical analysis, statistics, modeling and applied mathematics in general. A student who becomes entranced with mathematics, playing with such an easily formulated problem like Goldbach's conjecture, may then move into other fields of mathematics, perhaps knot theory, differential or algebraic geometry, statistics, or even diverge into some branch of engineering. And the fact is, I've had the need for results from diverse fields of mathematics in my career as a consulting mathematician.

You state that other tools exist that can serve the same purposes as this toolbox, i.e., Maple and Mathematica. While they do exist, and serve that purpose very well, not every individual has the money to pay for those tools. For example, I don't own the symbolic toolbox, nor a copy of Mathematica. For those individuals with an interest in the problems this toolbox can solve, they can do their work for no additional investment.

Finally, this toolbox can serve to teach others methods they can use to build their own toolboxes, a template of sorts. Alternatively, one could easily expand these tools to work with variable precision floating point arithmetic, or with fractions of arbitrary size. When enough people have the current matlab releases, I'll convert this toolbox to use the newer methods now found in matlab for working with classes and OOP code.

So I am truly sorry that this tool is of no interest to you. But that only says what limited horizons you live within. Just for once, expand those horizons. Learn some mathematics.

John

Kevin

Kevin (view profile)

I don't understand the utility of this function.

Sure its neat to have hundreds of digits for some esoteric project, but what is the goal. Give an example of the usage of this tool suite. Is it to prove a theorem? Is it to design a new compression algorithm? Is it to simulate a fluid flow problem? What is the goal?

In my line of work, 3 significant digits is ALL we can ever hope for. With 4, we would be virtually guaranteed contract renewal for the next 5 years. But 100s of digits? C'mon.

Please explain WHY calculating 100 factorial to the final place is of any utility.

BTW, mathematica, maple and other symbolic programs have been able to do this for years.

John D'Errico

John D'Errico (view profile)

In the first release, I had given new names to a few functions, like factorial, nchoosek, and prod. This allowed those functions to work on double inputs, yet still produce a vpi result. Petter makes a valid argument though, so I have changed that.

In this release, functions that already exist in matlab have been overloaded (along with a huge number of other enhancements.)

Petter

Petter (view profile)

Indroducing new function names is inconvenient for the user, overloading existing ones is preferable.

Great work!

Husam Aldahiyat

vpi2english is the icing on the cake.

Updates

1.45

powermod repair

1.44

Flag as a toolbox

1.43

Fix Totient function, some minor doc changes

1.41

Speedup for plus and minus operations, also new min and max functions by request

MATLAB Release
MATLAB 7.5 (R2007b)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

VariablePrecisionIntegers/html/