[ccan] Mathematics Module

[Rusty Russell]

On Friday 06 February 2009 13:56:53 Rupinder Singh wrote:
> Hello,
> 
>   I think your original code was in fact optimal.  For example, it takes
> > about 7.1 seconds to figure out that 36028797018963971 is prime on my
> > 64 bit machine here.
> 
> The original routine (is_prime) is still there. I had expressed in my
> previous mail
> > but is_prime should be used for that due to its superior efficiency.
> 
> It takes less than a second to determine the primality of 36028797018963971
> by is_prime. is_prime2 is the implementation of Sieve of Eratosthenes to
> determine the prime numbers, not very efficient to determine primality, as
> you just said.

Yes, but it's not useful.  I've dropped it.

> Here are a few additions :
> 
> 1) nCr - A Memory and Time efficient routine to calculate n! / (n-r)! r! .

This is OK.

> 2) newton_raphson - Determination of root, Newton-Raphson method. The user
> has to overload the function fn(x) & fnd(x)  [its derivative] to implement
> his given function.
>
> 3) ivp - Solution to the Initial Value Problem dy/dx = fnd2(x,y) , Given
> y(x0) = y0, y(x1) = ?
> 
> 4) romberg_integral - Definite Integration within closed limits [a, b],
> Romberg's Method. fn(x) has to be overloaded.

These need to take the function(s) as arguments, rather than using globals.

Here's what I currently have:

===
#ifndef CCAN_NUMERIC_H
#define CCAN_NUMERIC_H

#include <stdbool.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>

/**
 * is_square - Is the number a perfect square?
 * @n: number to test
 *
 * Example:
 *	if (is_square(4))
 *		printf("%d is indeed a Perfect Square!", 4);
 */
static inline bool is_square(unsigned long n)
{
	unsigned long tmp = sqrt((double)n);

	return tmp*tmp == n;
}

/**
 * is_prime - Is the number prime ?
 * @n: number to test
 *
 * Example:
 *	if (is_prime(17))
 *		printf("%d is Prime!\n", 17);
 */
static inline bool is_prime(unsigned long n)
{
	unsigned long i;
	unsigned int skip3 = 1;

	if(n == 2 || n == 3)
		return true;

	if(n%2 == 0 || n%3 == 0)
		return false;

	for(i=5; i*i <= n; i+=(2 << (skip3 ^= 1)))
		if(n%i == 0)
			return false;

	return true;
}

/**
 * gcd - Returns the GCD of two numbers
 * @a: first number
 * @b: second number
 * Example:
 *	if (gcd(5,7)==1)
 *		printf("gcd of 5 & 7 is 1");
 */
static inline int gcd(long a,long b)
{
    if(b==0) return (a);
    return (gcd(b,a%b));
}

/**
 * nCr - Calculates nCr efficiently, n! / (n-r)!r!
 *
 * Example:
 *	printf("%lu", nCr(6, 2));
 *
 *  15
 */
unsigned long nCr(unsigned long n, unsigned long r);

/**
 * newton_raphson - Returns the root of the function f(x) by the method of 
 * @fn: function to test
 * @fnd: derivative of the function
 * @seed: seed point
 * @eps: x is considered a root if |fn(x)| < eps
 * @delta: 
 * @iter: number of iterations to perform
 * Example:
 *	printf("%lf", newton_raphson(2, 0.001, 0.01, 1000);
 *
 *  3
 */
double newton_raphson(double (*fn)(double x),
		      double (*fnd)(double x),
		      double seed, double eps, double delta,
		      unsigned int iter);

/**
 * ivp - Solves the Initial Value Problem of the type dy/dx = f(x,y), y(x0) = y0, y(x1) = ?
 * @fnd2: dy/dx
 * 
 * In the example : y = e^x^2 ; dy/dx = xy
 *
 * Example:
 *	printf("%lf", ivp(0, 1, 1, 1000);
 *
 *  1.6487
 */
double ivp(double (*fnd2)(double x, double y),
	   double x0, double y0, double x1, unsigned int iter);

/**
 * romberg_integral - Calculates definite integral in the closed interval [a,b]
 *
 * @a : lower limit of the interval
 * @b : upper limit of the interval (up to ROMBERG_MAXITER)
 * Example:
 *	printf("%lf", romberg_integral(1, 2));
 *
 *  163.33333
 */
double romberg_integral(double (*fn)(double x),
			double a, double b,
			unsigned int itermin, unsigned int itermax);


/**
 * ROBERG_MAXITER - the ceiling of iterations for romberg_integral()
 */
#ifndef ROMBERG_MAXITER
#define ROMBERG_MAXITER 100
#endif

#endif /* CCAN_NUMERIC_H */
===
#include "numeric.h"

static inline void div_by_gcd(unsigned long *a, unsigned long *b)
{
	unsigned long g=gcd(*a,*b);
	*a/=g;
	*b/=g;
}

unsigned long nCr(unsigned long n, unsigned long r)
{
	unsigned long numerator=1,denominator=1,toMul,toDiv,i;
	if (r > n/2) r = n-r; // use smaller r
	for (i=r ; i ; i--)
	{
		toMul = n-r+i;
		toDiv = i;
		div_by_gcd(&toMul,&toDiv); /* always divide before multiply */
		div_by_gcd(&numerator,&toDiv);
		div_by_gcd(&toMul,&denominator);
		numerator *= toMul;
		denominator *= toDiv;
	}

	return numerator/denominator;
}

double newton_raphson(double (*fn)(double x),
		      double (*fnd)(double x),
		      double seed, double eps, double delta,
		      unsigned int iter)
{
	unsigned int i;
	double x2 = 0, x1;

	if ( !fnd(seed)) {seed += eps;}

	x1 = seed - fn(seed)/fnd(seed);

	for (i=0;i<iter;i++)
	{
		if (fabs(fnd(x1)) < delta)
			return 0.0/0.0;

		else x2 = x1 - fn(x1)/fnd(x1);
		if (fabs(fn(x2)) < eps) return x2;
		else {seed=x1; x1=x2;}
	}

	return x2;
}

double ivp(double (*fnd2)(double x, double y),
	   double x0, double y0, double x1, unsigned int iter)
{
	unsigned int i;
	double yn, yp=y0;
	double h = (x1-x0)/iter;

	for (i=0;i<iter;i++)
	{
		yn = yp + h*fnd2(x0, yp);
		yp = yn;
		x0 += h;
	}

	return yp;
}

double romberg_integral(double (*fn)(double x),
			double a, double b,
			unsigned int itermin, unsigned int itermax)
{
	int i,j,n=0;
    double d,r,s,ta;
    double I[ROMBERG_MAXITER][ROMBERG_MAXITER];
    if (itermax>ROMBERG_MAXITER) itermax=ROMBERG_MAXITER;
    r = fn(a);
    ta = (r + fn(b))/2;
    d = b-a;
    I[0][0] = ta*d;
	do
	{
		n = n+1;
		d /= 2;
		s=ta;
		for (i=1; i < pow(2.0,n); i++)
			s += fn(a + d*i);

		I[0][n] = s*d;
		r=1;
		for (i=1; i<n+1; i++)
		{
			r=r*4;
			j=n-i;
			I[i][j] = (r*I[i-1][j+1] - I[i-1][j])/(r-1);
		}

		if (n >= itermax) return I[n][0];
	} while (n < itermin);

	return I[n][0];
}