A Simple Windows Application for Factorizing Integers

So, I spend so much time working and doing mathematics and walking the dogs that I don’t always spend much time on programming. But, over the last year, I have realized that everything that I need to make windows applications is already on my computer. I don’t need any fancy software: just a text editing program, and an open source compiler to compile from the command line.
So my first code that I am posting here is a C++ application for factorizing integers. It just runs in a window with no fancy GUI. But it factorizes numbers very nicely.

The link for the downloadable exe file is below the code.

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstdint>
#include <string>
//  #include <NTL/ZZ.h>  get this number theory stuff eventually

bool primep (int64_t N);

int main()
{
using namespace std;
std::string another ="N";
do
{
int64_t N;
int64_t divisorlist[1000];
int64_t primedivisorlist[100][3];
int64_t sumofdivisors =0;
int64_t phi=1;
int counter (1);
bool squarep (false);
cout<< "This application factors integers of up to seventeen digits."<<endl;
cout << "Enter an integer to be factored: " <<endl;
cin>>N;
cin.ignore(32767, '\n');

for (int m=1; m<= floor(sqrt(N));m++ )
{
if (N%m == 0)
{divisorlist[counter]=m;
counter++;}
}

for(int m=counter-1;m> 0;m--)
{
int64_t quotient = N/divisorlist[m];
if (divisorlist[m] != quotient)
{divisorlist[counter]=quotient;
counter++;
}}
// now extract primes
int bar=1;
for(int m =2; m<=counter-1;m++)
{
if (primep(divisorlist[m]))
{primedivisorlist[bar][1]=divisorlist[m];
bar++;}
}
// now find powers of primes

for (int m=1; m<=bar-1; m++)
{
int power=1;
int64_t x = round(pow(primedivisorlist[m][1], power));

while (N % x ==0
&& x<=N)
{
power++;
x = round(pow(primedivisorlist[m][1], power));
}
primedivisorlist[m][2]=power-1;
}
cout<<N<<" = ";
for(int m =1; m<=bar-1;m++)

{
cout<<primedivisorlist[m][1];
if (primedivisorlist[m][2]>1 )
cout<< "^"<< primedivisorlist[m][2];
if (m<bar-1) cout<<"*";
phi=phi*round(pow(primedivisorlist[m][1], (primedivisorlist[m][2]-1)))*(primedivisorlist[m][1]-1);
}
cout<<endl;

cout<< "Number of divisors: "<<counter-1<<endl;
cout<<"The divisors are: "<<endl;
for(int m =1; m<=counter-1;m++)
{
cout<<divisorlist[m]<<" ";
sumofdivisors = sumofdivisors+divisorlist[m];
}
cout<<endl;
cout<<"The sum of the aliquot divisors of "<<N<<" is "<<sumofdivisors-N<<". Hence, this number is ";
if (sumofdivisors >(2*N)) {cout<<"abundant. "<<endl;}
else if (sumofdivisors<(2*N)){cout<<"deficient. "<<endl;}
else {cout<<"perfect. "<<endl;}

cout<<"Phi, the number of integers less than and relatively prime to "<<N<<" is "<<phi<<endl;
cout<<N<<" is ";
if (primep (N))
cout<<" ";
else cout<<"not ";
cout<<"prime."<<endl;
squarep = (counter-1)%2;
cout<<N<<" is ";
if (squarep)
cout<<" ";
else cout<<"not ";
cout<<"a perfect square."<<endl;
cout<<"Another number? ";

std::getline(std::cin, another);
}
while ((another == "yes") || (another == "y"));

return 0;
}

bool primep (int64_t N)

//returns t if N is prime
{
for (int m=2; m<= floor(sqrt (N)); m++)
{
if (N%m ==0) return false;
}
return 1;
}

The code is very simple and self-explanatory.  It is not the most efficient and beautiful way to factorize a number, but it works well enough within the realm of integers that can be easily handled by the base C++ variable types.

Here is the link for the exe file:

Factorizer

Multi-Period Binomial Options

Here is my latest program for calculating the prices of options using a multi-period binomial tree. It works, but represents a bunch of strange methods to solve unexpected problems.
First, I decided to calculate the price at each node by finding the price at the terminal nodes, then use a formula, rather than calculate the prices recursively. This technique, as you will see, led to the problem of calculating binomial coefficients.

