code guessing, round #23 (completed)

-- Code by SoundOfSpouting#6980 (UID: 151149148639330304)


entry = bigdiv -- Pointfree!!


bigdiv :: [Double] -> [Double] -> ([Double], [Double])
bigdiv [] _ = ([], [])
bigdiv (head : tail) b = do
    let (div, temp_mod) = bigdiv tail b
    let (x, mod) = smalldiv (head : temp_mod) b
    (x : div, mod)


smalldiv :: [Double] -> [Double] -> (Double, [Double])
smalldiv [] _ = (0, [])
smalldiv x [] = smalldiv x [0] -- Spot the zero division
smalldiv [a] [b] = (a/b, [])
smalldiv (a : ta) (b : tb) = do
    let (q, mod) = smalldiv ta tb
    (q, a - q * b : mod)

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

#include "polydivide.h"

// Prints a double array terminated by INFINITY
void print_array(double const *a) {
  putchar('{');
  for (; !isinf(*a); ++a) { printf(" %.2f", *a); }
  printf(" }");
}

// Compares using Total Squared Error; passes if error < tolerance
bool test(double const *a, double const *b, double const *expected_q,
          double const *expected_r, double const tolerance) {

  printf("\nTEST: a = ");
  print_array(a);

  printf(", b = ");
  print_array(b);
  putchar('\n');

  double *r;
  double *q = entry(a, b, &r);

  bool success = true;

  double error = 0;
  for (double const *qi = q, *eqi = expected_q; !isinf(*eqi); ++qi, ++eqi) {
    double const residual = *qi - *eqi;
    error += residual * residual;
  }
  if (error < tolerance) {
    printf("Q SUCCESS (error: %f): Q = ", error);
    print_array(q);
    putchar('\n');
  } else {
    printf("Q FAILED (error: %f): \n  EXPECTED: ", error);
    print_array(expected_q);
    printf("\n  GOT: ");
    print_array(q);
    putchar('\n');
    success = false;
  }

  error = 0;
  for (double const *ri = r, *eri = expected_r; !isinf(*eri); ++ri, ++eri) {
    error += (*ri - *eri) * (*ri - *eri);
  }
  if (error < tolerance) {
    printf("R SUCCESS (error: %f): R = ", error);
    print_array(r);
    putchar('\n');
  } else {
    printf("R FAILED (error: %f): \n  EXPECTED: ", error);
    print_array(expected_r);
    printf("\n  GOT: ");
    print_array(r);
    putchar('\n');
    success = false;
  }

  free(r);
  free(q);

  return success;
}

struct test_case {
  double const *a;
  double const *b;
  double const *expected_q;
  double const *expected_r;
};

int main(int argc, char const *argv[]) {
  double const tolerance = .000000001;

  struct test_case const test_suite[] = {
      // Basic tests
      {(double const[]){-42, 0, -12, 1, INFINITY},
       (double const[]){-3, 1, INFINITY},
       (double const[]){-27, -9, 1, INFINITY},
       (double const[]){-123, INFINITY}},

      {(double const[]){-21, -5, 4, INFINITY},
       (double const[]){-3, 1, INFINITY}, (double const[]){7, 4, INFINITY},
       (double const[]){0, INFINITY}},

      {(double const[]){12, 30, 40, INFINITY},
       (double const[]){1, 3, 4, INFINITY}, (double const[]){10, INFINITY},
       (double const[]){2, 0, INFINITY}},

      // Testing a case where the remainder isn't degree 0
      {(double const[]){-13, 2, -9, 0, 4, 2, INFINITY},
       (double const[]){1, -2, 1, INFINITY},
       (double const[]){11, 14, 8, 2, INFINITY},
       (double const[]){-24, 10, INFINITY}},

      // Testing a case that actually gives fractions
      {(double const[]){0, 0, 0, 1, INFINITY},
       (double const[]){1, -2, 3, INFINITY},
       (double const[]){2.0 / 9.0, 1.0 / 3.0, INFINITY},
       (double const[]){-2.0 / 9.0, 1.0 / 9.0, INFINITY}},

      // Testing a numerator with a lower degree than the denominator
      {(double const[]){1, 2, INFINITY},
       (double const[]){-3, 1, 2, 1, INFINITY}, (double const[]){INFINITY},
       (double const[]){1, INFINITY}},

      // Testing empty numerator
      {(double const[]){INFINITY}, (double const[]){-3, 1, INFINITY},
       (double const[]){INFINITY}, (double const[]){INFINITY}}};

  int const num_tests = sizeof test_suite / sizeof *test_suite;
  int num_success = 0;

  for (int test_index = 0; test_index < num_tests; ++test_index) {
    struct test_case const test_case = test_suite[test_index];

    num_success += test(test_case.a, test_case.b, test_case.expected_q,
                        test_case.expected_r, tolerance);
  }

  printf("\n%d of %d tests passed\n", num_success, num_tests);

  return 0;
}

