Arithmetic with integers of fully arbitrary size. Arrays and vectors of vpi numbers are supported.
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!
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 
Inspired by: Multiple Precision Toolbox for MATLAB
Inspired: Multipleroot polynomial solved by partial fraction expansion, nextprime, sumsqint, num2vpi  Converts input exactly into vpi, HPF  a big decimal class, The Computation of Pi by Archimedes, Fractions Toolbox, Number to Myriad, Words to Number, Number to Words, GEM Library
Lukas (view profile)
Nice! I needed to create iszero function for nonvpi variables to get going. But I assume this is due to a missing Symbolic toolbox in my case.
function [b] = iszero(n)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
b = boolean(n == 0);
end
Otherwise I was getting error:
>> b = vpi('5')
>> b * 2
Undefined function 'iszero' for input arguments of type 'double'.
Error in .* (line 49)
if iszero(INT1)  iszero(INT2)
Error in * (line 27)
INT = INT1.*INT2;
Paulo Fonte (view profile)
The complement
function hex=vpi2hex(INT,n)
HEX=vpi2base(INT,16);
if nargin==2
HEX=[zeros(1,nnumel(HEX)) HEX];
end
X='0123456789ABCDEF';
hex=X(HEX+1);
Paulo Fonte (view profile)
A little helper to convert large hex strings to vpi. (Or maybe it's already there and I couldn't find how to use?)
function INT=hex2vpi(hex)
X=NaN(1,55);
X(1:10)=0:9;
X(18:23)=10:15;
X(50:55)=10:15;
I=hex'0'+1;
INT=base2vpi(X(I),16);
pulengmoth (view profile)
Mohsin Shah (view profile)
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 (view profile)
Paul (view profile)
Jan Sieber (view profile)
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 (view profile)
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 (view profile)
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 (view profile)
Great work John. Keep it up. Is there any way I can convert a vpi number to a character string?
aiman saad (view profile)
you are amazing , thank you very much, how did you do that????
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 (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 (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 (view profile)
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.
selcuk caglar (view profile)
Can I use it in simulink?
Afefa Asiri (view profile)
Yu Ang Tan (view profile)
Well done! Thanks for sharing this awesome package!
Guangxi Liu (view profile)
Ander Biguri (view profile)
Poorren (view profile)
good work！
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 (view profile)
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?
Edgar Zakarian (view profile)
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 192bit size. thanks.
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.
Seth Wagenman (view profile)
Is this free to download? R2014a which I am using does not have vpi apparently...
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 rewrite 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.
nicolas rosselot (view profile)
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!
Edgar Zakarian (view profile)
Edgar Zakarian (view profile)
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 (view profile)
N = vpi('123456789012345678901234567890');
vpi2bin(N)
ans =
1100011101110100100001111111101101100001101110011111000001110111001001110001111110000101011010010
I suppose in hindsight, I should have used dec2bin as the name for this.
Edgar Zakarian (view profile)
Hi,
How can I convert a vpi number to binary ?
Thanks.
Giuseppe Cardillo (view profile)
Dear John,
I saw that you use the MillerRabin 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 MillerRabin 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 (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 (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?
JeanLuc (view profile)
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;
JeanLuc
Tétef (view profile)
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.
Edgar Zakarian (view profile)
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 (view profile)
Remarkable Sir
JeanLuc (view profile)
Hey John: thanks for implementing min and max in the latest version!
ShiuanTzuo Shen (view profile)
Nice submission!
Suggest the implementation of generating random numbers in Z_p^*, where p is a vpi number. Thanks!
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.
Fogato Abbestia (view profile)
how do i install it on windows ? i have matlab 2008b 64bit, thx.
John D'Errico (view profile)
Oops. Fixed totient(1) & uploaded.
the cyclist (view profile)
Great submission. I use it rarely, but it is a lifesaver 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.
Thomas Reasoner (view profile)
I used it for prime number factorization. Works great!
Seth Kosowsky (view profile)
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 ndigit numbers in basen 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 (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 (view profile)
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 (view profile)
Juliano (view profile)
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 (view profile)
thank you John!
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 (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 (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:ii1
if counts>0
break
end
temps=(ii*jj)*(ii+jj)*(iijj)/11; % here the m.
if round(temps)==temps
if~isempty(find(test256==mod(temps,256)))
if mod(temps,3)<2 &&(mod(temps,5)<2mod(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=xxxround(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 (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 (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, √n97 are all rational numbers ? I believe this would be a good verification method for the toolbox's calculating capability.
Michaela (view profile)
Wow ... I design floating point hardware, and this is just what the doctor ordered. Big Thanks xxx!
Tristan (view profile)
Very rad ! Worked right out of the box for what i needed, thanks!
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,1e10001]. 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 rewrite.)
Michael (view profile)
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 100200 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 perfectlyif 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 (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 rewriting VPI, but even then I will not reasonably expect more than about 101 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 (view profile)
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 (view profile)
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 resampling approach.
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:
>> 110^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^531, you could do a similar resampling scheme on vpi numbers using the randint function to generate the samples. unique is defined for vpi.
Michael (view profile)
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 (view profile)
Oh my God I love you.
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 (view profile)
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 (view profile)
Thanks john...ya you are right here after i will go through the helps as much as i can...
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 (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 nonpositive numbers.
Nanthini (view profile)
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 (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 (view profile)
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
Ahmet Karakucuk (view profile)
This is something truly great. Thank you for all your effort you've put into this great "product".
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 email 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 (view profile)
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 (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 (view profile)
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 (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.
JeanLuc (view profile)
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,
JeanLuc
JeanLuc (view profile)
Nanthini (view profile)
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 devc++ compiler...while running that function it asking to mention the compiler name used. i gave it as devc++ compiler. but for that it shows the following error message:
mex('C:\Users\Nanthini\Documents\MATLAB\num2strexact.c')
Select a compiler:
[0] None
Compiler: DevC++ IDE
Please select from 00
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^321 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/15C27B9/?solution=15C27B9
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 nonC 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 (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/22239num2strexactexactversionofnum2str
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 (view profile)
2960000000000000100663296, can you tell me the number of bits occupied for the above mentioned number...
Thanks in advance...
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^531.
>> 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 (view profile)
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 (view profile)
thank you so much.... i got the output...:)
Nanthini (view profile)
sorry i am using MATLAB R2011a mistakenly
i typed as R2021a...
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 mfiles 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 nono, 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 reinstall, this time following the instructions provided.
Nanthini (view profile)
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 (view profile)
M = vpi(2)^100;
M =
1267650600228229401496703205376
Modular operations work as you would expect.
mod((M+1)^2631,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 (view profile)
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 (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 (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^531 as a double.
>> 2^531
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 pseudorandomness 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 (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 (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 (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^531. 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 (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 (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 (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 (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 (view profile)
@Sarla
Use isprime(vpi(p))?
It's probably faster than factoring the whole thing. It return 1 when the number is prime.
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 (view profile)
Gabriel (view profile)
_
Lesson learned... Thank you =)
Is there a limit to the exponent in powermod?
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 (view profile)
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 (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 (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 (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 (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 (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 (view profile)
I want to factorize very large no
Example:
a = vpi(2);
b = a^2067;
c = b7;
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 (view profile)
Very Very Nice Job
Thanks for This.
Sarla (view profile)
okay. thank you anyway!
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 (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 (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 (view profile)
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.
Gustaf Lindberg (view profile)
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 (view profile)
This is quite possibly the best submission I have downloaded from the File Exchange. Great work John!
Harriet Fell (view profile)
Rob Campbell (view profile)
John D'Errico (view profile)
Thanks Ian. I've sent in an update to fix the legendresymbol bug. John
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 (view profile)
Guillaume (view profile)
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 (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 supercomputer in the world, even on a cloud of the 1000 fastest supercomputers 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
GennaroAlphonse Ramón Rodríguez (view profile)
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 (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.
Derek O'Connor (view profile)
John,
The demo_vpi.m needs your consolidator function:
http://www.mathworks.com/matlabcentral/fileexchange/8354consolidator
Great package.
Derek O'Connor
Khaled Hamed (view profile)
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 (view profile)
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^2y^6121*y^42)+550*y^8
z = 1.5112e+023
>> x=sym(77617);y=sym(33096);z=33375*y^6+100*x^2*(11*x^2*y^2y^6121*y^42)+550*y^8
z = 200
>>x=vpa(77617);y=vpa(33096);z=100*x^2*(11*x^2*y^2y^6121*y^42)+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^2y^6121*y^42)+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 (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
Gwendolyn Fischer (view profile)
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 (view profile)
Ben  I've submitted an update, that now handles the pathological cases for plus and times.
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 (view profile)
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 (view profile)
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 (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 (view profile)
John D'Errico (view profile)
I just uploaded a new release, fixing the bug in vpi2english.
Andreas Bonelli (view profile)
Andreas Bonelli (view profile)
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 (view profile)
Sorry. I'll fix the demo.
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 (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 (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 (view profile)
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.
Derek O'Connor (view profile)
Dear Alain,
Thanks for your reply and correct answer mean(x) = 4.661273948537807e005, (or
sum(x) = 7680*4.661273948537807e005 = 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 1norm 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 illconditioned.
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.579858392477036e001
Here is the result using Rump's Intlab 5.5 set to 30digit 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 26digit
precision to obtain the exact sum.
Intlab is here : http://www.ti3.tuharburg.de/rump/intlab/
Regards,
Derek O'Connor
Alain Barraud (view profile)
vpi is an excellent toolbox, thanks a lot.
I shoud add a comment to derek O Coonor question. I have found 4.661273948537807e005 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 (view profile)
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 (view profile)
Orlando Rodríguez (view profile)
What a marvelous tool! Thank you John!
Derek O'Connor (view profile)
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 (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 (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 (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 (view profile)
Indroducing new function names is inconvenient for the user, overloading existing ones is preferable.
Sebastian Randel (view profile)
Great work!
Husam Aldahiyat (view profile)
vpi2english is the icing on the cake.