The thing that I am most satisfied with is that I was able to get the program to output to a text file, which can be opened in Excel to visualize the data.  The only thing that the program does not do is to do the same calculation with the Black-Scholes formula, and visualize how the binomial calculation approaches the Black-Scholes result.

First, take a look at my output:

Binomial Pricing

We can see how the results approach a limit with each iteration.  The pricing for periods 34 and 35 is starting to get weird because of limits of integer sizes and floating point number sizes.

Here is a chart of the above data:

Binomial Option Tree

Finally, here is the code:

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <iomanip>
#include <fstream>

//Calculates price of European call or put option, and portfolio required to duplicate option, using binomial tree

int nchoosek(int n,int k);
//choice holds the binomial tree so that it doesn’t need to be recalculated
unsigned long long choice [50][50]={};

int main() {
using namespace std;
float S0, r,div, K, h,t,sigma, u, d,pStar, uS,dS,Cu,Cd,delta, bondValue;
double optionPrice;
double prices [50];
int n;
string callPut;
int putTrue;
//fill in nchoosek tree
for (int m=1; m<50;m++)
{ choice [m][0]=1;
choice [m][m]=1;}
for (int m=2;m<(50);m++)
{ for (int l=1;l<m;l++){
choice[m][l]=choice[m-1][l-1]+choice[m-1][l];
}
}
//get data from user
cout<<“Calculate the price of a European call or put, using a one period binary tree, then find rate with n period tree.”<<endl;
cout<<“Current Stock Price?”;
cin>>S0;
cout<<“Continuously Compounded Risk Free Rate?”;
cin>>r;
cout<<“Continuous Dividend Rate?”;
cin>>div;
cout<<“Volatility?”;
cin>>sigma;
cout<<“Strike Price?”;
cin>>K;
cout<<“Time to Expiration?”;
cin>>t;
cout<<“Number of Periods? (less than 50 please)”;
cin>>n;
cout<<“Call or Put?”;
cin>>callPut;
if ((callPut == “put”)||(callPut==”p”))
putTrue=-1;
else putTrue=1;
//calculate likely upward and downward motions
u=exp((r-div)*t+sigma*sqrt(t));
d=exp((r-div)*t-sigma*sqrt(t));
uS=u*S0;
dS=d*S0;
cout<<“Upward Move “<<uS<<endl;
cout<<“Downward Move “<<dS<<endl;
//calculate payoff of call, given motions above
//the variable putTrue turns Us-K into K-Us.  Cool.
Cu=fmax(0, putTrue*(uS-K));
Cd=fmax(0, putTrue*(dS-K));
//calculate Delta, Beta, and call option price
delta=exp(-div * t)*(Cu-Cd)/(S0*(u-d));
bondValue=exp(-r * t)*((u*Cd)-(d*Cu))/(u-d);
optionPrice=delta*S0+bondValue;
cout<<“To duplicate the option, buy “<<delta<<” shares=”” of=”” the=”” stock=”” for=”” “<<delta*s0<<“=”” dollars=”” and=”” borrow=”” “<<-1=”” *=”” bondvalue<<“=”” at=”” rate=”” “<<r<<“=”” a=”” total=”” price=”” “<<optionprice<<<span=”” class=”hiddenSpellError” pre=”” data-mce-bogus=”1″>endl;
cout<<“Price of Option “<<optionPrice<<endl;
cout<<“Prices above 30 levels deep may accumulate errors wildly”<<endl;

//Multi Period tree
int j=1;
while(j<(n+1)){
h=t/j;
u=exp((r-div)*h+sigma*sqrt(h));
d=exp((r-div)*h-sigma*sqrt(h));
pStar=(exp((r-div)*h)-d)/(u-d);
optionPrice=0;

for (int m=0; m<(j+1); m++) {
//This is the big doozy.
//Instead of computing the price at each interior node, only the terminal nodes are calculated
//and the price is found using the binomial theorem
optionPrice=optionPrice+nchoosek(j,m)*pow(pStar,j-m)*pow((1-pStar),m)*fmax(0,(putTrue* ((S0*pow(u, j-m)*pow(d,m))-K)));
}
optionPrice=optionPrice*exp(-r*t);
prices [j]=optionPrice;
cout<<“Periods: “<<j<<”  “;
cout<<“Option Price:  “<<optionPrice<<endl;
++j;  }

//output prices to text file
//open the text file in excel
//data is delimited by commas
ofstream outf(“Results.txt”);
if (!outf)
{  cerr << “File Error”<<endl;
exit(1); }
outf<<“Depth of Tree”<<“,”<<“Price”<<endl;
j=1;
while(j<(n+1)){
outf << j <<“,”<< prices[j]<< endl;
++j;  }
return 0;}