#ifndef POLYDIVIDE_H
#define POLYDIVIDE_H

#include <math.h>
#include <stdlib.h>
#include <string.h>

// Takes two polynomials a, b, ordered in order of increasing power, terminated
// by INFINITY and a pointer to a pre-allocated remainder polynomial r.
// Allocates memory to q and r
double *entry(double const *a, double const *b, double **r) {
  // We first need to find the end of the polynomial
  // This is because polynomial division starts at the highest power
  // This also lets us compute the length of the polynomials,
  // important because we to compute the length of the quotient
  // to allocate memory for it.
  double const *a_end = a, *b_end = b;

  while (!isinf(*a_end)) { ++a_end; }
  while (!isinf(*b_end)) { ++b_end; }

  // The number of coeffients (excludes INFINITY terminator)
  // This is one more than the degree of the polynomial
  size_t const a_len = a_end - a;
  size_t const b_len = b_end - b;

  // Handle the numerator having lower degree than denominator
  if (a_len < b_len) {
    double *const q = malloc(1 * sizeof *q);
    *q = INFINITY;

    // Just create a copy of the numerator as the remainder
    *r = malloc((a_len + 1) * sizeof **r);
    memcpy(*r, a, (a_len + 1) * sizeof **r);

    return q;
  }

  size_t const q_len = a_len - b_len + 1;

  double *const q = malloc((q_len + 1) * sizeof *q);
  double *q_end = q + q_len;

  // Terminate the quotient polynomial
  *q_end = INFINITY;
  --q_end;

  // Move back to the coefficient before INFINITY
  --a_end;
  --b_end;

  // temporary values for computation
  double *temps = malloc(b_len * sizeof *temps);
  double *const temps_end = temps + b_len - 1;
  memcpy(temps, a_end - b_len + 1, b_len * sizeof *temps);

  // a_end - a + 1 = number of coefficients left in a
  while (a_end - a + 1 >= b_len) {
    double const a_lead = *temps_end;
    double const b_lead = *b_end;

    double const lead_factor = a_lead / b_lead;

    *q_end = lead_factor;

    // Recalculate temporaries
    // Basically subtracting b * lead_factor from temps,
    // shifting over by one,
    // and then getting the next coefficient from a
    double *tptr;
    double const *bptr;
    for (tptr = temps_end, bptr = b_end; tptr != temps; --tptr, --bptr) {
      *tptr = *(tptr - 1) - lead_factor * *(bptr - 1);
    }

    // Last iteration
    // Since it's the last iteration, there are no more coefficients in a
    if (a_end - a + 1 == b_len) {
      break;
    }

    *temps = *(a_end - b_len);

    --a_end;
    --q_end;
  }

  // Temps now basically contains the remainder, just needs some touching up
  // meaning shifting it over by one and adding the INFINITY terminator
  memmove(temps, temps + 1, (b_len - 1) * sizeof *temps);
  *temps_end = INFINITY;

  *r = temps;
  return q;
}

#endif

entry=:[: (}: ; {:) ([ (] -/@,:&}. (* {:)) ] , %&{.~)^:(>:@-~&#)&.|.~

# bogopolynomialdivide
import struct
import random

def entry(a, b):
	"inefficient, also doesnt work in some cases because floating point but that isnt my fault"
	def l(p,x):
		if len(p)==1:return p*x
		else:return x*(p+l(p[:-1],x))
	while True:	
		p,q=[],[]
		for i in range(len(a)-len(b)):q+=struct.unpack('f',random.randbytes(4))
		for i in range(random.randint(len(a))):p+=struct.unpack('f',random.randbytes(4))
		e=1
		for i in range(500):
			x = struct.unpack('f',random.randbytes(4))
			if l(a,x)!=l(b,x)*l(p,x)+l(q,x):e=0
		if e>0.3409:return p,q