Code covered by the BSD License  

Highlights from
Variable Precision Integer Arithmetic

4.85

4.8 | 44 ratings Rate this file 104 Downloads (last 30 days) File Size: 2.86 MB File ID: #22725

Variable Precision Integer Arithmetic

by

 

19 Jan 2009 (Updated )

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

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

| Watch this File

File Information
Description

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!

Acknowledgements

Multiple Precision Toolbox For Matlab inspired this file.

This file inspired The Computation Of Pi By Archimedes, Multiple Root Polynomial Solved By Partial Fraction Expansion, Nextprime, Sumsqint, Num2vpi Converts Input Exactly Into Vpi, System Of Linear Congruences, Fractions Toolbox, and Hpf A Big Decimal Class.

Required Products MATLAB
MATLAB release MATLAB 7.5 (R2007b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (141)
05 Sep 2014 John D'Errico

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.

05 Sep 2014 Seth Wagenman

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

29 Aug 2014 John D'Errico

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.

29 Aug 2014 nicolas rosselot

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!

24 May 2014 Edgar Zakarian  
22 May 2014 Edgar Zakarian

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.

21 May 2014 John D'Errico

N = vpi('123456789012345678901234567890');

vpi2bin(N)
ans =
1100011101110100100001111111101101100001101110011111000001110111001001110001111110000101011010010

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

21 May 2014 Edgar Zakarian

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

15 May 2014 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

29 Apr 2014 John D'Errico

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.)

28 Apr 2014 Suma

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?

15 Feb 2014 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

04 Feb 2014 Tétef  
15 Jan 2014 John D'Errico

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.

14 Jan 2014 Edgar Zakarian

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.

02 Dec 2013 Cavin Dsouza

Remarkable Sir

21 Oct 2013 Jean-Luc

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

07 Oct 2013 Shiuan-Tzuo Shen

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

19 Aug 2013 John D'Errico

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.

19 Aug 2013 Fogato Abbestia

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

31 Jul 2013 John D'Errico

Oops. Fixed totient(1) & uploaded.

31 Jul 2013 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.

30 Apr 2013 Thomas Reasoner

I used it for prime number factorization. Works great!

16 Dec 2012 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

16 Dec 2012 John D'Errico

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.

16 Dec 2012 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.

13 Oct 2012 kevin  
20 Jun 2012 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.

19 Mar 2012 kevin

thank you John!

19 Mar 2012 John D'Errico

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.

19 Mar 2012 kevin

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 .

19 Mar 2012 kevin

% 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

18 Mar 2012 John D'Errico

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.

18 Mar 2012 kevin

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.

03 Mar 2012 Michaela

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

27 Feb 2012 Tristan

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

02 Feb 2012 John D'Errico

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.)

02 Feb 2012 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

02 Feb 2012 John D'Errico

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.

02 Feb 2012 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.

01 Feb 2012 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.

31 Jan 2012 John D'Errico

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.

31 Jan 2012 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?

31 Jan 2012 Michael

Oh my God I love you.

02 Dec 2011 John D'Errico

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.

02 Dec 2011 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

16 Nov 2011 Nanthini

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

16 Nov 2011 John D'Errico

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

16 Nov 2011 John D'Errico

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.

16 Nov 2011 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...

08 Nov 2011 John D'Errico

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.

08 Nov 2011 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

11 Oct 2011 Ahmet Karakucuk

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

04 Oct 2011 John D'Errico

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

04 Oct 2011 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....

03 Oct 2011 John D'Errico

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.

03 Oct 2011 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....

29 Sep 2011 John D'Errico

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.

29 Sep 2011 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

29 Sep 2011 Jean-Luc  
27 Sep 2011 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....

24 Sep 2011 John D'Errico

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.

24 Sep 2011 Nanthini

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

23 Sep 2011 John D'Errico

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.

23 Sep 2011 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..

22 Sep 2011 Nanthini

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

22 Sep 2011 Nanthini

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

22 Sep 2011 John D'Errico

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.

22 Sep 2011 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....

03 Aug 2011 John D'Errico

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

03 Aug 2011 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

12 Jun 2011 Sarla

:) 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!!

12 Jun 2011 John D'Errico

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.

12 Jun 2011 Sarla

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.

11 Jun 2011 John D'Errico

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.

11 Jun 2011 John D'Errico

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.

11 Jun 2011 Sarla

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?

10 Jun 2011 Sarla

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?

06 Jun 2011 John D'Errico

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.

06 Jun 2011 John D'Errico

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.

06 Jun 2011 Gabriel

@Sarla

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

06 Jun 2011 Sarla

:) 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!

06 Jun 2011 Gabriel  
06 Jun 2011 Gabriel

-_-

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

06 Jun 2011 John D'Errico

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.

06 Jun 2011 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

05 Jun 2011 John D'Errico

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.

05 Jun 2011 Sarla

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

04 Jun 2011 John D'Errico

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.

04 Jun 2011 Sarla

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?

20 May 2011 John D'Errico

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

20 May 2011 Zakir

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?

15 May 2011 Zakir

Very Very Nice Job
Thanks for This.

17 Feb 2011 Sarla

okay. thank you anyway!

16 Feb 2011 John D'Errico

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.

16 Feb 2011 Sarla

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.

11 Feb 2011 Sarla

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!

24 Jan 2011 Evgeny Pr  
04 Jul 2010 John D'Errico

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.

04 Jul 2010 Gustaf Lindberg

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

14 Apr 2010 Joseph Kirk

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

07 Feb 2010 Harriet Fell  
27 Jan 2010 Rob Campbell  
26 Jan 2010 John D'Errico

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

26 Jan 2010 Ian

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.

20 Jan 2010 Nathan Greco  
12 Jan 2010 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 !

31 Dec 2009 John D'Errico

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

31 Dec 2009 GennaroAlphonse Ramón Rodríguez

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.

22 Nov 2009 John D'Errico

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.

22 Nov 2009 Derek O'Connor

John,

The demo_vpi.m needs your consolidator function:

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

Great package.

Derek O'Connor

06 Aug 2009 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.

29 Jul 2009 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

22 Jul 2009 John D'Errico

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

22 Jul 2009 Gwendolyn Fischer

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

14 Jul 2009 John D'Errico

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

14 Jul 2009 John D'Errico

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.

14 Jul 2009 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

14 Jul 2009 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.

19 May 2009 Neha P

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.

23 Apr 2009 Lissa  
15 Apr 2009 John D'Errico

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

15 Apr 2009 Andreas Bonelli  
15 Apr 2009 Andreas Bonelli

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);

02 Apr 2009 John D'Errico

Sorry. I'll fix the demo.

01 Apr 2009 Giampiero Campa

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

24 Mar 2009 John D'Errico

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.

24 Mar 2009 John D'Errico

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.

24 Mar 2009 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.

23 Mar 2009 Derek O'Connor

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

23 Mar 2009 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.

14 Mar 2009 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.

13 Mar 2009 Kenneth Eaton  
13 Mar 2009 Orlando Rodríguez

What a marvelous tool! Thank you John!

13 Mar 2009 Derek O'Connor

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

12 Mar 2009 John D'Errico

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

11 Mar 2009 Kevin

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.

11 Mar 2009 John D'Errico

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.)

21 Feb 2009 Petter

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

15 Feb 2009 Sebastian Randel

Great work!

11 Feb 2009 Husam Aldahiyat

vpi2english is the icing on the cake.

Updates
25 Jan 2012

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

31 Jul 2013

Fix Totient function, some minor doc changes

Contact us