int nchoosek(int n, int k)
{ if (n>49 ||k>49)return 0;
else { return choice[n][k]; }}

The heart of everything is this ugly beast:

optionPrice=optionPrice+nchoosek(j,m)*pow(pStar,j-m)*pow((1-pStar),m)*fmax(0,(putTrue* ((S0*pow(u, j-m)*pow(d,m))-K)));

This is better put as C=C+\binom j m p^{*^{j-m}}(1-p*)^m*max(0,put(S_ou^{j-m}d^m-K)

You can figure out why this works fairly easily.  I like the bit where I multiply (F_{0,T}-K) by -1 to get (K-F_{0,T}.

More than anything, this exercise shows how errors creep into these long calculations.  Looking at the graph, you can see how our numbers are already wandering away from the correct numbers after the 20th iteration.

A Simple C++ Program for Computing Binomial Option Prices

I have been advised that C++ and VBA are commonly requested languages for actuaries. So, while I am studying for exam MFE, I have been writing simple programs.  My native language is LISP.  I like its beauty, and its adaptability.  I am sure that C++ has its promoters, but I think that it is a shame that such an ugly programming language has become one of the most used.

Of course, I am still a novice to the language, so my code is still extra ugly.

Computes call price using binomial model:

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <iomanip>

//Calculates price of call option, and portfolio required to duplicate option, using binomial tree

int main() {
using namespace std;
    float S0, r,div, K, h, sigma, u, d, uS,dS,Cu,Cd,delta, bondValue, callPrice;
//get data
    cout<<“Calculate the price of a call, using a one period binary tree”<<endl;
    cout<<“Current Stock Price?”;
    cin>>S0;
    cout<<“Continuously Compounded Risk Free Rate?”;
    cin>>r;
    cout<<“Continuous Dividend Rate?”;
    cin>>div;
    cout<<“Strike Price?”;
    cin>>K;
    cout<<“Time to Expiration?”;
    cin>>h;
    cout<<“Volatility?”;
    cin>>sigma;
//calculate likely upward and downward motions
    u=exp((r-div)*h+sigma*sqrt(h));
    d=exp((r-div)*h-sigma*sqrt(h));
    uS=u*S0;
    dS=d*S0;
    cout<<“Upward Move “<<uS<<endl;
    cout<<“Downward Move “<<dS<<endl;
//calculate payoff of call, given motions above
    if (uS > K)
    Cu = uS-K;
    else Cu = 0;
    if (dS > K)
    Cd = dS-K;
    else Cd = 0;
//calculate Delta, Beta, and call option price
    delta=exp(-div * h)*(Cu-Cd)/(S0*(u-d));
    bondValue=exp(-r * h)*((u*Cd)-(d*Cu))/(u-d);
    callPrice=delta*S0+bondValue;
    cout<<“To duplicate the call, buy “<<delta<<” shares=”” of=”” the=”” stock=”” for=”” “<<delta*s0<<“=”” dollars=”” and=”” borrow=”” “<<-1=”” *=”” bondvalue<<“=”” at=”” rate=”” “<<r<<“=”” a=”” total=”” price=”” “<<callprice<<<span=”” class=”hiddenSpellError” pre=”” data-mce-bogus=”1″>endl;
    cout<<“Price of Call “<<callPrice<<endl;

//find the price of some other related calls
    cout<<“Some other Call prices: “<<endl;
    float step;
    step = 0.05*K;
    K= K-5*step;
       for (int z=0;z<11;z++){
       if (uS > K)
    Cu = uS-K;
    else Cu = 0;
    if (dS > K)
    Cd = dS-K;
    else Cd = 0;

    delta=exp(-div * h)*(Cu-Cd)/(S0*(u-d));
    bondValue=exp(-r * h)*((u*Cd)-(d*Cu))/(u-d);
    callPrice=delta*S0+bondValue;

    cout<<setw(8)<<k<<” “<<callprice<<<span=”” class=”hiddenSpellError” pre=”” data-mce-bogus=”1″>endl;
    K=K+step;
}
    return 0;
}
//Wow what an ugly program

Here is a sample run:

Binomial option pricing