#import "CurveModel.h"

#import <complex.h>

@implementation CurveModel

//  y^2 = x^3+a*x+b
//
//

// return positive y for a given x
- (double)eval:(double)x {
    return sqrt(x*x*x+a*x+b);
}

double complex ccubert(double complex x)
{
    return cpow(x, 1.0/3.0);
}
int cmpdbl(const void *pa, const void *pb)
{
    double a= *(double*)pa;
    double b= *(double*)pb;
    if (a<b) return -1;
    if (a>b) return  1;
    return 0;
}

- (int)findroots:(double*)r {
    double complex disc= ((a==0 && b>0) ? -1 : 1 )*csqrt(4.0*a*a*a + 27.0*b*b);

    double complex c3d9b= ccubert(sqrt(3)*disc-9.0*b);

    double complex c1= c3d9b/ccubert(18)-ccubert(2.0/3.0)*a/c3d9b;

    double complex onemis3= (1+I*sqrt(3));
    double complex onepis3= (1-I*sqrt(3));
    double complex c2= onepis3*a/ccubert(12)/c3d9b - onemis3*c3d9b/2/ccubert(18);
    double complex c3= onemis3*a/ccubert(12)/c3d9b - onepis3*c3d9b/2/ccubert(18);

    double EPSILON= 0.000001;
    int n=0;
    if (fabs(cimag(c1)) < EPSILON) { r[n] = creal(c1); n++; }
    if (fabs(cimag(c2)) < EPSILON) { r[n] = creal(c2); n++; }
    if (fabs(cimag(c3)) < EPSILON) { r[n] = creal(c3); n++; }

    qsort(r, n, sizeof(*r), cmpdbl);


    int j=0;
    for (int i=1 ; i<n ; i++)
    {
        if ((r[i]-r[i-1])>EPSILON)
            r[j++]= r[i-1];
    }
    r[j++]= r[n-1];
    n= j;

    return n;
}

-(void)addsection:(NSMutableArray*)part view:(CurveView*)view xmin:(double)xmin xmax:(double)xmax
{
    for (double x = xmax ; x>xmin ; x-=view->_xdelta) {
        double y= [self eval:x];
        if (y > view->_ymin && y < view->_ymax) {
            [part addObject:[NSValue valueWithCGPoint:CGPointMake(x,y)]];
        }
    }
    [part addObject:[NSValue valueWithCGPoint:CGPointMake(xmin,0)]];
    for (double x = xmin ; x < xmax ; x+=view->_xdelta) {
        double y= -[self eval:x];
        if (y > view->_ymin && y < view->_ymax) {
            [part addObject:[NSValue valueWithCGPoint:CGPointMake(x,y)]];
        }
    }
}

- (NSArray*)makeCurve:(CurveView*)view {
    NSMutableArray *curve= [NSMutableArray array];
   
    double roots[3];
    int nroots= [self findroots:roots];


    NSMutableArray *part=[NSMutableArray array];
    [curve addObject:part];
    [self addsection:part view:view xmin:roots[nroots-1] xmax:view->_xmax];

    if (nroots>1) {
        part=[NSMutableArray array];
        [curve addObject:part];
        [part addObject:[NSValue valueWithCGPoint:CGPointMake(roots[1],0)]];
        [self addsection:part view:view xmin:roots[0] xmax:roots[1]];
        [part addObject:[NSValue valueWithCGPoint:CGPointMake(roots[1],0)]];
    }

    printf("generated %d part curve\n", [curve count]);

    return curve; // [curve autorelease];
}

- (CGPoint)add:(CGPoint)p q:(CGPoint)q
{
    double s= (p.x==q.x && p.y==q.y) ? (3*p.x*p.x+a)/(2*p.y) : (p.y-q.y)/(p.x-q.x);

    double rx= s*s-(p.x+q.x);
    double ry= p.y+s*(rx-p.x);

	printf("(%5.2f,%5.2f) + (%5.2f,%5.2f) = (%5.2f,%5.2f)\n", p.x, p.y, q.x, q.y, rx, ry);

    return CGPointMake(rx, -ry);
}

// todo: return NaN when point is not near curve
- (CGPoint)projectoncurve:(CGPoint)p
{
    double y= [self eval:p.x];

    if (y<-100 || y>100) {
        return CGPointMake(NAN,NAN);
    }
    else {
        if (p.y<0) y=-y;
        return CGPointMake(p.x,y);
    }
}

- (CGPoint)scalarmultiply:(CGPoint)p c:(double)c
{
	CGPoint result=CGPointMake(0,0);
	CGPoint addend= p;
	unsigned int m= floor(c);
	while(m)
	{
		if (m&1) {
			result= [self add:result q:addend];
		}
		addend= [self add:addend q:addend];
		m>>=1;
	}
	printf("%d x (%5.2f,%5.2f) = (%5.2f,%5.2f)\n", m, p.x,p.y, result.x, result.y);
	return result;
}
@end
