REM --------------------------------------------------------------------------------------
REM
REM Inspired by Project Euler - http://projecteuler.net - I needed a library to perform
REM     calculations with Really Large numbers.
REM
REM This INCLUDE file implements a subset of the GMP library at http://gmplib.org
REM
REM Integer functions implemented:
REM  ADD$, SUBSTRACT$, MULTIPLE$, DIVIDE$, MODULO$, SQUARE$, POWER$,
REM   FACTORIAL$, FIBONACCI$, ISPRIME, NEXTPRIME$
REM
REM Float functions implemented:
REM  FADD$, FSUBSTRACT$, FMULTIPLE$, FDIVIDE$, FSQUARE$, FPOWER$, PRECISION$
REM
REM Other functions:
REM  COMPARE, FCOMPARE, PRECISION
REM
REM (c) Peter van Eerten, September 2009 - GPL.
REM     Credits to the newLisp GMP module.
REM
REM Version 0.1: -initial release
REM Version 0.2: -added MODULO$, ISPRIME$
REM Version 0.3: -floating numbers return correct format now
REM              -incoming arguments in float functions always converted to float
REM              -added PRECISION to define precision of float calculations
REM Version 0.4: -fixed limitation in precision definition
REM Version 0.5: -better library naming
REM Version 0.6: -added heuristic library search
REM Version 0.7: -added nth-root
REM Version 0.8: -added COMPARE and FCOMPARE
REM Version 0.9: -adjustments for including separate functions, added INIT
REM Version 0.10: -ISPRIME should return integer, second argument should be integer
REM Version 0.11: -simplified heuristic library search
REM Version 0.12: -add GCD$ (Greatest Common Divisor)
REM --------------------------------------------------------------------------------------

CALL INIT

REM --------------------------------------------------------------------------------------------------

SUB INIT

        LOCAL gmp_lib$
        LOCAL gmp_seq = -1

    REM Check on the availability of GMP
    TRAP LOCAL
        CATCH GOTO Redo_Import

        REM Get the library
        LABEL Redo_Import
                INCR gmp_seq
                IF gmp_seq = 20 THEN GOTO Print_Error
                gmp_lib$ = CONCAT$("libgmp.so.", STR$(gmp_seq))
                IMPORT "__gmpz_init(long)" FROM gmp_lib$ TYPE void

    REM -------------------------- IMPORT FROM GMP ----------------------------------

    REM Integer funcs
    IMPORT "__gmpz_get_str(char*, int, long)" FROM gmp_lib$ TYPE char*
    IMPORT "__gmpz_set_str(long, char*, int)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpz_add(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_sub(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_mul(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_cdiv_q(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_mod(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_pow_ui(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_sqrt(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_root(long, long, long)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpz_cmp(long, long)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpz_fib_ui(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_fac_ui(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_gcd(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_probab_prime_p(long, long)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpz_nextprime(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpz_cmp_si(long, long)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpz_fac_ui(long, long)" FROM gmp_lib$ TYPE void
    REM Float funcs
    IMPORT "__gmpf_init2(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_set_default_prec(long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_set_prec(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_get_str(char*, long, int, int, long)" FROM gmp_lib$ TYPE char*
    IMPORT "__gmpf_set_str(long, char*, int)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpf_add(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_sub(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_mul(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_div(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_pow_ui(long, long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_sqrt(long, long)" FROM gmp_lib$ TYPE void
    IMPORT "__gmpf_cmp(long, long)" FROM gmp_lib$ TYPE int
    IMPORT "__gmpf_cmp_si(long, long)" FROM gmp_lib$ TYPE int

    REM Determine how many digits are significant in case of floats
    CONST gmp_signify = 64

    REM Determine memory size in bytes for GMP variables
    CONST gmp_int_size = 16
    CONST gmp_flt_size = 32

    REM Global use to avoid memswapping all the time
    gmp_op1 = MEMORY(gmp_int_size)
    gmp_op2 = MEMORY(gmp_int_size)
    gmp_res = MEMORY(gmp_int_size)
    gmp_fp1 = MEMORY(gmp_flt_size)
    gmp_fp2 = MEMORY(gmp_flt_size)
    gmp_fes = MEMORY(gmp_flt_size)

    REM Initialize these integer variables for GMP
    __gmpz_init(gmp_op1)
    __gmpz_init(gmp_op2)
    __gmpz_init(gmp_res)
    __gmpf_init2(gmp_fp1, gmp_signify)
    __gmpf_init2(gmp_fp2, gmp_signify)
    __gmpf_init2(gmp_fes, gmp_signify)

        REM Reset error trapping
        CATCH RESET
        TRAP SYSTEM

        GOTO Exit_Sub

        REM Show error
        LABEL Print_Error
                PRINT "GMP library not found!"
                END

        LABEL Exit_Sub

END SUB

REM -------------------------- IMPLEMENTATION ------------------------------------

