User:Starrynte/Remez

A utility for finding polynomial approximations for a function of your choice.

Current feature(s):

• Any function of your choice
• Any degree polynomial of your choice
• Any interval of your choice (i.e. approximates function over [-1,1] )

Planned feature(s):

• Rational approximations

How to Use

The code is below. Simply change the top part labeled "your settings", and the program does the rest. (Note: Attempting to do "crazy" approximations (i.e. high degree, low interval) will result in erratic behavior, presumably due to double precision)

public class Remez {
	//---YOUR SETTINGS---
	static double f(double x){ //THE FUNCTION TO BE APPROXIMATED
		return Math.sin(x);
	}
	static final int DEG = 6; //DEGREE OF THE POLYNOMIAL - MUST BE POSITIVE
	static final double lowx = 0d, highx = Math.PI*2; //APPROXIMATES F(X) FROM LOWX TO HIGHX 
	
	
	//---MAIN PROGRAM CODE---
	
	
	static final double RANGE = highx - lowx;
	static final int COLUMNS = DEG + 3;
	static final int ROWS = DEG + 2;
	static final double[] extrema = new double[ROWS];
	static final double[][] matrix = new double[ROWS][COLUMNS];
	
	public static void main(String[] args) {
		for(int i=0;i<ROWS;i++) extrema[i]=0.5 * (-Math.cos(i*Math.PI/(DEG + 1)) + 1) * RANGE + lowx;		
		for(int i=0;i<20;i++) remez();
		for(int i=0;i<DEG;i++) System.out.print(matrix[i][COLUMNS-1] + " + x * (");		
		System.out.print(matrix[DEG][COLUMNS-1]);
		for(int i=0;i<DEG;i++) System.out.print(")");
		System.out.println("\nMax error: " + matrix[ROWS-1][COLUMNS-1]);
	}	
	
	private static void remez(){
		for(int i=0;i<ROWS;i++){			
			for(int j=0;j<DEG+1;j++) matrix[i][j]=Math.pow(extrema[i],j);			
			matrix[i][DEG+1]=(i&1)==0?1d:-1d;
			matrix[i][DEG+2]=f(extrema[i]);
		}
		rref(matrix);
		extrema[0]=lowx;
		int z=1;
		for(double x=lowx+(RANGE/1000);x<=highx && z<(DEG+2);x+=(RANGE/1000)){
			if(deverf(x) * deverf(x-(RANGE/1000)) < 0){
				extrema[z]=extrema(x-(RANGE/1000),x);
				z++;
			}
		}
	}
	
	private static void rref(double[][] m){
		int lead; double factor;
		for(int r=1;r<ROWS;r++){ //starting at second row
			lead=0;
			for(int i=0;i<r;i++,lead++){
				factor=m[r][lead]/m[i][lead];
				for(int c=lead;c<COLUMNS;c++){				
					m[r][c]-=factor*m[i][c];
				}
			}
			for(int c=COLUMNS-1;c>=lead;c--){
				m[r][c]/=m[r][lead];
			}
		}
		for(int r=ROWS-2;r>=0;r--){ //starting at second last row
			lead=COLUMNS-2;
			for(int i=ROWS-1;i>r;i--,lead--){
				factor=m[r][lead];
				for(int c=lead;c<COLUMNS;c++){
					m[r][c]-=factor*m[i][c];
				}
			}
		}
	}
	
	private static double deverf(double x){
		final double h = 0.001 * (1 + Math.abs(x));
		double dev = (f(x+h) - f(x-h))/(2*h);			
		for(int r=1;r<(ROWS-1);r++){ //constant term doesn't affect derivative
			dev -= r * matrix[r][COLUMNS-1] * Math.pow(x,r-1);
		}
		return dev;
	}
	
	private static double extrema(double lbound,double rbound){
		double mid;
		while(Math.abs(rbound-lbound)>0.0000000002){
			if(deverf(lbound)*deverf(mid=(rbound+lbound)*0.5) > 0) lbound=mid;
			else rbound=mid;
		}
		return (lbound+rbound)*0.5;
	}
	
}