/*
 *  Eukleides  version 1.0.3
 *  Copyright (c) Christian Obrecht 2000-2004
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <stdlib.h>
#include <math.h>
#include "types.h"
#include "geometry.h"
#include "parser.tab.h"

extern void warning (char *);

extern void yyerror (char *);

extern void save (void *);

int Snd_degree (double *x1, double *x2, double a, double b, double c)
{
    double d;
    if (a == 0)
	return 0;
    d = b * b - 4 * a * c;
    if (d < 0)
	return 0;
    if (d == 0)
      {
	  *x1 = *x2 = -.5 * b / a;
	  return 1;
      }
    *x1 = (-b - sqrt (d)) / (2 * a);
    *x2 = (-b + sqrt (d)) / (2 * a);
    return 2;
}

double Pow (double x, double y)
{
    if (x < 0)
	return pow (x, rint (y));
    return pow (x, y);
}

double Sin (double x)
{
    return sin (RAD (x));
}

double Cos (double x)
{
    return cos (RAD (x));
}

double Tan (double x)
{
    return tan (RAD (x));
}

double Asin (double x)
{
    return DEG (asin (x));
}

double Acos (double x)
{
    return DEG (acos (x));
}

double Atan (double x)
{
    return DEG (atan (x));
}

double PP_distance (_point * A, _point * B)
{
    return LENG (A->x - B->x, A->y - B->y);
}

double PL_distance (_point * A, _line * l)
{
    return fabs ((l->y - A->y) * Cos (l->angle) -
		 (l->x - A->x) * Sin (l->angle));
}

double V_length (_vector * v)
{
    return LENG (v->x, v->y);
}

double S_length (_segment * s)
{
    return LENG (s->x1 - s->x2, s->y1 - s->y2);
}

double eccentricity (_conic * C)
{
    switch (C->kind)
      {
      case PARABOLA:
	  return 1;
      case ELLIPSE:
	  return sqrt (1 - C->b * C->b / (C->a * C->a));
      case HYPERBOLA:
	  return sqrt (1 + C->b * C->b / (C->a * C->a));
      }
    return 1;
}

double height (_point * A, _point * B, _point * C)
{
    double a, b, c;
    a = LENG (B->x - C->x, B->y - C->y);
    b = LENG (A->x - C->x, A->y - C->y);
    c = LENG (B->x - A->x, B->y - A->y);
    return sqrt ((a + b + c) * (a + b - c) * (a + c - b) * (b + c - a))
		/ (2 * b);
}

double angle (double x, double y)
{
    if ((ZERO (x)) && (ZERO (y)))
      {
	  warning ("undefined angle.");
	  return 0;
      }
    if (y >= 0)
      {
	  if (ZERO (x))
	      return 90;
	  else if (x > 0)
	      return DEG (atan (y / x));
	  else
	      return 180 + DEG (atan (y / x));
      }
    else
      {
	  if (ZERO (x))
	      return -90;
	  else if (x > 0)
	      return DEG (atan (y / x));
	  else
	      return -180 + DEG (atan (y / x));
      }
}

double format (double a)
{
    if (a > 180)
	return a - 360;
    else if (a <= -180)
	return a + 360;
    else
	return a;
}

double V_angle (_vector * v)
{
    return angle (v->x, v->y);
}

double S_angle (_segment * s)
{
    return angle (s->x2 - s->x1, s->y2 - s->y1);
}

double Co_angle (_conic * C)
{
    if (C->kind == PARABOLA)
	return format (C->d + 90);
    else
	return C->d;
}

double VV_angle (_vector * v1, _vector * v2)
{
    return format (V_angle (v2) - V_angle (v1));
}

double T_angle (_point * A, _point * B, _point * C)
{
    return format (angle (C->x - B->x, C->y - B->y) -
		   angle (A->x - B->x, A->y - B->y));
}

double C_argument (_circle * C, _point * A)
{
    return angle (A->x - C->x, A->y - C->y);
}

double Co_argument (_conic * C, _point * A)
{
    double c0, s0, x, y, z;

    c0 = Cos (C->d);
    s0 = Sin (C->d);
    x = (A->x - C->x) * c0 + (A->y - C->y) * s0;
    y = -(A->x - C->x) * s0 + (A->y - C->y) * c0;

    switch (C->kind)
      {
      case PARABOLA:
	  return x;
      case ELLIPSE:
	  if (ZERO (x))
	      return (y < 0) ? (-M_PI_2) : M_PI_2;
	  z = atan (C->a * y / (C->b * x));
	  if (x > 0)
	      return z;
	  else if (y < 0)
	      return -M_PI + z;
	  else
	      return M_PI + z;
      case HYPERBOLA:
	  z = atan (y / C->b);
	  if (x > 0)
	      return z;
	  else
	      return z + M_PI;
      }
    return 0;
}

double scalar (_vector * v1, _vector * v2)
{
    return v1->x * v2->x + v1->y * v2->y;
}

void P_assignment (symrec * s, _point * A)
{
    s->type = POINT;
    s->object.point = A;
}

_point *P_new (double x, double y)
{
    _point *ptr;

    ptr = (_point *) (malloc (sizeof (_point)));
    ptr->x = x;
    ptr->y = y;
    save ((void *) ptr);
    return ptr;
}

_point *P_polar (double r, double a)
{
    return P_new (r * Cos (a), r * Sin (a));
}

_point *P_line (_line * l, double d)
{
    return P_new (l->x + d * Cos (l->angle), l->y + d * Sin (l->angle));
}

_point *P_segment (_segment * s, double d)
{
    double a;
    a = S_angle (s);
    return P_new (s->x1 + d * Cos (a), s->y1 + d * Sin (a));
}

_point *P_circle (_circle * c, double a)
{
    return P_new (c->x + c->radius * Cos (a), c->y + c->radius * Sin (a));
}

_point *P_conic (_conic * C, double a)
{
    double c0, s0, A, B, x, y;

    c0 = Cos (C->d);
    s0 = Sin (C->d);
    A = C->a;
    B = C->b;

    switch (C->kind)
      {
      case PARABOLA:
	  x = a;
	  y = 2 * A * a * a;
	  break;
      case ELLIPSE:
	  x = A * cos (a);
	  y = B * sin (a);
	  break;
      case HYPERBOLA:
	  if (ZERO (cos (a)))
	      yyerror ("undefined point.");
	  else
	    {
		x = -A / cos (a);
		y = B * tan (a);
	    }
      }
    return P_new (C->x + c0 * x - s0 * y, C->y + s0 * x + c0 * y);
}

_point *P_translation (_point * A, _vector * v)
{
    return P_new (A->x + v->x, A->y + v->y);
}

_point *P_reflection (_point * A, _line * l)
{
    double a, le;

    a = format (2 * l->angle - angle (A->x - l->x, A->y - l->y));
    le = LENG (A->x - l->x, A->y - l->y);
    return P_new (l->x + le * Cos (a), l->y + le * Sin (a));
}

_point *P_rotation (_point * A, _point * O, double a)
{
    double c, s, x, y;
    c = Cos (a);
    s = Sin (a);
    x = A->x - O->x;
    y = A->y - O->y;
    return P_new (O->x + c * x - s * y, O->y + s * x + c * y);
}

_point *P_L_projection (_point * A, _line * l)
{
    _line *tmp;

    tmp = L_L_perpendicular (l, A);
    return P_LL_intersection (l, tmp);
}

_point *P_LL_projection (_point * A, _line * l1, _line * l2)
{
    _line *tmp;

    tmp = L_L_parallel (l2, A);
    return P_LL_intersection (l1, tmp);
}

_point *P_homothecy (_point * A, _point * O, double k)
{
    return P_new (O->x + k * (A->x - O->x), O->y + k * (A->y - O->y));
}

_point *P_2_barycenter (_point * A, double al, _point * B, double be)
{
    double s;
    s = al + be;
    if (ZERO (s))
      {
	  yyerror ("undefined barycenter.");
	  return NULL;
      }
    return P_new ((al * A->x + be * B->x) / s, (al * A->y + be * B->y) / s);
}

_point *P_3_barycenter (_point * A, double al, _point * B, double be,
			_point * C, double ga)
{
    double s;
    s = al + be + ga;
    if (ZERO (s))
      {
	  yyerror ("undefined barycenter.");
	  return NULL;
      }
    return P_new ((al * A->x + be * B->x + ga * C->x) / s,
		  (al * A->y + be * B->y + ga * C->y) / s);
}

_point *P_4_barycenter (_point * A, double al, _point * B, double be,
			_point * C, double ga, _point * D, double de)
{
    double s;
    s = al + be + ga + de;
    if (ZERO (s))
      {
	  yyerror ("undefined barycenter.");
	  return NULL;
      }
    return P_new ((al * A->x + be * B->x + ga * C->x + de * D->x) / s,
		  (al * A->y + be * B->y + ga * C->y + de * D->y) / s);
}

_point *P_LL_intersection (_line * l1, _line * l2)
{
    double c1, c2, s1, s2, d1, d2, d;
    c1 = Cos (l1->angle);
    c2 = Cos (l2->angle);
    s1 = Sin (l1->angle);
    s2 = Sin (l2->angle);
    d = c1 * s2 - c2 * s1;
    if (ZERO (d))
      {
	  yyerror ("parallel lines.");
	  return NULL;
      }
    d1 = s1 * l1->x - c1 * l1->y;
    d2 = s2 * l2->x - c2 * l2->y;
    return P_new ((c1 * d2 - c2 * d1) / d, (s1 * d2 - s2 * d1) / d);
}

_point *P_abscissa (_line * l, double x)
{
    if (ZERO (Cos (l->angle)))
      {
	  yyerror ("undefined point.");
	  return NULL;
      }
    return P_new (x, Tan (l->angle) * (x - l->x) + l->y);
}

_point *P_ordinate (_line * l, double y)
{
    if (ZERO (Sin (l->angle)))
      {
	  yyerror ("undefined point.");
	  return NULL;
      }
    return P_new ((y - l->y) / Tan (l->angle) + l->x, y);
}

_point *P_midpoint (_segment * s)
{
    return P_new ((s->x1 + s->x2) / 2, (s->y1 + s->y2) / 2);
}

_point *P_begin (_segment * s)
{
    return P_new (s->x1, s->y1);
}

_point *P_end (_segment * s)
{
    return P_new (s->x2, s->y2);
}

_point *P_center (_circle * c)
{
    return P_new (c->x, c->y);
}

_point *P_Co_center (_conic * C)
{
    return P_new (C->x, C->y);
}

_point *P_orthocenter (_point * A, _point * B, _point * C)
{
    _line *tmp1, *tmp2;

    tmp1 = L_altitude (A, B, C);
    tmp2 = L_altitude (B, A, C);
    return P_LL_intersection (tmp1, tmp2);
}

void V_assignment (symrec * s, _vector * v)
{
    s->type = VECTOR;
    s->object.vector = v;
}

_vector *V_NN_new (double x, double y)
{
    _vector *ptr;
    ptr = (_vector *) (malloc (sizeof (_vector)));
    ptr->x = x;
    ptr->y = y;
    save ((void *) ptr);
    return ptr;
}

_vector *V_NA_new (double l, double a)
{
    return V_NN_new (l * Cos (a), l * Sin (a));
}

_vector *V_PP_new (_point * A, _point * B)
{
    return V_NN_new (B->x - A->x, B->y - A->y);
}

_vector *V_L_new (_line * l)
{
    return V_NN_new (Cos (l->angle), Sin (l->angle));
}

_vector *V_S_new (_segment * s)
{
    return V_NN_new (s->x2 - s->x1, s->y2 - s->y1);
}

_vector *V_rotation (_vector * u, double a)
{
    return V_NN_new (Cos (a) * u->x - Sin (a) * u->y,
		     Sin (a) * u->x + Cos (a) * u->y);
}

_vector *V_add (_vector * u, _vector * v)
{
    return V_NN_new (u->x + v->x, u->y + v->y);
}

_vector *V_sub (_vector * u, _vector * v)
{
    return V_NN_new (u->x - v->x, u->y - v->y);
}

_vector *V_mult (double k, _vector * u)
{
    return V_NN_new (k * u->x, k * u->y);
}

void L_assignment (symrec * s, _line * l)
{
    s->type = LINE;
    s->object.line = l;
}

_line *L_NNA_new (double x, double y, double a)
{
    _line *ptr;
    ptr = (_line *) (malloc (sizeof (_line)));
    ptr->x = x;
    ptr->y = y;
    ptr->angle = a;
    save ((void *) ptr);
    return ptr;
}

_line *L_PA_new (_point * A, double a)
{
    return L_NNA_new (A->x, A->y, a);
}

_line *L_PP_new (_point * A, _point * B)
{
    return L_PA_new (A, angle (B->x - A->x, B->y - A->y));
}

_line *L_PV_new (_point * A, _vector * v)
{
    return L_PA_new (A, V_angle (v));
}

_line *L_S_new (_segment * s)
{
    return L_NNA_new (s->x1, s->y1, S_angle (s));
}

_line *L_CA_new (_circle * c, double a)
{
    return L_NNA_new (c->x + c->radius * Cos (a), c->y + c->radius * Sin (a),
		      format (a + 90));
}

_line *L_Co_new (_conic * C, double a)
{
    double c0, s0, A, B, x, y, d;

    c0 = Cos (C->d);
    s0 = Sin (C->d);
    A = C->a;
    B = C->b;

    switch (C->kind)
      {
      case PARABOLA:
	  x = a;
	  y = 2 * A * a * a;
	  d = angle (1, 4 * A * a);
	  break;
      case ELLIPSE:
	  x = A * cos (a);
	  y = B * sin (a);
	  d = angle (-A * sin (a), B * cos (a));
	  break;
      case HYPERBOLA:
	  if (ZERO (cos (a)))
	      x = y = 0;
	  else
	    {
		x = -A / cos (a);
		y = B * tan (a);
	    }
	  d = angle (-A * sin (a), B);
      }
    return L_NNA_new (C->x + c0 * x - s0 * y, C->y + s0 * x + c0 * y,
		      C->d + d);
}

_line *L_translation (_line * l, _vector * v)
{
    return L_NNA_new (l->x + v->x, l->y + v->y, l->angle);
}

_line *L_reflection (_line * l, _line * l1)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (l->x, l->y);
    tmp2 = P_reflection (tmp1, l1);
    return L_PA_new (tmp2, format (2 * l1->angle - l->angle));
}

_line *L_rotation (_line * l, _point * O, double a)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (l->x, l->y);
    tmp2 = P_rotation (tmp1, O, a);
    return L_PA_new (tmp2, format (l->angle + a));
}

_line *L_homothecy (_line * l, _point * O, double k)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (l->x, l->y);
    tmp2 = P_homothecy (tmp1, O, k);
    return L_PA_new (tmp2, l->angle);
}

_line *L_L_parallel (_line * l, _point * A)
{
    return L_PA_new (A, l->angle);
}

_line *L_S_parallel (_segment * s, _point * A)
{
    return L_PA_new (A, S_angle (s));
}

_line *L_L_perpendicular (_line * l, _point * A)
{
    return L_PA_new (A, format (l->angle + 90));
}

_line *L_S_perpendicular (_segment * s, _point * A)
{
    return L_PA_new (A, format (S_angle (s) + 90));
}

_line *L_S_bisector (_segment * s)
{
    return L_NNA_new ((s->x1 + s->x2) / 2, (s->y1 + s->y2) / 2,
		      format (S_angle (s) + 90));
}

_line *L_P_bisector (_point * A, _point * B, _point * C)
{
    _line *ptr;
    ptr = L_PP_new (B, A);
    ptr->angle = format (ptr->angle + T_angle (A, B, C) / 2);
    return ptr;
}

_line *L_LL_bisector (_line * l1, _line * l2)
{
    _point *A;
    double a;

    if (ZERO (format (l1->angle - l2->angle)))
	return L_NNA_new ((l1->x + l2->x) / 2, (l1->y + l2->y) / 2,
			  l1->angle);
    A = P_LL_intersection (l1, l2);
    if (abs (format (l1->angle - l2->angle)) > 90)
	a = format ((l1->angle + l2->angle) / 2 + 90);
    else
	a = format ((l1->angle + l2->angle) / 2);
    return L_PA_new (A, a);
}

_line *L_altitude (_point * A, _point * B, _point * C)
{
    return L_PA_new (A, format (angle (B->x - C->x, B->y - C->y) + 90));
}

_line *L_median (_point * A, _point * B, _point * C)
{
    _point *tmp;
    tmp = P_2_barycenter (B, 1, C, 1);
    return L_PP_new (A, tmp);
}

_line *L_invert (_line * l)
{
    return L_NNA_new (l->x, l->y, format (l->angle - 180));
}

void S_assignment (symrec * s, _segment * sg)
{
    s->type = SEGMENT;
    s->object.segment = sg;
}

_segment *S_N_new (double x1, double y1, double x2, double y2)
{
    _segment *ptr;
    ptr = (_segment *) (malloc (sizeof (_segment)));
    ptr->x1 = x1;
    ptr->y1 = y1;
    ptr->x2 = x2;
    ptr->y2 = y2;
    save ((void *) ptr);
    return ptr;
}

_segment *S_P_new (_point * A, _point * B)
{
    return S_N_new (A->x, A->y, B->x, B->y);
}

_segment *S_V_new (_point * A, _vector * v)
{
    return S_N_new (A->x, A->y, A->x + v->x, A->y + v->y);
}

_segment *S_NA_new (_point * A, double l, double a)
{
    return S_N_new (A->x, A->y, A->x + l * Cos (a), A->y + l * Sin (a));
}

_segment *S_CA_new (_circle * c, double a)
{
    return S_N_new (c->x, c->y, c->x + c->radius * Cos (a),
		    c->y + c->radius * Sin (a));
}

_segment *S_translation (_segment * s, _vector * v)
{
    return S_N_new (s->x1 + v->x, s->y1 + v->y, s->x2 + v->x, s->y2 + v->y);
}

_segment *S_reflection (_segment * s, _line * l)
{
    _point *tmp1, *tmp2, *tmp3, *tmp4;

    tmp1 = P_new (s->x1, s->y1);
    tmp2 = P_new (s->x2, s->y2);
    tmp3 = P_reflection (tmp1, l);
    tmp4 = P_reflection (tmp2, l);
    return S_P_new (tmp3, tmp4);
}

_segment *S_rotation (_segment * s, _point * O, double a)
{
    _point *tmp1, *tmp2, *tmp3, *tmp4;

    tmp1 = P_new (s->x1, s->y1);
    tmp2 = P_new (s->x2, s->y2);
    tmp3 = P_rotation (tmp1, O, a);
    tmp4 = P_rotation (tmp2, O, a);
    return S_P_new (tmp3, tmp4);
}

_segment *S_homothecy (_segment * s, _point * O, double k)
{
    _point *tmp1, *tmp2, *tmp3, *tmp4;

    tmp1 = P_new (s->x1, s->y1);
    tmp2 = P_new (s->x2, s->y2);
    tmp3 = P_homothecy (tmp1, O, k);
    tmp4 = P_homothecy (tmp2, O, k);
    return S_P_new (tmp3, tmp4);
}

_segment *S_invert (_segment * s)
{
    return S_N_new (s->x2, s->y2, s->x1, s->y1);
}

void C_assignment (symrec * s, _circle * c)
{
    s->type = CIRCLE;
    s->object.circle = c;
}

_circle *C_NNN_new (double x, double y, double r)
{
    _circle *ptr;
    ptr = (_circle *) (malloc (sizeof (_circle)));
    ptr->x = x;
    ptr->y = y;
    ptr->radius = r;
    save ((void *) ptr);
    return ptr;
}

_circle *C_PP_new (_point * A, _point * B)
{
    return C_NNN_new ((A->x + B->x) / 2, (A->y + B->y) / 2,
		      PP_distance (A, B) / 2);
}

_circle *C_PPP_new (_point * A, _point * B, _point * C)
{
    _line *tmp1, *tmp2;
    _point *tmp3;

    tmp1 =
	L_NNA_new ((A->x + B->x) / 2, (A->y + B->y) / 2,
		   format (angle (A->x - B->x, A->y - B->y) + 90));
    tmp2 =
	L_NNA_new ((C->x + B->x) / 2, (C->y + B->y) / 2,
		   format (angle (C->x - B->x, C->y - B->y) + 90));
    tmp3 = P_LL_intersection (tmp1, tmp2);
    return C_NNN_new (tmp3->x, tmp3->y,
		      LENG (A->x - tmp3->x, A->y - tmp3->y));
}

_circle *C_PN_new (_point * A, double r)
{
    return C_NNN_new (A->x, A->y, r);
}

_circle *C_translation (_circle * c, _vector * v)
{
    return C_NNN_new (c->x + v->x, c->y + v->y, c->radius);
}

_circle *C_reflection (_circle * c, _line * l)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (c->x, c->y);
    tmp2 = P_reflection (tmp1, l);
    return C_NNN_new (tmp2->x, tmp2->y, c->radius);
}

_circle *C_rotation (_circle * c, _point * O, double a)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (c->x, c->y);
    tmp2 = P_rotation (tmp1, O, a);
    return C_NNN_new (tmp2->x, tmp2->y, c->radius);
}

_circle *C_homothecy (_circle * c, _point * O, double k)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (c->x, c->y);
    tmp2 = P_homothecy (tmp1, O, k);
    return C_NNN_new (tmp2->x, tmp2->y, k * c->radius);
}

_circle *C_incircle (_point * A, _point * B, _point * C)
{
    _line *tmp1, *tmp2;
    _point *tmp3;
    double a, b, c;

    tmp1 = L_P_bisector (A, B, C);
    tmp2 = L_P_bisector (B, C, A);
    tmp3 = P_LL_intersection (tmp1, tmp2);
    a = LENG (B->x - C->x, B->y - C->y);
    b = LENG (A->x - C->x, A->y - C->y);
    c = LENG (B->x - A->x, B->y - A->y);
    return C_NNN_new (tmp3->x, tmp3->y,
		      sqrt ((a + b - c) * (a + c - b) * (b + c - a)
			    / (a + b + c)) / 2);
}


void Co_assignment (symrec * s, _conic * C)
{
    s->type = CONIC;
    s->object.conic = C;
}

_conic *Co_N_new (double x, double y, double d, double a, double b, int k)
{
    _conic *ptr;
    ptr = (_conic *) (malloc (sizeof (_conic)));
    ptr->x = x;
    ptr->y = y;
    ptr->d = d;
    ptr->a = a;
    ptr->b = b;
    ptr->kind = k;
    save ((void *) ptr);
    return ptr;
}

_conic *Pa_new (_point * A, _line * l)
{
    _point *tmp1;
    _vector *tmp2;
    if (ZERO (PL_distance (A, l)))
	yyerror ("undefined conic curve.");
    tmp1 = P_L_projection (A, l);
    tmp2 = V_PP_new (tmp1, A);
    return Co_N_new ((A->x + tmp1->x) / 2, (A->y + tmp1->y) / 2,
		     format (V_angle (tmp2) - 90), .25 / V_length (tmp2), 0,
		     PARABOLA);
}

_conic *Co_new (_point * A, _line * l, double e)
{
    _point *tmp1;
    _vector *tmp2;
    double d, f, g;

    if (ZERO (PL_distance (A, l)))
	yyerror ("undefined conic curve.");
    if (ZERO (fabs (e) - 1))
	return Pa_new (A, l);
    tmp1 = P_L_projection (A, l);
    tmp2 = V_PP_new (A, tmp1);
    d = V_angle (tmp2);
    g = V_length (tmp2);
    f = e * e / (1 - e * e);

    return Co_N_new (A->x - Cos (d) * f * g, A->y - Sin (d) * f * g, d,
		     sqrt (f * g * g / (1 - e * e)), sqrt (fabs (f * g * g)),
		     (fabs (e) < 1) ? ELLIPSE : HYPERBOLA);
}

_conic *Co_PP_new (_point * A, _point * B, double a)
{
    _segment *tmp1;
    _point *tmp2;
    double b, c;
    int k;

    tmp1 = S_P_new (A, B);
    tmp2 = P_midpoint (tmp1);
    c = S_length (tmp1) / 2;

    if (ZERO (a))
	return Co_N_new (0, 0, 0, 0, 0, -1);
    if (c < fabs (a))
      {
	  b = sqrt (a * a - c * c);
	  k = ELLIPSE;
      }
    else
      {
	  b = sqrt (c * c - a * a);
	  k = HYPERBOLA;
      }
    return Co_N_new (tmp2->x, tmp2->y, S_angle (tmp1), fabs (a), b, k);
}

_conic *Co_translation (_conic * C, _vector * v)
{
    return Co_N_new (C->x + v->x, C->y + v->y, C->d, C->a, C->b, C->kind);
}

_conic *Co_reflection (_conic * C, _line * l)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (C->x, C->y);
    tmp2 = P_reflection (tmp1, l);
    return Co_N_new (tmp2->x, tmp2->y, format (2 * l->angle - C->d), C->a,
		     C->b, C->kind);
}

_conic *Co_rotation (_conic * C, _point * O, double a)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (C->x, C->y);
    tmp2 = P_rotation (tmp1, O, a);
    return Co_N_new (tmp2->x, tmp2->y, format (a + C->d), C->a, C->b,
		     C->kind);
}

_conic *Co_homothecy (_conic * C, _point * O, double k)
{
    _point *tmp1, *tmp2;

    tmp1 = P_new (C->x, C->y);
    tmp2 = P_homothecy (tmp1, O, k);
    return Co_N_new (tmp2->x, tmp2->y, C->d, k * C->a, k * C->b, C->kind);
}


void LC_intersection (symrec * s1, symrec * s2, _line * l, _circle * c)
{
    double r, d, a, h, p;

    r = c->radius;
    d = LENG (c->x - l->x, c->y - l->y);
    a = format (angle (c->x - l->x, c->y - l->y) - l->angle);
    h = d * Sin (a);
    if (h > r)
	yyerror ("undefined points.");
    p = sqrt (r * r - h * h);
    P_assignment (s1, P_line (l, d * Cos (a) - p));
    P_assignment (s2, P_line (l, d * Cos (a) + p));
}

void CC_intersection (symrec * s1, symrec * s2, _circle * c1, _circle * c2)
{
    double d, r1, r2, h, p, a;

    d = LENG (c1->x - c2->x, c1->y - c2->y);
    r1 = c1->radius;
    r2 = c2->radius;
    if ((d > r1 + r2) || (d < fabs (r1 - r2)))
	yyerror ("undefined points.");
    h = sqrt ((r1 + r2 + d) * (r1 + d - r2) * (r1 + r2 - d) * (r2 + d - r1))
	     / (2 * d);
    p = sqrt (r1 * r1 - h * h) * (r2 * r2 > r1 * r1 + d * d ? -1 : 1);
    a = angle (c2->x - c1->x, c2->y - c1->y);
    P_assignment (s1,
		  P_new (c1->x + p * Cos (a) - h * Sin (a),
			 c1->y + p * Sin (a) + h * Cos (a)));
    P_assignment (s2,
		  P_new (c1->x + p * Cos (a) + h * Sin (a),
			 c1->y + p * Sin (a) - h * Cos (a)));
}

void LCo_intersection (symrec * s1, symrec * s2, _line * l, _conic * C)
{
    double c0, s0, dx, dy, a, b, c, A, B, a2, b2, c2, A2, B2, x1, y1, x2, y2;

    c0 = Cos (C->d);
    s0 = Sin (C->d);
    dx = l->x - C->x;
    dy = l->y - C->y;
    a = -Sin (l->angle - C->d);
    b = Cos (l->angle - C->d);
    c = (-c0 * dx - s0 * dy) * a + (s0 * dx - c0 * dy) * b;
    A = C->a;
    B = C->b;
    a2 = a * a;
    b2 = b * b;
    c2 = c * c;
    A2 = A * A;
    B2 = B * B;

    switch (C->kind)
      {
      case PARABOLA:
	  if (ZERO (b))
	    {
		x1 = x2 = -c / a;
		y1 = y2 = 2 * A * c2 / a2;
	    }
	  else
	    {
		if (Snd_degree (&x1, &x2, 2 * b * A, a, c))
		  {
		      y1 = -(a * x1 + c) / b;
		      y2 = -(a * x2 + c) / b;
		  }
		else
		    yyerror ("undefined points.");
	    }
	  break;
      case ELLIPSE:
	  if (ZERO (b))
	    {
		x1 = x2 = -c / a;
		y1 = B * sqrt (1 - x1 * x1 / A2);
		y2 = -y1;
	    }
	  else
	    {
		if (Snd_degree
		    (&x1, &x2, a2 * A2 + b2 * B2, 2 * a * c * A2,
		     A2 * (c2 - b2 * B2)))
		  {
		      y1 = -(a * x1 + c) / b;
		      y2 = -(a * x2 + c) / b;
		  }
		else
		    yyerror ("undefined points.");
	    }
	  break;
      case HYPERBOLA:
	  if (ZERO (b))
	    {
		x1 = x2 = -c / a;
		y1 = B * sqrt (x1 * x1 / A2 - 1);
		y2 = -y1;

	    }
	  else
	    {
		if (Snd_degree
		    (&x1, &x2, -a2 * A2 + b2 * B2, -2 * a * c * A2,
		     -A2 * (c2 + b2 * B2)))
		  {
		      y1 = -(a * x1 + c) / b;
		      y2 = -(a * x2 + c) / b;
		  }
		else
		    yyerror ("undefined points.");
	    }
      }
    P_assignment (s1,
		  P_new (C->x + c0 * x1 - s0 * y1, C->y + s0 * x1 + c0 * y1));
    P_assignment (s2,
		  P_new (C->x + c0 * x2 - s0 * y2, C->y + s0 * x2 + c0 * y2));
}

void vertices (symrec * s1, symrec * s2, _conic * C)
{
    _point *tmp;

    if (C->kind == PARABOLA)
      {
	  tmp = P_new (C->x, C->y);
	  P_assignment (s1, tmp);
	  P_assignment (s2, tmp);
      }
    else
      {
	  P_assignment (s1,
			P_new (C->x + C->a * Cos (C->d),
			       C->y + C->a * Sin (C->d)));
	  P_assignment (s2,
			P_new (C->x - C->a * Cos (C->d),
			       C->y - C->a * Sin (C->d)));
      }
}

void foci (symrec * s1, symrec * s2, _conic * C)
{
    _point *tmp;
    double c;

    if (C->kind == PARABOLA)
      {
	  tmp =
	      P_new (C->x - .5 * C->a * Sin (C->d),
		     C->y + .5 * C->a * Cos (C->d));
	  P_assignment (s1, tmp);
	  P_assignment (s2, tmp);
      }
    else
      {
	  if (C->kind == ELLIPSE)
	      c = sqrt (C->a * C->a - C->b * C->b);
	  else
	      c = sqrt (C->a * C->a + C->b * C->b);
	  P_assignment (s1,
			P_new (C->x + c * Cos (C->d), C->y + c * Sin (C->d)));
	  P_assignment (s2,
			P_new (C->x - c * Cos (C->d), C->y - c * Sin (C->d)));
      }
}

void
T_assignment (symrec * s, double x, double y, 
	      double h, double p, double l, double a)
{
    P_assignment (s,
		  P_new (x + l * (p * Cos (a) - h * Sin (a)),
			 y + l * (h * Cos (a) + p * Sin (a))));
}

void
T_assignments (symrec * s1, symrec * s2, double x, double y, double b,
	      double h, double p, double l, double a)
{
    P_assignment (s1, P_new (x + l * b * Cos (a), y + l * b * Sin (a)));
    P_assignment (s2,
		  P_new (x + l * (p * Cos (a) - h * Sin (a)),
			 y + l * (h * Cos (a) + p * Sin (a))));
}

void T_check (symrec * s1, symrec * s2)
{
    if (s1->type != POINT) {
	P_assignment (s1, P_new (0, 0));
	P_assignment (s2, P_new (6, 0));
    }
    if (s2->type != POINT)
	P_assignment (s2, P_new (s1->object.point->x+6, s1->object.point->y));
}
  
void T_scalenes (symrec * s1, symrec * s2, symrec * s3, double l, double a)
{
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    T_assignments (s2, s3, s1->object.point->x, s1->object.point->y, 1,
		  .61237244, .375, l, a);
}

void
T_3N_triangle (symrec * s1, symrec * s2, symrec * s3, double l1, double l2,
	      double l3, double a)
{
    double h, p;
    
    if (ZERO (l1) || (l1 > l2 + l3) || (l2 > l1 + l3) || (l3 > l1 + l2))
	yyerror ("invalid length.");
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    h = sqrt ((l1 + l2 + l3) * (l1 + l3 - l2) * (l1 + l2 - l3) * (l2 + l3 - l1))
	 / (2 * l1);
    p = sqrt (fabs(l3 * l3 - h * h))
	 * ((l2 * l2 > (l1 * l1 + l3 * l3)) ? -1 : 1);
    T_assignments (s2, s3, s1->object.point->x, s1->object.point->y, l1, h, p,
		  1, a);
}

void
T_2N_triangle (symrec * s1, symrec * s2, symrec * s3, double l2, double l3)
{
    double a, l1, h, p;
   
    T_check (s1, s2);
    l1 = PP_distance(s1->object.point,s2->object.point);
    if (ZERO (l1) || (l1 > l2 + l3) || (l2 > l1 + l3) || (l3 > l1 + l2))
	yyerror ("invalid length.");
    a = angle (s2->object.point->x-s1->object.point->x,
	       s2->object.point->y-s1->object.point->y);
    h = sqrt ((l1 + l2 + l3) * (l1 + l3 - l2) * (l1 + l2 - l3) * (l2 + l3 - l1))
	 / (2 * l1);
    p = sqrt (fabs(l3 * l3 - h * h))
	 * ((l2 * l2 > (l1 * l1 + l3 * l3)) ? -1 : 1);
    T_assignment (s3, s1->object.point->x, s1->object.point->y, h, p, 1, a);
}

void
T_2A_triangle (symrec * s1, symrec * s2, symrec * s3, double l, double a1,
	      double a2, double a)
{
    _line *tmp1, *tmp2;

    if (ZERO (Sin (a1 + a2)))
	yyerror ("invalid angle.");
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    tmp1 = L_PA_new (s1->object.point, a1 + a);
    P_assignment (s2,
		  P_new (s1->object.point->x + l * Cos (a),
			 s1->object.point->y + l * Sin (a)));
    tmp2 = L_PA_new (s2->object.point, format (180 - a2 + a));
    P_assignment (s3, P_LL_intersection (tmp1, tmp2));
}

void
T_A_triangle (symrec * s1, symrec * s2, symrec * s3, double a1, double a2)
{
    _line *tmp1, *tmp2;
    
    double a;
    if (ZERO (Sin (a1 + a2)))
	yyerror ("invalid angle.");
    T_check (s1, s2);
    a = angle (s2->object.point->x-s1->object.point->x,
	       s2->object.point->y-s1->object.point->y);    
    tmp1 = L_PA_new (s1->object.point, a1 + a);
    tmp2 = L_PA_new (s2->object.point, format (180 - a2 + a));
    P_assignment (s3, P_LL_intersection (tmp1, tmp2));
}

void
T_2N_right (symrec * s1, symrec * s2, symrec * s3, double l1, double l2,
	    double a)
{
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    T_assignments (s2, s3, s1->object.point->x, s1->object.point->y, l1, l2, 0,
		  1, a);
}

void
T_N_right (symrec * s1, symrec * s2, symrec * s3, double l2)
{
    double a;
    T_check (s1, s2);
    a = angle (s2->object.point->x-s1->object.point->x,
	       s2->object.point->y-s1->object.point->y);
    T_assignment (s3, s1->object.point->x, s1->object.point->y, l2, 0,
		  1, a);
}

void
T_NA_isosceles (symrec * s1, symrec * s2, symrec * s3, double l1, double a1,
	       double a)
{
    if (ZERO (Cos (a1)))
	yyerror ("invalid angle.");
    else {
	if (s1->type != POINT)
	    P_assignment (s1, P_new (0, 0));
	T_assignments (s2, s3, s1->object.point->x, s1->object.point->y, l1,
		      l1 * Tan (a1) / 2, l1 / 2, 1, a);
    }
}

void
T_A_isosceles (symrec * s1, symrec * s2, symrec * s3, double a1)
{
    double a, l1;

    if (ZERO (Cos (a1)))
	yyerror ("invalid angle.");
    else {
	T_check (s1, s2);
	a = angle (s2->object.point->x-s1->object.point->x,
		   s2->object.point->y-s1->object.point->y);
	l1 = PP_distance(s1->object.point,s2->object.point);
	T_assignment (s3, s1->object.point->x, s1->object.point->y,
		      l1 * Tan (a1) / 2, l1 / 2, 1, a);
    }
}

void
T_2N_isosceles (symrec * s1, symrec * s2, symrec * s3, double l1, double l2,
	       double a)
{
    if (l1 > 2 * l2)
	yyerror ("invalid length.");
    else {
	if (s1->type != POINT)
	    P_assignment (s1, P_new (0, 0));
	T_assignments (s2, s3, s1->object.point->x, s1->object.point->y, l1,
		      sqrt (l2 * l2 - l1 * l1 / 4), l1 / 2, 1, a);
    }
}

void
T_N_isosceles (symrec * s1, symrec * s2, symrec * s3, double l2)
{
    double a, l1;

    if (l1 > 2 * l2)
	yyerror ("invalid length.");
    else {
	T_check (s1, s2);
	a = angle (s2->object.point->x-s1->object.point->x,
		   s2->object.point->y-s1->object.point->y);
	l1 = PP_distance(s1->object.point,s2->object.point);
	T_assignment (s3, s1->object.point->x, s1->object.point->y,
		      sqrt (l2 * l2 - l1 * l1 / 4), l1 / 2, 1, a);
    }
}

void T_E_equilateral (symrec * s1, symrec * s2, symrec * s3)
{
    double l, a;
    a = angle (s2->object.point->x-s1->object.point->x,
	       s2->object.point->y-s1->object.point->y);
    l = PP_distance(s1->object.point,s2->object.point);
    T_assignment (s3, s1->object.point->x, s1->object.point->y,
	      .866025, .5, l, a);
}

void T_equilateral (symrec * s1, symrec * s2, symrec * s3, double l, double a)
{
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    T_assignments (s2, s3, s1->object.point->x, s1->object.point->y, 1,
	      .866025, .5, l, a);
}

void
Q_E_parallelogram (symrec * s1, symrec * s2, symrec * s3, symrec * s4)
{
    if (s1->type != POINT) {
	P_assignment (s1, P_new (0, 0));
	P_assignment (s2, P_new (5, 0));
	P_assignment (s3, P_new (6.0352, 3.8637));
    }
    if (s2->type != POINT)
	P_assignment (s2, P_new (s1->object.point->x+5, s1->object.point->y));
    if (s2->type != POINT || s3->type != POINT)
	P_assignment (s3, s2->object.point);
    P_assignment (s4, P_new (
	s1->object.point->x + (s3->object.point->x - s2->object.point->x),
	s1->object.point->y + (s3->object.point->y - s2->object.point->y)));
}

void
Q_S_parallelogram (symrec * s1, symrec * s2, symrec * s3, symrec * s4,
		   double l2, double a1)
{
    double l1, a, xu, xv, yu, yv;
    if (s1->type != POINT) {
	P_assignment (s1, P_new (0, 0));
	P_assignment (s2, P_new (5, 0));
    }
    if (s2->type != POINT)
	P_assignment (s2, P_new (s1->object.point->x+5, s1->object.point->y));
    a = angle (s2->object.point->x-s1->object.point->x,
	       s2->object.point->y-s1->object.point->y);
    l1 = PP_distance(s1->object.point,s2->object.point);
    xu = l1 * Cos (a);
    xv = l2 * Cos (a + a1);
    yu = l1 * Sin (a);
    yv = l2 * Sin (a + a1);
    P_assignment (s3,
		  P_new (s1->object.point->x + xu + xv,
			 s1->object.point->y + yu + yv));
    P_assignment (s4,
		  P_new (s1->object.point->x + xv, s1->object.point->y + yv));
}

void
Q_parallelogram (symrec * s1, symrec * s2, symrec * s3, symrec * s4,
		 double l1, double l2, double a1, double a)
{
    double xu, xv, yu, yv;

    xu = l1 * Cos (a);
    xv = l2 * Cos (a + a1);
    yu = l1 * Sin (a);
    yv = l2 * Sin (a + a1);
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    P_assignment (s2,
		  P_new (s1->object.point->x + xu, s1->object.point->y + yu));
    P_assignment (s3,
		  P_new (s1->object.point->x + xu + xv,
			 s1->object.point->y + yu + yv));
    P_assignment (s4,
		  P_new (s1->object.point->x + xv, s1->object.point->y + yv));
}

void
Q_V_parallelogram (symrec * s1, symrec * s2, symrec * s3, symrec * s4,
		   _vector * u, _vector * v)
{
    if (s1->type != POINT)
	P_assignment (s1, P_new (0, 0));
    P_assignment (s2,
		  P_new (s1->object.point->x + u->x,
			 s1->object.point->y + u->y));
    P_assignment (s3,
		  P_new (s1->object.point->x + u->x + v->x,
			 s1->object.point->y + u->y + v->y));
    P_assignment (s4,
		  P_new (s1->object.point->x + v->x,
			 s1->object.point->y + v->y));
}

void
Q_E_rectangle (symrec * s1, symrec * s2, symrec * s3, symrec * s4)
{
    if (s1->type != POINT) {
	P_assignment (s1, P_new (0, 0));
	P_assignment (s2, P_new (6, 0));
	P_assignment (s3, P_new (6, 3.7082));
    }
    if (s2->type != POINT)
	P_assignment (s2, P_new (s1->object.point->x+6, s1->object.point->y));
    if (s3->type != POINT)
	P_assignment (s3, s2->object.point);
    P_assignment (s4, P_new (
	s1->object.point->x + (s3->object.point->x - s2->object.point->x),
	s1->object.point->y + (s3->object.point->y - s2->object.point->y)));
}

void
Q_S_rectangle (symrec * s1, symrec * s2, symrec * s3, symrec * s4,
		   double l2)
{
    double l1, a, xu, xv, yu, yv;
    if (s1->type != POINT) {
	P_assignment (s1, P_new (0, 0));
	P_assignment (s2, P_new (6, 0));
    }
    if (s2->type != POINT)
	P_assignment (s2, P_new (s1->object.point->x+6, s1->object.point->y));
    a = angle (s2->object.point->x-s1->object.point->x,
	       s2->object.point->y-s1->object.point->y);
    l1 = PP_distance(s1->object.point,s2->object.point);
    xu = l1 * Cos (a);
    xv = -l2 * Sin (a);
    yu = l1 * Sin (a);
    yv = l2 * Cos (a);
    P_assignment (s3,
		  P_new (s1->object.point->x + xu + xv,
			 s1->object.point->y + yu + yv));
    P_assignment (s4,
		  P_new (s1->object.point->x + xv, s1->object.point->y + yv));
}

void
Q_square (symrec * s1, symrec * s2, symrec * s3, symrec * s4)
{
    if (s1->type != POINT) {
	P_assignment (s1, P_new (0, 0));
	P_assignment (s2, P_new (4, 0));
    }
    if (s2->type != POINT)
	P_assignment (s2, P_new (s1->object.point->x+4, s1->object.point->y));
    Q_S_rectangle (s1, s2, s3, s4,
	    PP_distance(s1->object.point,s2->object.point));
}

void pentagon (symrec * s1, symrec * s2, symrec * s3, symrec * s4,
	       symrec * s5, _point * O, double r, double a)
{
    P_assignment (s1, P_new (O->x + r * Cos (a), O->y + r * Sin (a)));
    P_assignment (s2,
		  P_new (O->x + r * Cos (a + 72), O->y + r * Sin (a + 72)));
    P_assignment (s3,
		  P_new (O->x + r * Cos (a + 144), O->y + r * Sin (a + 144)));
    P_assignment (s4,
		  P_new (O->x + r * Cos (a + 216), O->y + r * Sin (a + 216)));
    P_assignment (s5,
		  P_new (O->x + r * Cos (a + 288), O->y + r * Sin (a + 288)));
}

void hexagon (symrec * s1, symrec * s2, symrec * s3, symrec * s4, symrec * s5,
	      symrec * s6, _point * O, double r, double a)
{
    P_assignment (s1, P_new (O->x + r * Cos (a), O->y + r * Sin (a)));
    P_assignment (s2,
		  P_new (O->x + r * Cos (a + 60), O->y + r * Sin (a + 60)));
    P_assignment (s3,
		  P_new (O->x + r * Cos (a + 120), O->y + r * Sin (a + 120)));
    P_assignment (s4,
		  P_new (O->x + r * Cos (a + 180), O->y + r * Sin (a + 180)));
    P_assignment (s5,
		  P_new (O->x + r * Cos (a + 240), O->y + r * Sin (a + 240)));
    P_assignment (s6,
		  P_new (O->x + r * Cos (a + 300), O->y + r * Sin (a + 300)));
}