FUNCTION ADD$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

    __gmpz_add(gmp_res, gmp_op1, gmp_op2)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION SUBSTRACT$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

    __gmpz_sub(gmp_res, gmp_op1, gmp_op2)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION MULTIPLY$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

    __gmpz_mul(gmp_res, gmp_op1, gmp_op2)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION DIVIDE$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

    IF ISFALSE(__gmpz_cmp_si(gmp_op2, 0)) THEN
        PRINT "ERROR: division by zero in GMP module!"
        END
    END IF
    __gmpz_cdiv_q(gmp_res, gmp_op1, gmp_op2)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION MODULO$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

    IF ISFALSE(__gmpz_cmp_si(gmp_op2, 0)) THEN
        PRINT "ERROR: division by zero in GMP module!"
        END
    END IF
    __gmpz_mod(gmp_res, gmp_op1, gmp_op2)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION POWER$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)

    __gmpz_pow_ui(gmp_res, gmp_op1, VAL(op2$))

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION SQUARE$(STRING op1$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)

    __gmpz_sqrt(gmp_res, gmp_op1)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION ROOT$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)

    __gmpz_root(gmp_res, gmp_op1, VAL(op2$))

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION COMPARE(STRING op1$, STRING op2$)

    LOCAL result

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

    result = __gmpz_cmp(gmp_op1, gmp_op2)

    RETURN result

END FUNCTION

FUNCTION GCD$(STRING op1$, STRING op2$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)
    __gmpz_set_str(gmp_op2, op2$, 10)

        __gmpz_gcd(gmp_res, gmp_op1, gmp_op2)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

SUB PRECISION(NUMBER n)

    __gmpf_set_prec(gmp_fp1, n)
    __gmpf_set_prec(gmp_fp2, n)
    __gmpf_set_prec(gmp_fes, n)

END SUB

