Generating a power-law distribution in C and testing it with python

I know that, given a rng which generates random numbers uniformly distributed, a way to obtain power-like data is, following Wolfram Mathworld the following: let y be a random variable uniformly distributed in (0,1) and x another random variable distributed as P(x) = C*x**n (for x in (xmin,xmax)). We have that

x=[ (xmax**(n+1) - xmin**(n-1))y+xmin**(n+1)  ]**(1/(n+1))

So i made this program in C that generates 50k numbers from 1 to 100 that should be distributed as x^(-2) and prints the frequency of the outcomes on a file DATA.txt:

void random_powerlike(int *k, int dim,  double degree, int xmin, int xmax, unsigned int *seed)
{
int i; 
double aux;
for(i=0; i<dim; i++)
    {
    aux=(powq(xmax, degree +1 ) - powq(xmin, degree +1 ))*((double)rand_r(seed)/RAND_MAX)+ powq(xmin, degree +1);

    k[i]=(int) powq(aux, 1/(degree+1));

    }
}

int main()
{
    unsigned int seed = 1934123471792583;

    FILE *tmp; 
    char  stringa[50];
    sprintf(stringa, "Data.txt");
    tmp=fopen(stringa, "w");

    int dim=50000;
    int *k;
    k= (int *) malloc(dim*sizeof(int));
    int degree=-2;
    int freq[100];  

    random_powerlike(k,dim, degree, 1,100,&seed);
    fprintf(tmp, "#degree = %d  x=[%d,%d]\n",degree,1,100);
    for(int j=0; j< 100;j++)
    {   
        freq[j]=0;
        for(int i = 0; i< dim; ++i)
        {
            if(k[i]==j+1)
            freq[j]++;
        }
        fprintf(tmp, "%d    %d\n", j+1, freq[j]);
    }
    fflush(tmp);
    fclose(tmp);

return 0;
}

I decided to fit these numbers with pylab, to see if the best power-law to fit them is something as a*x**b, with b = -2. I wrote this program in python:

import numpy
from scipy.optimize import curve_fit
import pylab

num, freq = pylab.loadtxt("Data.txt", unpack=True)
freq=freq/freq[0]

def funzione(num, a,b):
    return a*num**(b)

pars, covm =  curve_fit(funzione, num, freq, absolute_sigma=True)
xx=numpy.linspace(1, 99)
pylab.plot(xx, funzione(xx, pars[0],pars[1]), color='red')
pylab.errorbar(num, freq, linestyle='', marker='.',color='black')
pylab.show()
print pars

The problem is that when i fit the data, I obtain an exponent value of ~-1.65.

I think that I made a mistake somewhere, but I can't figure it out where.