FUNCTION FADD$(STRING fp1$, STRING fp2$)

    LOCAL pointer TYPE char*
    LOCAL result$
    LOCAL exp

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF
    IF ISFALSE(INSTR(fp2$, ".")) THEN
        fp2$ = CONCAT$(fp2$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)
    __gmpf_set_str(gmp_fp2, fp2$, 10)

    __gmpf_add(gmp_fes, gmp_fp1, gmp_fp2)

    pointer = __gmpf_get_str(0, ADDRESS(exp), 10, 0, gmp_fes)

    IF ISFALSE(exp) THEN
        result$ = CONCAT$("0.", pointer)
    ELIF exp < 0 THEN
        result$ = CONCAT$("0.", FILL$(ABS(exp), 48), pointer)
    ELSE
        IF LEN(pointer) > exp THEN
            result$ = CONCAT$(LEFT$(pointer, exp), ".", MID$(pointer, exp+1))
        ELIF LEN(pointer) EQ exp THEN
            result$ = CONCAT$(pointer, ".0")
        ELSE
            result$ = CONCAT$(pointer, FILL$(exp-LEN(pointer), 48), ".0")
        ENDIF
    ENDIF

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FSUBSTRACT$(STRING fp1$, STRING fp2$)

    LOCAL pointer TYPE char*
    LOCAL result$
    LOCAL exp

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF
    IF ISFALSE(INSTR(fp2$, ".")) THEN
        fp2$ = CONCAT$(fp2$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)
    __gmpf_set_str(gmp_fp2, fp2$, 10)

    __gmpf_sub(gmp_fes, gmp_fp1, gmp_fp2)

    pointer = __gmpf_get_str(0, ADDRESS(exp), 10, 0, gmp_fes)

    IF ISFALSE(exp) THEN
        result$ = CONCAT$("0.", pointer)
    ELIF exp < 0 THEN
        result$ = CONCAT$("0.", FILL$(ABS(exp), 48), pointer)
    ELSE
        IF LEN(pointer) > exp THEN
            result$ = CONCAT$(LEFT$(pointer, exp), ".", MID$(pointer, exp+1))
        ELIF LEN(pointer) EQ exp THEN
            result$ = CONCAT$(pointer, ".0")
        ELSE
            result$ = CONCAT$(pointer, FILL$(exp-LEN(pointer), 48), ".0")
        ENDIF
    ENDIF

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FMULTIPLY$(STRING fp1$, STRING fp2$)

    LOCAL pointer TYPE char*
    LOCAL result$
    LOCAL exp

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF
    IF ISFALSE(INSTR(fp2$, ".")) THEN
        fp2$ = CONCAT$(fp2$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)
    __gmpf_set_str(gmp_fp2, fp2$, 10)

    __gmpf_mul(gmp_fes, gmp_fp1, gmp_fp2)

    pointer = __gmpf_get_str(0, ADDRESS(exp), 10, 0, gmp_fes)

    IF ISFALSE(exp) THEN
        result$ = CONCAT$("0.", pointer)
    ELIF exp < 0 THEN
        result$ = CONCAT$("0.", FILL$(ABS(exp), 48), pointer)
    ELSE
        IF LEN(pointer) > exp THEN
            result$ = CONCAT$(LEFT$(pointer, exp), ".", MID$(pointer, exp+1))
        ELIF LEN(pointer) EQ exp THEN
            result$ = CONCAT$(pointer, ".0")
        ELSE
            result$ = CONCAT$(pointer, FILL$(exp-LEN(pointer), 48), ".0")
        ENDIF
    ENDIF

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FDIVIDE$(STRING fp1$, STRING fp2$)

    LOCAL pointer TYPE char*
    LOCAL result$
    LOCAL exp

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF
    IF ISFALSE(INSTR(fp2$, ".")) THEN
        fp2$ = CONCAT$(fp2$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)
    __gmpf_set_str(gmp_fp2, fp2$, 10)

    IF ISFALSE(__gmpf_cmp_si(gmp_fp2, 0)) THEN
        PRINT "ERROR: division by zero in GMP module!"
        END
    END IF
    __gmpf_div(gmp_fes, gmp_fp1, gmp_fp2)

    pointer = __gmpf_get_str(0, ADDRESS(exp), 10, 0, gmp_fes)

    IF ISFALSE(exp) THEN
        result$ = CONCAT$("0.", pointer)
    ELIF exp < 0 THEN
        result$ = CONCAT$("0.", FILL$(ABS(exp), 48), pointer)
    ELSE
        IF LEN(pointer) > exp THEN
            result$ = CONCAT$(LEFT$(pointer, exp), ".", MID$(pointer, exp+1))
        ELIF LEN(pointer) EQ exp THEN
            result$ = CONCAT$(pointer, ".0")
        ELSE
            result$ = CONCAT$(pointer, FILL$(exp-LEN(pointer), 48), ".0")
        ENDIF
    ENDIF

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FPOWER$(STRING fp1$, STRING fp2$)

    LOCAL pointer TYPE char*
    LOCAL result$
    LOCAL exp

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)

    __gmpf_pow_ui(gmp_fes, gmp_fp1, VAL(fp2$))

    pointer = __gmpf_get_str(0, ADDRESS(exp), 10, 0, gmp_fes)

    IF ISFALSE(exp) THEN
        result$ = CONCAT$("0.", pointer)
    ELIF exp < 0 THEN
        result$ = CONCAT$("0.", FILL$(ABS(exp), 48), pointer)
    ELSE
        IF LEN(pointer) > exp THEN
            result$ = CONCAT$(LEFT$(pointer, exp), ".", MID$(pointer, exp+1))
        ELIF LEN(pointer) EQ exp THEN
            result$ = CONCAT$(pointer, ".0")
        ELSE
            result$ = CONCAT$(pointer, FILL$(exp-LEN(pointer), 48), ".0")
        ENDIF
    ENDIF

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FSQUARE$(STRING fp1$)

    LOCAL pointer TYPE char*
    LOCAL result$
    LOCAL exp

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)

    __gmpf_sqrt(gmp_fes, gmp_fp1)

    pointer = __gmpf_get_str(0, ADDRESS(exp), 10, 0, gmp_fes)

    IF ISFALSE(exp) THEN
        result$ = CONCAT$("0.", pointer)
    ELIF exp < 0 THEN
        result$ = CONCAT$("0.", FILL$(ABS(exp), 48), pointer)
    ELSE
        IF LEN(pointer) > exp THEN
            result$ = CONCAT$(LEFT$(pointer, exp), ".", MID$(pointer, exp+1))
        ELIF LEN(pointer) EQ exp THEN
            result$ = CONCAT$(pointer, ".0")
        ELSE
            result$ = CONCAT$(pointer, FILL$(exp-LEN(pointer), 48), ".0")
        ENDIF
    ENDIF

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FCOMPARE(STRING fp1$, STRING fp2$)

    LOCAL result

    IF ISFALSE(INSTR(fp1$, ".")) THEN
        fp1$ = CONCAT$(fp1$, ".0")
    ENDIF
    IF ISFALSE(INSTR(fp2$, ".")) THEN
        fp2$ = CONCAT$(fp2$, ".0")
    ENDIF

    __gmpf_set_str(gmp_fp1, fp1$, 10)
    __gmpf_set_str(gmp_fp2, fp2$, 10)

    result = __gmpf_cmp(gmp_fp1, gmp_fp2)

    RETURN result

END FUNCTION

FUNCTION FIBONACCI$(STRING op1$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_fib_ui(gmp_res, VAL(op1$))

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION FACTORIAL$(STRING op1$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_fac_ui(gmp_res, VAL(op1$))

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION

FUNCTION ISPRIME(STRING op1$, NUMBER op2)

    LOCAL result

    __gmpz_set_str(gmp_op1, op1$, 10)

    result = __gmpz_probab_prime_p(gmp_op1, op2)

    RETURN result

END FUNCTION

FUNCTION NEXTPRIME$(STRING op1$)

    LOCAL pointer TYPE char*
    LOCAL result$

    __gmpz_set_str(gmp_op1, op1$, 10)

    __gmpz_nextprime(gmp_res, gmp_op1)

    pointer = __gmpz_get_str(0, 10, gmp_res)
    result$ = pointer

    FREE pointer

    RETURN result$

END FUNCTION