Changeset 0c2a2e in git


Ignore:
Timestamp:
Sep 27, 2023, 1:26:41 PM (8 months ago)
Author:
slap <slaplagne@…>
Branches:
(u'spielwiese', 'd08f5f0bb3329b8ca19f23b74cb1473686415c3a')
Children:
3d86bb7cdede7aa763e50a9436afbeb714fc26c2
Parents:
8bdde5e4b92ee38ec82a194d2a8c7f52517dd2e2
Message:
old optimization removed

The previous implementation of the optimization algorithm is removed.
Location:
Singular/LIB
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/.singularhistory

    r8bdde5 r0c2a2e  
    25932593timer - t;
    25942594list l1 = integralBasis(f, 2, "atOrigin");
    2595 
    25962595l1
    25972596;
     2597LIB"integralbasis.lib";
     2598ring r = 0,(x,y),dp;
     2599poly f = (y^5 + y^4*x^7 + 2*x^8) * (y^3 + 7*x^4) * (y^7 + 2*x^12) * (y^11 + 2*x^18) + y^30;
     2600int t = timer;
     2601list l1 = integralBasis(f, 2, "atOrigin");
     2602LIB"integralbasis.lib";
     2603ring r = 0,(x,y),(dp, L(100000));
     2604poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52) + y^36;
     2605int t = timer;
     2606list l1 = integralBasis(f, 2);
     2607ibSeg[1]
     2608segmentSlope
     2609out
     2610LIB"integralbasis.lib";
     2611ring r = 0,(x,y),(dp, L(100000));
     2612poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     2613int t = timer;
     2614list l1 = integralBasis(f, 2);
     2615LIB"integralbasis.lib";
     2616ring r = 0,(x,y),(dp, L(100000));
     2617poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     2618int t = timer;
     2619list l1 = integralBasis(f, 2, "atOrigin");
     2620timer - t;
     2621LIB"integralbasis.lib";
     2622ring r = 0,(x,y),(dp, L(100000));
     2623poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     2624int t = timer;
     2625list l1 = integralBasis(f, 2, "atOrigin");
     2626timer - t;
     2627LIB"integralbasis.lib";
     2628ring r = 0,(x,y),(dp, L(100000));
     2629poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     2630int t = timer;
     2631list l1 = integralBasis(f, 2, "atOrigin");
     2632timer - t;
     2633l1;
     2634t = timer;
     2635list l2 = integralBasis(f, 2, "noOpti");
     2636timer - t;
     2637LIB"integralbasis.lib";
     2638ring r = 0,(x,y),(dp, L(100000));
     2639poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     2640int t = timer;
     2641list l1 = integralBasis(f, 2, "atOrigin");
     2642timer - t;
     2643t = timer;
     2644list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2645timer - t;
     2646LIB"integralbasis.lib";
     2647ring r = 0,(x,y),(dp, L(100000));
     2648poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2649int t = timer;
     2650list l1 = integralBasis(f, 2, "atOrigin");
     2651timer - t;
     2652t = timer;
     2653list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2654timer - t;
     2655reduce(l1[1], groebner(l2[1]));
     2656reduce(l2[1], groebner(l1[1]));
     2657with(algcurves);
     2658f := (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2659intBasis := integral_basis(f, x, y, {x});
     2660LIB"integralbasis.lib";
     2661ring r = 0,(x,y),(dp, L(100000));
     2662poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2663int t = timer;
     2664list l1 = integralBasis(f, 2, "atOrigin");
     2665l1;
     2666LIB"integralbasis.lib";
     2667ring r = 0,(x,y),(dp, L(100000));
     2668poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2669int t = timer;
     2670list l1 = integralBasis(f, 2, "atOrigin");
     2671timer - t;
     2672t = timer;
     2673list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2674timer - t;
     2675[200~reduce(l1[1], groebner(l2[1]));
     2676reduce(l2[1], groebner(l1[1]));
     2677LIB"integralbasis.lib";
     2678ring r = 0,(x,y),(dp, L(100000));
     2679poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2680int t = timer;
     2681list l1 = integralBasis(f, 2, "atOrigin");
     2682timer - t;
     2683t = timer;
     2684list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2685timer - t;
     2686reduce(l1[1], groebner(l2[1]));
     2687reduce(l2[1], groebner(l1[1]));
     2688LIB"integralbasis.lib";
     2689ring r = 0,(x,y),(dp, L(100000));
     2690poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2691int t = timer;
     2692list l1 = integralBasis(f, 2, "atOrigin");
     2693timer - t;
     2694t = timer;
     2695list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2696timer - t;
     2697reduce(l1[1], groebner(l2[1]));
     2698reduce(l2[1], groebner(l1[1]));
     2699LIB"integralbasis.lib";
     2700ring r = 0,(x,y),(dp, L(100000));
     2701poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2702int t = timer;
     2703list l1 = integralBasis(f, 2, "atOrigin");
     2704timer - t;
     2705t = timer;
     2706list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2707timer - t;
     2708reduce(l1[1], groebner(l2[1]));
     2709reduce(l2[1], groebner(l1[1]));
     2710LIB"integralbasis.lib";
     2711ring r = 0,(x,y),(dp, L(100000));
     2712poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2) + y^25;
     2713int t = timer;
     2714list l1 = integralBasis(f, 2, "atOrigin");
     2715l1;
     2716timer - t;
     2717t = timer;
     2718list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2719timer - t;
     2720reduce(l1[1], groebner(l2[1]));
     2721reduce(l2[1], groebner(l1[1]));
     2722LIB"integralbasis.lib";
     2723ring r = 0,(x,y),(dp, L(100000));
     2724poly f = (y^5 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^34;
     2725int t = timer;
     2726list l1 = integralBasis(f, 2, "atOrigin");
     2727timer - t;
     2728t = timer;
     2729list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2730timer - t;
     2731reduce(l1[1], groebner(l2[1]));
     2732reduce(l2[1], groebner(l1[1]));
     2733LIB"integralbasis.lib";
     2734ring r = 0,(x,y),(dp, L(100000));
     2735poly f = (y^15 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^44;
     2736int t = timer;
     2737list l1 = integralBasis(f, 2, "atOrigin");
     2738timer - t;
     2739t = timer;
     2740list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2741timer - t;
     2742reduce(l1[1], groebner(l2[1]));
     2743reduce(l2[1], groebner(l1[1]));
     2744LIB"integralbasis.lib";
     2745ring r = 0,(x,y),(dp, L(100000));
     2746poly f = (y^15 + 2*x^7)*(y^15 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^60;
     2747int t = timer;
     2748list l1 = integralBasis(f, 2, "atOrigin");
     2749timer - t;
     2750t = timer;
     2751list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2752timer - t;
     2753LIB"integralbasis.lib";
     2754ring r = 0,(x,y),(dp, L(100000));
     2755poly f = (y^15 + 2*x^7)*(y^15 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^60;
     2756int t = timer;
     2757list l1 = integralBasis(f, 2, "atOrigin");
     2758timer - t;
     2759l1;
     2760LIB"integralbasis.lib";
     2761ring r = 0, (x,y), dp;
     2762poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2763int t = timer;
     2764list l1 = integralBasis(f,2, "atOrigin");
     2765LIB"integralbasis.lib";
     2766ring r = 0, (x,y), dp;
     2767poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2768int t = timer;
     2769list l1 = integralBasis(f,2, "atOrigin");
     2770LIB"integralbasis.lib";
     2771ring r = 0, (x,y), dp;
     2772poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2773int t = timer;
     2774list l1 = integralBasis(f,2, "atOrigin");
     2775LIB"integralbasis.lib";
     2776ring r = 0, (x,y), dp;
     2777poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2778int t = timer;
     2779list l1 = integralBasis(f,2, "atOrigin");
     2780basering
     2781fFrac
     2782LIB"integralbasis.lib";
     2783ring r = 0, (x,y), dp;
     2784poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2785int t = timer;
     2786list l1 = integralBasis(f,2, "atOrigin");
     2787fNew
     2788sD
     2789setring R
     2790fFrac
     2791printlevel = 5;
     2792LIB"integralbasis.lib";
     2793ring r = 0, (x,y), dp;
     2794poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2795int t = timer;
     2796list l1 = integralBasis(f,2, "atOrigin");
     2797printlevel = 5;
     2798LIB"integralbasis.lib";
     2799ring r = 0, (x,y), dp;
     2800poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2801int t = timer;
     2802list l1 = integralBasis(f,2, "atOrigin");
     2803classes
     2804def S = classes[1][1];
     2805setring S
     2806PE
     2807basering
     2808LIB"integralbasis.lib";
     2809ring r = 0, (x,y), dp;
     2810poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2811int t = timer;
     2812list l1 = integralBasis(f,2, "atOrigin");
     2813LIB"integralbasis.lib";
     2814ring r = 0, (x,y), dp;
     2815poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2816int t = timer;
     2817list l1 = integralBasis(f,2, "atOrigin");
     2818classes
     2819def S = classes[1][1];
     2820setring S;
     2821PE
     2822setring R;
     2823f
     2824list l = puiseux(f, -1, 1);
     2825l;
     2826with(algcurves);
     2827f := (y^7 + x^4) * (y^7 + y^5*x^3 + x^4)+y^16;
     2828puiseux(f,x=0,y,0);
     2829LIB"integralbasis.lib";
     2830ring r = 0, (x,y), dp;
     2831poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2832int t = timer;
     2833list l1 = integralBasis(f,2, "atOrigin");
     2834LIB"integralbasis.lib";
     2835ring r = 0, (x,y), dp;
     2836poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2837int t = timer;
     2838list l1 = integralBasis(f,2, "atOrigin");
     2839fFr
     2840den
     2841basering
     2842PE
     2843LIB"integralbasis.lib";
     2844ring r = 0, (x,y), dp;
     2845poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2846int t = timer;
     2847list l1 = integralBasis(f,2, "atOrigin");
     2848setring S
     2849PE
     2850setring R
     2851fFr
     2852dMP
     2853pardeg(minpoly);
     2854LIB"integralbasis.lib";
     2855ring r = 0, (x,y), dp;
     2856poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2857int t = timer;
     2858list l1 = integralBasis(f,2, "atOrigin");
     2859d
     2860setring S;
     2861PE
     2862LIB"integralbasis.lib";
     2863ring r = 0, (x,y), dp;
     2864poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2865int t = timer;
     2866list l1 = integralBasis(f,2, "atOrigin");
     2867fFr
     2868LIB"integralbasis.lib";
     2869ring r = 0, (x,y), dp;
     2870poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2871int t = timer;
     2872list l1 = integralBasis(f,2, "atOrigin");
     2873fF
     2874basering;
     2875den
     2876d
     2877buildPolyGroundXRoot(fFr, den, d);
     2878LIB"integralbasis.lib";
     2879ring r = 0, (x,y), dp;
     2880poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2881int t = timer;
     2882list l1 = integralBasis(f,2, "atOrigin");
     2883buildPolyGroundXRoot(fF, den, d);
     2884LIB"integralbasis.lib";
     2885ring r = 0, (x,y), dp;
     2886poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2887int t = timer;
     2888list l1 = integralBasis(f,2, "atOrigin");
     2889LIB"integralbasis.lib";
     2890ring r = 0, (x,y), dp;
     2891poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2892int t = timer;
     2893list l1 = integralBasis(f,2, "atOrigin");
     2894LIB"integralbasis.lib";
     2895ring r = 0, (x,y), dp;
     2896poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2897int t = timer;
     2898list l1 = integralBasis(f,2, "atOrigin");
     2899LIB"integralbasis.lib";
     2900ring r = 0, (x,y), dp;
     2901poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2902int t = timer;
     2903list l1 = integralBasis(f,2, "atOrigin");
     2904LIB"integralbasis.lib";
     2905ring r = 0, (x,y), dp;
     2906poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2907int t = timer;
     2908list l1 = integralBasis(f,2, "atOrigin");
     2909l1;
     2910LIB"integralbasis.lib";
     2911ring r = 0, (x,y), dp;
     2912poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2913int t = timer;
     2914list l1 = integralBasis(f, 2, "atOrigin");
     2915timer - t;
     2916t = timer;
     2917list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2918timer - t;
     2919reduce(l1[1], groebner(l2[1]));
     2920reduce(l2[1], groebner(l1[1]));
     2921LIB"integralbasis.lib";
     2922ring r = 0, (x,y), dp;
     2923poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2924int t = timer;
     2925list l1 = integralBasis(f, 2, "atOrigin");
     2926LIB"integralbasis.lib";
     2927ring r = 0, (x,y), dp;
     2928poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2929int t = timer;
     2930list l1 = integralBasis(f, 2, "atOrigin");
     2931LIB"integralbasis.lib";
     2932ring r = 0, (x,y), dp;
     2933poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2934int t = timer;
     2935list l1 = integralBasis(f, 2, "atOrigin");
     2936LIB"integralbasis.lib";
     2937ring r = 0, (x,y), dp;
     2938poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     2939int t = timer;
     2940list l1 = integralBasis(f, 2, "atOrigin");
     2941printlevel = 5;
     2942LIB"integralbasis.lib";
     2943ring r = 0,(x,y),(dp, L(100000));
     2944poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     2945int t = timer;
     2946list l1 = integralBasis(f, 2, "atOrigin");
     2947timer - t;
     2948t = timer;
     2949list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2950timer - t;
     2951reduce(l1[1], groebner(l2[1]));
     2952reduce(l2[1], groebner(l1[1]));
     2953LIB"integralbasis.lib";
     2954ring r = 0,(x,y),(dp, L(100000));
     2955poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52) + y^100;
     2956int t = timer;
     2957list l1 = integralBasis(f, 2, "atOrigin");
     2958timer - t;
     2959t = timer;
     2960list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2961timer - t;
     2962reduce(l1[1], groebner(l2[1]));
     2963reduce(l2[1], groebner(l1[1]));
     2964[200~LIB"integralbasis.lib";
     2965ring r = 0,(x,y),(dp, L(100000));
     2966poly f = (y^15 + 2*x^7)*(y^15 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^60;
     2967int t = timer;
     2968list l1 = integralBasis(f, 2, "atOrigin");
     2969LIB"integralbasis.lib";
     2970ring r = 0,(x,y),(dp, L(100000));
     2971poly f = (y^15 + 2*x^7)*(y^15 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^60;
     2972int t = timer;
     2973list l1 = integralBasis(f, 2, "atOrigin");
     2974LIB"integralbasis.lib";
     2975ring r = 0,(x,y),(dp, L(100000));
     2976poly f = (y^15 + 2*x^7)*(y^15 + 2*x^8)*(y^9 + 7*x^2)*(y^7+3*x^2)*(y^11+7*x6*y6+x^12) + y^60;
     2977int t = timer;
     2978list l1 = integralBasis(f, 2, "atOrigin");
     2979timer - t;
     2980t = timer;
     2981list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2982timer - t;
     2983LIB"integralbasis.lib";
     2984ring r = 0, (x,y), dp;
     2985poly f = (y7 + x4) * (y7 + y5x3 + x4) + y30;
     2986int t = timer;
     2987list l1 = integralBasis(f, 2, "atOrigin");
     2988timer - t;
     2989t = timer;
     2990list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     2991timer - t;
     2992LIB"integralbasis.lib";
     2993ring r = 0, (x,y), dp;
     2994poly f = (y7 + x4) * (y7 + y5x3 + x4) + y30;
     2995int t = timer;
     2996list l1 = integralBasis(f, 2, "atOrigin");
     2997timer - t;
     2998t = timer;
     2999list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3000timer - t;
     3001LIB"integralbasis.lib";
     3002ring r = 0, (x,y), dp;
     3003poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     3004int t = timer;
     3005list l1 = integralBasis(f, 2, "atOrigin");
     3006timer - t;
     3007t = timer;
     3008list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3009timer - t;
     3010reduce(l1[1], groebner(l2[1]));
     3011reduce(l2[1], groebner(l1[1]));
     3012LIB"integralbasis.lib";
     3013ring r = 0,(x,y),dp;
     3014int k = 40; int d = 50;
     3015poly f = y^2 + x^(k+1) + y^d;
     3016int t = timer;
     3017list l1 = integralBasis(f, 2, "atOrigin");
     3018timer - t;
     3019t = timer;
     3020list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3021timer - t;
     3022reduce(l1[1], groebner(l2[1]));
     3023reduce(l2[1], groebner(l1[1]));
     3024LIB"integralbasis.lib";
     3025ring r = 0,(x,y),dp;
     3026int k = 400; int d = 500;
     3027poly f = y^2 + x^(k+1) + y^d;
     3028int t = timer;
     3029list l1 = integralBasis(f, 2, "atOrigin");
     3030timer - t;
     3031t = timer;
     3032list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3033timer - t;
     3034reduce(l1[1], groebner(l2[1]));
     3035reduce(l2[1], groebner(l1[1]));
     3036LIB"integralbasis.lib";
     3037ring r = 0,(x,y),dp;
     3038int k = 40; int d = 50;
     3039poly f = x*(x^(k-1) + y^2) + y^d;
     3040list l1 = integralBasis(f, 2, "atOrigin");
     3041timer - t;
     3042t = timer;
     3043list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3044timer - t;
     3045reduce(l1[1], groebner(l2[1]));
     3046reduce(l2[1], groebner(l1[1]));
     3047int t = timer;
     3048list l1 = integralBasis(f, 2, "atOrigin");
     3049timer - t;
     3050t = timer;
     3051list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3052timer - t;
     3053reduce(l1[1], groebner(l2[1]));
     3054reduce(l2[1], groebner(l1[1]));
     3055LIB"integralbasis.lib";
     3056ring r = 0,(x,y),dp;
     3057int k = 400; int d = 500;
     3058poly f = x*(x^(k-1) + y^2) + y^d;
     3059int t = timer;
     3060list l1 = integralBasis(f, 2, "atOrigin");
     3061timer - t;
     3062t = timer;
     3063list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3064timer - t;
     3065LIB"integralbasis.lib";
     3066ring r = 0,(x,y),dp;
     3067int k = 400; int d = 500;
     3068poly f = x*(x^(k-1) + y^2) + y^d;
     3069int t = timer;
     3070list l1 = integralBasis(f, 2, "atOrigin");
     3071timer - t;
     3072t = timer;
     3073list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3074timer - t;
     3075LIB"integralbasis.lib";
     3076ring r = 0,(x,y),dp;
     3077int k = 400; int d = 500;
     3078poly f = x*(x^(k-1) + y^2) + y^d;
     3079int t = timer;
     3080list l1 = integralBasis(f, 2, "atOrigin");
     3081timer - t;
     3082t = timer;
     3083list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3084timer - t;
     3085
     3086LIB"integralbasis.lib";
     3087ring r = 0,(x,y),dp;
     3088int k = 400; int d = 500;
     3089poly f = x*(x^(k-1) + y^2) + y^d;
     3090int t = timer;
     3091list l1 = integralBasis(f, 2, "atOrigin");
     3092timer - t;
     3093t = timer;
     3094list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3095timer - t;
     3096LIB"integralbasis.lib";
     3097ring r = 0,(x,y),dp;
     3098int k = 400; int d = 500;
     3099poly f = x*(x^(k-1) + y^2) + y^d;
     3100int t = timer;
     3101list l1 = integralBasis(f, 2, "atOrigin");
     3102timer - t;
     3103t = timer;
     3104list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3105timer - t;
     3106LIB"integralbasis.lib";
     3107ring r = 0,(x,y),dp;
     3108int k = 400; int d = 500;
     3109poly f = x*(x^(k-1) + y^2) + y^d;
     3110int t = timer;
     3111list l1 = integralBasis(f, 2, "atOrigin");
     3112timer - t;
     3113t = timer;
     3114list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3115timer - t;
     3116min(10,20);
     3117LIB"integralbasis.lib";
     3118ring r = 0,(x,y),dp;
     3119int k = 400; int d = 500;
     3120poly f = x*(x^(k-1) + y^2) + y^d;
     3121int t = timer;
     3122list l1 = integralBasis(f, 2, "atOrigin");
     3123timer - t;
     3124t = timer;
     3125list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3126timer - t;
     3127LIB"integralbasis.lib";
     3128ring r = 0,(x,y),dp;
     3129int k = 400; int d = 500;
     3130poly f = x*(x^(k-1) + y^2) + y^d;
     3131int t = timer;
     3132list l1 = integralBasis(f, 2, "atOrigin");
     3133timer - t;
     3134t = timer;
     3135list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3136timer - t;
     3137LIB"integralbasis.lib";
     3138ring r = 0,(x,y),dp;
     3139int k = 400; int d = 500;
     3140poly f = x*(x^(k-1) + y^2) + y^d;
     3141int t = timer;
     3142list l1 = integralBasis(f, 2, "atOrigin");
     3143timer - t;
     3144t = timer;
     3145list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3146timer - t;
     3147reduce(l1[1], groebner(l2[1]));
     3148reduce(l2[1], groebner(l1[1]));
     3149LIB"integralbasis.lib";
     3150ring r = 0,(x,y),dp;
     3151int k = 400; int d = 500;
     3152poly f = y^2 + x^(k+1) + y^d;
     3153int t = timer;
     3154list l1 = integralBasis(f, 2, "atOrigin");
     3155timer - t;
     3156t = timer;
     3157list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3158timer - t;
     3159reduce(l1[1], groebner(l2[1]));
     3160reduce(l2[1], groebner(l1[1]));
     3161LIB"integralbasis.lib";
     3162ring r = 0,(x,y),dp;
     3163int k = 400; int d = 500;
     3164poly f = x*(x^(k-1) + y^2) + y^d;
     3165int t = timer;
     3166list l1 = integralBasis(f, 2, "atOrigin");
     3167timer - t;
     3168t = timer;
     3169list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3170timer - t;
     3171LIB"integralbasis.lib";
     3172ring r = 0, (x,y), dp;
     3173int k = 5;
     3174int d = 10;
     3175poly f = polyDK(d, k, 1231);
     3176int t = timer;
     3177list l1 = integralBasis(f, 2, "atOrigin");
     3178timer - t;
     3179t = timer;
     3180list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3181timer - t;
     3182reduce(l1[1], groebner(l2[1]));
     3183reduce(l2[1], groebner(l1[1]));
     3184l1;
     3185l2;
     3186f;
     3187LIB"integralbasis.lib";
     3188ring r = 0, (x,y), dp;
     3189int k = 5;
     3190int d = 10;
     3191poly f = polyDK(d, k, 1231);
     3192int t = timer;
     3193list l1 = integralBasis(f, 2, "atOrigin");
     3194l1;
     3195f;
     3196deg(f);
     3197subst(f,x,0);
     3198puiseux(f,-1,1);
     3199def S = puiseux(f,-1,1);
     3200setring S;
     3201def S = puiseux(f,-1,1)[1];
     3202setring S;
     3203PE;
     3204basering;
     3205LIB"integralbasis.lib";
     3206ring r = 0, (x,y), dp;
     3207int k = 5;
     3208int d = 10;
     3209poly f = polyDK(d, k, 1231);
     3210int t = timer;
     3211list l1 = integralBasis(f, 2, "atOrigin");
     3212l1;
     3213LIB"integralbasis.lib";
     3214ring r = 0, (x,y), dp;
     3215int k = 5;
     3216int d = 10;
     3217poly f = polyDK(d, k, 1231);
     3218int t = timer;
     3219list l1 = integralBasis(f, 2, "atOrigin");
     3220LIB"integralbasis.lib";
     3221ring r = 0, (x,y), dp;
     3222int k = 5;
     3223int d = 10;
     3224poly f = polyDK(d, k, 1231);
     3225int t = timer;
     3226list l1 = integralBasis(f, 2, "atOrigin");
     3227l1;
     3228ring r = 0,(x,y),dp;
     3229int k = 400; int d = 500;
     3230poly f = x*(x^(k-1) + y^2) + y^d;
     3231int t = timer;
     3232list l1 = integralBasis(f, 2, "atOrigin");
     3233timer - t;
     3234t = timer;
     3235list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3236LIB"integralbasis.lib";
     3237ring r = 0,(x,y),dp;
     3238int k = 400; int d = 500;
     3239poly f = x*(x^(k-1) + y^2) + y^d;
     3240int t = timer;
     3241list l1 = integralBasis(f, 2, "atOrigin");
     3242l1;
     3243timer - t;
     3244t = timer;
     3245list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3246timer - t;
     3247l2;
     3248LIB"integralbasis.lib";
     3249ring r = 0,(x,y),dp;
     3250int k = 400; int d = 500;
     3251poly f = x*(x^(k-1) + y^2) + y^d;
     3252int t = timer;
     3253list l1 = integralBasis(f, 2, "atOrigin");
     3254l1;
     3255LIB"integralbasis.lib";
     3256ring r = 0,(x,y),dp;
     3257int k = 400; int d = 500;
     3258poly f = x*(x^(k-1) + y^2) + y^d;
     3259int t = timer;
     3260list l1 = integralBasis(f, 2, "atOrigin");
     3261timer - t;
     3262t = timer;
     3263list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3264timer - t;
     3265LIB"integralbasis.lib";
     3266ring r = 0,(x,y),dp;
     3267int k = 400; int d = 500;
     3268poly f = x*(x^(k-1) + y^2) + y^d;
     3269int t = timer;
     3270list l1 = integralBasis(f, 2, "atOrigin");
     3271timer - t;
     3272t = timer;
     3273list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3274timer - t;
     3275LIB"integralbasis.lib";
     3276ring r = 0, (x,y), dp;
     3277int k = 5;
     3278int d = 10;
     3279poly f = polyDK(d, k, 1231);
     3280int t = timer;
     3281list l1 = integralBasis(f, 2, "atOrigin");
     3282timer - t;
     3283t = timer;
     3284list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3285timer - t;
     3286reduce(l1[1], groebner(l2[1]));
     3287reduce(l2[1], groebner(l1[1]));
     3288LIB"integralbasis.lib";
     3289ring r = 0, (x,y), dp;
     3290int k = 5;
     3291int d = 10;
     3292poly f = polyDK(d, k, 1231);
     3293int t = timer;
     3294list l1 = integralBasis(f, 2, "atOrigin");
     3295intExp;
     3296ib
     3297LIB"integralbasis.lib";
     3298ring r = 0, (x,y), dp;
     3299int k = 5;
     3300int d = 10;
     3301poly f = polyDK(d, k, 1231);
     3302int t = timer;
     3303list l1 = integralBasis(f, 2, "atOrigin");
     3304timer - t;
     3305t = timer;
     3306list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3307timer - t;
     3308reduce(l1[1], groebner(l2[1]));
     3309reduce(l2[1], groebner(l1[1]));
     3310LIB"integralbasis.lib";
     3311ring r = 0,(x,y),dp;
     3312int k = 400; int d = 500;
     3313poly f = x*(x^(k-1) + y^2) + y^d;
     3314int t = timer;
     3315list l1 = integralBasis(f, 2, "atOrigin");
     3316timer - t;
     3317t = timer;
     3318list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3319timer - t;
     3320reduce(l1[1], groebner(l2[1]));
     3321reduce(l2[1], groebner(l1[1]));
     3322LIB"integralbasis.lib";
     3323ring r = 0,(x,y),dp;
     3324int k = 400; int d = 500;
     3325poly f = y^2 + x^(k+1) + y^d;
     3326int t = timer;
     3327list l1 = integralBasis(f, 2, "atOrigin");
     3328timer - t;
     3329t = timer;
     3330list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3331timer - t;
     3332reduce(l1[1], groebner(l2[1]));
     3333reduce(l2[1], groebner(l1[1]));
     3334printlevel = 5;
     3335LIB"integralbasis.lib";
     3336ring r = 0, (x,y), dp;
     3337poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     3338int t = timer;
     3339list l1 = integralBasis(f, 2, "atOrigin");
     3340timer - t;
     3341t = timer;
     3342list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3343timer - t;
     3344reduce(l1[1], groebner(l2[1]));
     3345reduce(l2[1], groebner(l1[1]));
     3346printlevel = 5;
     3347LIB"integralbasis.lib";
     3348ring r = 0, (x,y), dp;
     3349poly f = (y7 + x4) * (y7 + y5x3 + x4) + y30;
     3350int t = timer;
     3351list l1 = integralBasis(f, 2, "atOrigin");
     3352timer - t;
     3353t = timer;
     3354list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3355timer - t;
     3356reduce(l1[1], groebner(l2[1]));
     3357reduce(l2[1], groebner(l1[1]));
     3358LIB"integralbasis.lib";
     3359ring r = 0,(x,y),dp;
     3360int k = 6;
     3361poly z = 2*x - y + 1;
     3362poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3363f = monic(f);
     3364int t = timer;
     3365list l1 = integralBasis(f, 2, "atOrigin");
     3366timer - t;
     3367t = timer;
     3368list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3369timer - t;
     3370reduce(l1[1], groebner(l2[1]));
     3371reduce(l2[1], groebner(l1[1]));
     3372reduce(l1[1], groebner(l2[1]));
     3373reduce(l2[1], groebner(l1[1]));
     3374printlevel = 4;
     3375LIB"integralbasis.lib";
     3376ring r = 0,(x,y),dp;
     3377int k = 11;
     3378poly z = x - 2*y + 1;
     3379poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3380f = monic(f);
     3381list l = integralBasis(f, 2, "nonModular");
     3382LIB"integralbasis.lib";
     3383ring r = 0,(x,y),dp;
     3384poly f = x^15-21*x^14+8*x^13*y-6*x^13-16*x^12*y+20*x^11*y^2-x^12+8*x^11*y-36*x^10*y^2+24*x^9*y^3+4*x^9*y^2-16*x^8*y^3+26*x^7*y^4-6*x^6*y^4+8*x^5*y^5+4*x^3*y^6-y^8;
     3385f = f *(-1);
     3386list l = integralBasis(f, 2, "nonModular");
     3387LIB"integralbasis.lib";
     3388ring r = 0,(x,y),dp;
     3389poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3390list l = integralBasis(f, 2, "nonModular");
     3391LIB"integralbasis.lib";
     3392ring r = 0,(x,y),dp;
     3393poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3394list l = integralBasis(f, 2, "nonModular");
     3395int t = timer;
     3396list l1 = integralBasis(f, 2, "atOrigin");
     3397timer - t;
     3398t = timer;
     3399list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3400timer - t;
     3401reduce(l1[1], groebner(l2[1]));
     3402reduce(l2[1], groebner(l1[1]));
     3403LIB"integralbasis.lib";
     3404ring r = 0,(x,y),dp;
     3405poly f = (y^5 + y^4*x^7 + 2*x^8) * (y^3 + 7*x^4) * (y^7 + 2*x^12) * (y^11 + 2*x^18) + y^30;
     3406int t = timer;
     3407list l1 = integralBasis(f, 2, "atOrigin");
     3408timer - t;
     3409t = timer;
     3410list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3411timer - t;
     3412LIB"integralbasis.lib";
     3413ring r = 0,(x,y),(dp, L(100000));
     3414poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52) + y^36;
     3415int t = timer;
     3416list l1 = integralBasis(f, 2);
     3417timer - t;
     3418t = timer;
     3419list l2 = integralBasis(f, 2, "noOpti");
     3420timer - t;
     3421LIB"integralbasis.lib";
     3422ring r = 0,(x,y),dp;
     3423poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     3424list l = integralBasis(f, 2, "atOrigin");
     3425int t = timer;
     3426list l1 = integralBasis(f, 2);
     3427timer - t;
     3428t = timer;
     3429list l2 = integralBasis(f, 2, "noOpti");
     3430timer - t;
     3431LIB"integralbasis.lib";
     3432ring r = 0,(x,y),dp;
     3433poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     3434list l = integralBasis(f, 2, "atOrigin");
     3435int t = timer;
     3436list l1 = integralBasis(f, 2);
     3437timer - t;
     3438t = timer;
     3439list l2 = integralBasis(f, 2, "noOpti");
     3440timer - t;
     3441LIB"integralbasis.lib";
     3442ring r = 0,(x,y),dp;
     3443poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     3444list l = integralBasis(f, 2, "atOrigin");
     3445LIB"integralbasis.lib";
     3446ring r = 0,(x,y),dp;
     3447int k = 11;
     3448poly z = x - 2*y + 1;
     3449poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3450f = monic(f);
     3451list l = integralBasis(f, 2, "nonModular");
     3452LIB"integralbasis.lib";
     3453ring r = 0,(x,y),dp;
     3454int k = 11;
     3455poly z = x - 2*y + 1;
     3456poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3457f = monic(f);
     3458list l = integralBasis(f, 2, "nonModular");
     3459LIB"integralbasis.lib";
     3460ring r = 0,(x,y),dp;
     3461int k = 11;
     3462poly z = x - 2*y + 1;
     3463poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3464f = monic(f);
     3465list l = integralBasis(f, 2, "nonModular");
     3466list l = integralBasis(f, 2, "nonModular");
     3467LIB"integralbasis.lib";
     3468ring r = 0,(x,y),dp;
     3469int k = 7;
     3470poly z = x - 2*y + 1;
     3471poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3472f = monic(f);
     3473list l = integralBasis(f, 2, "nonModular");
     3474int t = timer;
     3475list l1 = integralBasis(f, 2, "atOrigin");
     3476timer - t;
     3477t = timer;
     3478list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3479timer - t;
     3480reduce(l1[1], groebner(l2[1]));
     3481reduce(l2[1], groebner(l1[1]));
     3482printlevel = 4;
     3483LIB"integralbasis.lib";
     3484ring r = 0,(x,y),dp;
     3485int k = 7;
     3486poly z = x - 2*y + 1;
     3487poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3488f = monic(f);
     3489int t = timer;
     3490list l1 = integralBasis(f, 2, "atOrigin");
     3491timer - t;
     3492t = timer;
     3493list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3494timer - t;
     3495reduce(l1[1], groebner(l2[1]));
     3496reduce(l2[1], groebner(l1[1]));
     3497LIB"integralbasis.lib";
     3498ring r = 0,(x,y),dp;
     3499int k = 7;
     3500poly z = x - 2*y + 1;
     3501poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3502f = monic(f);
     3503int t = timer;
     3504list l1 = integralBasis(f, 2);
     3505timer - t;
     3506t = timer;
     3507list l2 = integralBasis(f, 2, "noOpti");
     3508timer - t;
     3509reduce(l1[1], groebner(l2[1]));
     3510reduce(l2[1], groebner(l1[1]));
     3511LIB"integralbasis.lib";
     3512ring r = 0,(x,y),dp;
     3513int k = 7;
     3514poly z = x - 2*y + 1;
     3515poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3516f = monic(f);
     3517int t = timer;
     3518list l1 = integralBasis(f, 2);
     3519timer - t;
     3520t = timer;
     3521list l2 = integralBasis(f, 2, "noOpti");
     3522timer - t;
     3523l1;
     3524l2;
     3525l1;
     3526l2;
     3527reduce(l1[1], groebner(l2[1]));
     3528LIB"integralbasis.lib";
     3529ring r = 0,(x,y),dp;
     3530int k = 7;
     3531poly z = x - 2*y + 1;
     3532poly f = x^(2*k) + y^(2*k) + z^(2*k) + 2 * (x^k * z^k - x^k * y^k + y^k * z^k);
     3533f = monic(f);
     3534int t = timer;
     3535list l1 = integralBasis(f, 2);
     3536timer - t;
     3537t = timer;
     3538list l2 = integralBasis(f, 2, "noOpti");
     3539timer - t;
     3540reduce(l2[1], groebner(l1[1]));
     3541LIB"integralbasis.lib";
     3542ring r = 0,(x,y),dp;
     3543int k = 4;
     3544poly z = 2*x - y + 1;
     3545poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3546f = monic(f);
     3547int t = timer;
     3548list l1 = integralBasis(f, 2, "atOrigin");
     3549timer - t;
     3550t = timer;
     3551list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3552timer - t;
     3553reduce(l1[1], groebner(l2[1]));
     3554reduce(l2[1], groebner(l1[1]));
     3555LIB"integralbasis.lib";
     3556ring r = 0,(x,y),dp;
     3557int k = 4;
     3558poly z = 2*x - y + 1;
     3559poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3560f = monic(f);
     3561int t = timer;
     3562list l1 = integralBasis(f, 2, "atOrigin");
     3563LIB"integralbasis.lib";
     3564ring r = 0,(x,y),dp;
     3565int k = 4;
     3566poly z = 2*x - y + 1;
     3567poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3568f = monic(f);
     3569int t = timer;
     3570list l1 = integralBasis(f, 2);
     3571LIB"integralbasis.lib";
     3572ring r = 0,(x,y),dp;
     3573int k = 4;
     3574poly z = 2*x - y + 1;
     3575poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3576f = monic(f);
     3577int t = timer;
     3578list l1 = integralBasis(f, 2);
     3579timer - t;
     3580t = timer;
     3581list l2 = integralBasis(f, 2, "noOpti");
     3582timer - t;
     3583reduce(l1[1], groebner(l2[1]));
     3584reduce(l2[1], groebner(l1[1]));
     3585l1;
     3586l2;
     3587l1;
     3588l2;
     3589LIB"integralbasis.lib";
     3590ring r = 0,(x,y),dp;
     3591int k = 4;
     3592poly z = 2*x - y + 1;
     3593poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3594f = monic(f);
     3595int t = timer;
     3596list l1 = integralBasis(f, 2);
     3597basering
     3598basering
     3599LIB"integralbasis.lib";
     3600ring r = 0,(x,y),dp;
     3601int k = 4;
     3602poly z = 2*x - y + 1;
     3603poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3604f = monic(f);
     3605int t = timer;
     3606list l1 = integralBasis(f, 2);
     3607minpoly
     3608pardeg(minpoly)
     3609pardeg(minpoly)
     3610basering
     3611LIB"integralbasis.lib";
     3612ring r = 0,(x,y),dp;
     3613int k = 4;
     3614poly z = 2*x - y + 1;
     3615poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3616f = monic(f);
     3617int t = timer;
     3618list l1 = integralBasis(f, 2);
     3619timer - t;
     3620t = timer;
     3621list l2 = integralBasis(f, 2, "noOpti");
     3622timer - t;
     3623reduce(l1[1], groebner(l2[1]));
     3624reduce(l2[1], groebner(l1[1]));
     3625LIB"integralbasis.lib";
     3626ring r = 0,(x,y),dp;
     3627int k = 4;
     3628poly z = 2*x - y + 1;
     3629poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3630f = monic(f);
     3631int t = timer;
     3632list l1 = integralBasis(f, 2);
     3633timer - t;
     3634t = timer;
     3635list l2 = integralBasis(f, 2, "noOpti");
     3636timer - t;
     3637reduce(l1[1], groebner(l2[1]));
     3638reduce(l2[1], groebner(l1[1]));
     3639LIB"integralbasis.lib";
     3640ring r = 0,(x,y),dp;
     3641int k = 4;
     3642poly z = 2*x - y + 1;
     3643poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3644f = monic(f);
     3645int t = timer;
     3646list l1 = integralBasis(f, 2);
     3647timer - t;
     3648t = timer;
     3649list l2 = integralBasis(f, 2, "noOpti");
     3650timer - t;
     3651reduce(l1[1], groebner(l2[1]));
     3652reduce(l2[1], groebner(l1[1]));
     3653LIB"integralbasis.lib";
     3654ring r = 0,(x,y),dp;
     3655int k = 6;
     3656poly z = 2*x - y + 1;
     3657poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3658f = monic(f);
     3659int t = timer;
     3660list l1 = integralBasis(f, 2);
     3661timer - t;
     3662t = timer;
     3663list l2 = integralBasis(f, 2, "noOpti");
     3664timer - t;
     3665reduce(l1[1], groebner(l2[1]));
     3666reduce(l2[1], groebner(l1[1]));
     3667LIB"integralbasis.lib";
     3668ring r = 0,(x,y),dp;
     3669int k = 6;
     3670poly z = 2*x - y + 1;
     3671poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3672f = monic(f);
     3673int t = timer;
     3674list l1 = integralBasis(f, 2);
     3675timer - t;
     3676t = timer;
     3677list l2 = integralBasis(f, 2, "noOpti");
     3678timer - t;
     3679reduce(l1[1], groebner(l2[1]));
     3680reduce(l2[1], groebner(l1[1]));
     3681reduce(l1[1], groebner(l2[1]));
     3682reduce(l2[1], groebner(l1[1]));
     3683LIB"integralbasis.lib";
     3684ring r = 0,(x,y),dp;
     3685int k = 6;
     3686poly z = 2*x - y + 1;
     3687poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3688f = monic(f);
     3689int t = timer;
     3690list l1 = integralBasis(f, 2);
     3691timer - t;
     3692t = timer;
     3693list l2 = integralBasis(f, 2, "noOpti");
     3694timer - t;
     3695basering
     3696PE
     3697minPolys
     3698composePolys(minPolys)
     3699poly mpX = subst(mp, var(2), var(1));
     3700LIB"integralbasis.lib";
     3701ring r = 0,(x,y),dp;
     3702int k = 6;
     3703poly z = 2*x - y + 1;
     3704poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3705f = monic(f);
     3706int t = timer;
     3707list l1 = integralBasis(f, 2);
     3708timer - t;
     3709t = timer;
     3710list l2 = integralBasis(f, 2, "noOpti");
     3711timer - t;
     3712mp
     3713LIB"integralbasis.lib";
     3714ring r = 0,(x,y),dp;
     3715int k = 6;
     3716poly z = 2*x - y + 1;
     3717poly f = (x^(k+1) + y^(k+1) + z^(k+1))^2 - 4*(x^(k+1)*y^(k+1) + y^(k+1)*z^(k+1) + z^(k+1)*x^(k+1));
     3718f = monic(f);
     3719int t = timer;
     3720list l1 = integralBasis(f, 2);
     3721timer - t;
     3722t = timer;
     3723list l2 = integralBasis(f, 2, "noOpti");
     3724timer - t;
     3725            poly mp = composePolys(minPolys);
     3726            poly mpX = subst(mp, var(2), var(1));
     3727fF
     3728            poly fFB = buildPolyFracNew(fF, mpX);
     3729fFB
     3730PE
     3731d
     3732LIB"integralbasis.lib";
     3733ring r = 0,(x,y),(dp, L(100000));
     3734poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52) + y^100;
     3735int t = timer;
     3736list l1 = integralBasis(f, 2, "atOrigin");
     3737timer - t;
     3738t = timer;
     3739list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3740timer - t;
     3741reduce(l1[1], groebner(l2[1]));
     3742reduce(l2[1], groebner(l1[1]));
     3743size(reduce(l1[1], groebner(l2[1])));
     3744size(reduce(l2[1], groebner(l1[1])));
     3745LIB"integralbasis.lib";
     3746ring r = 0,(x,y),(dp, L(100000));
     3747poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52)*(y^7+3*x^2) + y^100;
     3748int t = timer;
     3749list l1 = integralBasis(f, 2, "atOrigin");
     3750timer - t;
     3751t = timer;
     3752list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3753timer - t;
     3754size(reduce(l1[1], groebner(l2[1])));
     3755size(reduce(l2[1], groebner(l1[1])));
     3756LIB"integralbasis.lib";
     3757ring r = 0, (x,y), dp;
     3758poly f = (y7 + x4) * (y7 + y5x3 + x4) + y30;
     3759int t = timer;
     3760list l1 = integralBasis(f, 2, "atOrigin");
     3761timer - t;
     3762t = timer;
     3763list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3764timer - t;
     3765size(reduce(l1[1], groebner(l2[1])));
     3766size(reduce(l2[1], groebner(l1[1])));
     3767LIB"integralbasis.lib";
     3768ring r = 0, (x,y), dp;
     3769poly f = (y7 + x4) * (y7 + y5x3 + x4) + y30;
     3770int t = timer;
     3771list l1 = integralBasis(f, 2, "atOrigin");
     3772timer - t;
     3773t = timer;
     3774list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3775timer - t;
     3776LIB"integralbasis.lib";
     3777ring r = 0, (x,y), dp;
     3778poly f = (y7 + x4) * (y7 + y5x3 + x4) + y30;
     3779int t = timer;
     3780list l1 = integralBasis(f, 2, "atOrigin");
     3781timer - t;
     3782t = timer;
     3783list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3784timer - t;
     3785size(reduce(l1[1], groebner(l2[1])));
     3786size(reduce(l2[1], groebner(l1[1])));
     3787LIB"integralbasis.lib";
     3788ring r = 0, (x,y), dp;
     3789poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     3790int t = timer;
     3791list l1 = integralBasis(f, 2, "atOrigin");
     3792timer - t;
     3793t = timer;
     3794list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3795timer - t;
     3796size(reduce(l1[1], groebner(l2[1])));
     3797size(reduce(l2[1], groebner(l1[1])));
     3798LIB"integralbasis.lib";
     3799ring r = 0, (x,y), dp;
     3800poly f = (y7 + x4) * (y7 + y5x3 + x4) + y16;
     3801int t = timer;
     3802list l1 = integralBasis(f, 2, "atOrigin");
     3803timer - t;
     3804t = timer;
     3805list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3806timer - t;
     3807size(reduce(l1[1], groebner(l2[1])));
     3808size(reduce(l2[1], groebner(l1[1])));
     3809LIB"integralbasis.lib";
     3810ring r = 0,(x,y),dp;
     3811int k = 400; int d = 500;
     3812poly f = y^2 + x^(k+1) + y^d;
     3813int t = timer;
     3814list l1 = integralBasis(f, 2, "atOrigin");
     3815timer - t;
     3816t = timer;
     3817list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3818timer - t;
     3819size(reduce(l1[1], groebner(l2[1])));
     3820size(reduce(l2[1], groebner(l1[1])));
     3821LIB"integralbasis.lib";
     3822ring r = 0,(x,y),dp;
     3823int k = 400; int d = 500;
     3824poly f = x*(x^(k-1) + y^2) + y^d;
     3825int t = timer;
     3826list l1 = integralBasis(f, 2, "atOrigin");
     3827timer - t;
     3828t = timer;
     3829list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3830timer - t;
     3831size(reduce(l1[1], groebner(l2[1])));
     3832size(reduce(l2[1], groebner(l1[1])));
     3833LIB"integralbasis.lib";
     3834ring r = 0,(x,y),dp;
     3835int k = 400; int d = 500;
     3836poly f = x*(x^(k-1) + y^2) + y^d;
     3837int t = timer;
     3838list l1 = integralBasis(f, 2, "atOrigin");
     3839timer - t;
     3840t = timer;
     3841list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3842timer - t;
     3843size(reduce(l1[1], groebner(l2[1])));
     3844size(reduce(l2[1], groebner(l1[1])));
     3845LIB"integralbasis.lib";
     3846ring r = 0, (x,y), dp;
     3847int k = 5;
     3848int d = 10;
     3849poly f = polyDK(d, k, 1231);
     3850int t = timer;
     3851list l1 = integralBasis(f, 2, "atOrigin");
     3852timer - t;
     3853t = timer;
     3854list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3855timer - t;
     3856size(reduce(l1[1], groebner(l2[1])));
     3857size(reduce(l2[1], groebner(l1[1])));
     3858LIB"integralbasis.lib";
     3859ring r = 0, (x,y), dp;
     3860int k = 20;
     3861int d = 30;
     3862poly f = polyDK(d, k, 1231);
     3863int t = timer;
     3864list l1 = integralBasis(f, 2, "atOrigin");
     3865timer - t;
     3866t = timer;
     3867list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3868timer - t;
     3869size(reduce(l1[1], groebner(l2[1])));
     3870size(reduce(l2[1], groebner(l1[1])));
     3871LIB"integralbasis.lib";
     3872ring r = 0, (x,y), dp;
     3873int k = 50;
     3874int d = 60;
     3875poly f = polyDK(d, k, 1231);
     3876int t = timer;
     3877list l1 = integralBasis(f, 2, "atOrigin");
     3878timer - t;
     3879t = timer;
     3880list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3881timer - t;
     3882size(reduce(l1[1], groebner(l2[1])));
     3883size(reduce(l2[1], groebner(l1[1])));
     3884ring r = 0, (x,y), dp;
     3885int k = 50;
     3886int d = 70;
     3887poly f = polyDK(d, k, 1231);
     3888int t = timer;
     3889list l1 = integralBasis(f, 2, "atOrigin");
     3890timer - t;
     3891t = timer;
     3892list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3893timer - t;
     3894size(reduce(l1[1], groebner(l2[1])));
     3895size(reduce(l2[1], groebner(l1[1])));
     3896ring r = 0, (x,y), dp;
     3897int k = 50;
     3898int d = 70;
     3899poly f = polyDK(d, k, 1231);
     3900int t = timer;
     3901list l1 = integralBasis(f, 2, "atOrigin");
     3902timer - t;
     3903t = timer;
     3904list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3905timer - t;
     3906size(reduce(l1[1], groebner(l2[1])));
     3907size(reduce(l2[1], groebner(l1[1])));
     3908LIB"integralbasis.lib";
     3909ring r = 0, (x,y), dp;
     3910int k = 50;
     3911int d = 70;
     3912poly f = polyDK(d, k, 1231);
     3913int t = timer;
     3914list l1 = integralBasis(f, 2, "atOrigin");
     3915timer - t;
     3916t = timer;
     3917list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3918timer - t;
     3919size(reduce(l1[1], groebner(l2[1])));
     3920size(reduce(l2[1], groebner(l1[1])));
     3921LIB"integralbasis.lib";
     3922ring r = 0, (x,y), dp;
     3923int k = 50;
     3924int d = 70;
     3925poly f = polyDK(d, k, 1231);
     3926int t = timer;
     3927list l1 = integralBasis(f, 2, "atOrigin");
     3928timer - t;
     3929t = timer;
     3930list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3931timer - t;
     3932size(reduce(l1[1], groebner(l2[1])));
     3933size(reduce(l2[1], groebner(l1[1])));
     3934LIB"integralbasis.lib";
     3935ring r = 0,(x,y),dp;
     3936poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3937int t = timer;
     3938list l1 = integralBasis(f, 2, "atOrigin");
     3939timer - t;
     3940t = timer;
     3941list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3942timer - t;
     3943reduce(l1[1], groebner(l2[1]));
     3944reduce(l2[1], groebner(l1[1]));
     3945LIB"integralbasis.lib";
     3946ring r = 0,(x,y),dp;
     3947poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3948int t = timer;
     3949list l1 = integralBasis(f, 2, "atOrigin");
     3950timer - t;
     3951t = timer;
     3952list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3953timer - t;
     3954reduce(l1[1], groebner(l2[1]));
     3955reduce(l2[1], groebner(l1[1]));
     3956LIB"integralbasis.lib";
     3957ring r = 0,(x,y),dp;
     3958poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3959int t = timer;
     3960list l1 = integralBasis(f, 2, "atOrigin");
     3961timer - t;
     3962t = timer;
     3963list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3964timer - t;
     3965size(reduce(l1[1], groebner(l2[1])));
     3966size(reduce(l2[1], groebner(l1[1])));
     3967LIB"integralbasis.lib";
     3968ring r = 0,(x,y),dp;
     3969poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3970int t = timer;
     3971list l1 = integralBasis(f, 2, "atOrigin");
     3972timer - t;
     3973t = timer;
     3974list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3975timer - t;
     3976size(reduce(l1[1], groebner(l2[1])));
     3977size(reduce(l2[1], groebner(l1[1])));
     3978LIB"integralbasis.lib";
     3979ring r = 0,(x,y),dp;
     3980poly f = (y^4 + 2*x^3*y^2 + x^6 + x^5*y)^3 + x^11*y^11;
     3981int t = timer;
     3982list l1 = integralBasis(f, 2, "atOrigin");
     3983timer - t;
     3984t = timer;
     3985list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3986timer - t;
     3987size(reduce(l1[1], groebner(l2[1])));
     3988size(reduce(l2[1], groebner(l1[1])));
     3989LIB"integralbasis.lib";
     3990ring r = 0,(x,y),dp;
     3991poly f = (y^5 + y^4*x^7 + 2*x^8) * (y^3 + 7*x^4) * (y^7 + 2*x^12) * (y^11 + 2*x^18) + y^30;
     3992int t = timer;
     3993list l1 = integralBasis(f, 2, "atOrigin");
     3994timer - t;
     3995t = timer;
     3996list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     3997timer - t;
     3998size(reduce(l1[1], groebner(l2[1])));
     3999size(reduce(l2[1], groebner(l1[1])));
     4000LIB"integralbasis.lib";
     4001ring r = 0,(x,y),(dp, L(100000));
     4002poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52) + y^36;
     4003int t = timer;
     4004list l1 = integralBasis(f, 2);
     4005timer - t;
     4006t = timer;
     4007list l2 = integralBasis(f, 2, "noOpti");
     4008timer - t;
     4009LIB"integralbasis.lib";
     4010ring r = 0,(x,y),(dp, L(100000));
     4011poly f = (y^15 + 2*x^38)*(y^19 + 7*x^52) + y^36;
     4012int t = timer;
     4013list l1 = integralBasis(f, 2);
     4014timer - t;
     4015t = timer;
     4016list l2 = integralBasis(f, 2, "noOpti");
     4017timer - t;
     4018size(reduce(l1[1], groebner(l2[1])));
     4019size(reduce(l2[1], groebner(l1[1])));
     4020LIB"integralbasis.lib";
     4021ring r = 0,(x,y),dp;
     4022poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4023list l = integralBasis(f, 2, "atOrigin");
     4024int t = timer;
     4025list l1 = integralBasis(f, 2);
     4026timer - t;
     4027t = timer;
     4028list l2 = integralBasis(f, 2, "noOpti");
     4029timer - t;
     4030size(reduce(l1[1], groebner(l2[1])));
     4031size(reduce(l2[1], groebner(l1[1])));
     4032printlevel = 4;
     4033LIB"integralbasis.lib";
     4034ring r = 0,(x,y),dp;
     4035poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4036list l = integralBasis(f, 2, "atOrigin");
     4037int t = timer;
     4038list l1 = integralBasis(f, 2);
     4039timer - t;
     4040t = timer;
     4041list l2 = integralBasis(f, 2, "noOpti");
     4042timer - t;
     4043printlevel = 4;
     4044LIB"integralbasis.lib";
     4045ring r = 0,(x,y),dp;
     4046poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4047list l = integralBasis(f, 2, "atOrigin");
     4048int t = timer;
     4049list l1 = integralBasis(f, 2);
     4050timer - t;
     4051t = timer;
     4052list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     4053timer - t;
     4054size(reduce(l1[1], groebner(l2[1])));
     4055size(reduce(l2[1], groebner(l1[1])));
     4056printlevel = 4;
     4057LIB"integralbasis.lib";
     4058ring r = 0,(x,y),dp;
     4059poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4060int t = timer;
     4061list l1 = integralBasis(f, 2);
     4062timer - t;
     4063t = timer;
     4064list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     4065timer - t;
     4066LIB"integralbasis.lib";
     4067ring r = 0,(x,y),dp;
     4068poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4069int t = timer;
     4070list l1 = integralBasis(f, 2, "atOrigin");
     4071timer - t;
     4072t = timer;
     4073list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     4074timer - t;
     4075size(reduce(l1[1], groebner(l2[1])));
     4076size(reduce(l2[1], groebner(l1[1])));
     4077printlevel = 4;
     4078LIB"integralbasis.lib";
     4079ring r = 0,(x,y),dp;
     4080poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4081int t = timer;
     4082list l1 = integralBasis(f, 2, "atOrigin");
     4083timer - t;
     4084t = timer;
     4085list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     4086timer - t;
     4087"Verification...";
     4088size(reduce(l1[1], groebner(l2[1])));
     4089size(reduce(l2[1], groebner(l1[1])));
     4090LIB"integralbasis.lib";
     4091ring r = 0,(x,y),dp;
     4092poly f = (y^35 + y^34*x^7 + 2*x^38) * (y^33 + 7*x^44) * (y^37 + 2*x^52) + y^110;
     4093int t = timer;
     4094list l1 = integralBasis(f, 2, "atOrigin");
     4095timer - t;
     4096t = timer;
     4097list l2 = integralBasis(f, 2, "atOrigin", "noOpti");
     4098timer - t;
     4099"Verification...";
     4100size(reduce(l1[1], groebner(l2[1])));
     4101size(reduce(l2[1], groebner(l1[1])));
     4102ring r = 0, (x,y,z), dp;
     4103ideal I = x5, y5, z5, xz4, xyz3, xy3z, x3yz;
     4104hilb(groebner(I));
     4105I;
     4106hilb(I);
     4107ideal J = x5, y5, z5;
     4108hilb(groebner(J));
     4109ideal I = x3,y3,z3;
     4110hilb(groebner(I));
     4111ring R=32003,(x,y,z),dp;
     4112  ideal i=x2,y2,z2;
     4113  ideal s=std(i);
     4114  hilb(s);
     4115ring R=32003,(x,y,z),dp;
     4116  ideal i=x2,y2,z2;
     4117  ideal s=std(i);
     4118  hilb(s,1);
     4119ring R=32003,(x,y,z),dp;
     4120  ideal i=x2,y2,z2;
     4121  ideal s=std(i);
     4122  hilb(s,2);
     4123hilb(groebner(I));
     4124ring r = 0, (x,y,z), dp;
     4125ideal I = x5, y5, z5, xz4, xyz3, xy3z, x3yz;
     4126hilb(groebner(I));
     4127def H = hilb(groebner(I));
     4128ring r = 0, (x,y,z), dp;
     4129ideal I = x5, y5, z5, xz4, xyz3, xy3z, x3yz;
     4130hilb(groebner(I),2);
     4131ring r = 0, (x,y,z), dp;
     4132ideal I = x5, y5, z5;
     4133hilb(groebner(I),2);
     4134ideal I = x5, y5, z5, xz4, xyz3, xy3z, x3yz;
     4135hilb(groebner(I),2);
     4136ideal I = x5, y5, z5, x*z^4, y^2*z^3, x*y*z^3, x^3*y*z;
     4137hilb(groebner(I),2);
     4138ideal I = x5, y5, z5, x*z^4, x^4*y, y^4*z, x^2*y*z^2;
     4139hilb(groebner(I),2);
     4140ideal I = x5, y5, z5, x^2*y^3, x^3*z^2, x*y*z^3, x^2*y*z^2;
     4141hilb(groebner(I),2);
     4142ideal I = x5, y5, z5, x^2*y^3, x^3*z^2, x^2*y*z^2, x^2*y*z^2;
     4143hilb(groebner(I),2);
     4144ideal I = x5, y5, z5, x^2*y^3, x^3*z^2, x^3*y*z, x^2*y*z^2;
     4145hilb(groebner(I),2);
     4146LIB"integralbasis.lib";
     4147ring r = 0,(x,y),dp;
     4148poly f = x^(30)+x^2*y^2+y^(11);
     4149map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4150f=phi(f);
     4151int t = timer;
     4152list l1 = integralBasis(f, 2, "atOrigin");
     4153timer - t;
     4154LIB"integralbasis.lib";
     4155ring R = 0,(x,y),dp;
     4156poly f = x^(30)+x^2*y^2+y^(11);
     4157map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4158f=phi(f);
     4159int t = timer;
     4160list l1 = integralBasis(f, 2, "atOrigin");
     4161timer - t;
     4162poly f = x^(30)+x^2*y^2+y^(11);
     4163int t = timer;
     4164int t = timer;
     4165list l1 = integralBasis(f, 2, "atOrigin");
     4166timer - t;
     4167printlevel = 5;
     4168LIB"integralbasis.lib";
     4169ring R = 0,(x,y),dp;
     4170poly f = x^(30)+x^2*y^2+y^(11);
     4171map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4172poly g =phi(f);
     4173int t = timer;
     4174list l1 = integralBasis(f, 2, "atOrigin");
     4175timer - t;
     4176int t = timer;
     4177list l1 = integralBasis(g, 2, "atOrigin");
     4178timer - t;
     4179g;
     4180printlevel = 8;
     4181LIB"integralbasis.lib";
     4182ring R = 0,(x,y),dp;
     4183poly f = x^(30)+x^2*y^2+y^(11);
     4184map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4185poly g =phi(f);
     4186int t = timer;
     4187list l1 = integralBasis(g, 2, "atOrigin");
     4188timer - t;
     4189[200~printlevel = 8;
     4190LIB"integralbasis.lib";
     4191ring R = 0,(x,y),dp;
     4192poly f = x^(30)+x^2*y^2+y^(11);
     4193map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4194poly g =phi(f);
     4195printlevel = 8;
     4196LIB"integralbasis.lib";
     4197ring R = 0,(x,y),dp;
     4198poly f = x^(30)+x^2*y^2+y^(11);
     4199map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4200poly g =phi(f);
     4201puisex(f,1,-1);
     4202puiseux(f,1,-1);
     4203puiseux(g,1,-1);
     4204list l = puiseux(g,1,-1);
     4205def S = l[1];
     4206setring S;
     4207PE;
     4208basering
     4209;
     4210ring R = 0,(x,y),dp;
     4211poly f = x^(30)+x^2*y^2+y^(11);
     4212map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4213poly g =phi(f);
     4214g;
     4215f;
     4216puiseux(g, 1, -1);
     4217LIB"integralbasis.lib";
     4218ring R = 0,(x,y),dp;
     4219poly f = x^(30)+x^2*y^2+y^(11);
     4220map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4221poly g =phi(f);
     4222puiseux(g, 1, -1);
     4223puiseux(f, 1, -1);
     4224list l1 = puiseux(g, 1, -1);
     4225def S = l1[1];
     4226setring S;
     4227PE;
     4228setring r;
     4229setring R;
     4230l;
     4231l1;
     4232g;
     4233LIB"integralbasis.lib";
     4234ring R = 0,(x,y),dp;
     4235poly h = x31+x30y+x29y2+x28y3+x27y4+x26y5+x25y6+x24y7+x23y8+x22y9+x21y10+x20y11+x19y12+x18y13+x17y14+x16y15+x15y16+x14y17+x13y18+x12y19+x11y20+x10y21+x9y22+x8y23+x7y24+x6y25+x5y26+x4y27+x3y28+x2y29+xy30+y31+x30+x29y+x28y2+x27y3+x26y4+x25y5+x24y6+x23y7+x22y8+x21y9+x20y10+x19y11+x18y12+x17y13+x16y14+x15y15+x14y16+x13y17+x12y18+x11y19+x10y20+x9y21+x8y22+x7y23+x6y24+x5y25+x4y26+x3y27+x2y28+xy29+y30+x29+x28y+x27y2+x26y3+x25y4+x24y5+x23y6+x22y7+x21y8+x20y9+x19y10+x18y11+x17y12+x16y13+x15y14+x14y15+x13y16+x12y17+x11y18+x10y19+x9y20+x8y21+x7y22+x6y23+x5y24+x4y25+x3y26+x2y27+xy28+x28+x27y+x26y2+x25y3+x24y4+x23y5+x22y6+x21y7+x20y8+x19y9+x18y10+x17y11+x16y12+x15y13+x14y14+x13y15+x12y16+x11y17+x10y18+x9y19+x8y20+x7y21+x6y22+x5y23+x4y24+x3y25+x2y26+xy27+x27+x26y+x25y2+x24y3+x23y4+x22y5+x21y6+x20y7+x19y8+x18y9+x17y10+x16y11+x15y12+x14y13+x13y14+x12y15+x11y16+x10y17+x9y18+x8y19+x7y20+x6y21+x5y22+x4y23+x3y24+x2y25+xy26+x26+x25y+x24y2+x23y3+x22y4+x21y5+x20y6+x19y7+x18y8+x17y9+x16y10+x15y11+x14y12+x13y13+x12y14+x11y15+x10y16+x9y17+x8y18+x7y19+x6y20+x5y21+x4y22+x3y23+x2y24+xy25+x25+x24y+x23y2+x22y3+x21y4+x20y5+x19y6+x18y7+x17y8+x16y9+x15y10+x14y11+x13y12+x12y13+x11y14+x10y15+x9y16+x8y17+x7y18+x6y19+x5y20+x4y21+x3y22+x2y23+xy24+x24+x23y+x22y2+x21y3+x20y4+x19y5+x18y6+x17y7+x16y8+x15y9+x14y10+x13y11+x12y12+x11y13+x10y14+x9y15+x8y16+x7y17+x6y18+x5y19+x4y20+x3y21+x2y22+xy23+x23+x22y+x21y2+x20y3+x19y4+x18y5+x17y6+x16y7+x15y8+x14y9+x13y10+x12y11+x11y12+x10y13+x9y14+x8y15+x7y16+x6y17+x5y18+x4y19+x3y20+x2y21+xy22+x22+x21y+x20y2+x19y3+x18y4+x17y5+x16y6+x15y7+x14y8+x13y9+x12y10+x11y11+x10y12+x9y13+x8y14+x7y15+x6y16+x5y17+x4y18+x3y19+x2y20+xy21+x21+x20y+x19y2+x18y3+x17y4+x16y5+x15y6+x14y7+x13y8+x12y9+x11y10+x10y11+x9y12+x8y13+x7y14+x6y15+x5y16+x4y17+x3y18+x2y19+xy20+x20+x19y+x18y2+x17y3+x16y4+x15y5+x14y6+x13y7+x12y8+x11y9+x10y10+x9y11+x8y12+x7y13+x6y14+x5y15+x4y16+x3y17+x2y18+xy19+x19+x18y+x17y2+x16y3+x15y4+x14y5+x13y6+x12y7+x11y8+x10y9+x9y10+x8y11+x7y12+x6y13+x5y14+x4y15+x3y16+x2y17+xy18+x18+x17y+x16y2+x15y3+x14y4+x13y5+x12y6+x11y7+x10y8+x9y9+x8y10+x7y11+x6y12+x5y13+x4y14+x3y15+x2y16+xy17+x17+x16y+x15y2+x14y3+x13y4+x12y5+x11y6+x10y7+x9y8+x8y9+x7y10+x6y11+x5y12+x4y13+x3y14+x2y15+xy16+x16+x15y+x14y2+x13y3+x12y4+x11y5+x10y6+x9y7+x8y8+x7y9+x6y10+x5y11+x4y12+x3y13+x2y14+x15+x14y+x13y2+x12y3+x11y4+x10y5+x9y6+x8y7+x7y8+x6y9+x5y10+x4y11+x3y12+x2y13+x14+x13y+x12y2+x11y3+x10y4+x9y5+x8y6+x7y7+x6y8+x5y9+x4y10+x3y11+x2y12+x13+x12y+x11y2+x10y3+x9y4+x8y5+x7y6+x6y7+x5y8+x4y9+x3y10+x2y11+x12+x11y+x10y2+x9y3+x8y4+x7y5+x6y6+x5y7+x4y8+x3y9+x2y10+x11+x10y+x9y2+x8y3+x7y4+x6y5+x5y6+x4y7+x3y8+x2y9+x9y+x8y2+x7y3+x6y4+x5y5+x4y6+x3y7+x2y8+x8y+x7y2+x6y3+x5y4+x4y5+x3y6+x2y7+x7y+x6y2+x5y3+x4y4+x3y5+x2y6+x5y2+x4y3+x3y4+x2y5+x4y2+x3y3+x2y4+x3y2+x2y3+x2y2;
     4236poly h2 = phi(h);
     4237int t = timer;
     4238list l1 = integralBasis(h2, 2, "atOrigin");
     4239timer - t;
     4240LIB"integralbasis.lib";
     4241ring R = 0,(x,y),dp;
     4242poly h = x31+x30y+x29y2+x28y3+x27y4+x26y5+x25y6+x24y7+x23y8+x22y9+x21y10+x20y11+x19y12+x18y13+x17y14+x16y15+x15y16+x14y17+x13y18+x12y19+x11y20+x10y21+x9y22+x8y23+x7y24+x6y25+x5y26+x4y27+x3y28+x2y29+xy30+y31+x30+x29y+x28y2+x27y3+x26y4+x25y5+x24y6+x23y7+x22y8+x21y9+x20y10+x19y11+x18y12+x17y13+x16y14+x15y15+x14y16+x13y17+x12y18+x11y19+x10y20+x9y21+x8y22+x7y23+x6y24+x5y25+x4y26+x3y27+x2y28+xy29+y30+x29+x28y+x27y2+x26y3+x25y4+x24y5+x23y6+x22y7+x21y8+x20y9+x19y10+x18y11+x17y12+x16y13+x15y14+x14y15+x13y16+x12y17+x11y18+x10y19+x9y20+x8y21+x7y22+x6y23+x5y24+x4y25+x3y26+x2y27+xy28+x28+x27y+x26y2+x25y3+x24y4+x23y5+x22y6+x21y7+x20y8+x19y9+x18y10+x17y11+x16y12+x15y13+x14y14+x13y15+x12y16+x11y17+x10y18+x9y19+x8y20+x7y21+x6y22+x5y23+x4y24+x3y25+x2y26+xy27+x27+x26y+x25y2+x24y3+x23y4+x22y5+x21y6+x20y7+x19y8+x18y9+x17y10+x16y11+x15y12+x14y13+x13y14+x12y15+x11y16+x10y17+x9y18+x8y19+x7y20+x6y21+x5y22+x4y23+x3y24+x2y25+xy26+x26+x25y+x24y2+x23y3+x22y4+x21y5+x20y6+x19y7+x18y8+x17y9+x16y10+x15y11+x14y12+x13y13+x12y14+x11y15+x10y16+x9y17+x8y18+x7y19+x6y20+x5y21+x4y22+x3y23+x2y24+xy25+x25+x24y+x23y2+x22y3+x21y4+x20y5+x19y6+x18y7+x17y8+x16y9+x15y10+x14y11+x13y12+x12y13+x11y14+x10y15+x9y16+x8y17+x7y18+x6y19+x5y20+x4y21+x3y22+x2y23+xy24+x24+x23y+x22y2+x21y3+x20y4+x19y5+x18y6+x17y7+x16y8+x15y9+x14y10+x13y11+x12y12+x11y13+x10y14+x9y15+x8y16+x7y17+x6y18+x5y19+x4y20+x3y21+x2y22+xy23+x23+x22y+x21y2+x20y3+x19y4+x18y5+x17y6+x16y7+x15y8+x14y9+x13y10+x12y11+x11y12+x10y13+x9y14+x8y15+x7y16+x6y17+x5y18+x4y19+x3y20+x2y21+xy22+x22+x21y+x20y2+x19y3+x18y4+x17y5+x16y6+x15y7+x14y8+x13y9+x12y10+x11y11+x10y12+x9y13+x8y14+x7y15+x6y16+x5y17+x4y18+x3y19+x2y20+xy21+x21+x20y+x19y2+x18y3+x17y4+x16y5+x15y6+x14y7+x13y8+x12y9+x11y10+x10y11+x9y12+x8y13+x7y14+x6y15+x5y16+x4y17+x3y18+x2y19+xy20+x20+x19y+x18y2+x17y3+x16y4+x15y5+x14y6+x13y7+x12y8+x11y9+x10y10+x9y11+x8y12+x7y13+x6y14+x5y15+x4y16+x3y17+x2y18+xy19+x19+x18y+x17y2+x16y3+x15y4+x14y5+x13y6+x12y7+x11y8+x10y9+x9y10+x8y11+x7y12+x6y13+x5y14+x4y15+x3y16+x2y17+xy18+x18+x17y+x16y2+x15y3+x14y4+x13y5+x12y6+x11y7+x10y8+x9y9+x8y10+x7y11+x6y12+x5y13+x4y14+x3y15+x2y16+xy17+x17+x16y+x15y2+x14y3+x13y4+x12y5+x11y6+x10y7+x9y8+x8y9+x7y10+x6y11+x5y12+x4y13+x3y14+x2y15+xy16+x16+x15y+x14y2+x13y3+x12y4+x11y5+x10y6+x9y7+x8y8+x7y9+x6y10+x5y11+x4y12+x3y13+x2y14+x15+x14y+x13y2+x12y3+x11y4+x10y5+x9y6+x8y7+x7y8+x6y9+x5y10+x4y11+x3y12+x2y13+x14+x13y+x12y2+x11y3+x10y4+x9y5+x8y6+x7y7+x6y8+x5y9+x4y10+x3y11+x2y12+x13+x12y+x11y2+x10y3+x9y4+x8y5+x7y6+x6y7+x5y8+x4y9+x3y10+x2y11+x12+x11y+x10y2+x9y3+x8y4+x7y5+x6y6+x5y7+x4y8+x3y9+x2y10+x11+x10y+x9y2+x8y3+x7y4+x6y5+x5y6+x4y7+x3y8+x2y9+x9y+x8y2+x7y3+x6y4+x5y5+x4y6+x3y7+x2y8+x8y+x7y2+x6y3+x5y4+x4y5+x3y6+x2y7+x7y+x6y2+x5y3+x4y4+x3y5+x2y6+x5y2+x4y3+x3y4+x2y5+x4y2+x3y3+x2y4+x3y2+x2y3+x2y2;
     4243map phi = R,x+2*y+x^2+x*y+y^2+x^3,3*x+5*y+x^2*y+y^2;
     4244poly h2 = phi(h);
     4245int t = timer;
     4246list l1 = integralBasis(h2, 2, "atOrigin");
     4247timer - t;
     4248

  • Singular/LIB/integralbasis.lib

    r8bdde5 r0c2a2e  
    10991099
    11001100    // We compute the integral basis.
    1101     list ib = ibNonRatXY(fCh, px, py, locBasis);
     1101    list ib = ibNonRatXY(fCh, px, py, locBasis, optimize);
    11021102
    11031103    if(li > 0)
    … …  
    11661166
    11671167
    1168 ///////////////////////////////////////////////////////////////////////
    1169 // Computation of Puiseux segments and classes in each segment
    1170 // Finding classes should be easy, just looking at the slope of each segment
    1171 
    1172 // Returns a list of two lists
    1173 // In the first list, the polynomials of each segment
    1174 // In the second list, the classes of expansions corresponding to each segment
    1175 
    1176 proc in_array(list l, def elem)
     1168//////////////////////////////////////////////////////////////////////////////
     1169//
     1170//    OPTIMIZATION APPROACH TO COMPUTE THE INTEGRAL BASIS AS EXPLAINED IN THE
     1171//    PAPER by Boehm, Decker, Laplagne & Pfister
     1172//
     1173//    In most cases this is faster than locLocAlgorithm
     1174//
     1175/////////////////////////////////////////////////////////////////////////////
     1176
     1177// Auxiliary function
     1178static proc in_array(list l, def elem)
    11771179{
    11781180  int i;
    … …  
    11881190}
    11891191
     1192///////////////////////////////////////////////////////////////////////
     1193// Computation of Puiseux segments and classes in each segment
     1194//
    11901195// We compute the Puiseux segments using Hensel lifting.
    11911196// The segments are order according to the initial exponent from
    … …  
    11931198// The classes of Puiseux expansions are grouped matching the
    11941199// Puiseux segments.
     1200//
     1201// Returns a list of two lists
     1202// In the first list, the polynomials of each segment
     1203// In the second list, the classes of expansions corresponding to each segment
     1204
    11951205
    11961206proc getSegments(poly f, list classes, list slopes, int globOrder)
    … …  
    12341244    dbprint(dbg, "----Computing Puiseux Segments...");
    12351245    fSegment = henselSegments(fLoc, maxOrder);
     1246    "Number of segments: ", size(fSegment);
     1247   
     1248    list pRecomputed;
     1249    list nPoly;
    12361250   
    12371251    //"Check output of henselSegments";
    … …  
    12411255      // We compute the newton poly of each polynomial to find out the
    12421256      // corresponding slope
    1243       list l = newtonpoly(fSegment[i]);
    1244       slopeSegment = number(l[2][1]) / number(l[1][2]);
     1257      nPoly = newtonpoly(fSegment[i]);
     1258      slopeSegment = number(nPoly[2][1]) / number(nPoly[1][2]);
    12451259
    12461260      // We look for all classes with same slope as each segment
    … …  
    12481262      for(j = 1; j <= size(classes); j++)
    12491263      {
    1250         if(slopes[j][1] * l[1][2] == slopes[j][2] * l[2][1])  // same slope
     1264        if(slopes[j][1] * nPoly[1][2] == slopes[j][2] * nPoly[2][1])  // same slope
    12511265        {
    12521266          // If slopes coincide we add the class to the list
    … …  
    12541268        }
    12551269      } 
     1270     
     1271      // If the number of classes in the segments is larger than one,
     1272      // we need to recompute expansions up to degree maxOrder so that
     1273      // the factors can be computed correctly.
     1274      if(size(classesNew) > 1)
     1275      {
     1276        // puiseux(poly f, int maxDeg, int atOrigin)
     1277        list lRecomputed = puiseux(fSegment[i], globOrder, 1);
     1278        debug_log_intbas(3, "------Puiseux expansions recomputed: ");
     1279        classesNew = getClasses(lRecomputed);
     1280      }
    12561281     
    12571282      // We put everything together
    … …  
    13021327 
    13031328  // We get all the possible polynomials, including the full polys
     1329  //"Building factors from classes...";
    13041330  list bF = buildFactors(classes);
    13051331 
    … …  
    13181344      }
    13191345    }
    1320     // We compute the full polynomial up to the integrality exponent
    1321     pGround = buildPolyGroundXRootClass(classes[i], degExpand);
    1322     if(in_array(polySet, pGround) == 0)
    1323     {
    1324       polySet[size(polySet)+1] = pGround;
     1346   
     1347    // We compute the full polynomial up to the integrality exponent.
     1348    // When there is only one class we use the polynomial obtained from Hensel lifting.
     1349    // If there are more than 1 class in the segment we need to compute the
     1350    // polynomials from the Puiseux expansions (or use the blocks, but not implemented).
     1351    if(size(classes) == 1)
     1352    {
     1353      polySet[size(polySet)+1] = f;
     1354    } else
     1355    {
     1356      pGround = buildPolyGroundXRootClass(classes[i], degExpand);
     1357      if(in_array(polySet, pGround) == 0)
     1358      {
     1359        polySet[size(polySet)+1] = pGround;
     1360      }
    13251361    }
    13261362  }
    … …  
    13581394}
    13591395 
    1360 // For each degree up to the degree of f, with compute the best
     1396// For each degree up to the degree of f, we compute the best
    13611397// possible product of factors by checking exhaustively all possible
    13621398// combinations such that the product has the correct degree.
     1399// Used for the direct approach.
    13631400proc ibExhaustive(list polySet, intvec degPolySet, list ordExpAtPoly, int degF)
    13641401{
    … …  
    14861523proc optiBase(list ibSeg, list segmentSlope)
    14871524{
    1488   int i;
     1525  int i, k;
    14891526  int s = size(ibSeg);
    14901527 
    14911528  if(s == 1) {return(ibSeg[1]);}
    1492   // Implement the optimization algorithm...
    1493   // Use the code already implemented??
    14941529 
    1495   //for(i = s-1; i >= 1; i++)
    1496   //{
    1497   if(s == 2)
    1498   {
    1499     return(optiTwoSets(ibSeg[1], segmentSlope[1][2], ibSeg[2], segmentSlope[2][2]));
    1500   }
     1530  // We initialize the recursion step
     1531  list ibSegRec = ibSeg[s];
     1532  list segmentSlopeRec = segmentSlope[s]; // Should be fixed to contain only slopes
     1533 
     1534  for(k = s-1; k >= 1; k--)
     1535  {
     1536    ibSegRec = optiTwoSets(ibSeg[k], segmentSlope[k][2], ibSegRec, segmentSlopeRec[2]);
     1537    segmentSlopeRec = segmentSlope[k];
     1538  }
     1539  return(ibSegRec);
    15011540}
    15021541
    … …  
    15041543proc optiTwoSets(list ib1, number slope1, list ib2, number slope2)
    15051544{
    1506   "optiTwoSets begin...";
     1545  //"optiTwoSets begin...";
    15071546 
    15081547  int s1 = size(ib1)-1;  // max degree
    … …  
    15111550  intvec cBest;
    15121551  int d;
     1552 
     1553  //"s1 + s2 = ", s1 + s2;
    15131554 
    15141555  number o, o1, o2;
    … …  
    15191560  out[1] = list(0, intvec(0,0), 1);
    15201561  // Best element of degree d
     1562 
     1563  int m1, m2;
    15211564  for(d = 1; d <= s1 + s2; d++)
    15221565  {
    15231566    oBest = 0;
    1524     for(c1 = 0; c1 <= s1; c1++)
    1525     {
    1526       if(c1 <= d)
    1527       {
    1528         c2 = d - c1;
    1529         if(c2 <= s2)
     1567    m1 = max(0, d - s2);
     1568    m2 = min(s1, d);
     1569    for(c1 = m1; c1 <= m2; c1++)
     1570    {
     1571      c2 = d - c1;
     1572      if(c2 <= s2)
     1573      {
     1574        // The first component of ib1 contains the orders of the best polynomial of each degree
     1575        if(c1 < s1)
    15301576        {
    1531           // The first component of ib1 contains the orders of the best polynomial of each degree
    1532           if(c1 < s1)
     1577          o1 = c2 * slope1 + ib1[c1+1][1];  // The valuation of an expansion from the first group
     1578        }
     1579        if(c2 < s2)
     1580        {
     1581          o2 = c1 * slope1 + ib2[c2+1][1];  // The valuation of an expansion from the second group
     1582        }
     1583        if(c1 == s1)
     1584        {
     1585          o = o2;
     1586        } else
     1587        {
     1588          if(c2 == s2)
    15331589          {
    1534             o1 = c2 * slope1 + ib1[c1+1][1];  // The valuation of an expansion from the first group
    1535           }
    1536           if(c2 < s2)
    1537           {
    1538             o2 = c1 * slope1 + ib2[c2+1][1];  // The valuation of an expansion from the second group
    1539           }
    1540           if(c1 == s1)
    1541           {
    1542             o = o2;
     1590            o = o1;
    15431591          } else
    15441592          {
    1545             if(c2 == s2)
    1546             {
    1547               o = o1;
    1548             } else
    1549             {
    1550               if(o1 < o2) {o = o1;}
    1551               else {o = o2;}
    1552             }
    1553           }
    1554           //"c1 = ", c1, "c2 = ", c2, "s1 = ", s1, "s2 = ", s2, "o = ", o;
    1555          
    1556           if(o > oBest)
    1557           {
    1558             oBest = o;
    1559             cBest = intvec(c1, c2);
    1560             polyBest = ib1[c1+1][3] * ib2[c2+1][3];
     1593            if(o1 < o2) {o = o1;}
     1594            else {o = o2;}
    15611595          }
    15621596        }
    1563       }
    1564     }
     1597        //"c1 = ", c1, "c2 = ", c2, "s1 = ", s1, "s2 = ", s2, "o = ", o;
     1598       
     1599        if(o > oBest)
     1600        {
     1601          oBest = o;
     1602          cBest = intvec(c1, c2);
     1603        }
     1604      }
     1605    }
     1606    polyBest = ib1[cBest[1]+1][3] * ib2[cBest[2]+1][3];
    15651607    out[size(out) + 1] = list(oBest, cBest, polyBest);
    15661608  }
     1609  //"optiTwoSets finished...";
     1610 
    15671611  return(out);
    15681612}
    15691613   
    15701614 
    1571 
    1572 
    15731615
    15741616////////////////////////////////////////////////////////////////////////
    … …  
    16451687  list bestElem = vS[2];
    16461688
    1647   "MClass";
    1648   MClass;
     1689  //"MClass";
     1690  //MClass;
    16491691
    16501692  //classes;
     1693  //~;
    16511694
    16521695
    … …  
    16821725
    16831726  // Optimization algorithm from scratch
     1727 
     1728  if((pardeg(minpoly) > 0) and (optimize == 1))
     1729  {
     1730    "Optimization not implemented for conjugated singularities. Using chinese remainder algorithm for this component.";
     1731    optimize = 0;
     1732  }
     1733
    16841734  if(optimize == 1)
    16851735  {
     1736    "Using optimization algorithm";
    16861737    // 1) Compute Puiseux sets (we use segments for easier implementation)
    16871738    list pSegmentFull = getSegments(f, classes, slopes, degExpand);
    … …  
    17251776      IOut[k+1] = condu * var(2)^(k);
    17261777    }
    1727     for(k = 1; k <= size(ib) - 2; k++)
    1728     {
    1729       IOut[size(IOut)+1] = ib[k+1][3] * var(1)^(intExp - int(ib[k+1][1])) * unitPoly;
     1778    //~;
     1779    for(k = 1; k <= size(ib) - 1; k++)
     1780    {
     1781      IOut[size(IOut)+1] = ib[k][3] * var(1)^(intExp - int(ib[k][1])) * unitPoly;
    17301782    }
    17311783    "step 4 finished";
    … …  
    18291881
    18301882    debug_log_intbas(3, "------Time for ordFull: ", timer - tt);
    1831     "ordFull is wrong for us...";
    1832     //~;
    18331883   
    18341884    tt = timer;
    … …  
    18611911    for(i = 1; i <= nClasses; i++)
    18621912    {
    1863       mdmTemp = maxDegMerge(n, i, locBasis, d0, M, MSelf, md, classes, I2Lifted, bestElem, ordsFull, ordsBest, optimize);
     1913      mdmTemp = maxDegMerge(n, i, locBasis, d0, M, MSelf, md, classes, I2Lifted, bestElem, ordsFull, ordsBest);
    18641914      //"mdmTemp = ", mdmTemp;
    18651915      mdm = maxInt(mdm, mdmTemp);
    … …  
    20032053    tt = timer;
    20042054
    2005     if(optimize == 0)
    2006     {
    2007       for(i = 1; i <= nClasses; i++)
    2008       {
    2009         // Computing ~A(i)_new... (using modified chinese remainder)
    2010         IdOut[i] = locLocAlgorithm(n, i, locBasis, d0, M, MSelf, md, classes, I2LiftedFull, mdm, bestElem, ordsFull, ordsBest);
    2011         bases = bases + list(IdOut[i]);
    2012       }
    2013       debug_log_intbas(3, "------Time integral basis for all branches: ", timer - tt);
    2014 
    2015       if(nClasses > 1)
    2016       {
    2017         // We add some extra polynomials to the integral basis
    2018         // which will simplify the computations
    2019         list enzNum;
    2020         list enzDen;
    2021         poly elemNum;
    2022         number euNum;
    2023         int eu;
    2024         intvec classesInd = 1:nClasses;
    2025         number ordY = ordAtClassesPol(var(2), classes, classesInd);
    2026         for(k = 0; k < n - deg(I2Lifted[1], vy); k++)
    2027         {
    2028           elemNum = var(2)^k;
    2029           euNum = ordY * k;
    2030           eu = int(euNum);
    2031           enzNum = enzNum + list(elemNum);
    2032           enzDen = enzDen + list(var(1)^(eu));
    2033         }
    2034         list enzOut = comDen(enzNum, enzDen);
    2035         ideal enzId = enzOut[2];
    2036         for(k=1; k <=size(enzOut[1]); k++)
    2037         {
    2038           enzId[k+1] = enzOut[1][k];
    2039         }
    2040         list enzOutList = enzId, enzId[1];
    2041         bases = bases + list(enzOutList);
    2042       }
    2043 
    2044 
    2045       tt = timer;
    2046       if(size(bases) > 1)
    2047       {
    2048         // We compute the local integral basis at the origin
    2049         dbprint(dbg, "--Merging the integral bases for the branches...");
    2050         outNormal = mergeBases(bases);
    2051 
    2052         // Check the correct degree for reduction
    2053         //outNormal[1] = reduce(outNormal[1], groebner(x^mdm));
    2054 
    2055 
    2056         //dbprint(dbg, "----Computing the groebner basis...", outNormal[1]);
    2057         //dbprint(dbg, "----Denominator:", outNormal[2]);
    2058         outNormal[1] = groebner(outNormal[1]);
    2059         dbprint(dbg, "----Merging finished");
    2060         poly fLoc = 1;
    2061         tmp_var_mdm=std(x^mdm);
    2062         for(k = 2; k <= size(I2LiftedFull); k++)
    2063         {
    2064           fLoc = reduce(fLoc * reduce(I2LiftedFull[k], tmp_var_mdm), tmp_var_mdm);
    2065         }
    2066         int needGroeb = 1;
    2067         list intBas = normToInt(fLoc, outNormal[1], outNormal[2], needGroeb);
    2068       } else
    2069       {
    2070         list intBas;
    2071         intBas[1] = bases[1][1];
    2072         intBas[2] = bases[1][2];
    2073       }
    2074       debug_log_intbas(3, "------Time for merging bases: ", timer - tt);
    2075 
    2076 
    2077       // We add the unit outside the origin
    2078       for(k = 0; k < deg(I2LiftedFull[1], vy); k++)
    2079       {
    2080         IOut[k+1] = intBas[2]*var(2)^(k);
    2081       }
    2082       for(k = 1; k <= size(intBas[1]); k++)
    2083       {
    2084         IOut[size(IOut)+1] = intBas[1][k]  * I2Lifted[1];
    2085       }
    2086     } else {
    2087       // Computation of the smallest slope
    2088       list sl = getSlope(slopes);
    2089       slN = sl[1];
    2090       slD = sl[2];
    2091       bestSl = sl[3];
    2092       intvec vxy = (slD, slN);
    2093 
    2094       IOut = optimAlgorithm(n, locBasis, d0, d1, vxy, M, MSelf, md, classes, I2Lifted, bestElem, optimize);
     2055    for(i = 1; i <= nClasses; i++)
     2056    {
     2057      // Computing ~A(i)_new... (using modified chinese remainder)
     2058      IdOut[i] = locLocAlgorithm(n, i, locBasis, d0, M, MSelf, md, classes, I2LiftedFull, mdm, bestElem, ordsFull, ordsBest);
     2059      bases = bases + list(IdOut[i]);
     2060    }
     2061    debug_log_intbas(3, "------Time integral basis for all branches: ", timer - tt);
     2062
     2063    if(nClasses > 1)
     2064    {
     2065      // We add some extra polynomials to the integral basis
     2066      // which will simplify the computations
     2067      list enzNum;
     2068      list enzDen;
     2069      poly elemNum;
     2070      number euNum;
     2071      int eu;
     2072      intvec classesInd = 1:nClasses;
     2073      number ordY = ordAtClassesPol(var(2), classes, classesInd);
     2074      for(k = 0; k < n - deg(I2Lifted[1], vy); k++)
     2075      {
     2076        elemNum = var(2)^k;
     2077        euNum = ordY * k;
     2078        eu = int(euNum);
     2079        enzNum = enzNum + list(elemNum);
     2080        enzDen = enzDen + list(var(1)^(eu));
     2081      }
     2082      list enzOut = comDen(enzNum, enzDen);
     2083      ideal enzId = enzOut[2];
     2084      for(k=1; k <=size(enzOut[1]); k++)
     2085      {
     2086        enzId[k+1] = enzOut[1][k];
     2087      }
     2088      list enzOutList = enzId, enzId[1];
     2089      bases = bases + list(enzOutList);
     2090    }
     2091
     2092
     2093    tt = timer;
     2094    if(size(bases) > 1)
     2095    {
     2096      // We compute the local integral basis at the origin
     2097      dbprint(dbg, "--Merging the integral bases for the branches...");
     2098      outNormal = mergeBases(bases);
     2099
     2100      // Check the correct degree for reduction
     2101      //outNormal[1] = reduce(outNormal[1], groebner(x^mdm));
     2102
     2103
     2104      //dbprint(dbg, "----Computing the groebner basis...", outNormal[1]);
     2105      //dbprint(dbg, "----Denominator:", outNormal[2]);
     2106      outNormal[1] = groebner(outNormal[1]);
     2107      dbprint(dbg, "----Merging finished");
     2108      poly fLoc = 1;
     2109      tmp_var_mdm=std(x^mdm);
     2110      for(k = 2; k <= size(I2LiftedFull); k++)
     2111      {
     2112        fLoc = reduce(fLoc * reduce(I2LiftedFull[k], tmp_var_mdm), tmp_var_mdm);
     2113      }
     2114      int needGroeb = 1;
     2115      list intBas = normToInt(fLoc, outNormal[1], outNormal[2], needGroeb);
     2116    } else
     2117    {
     2118      list intBas;
     2119      intBas[1] = bases[1][1];
     2120      intBas[2] = bases[1][2];
     2121    }
     2122    debug_log_intbas(3, "------Time for merging bases: ", timer - tt);
     2123
     2124
     2125    // We add the unit outside the origin
     2126    for(k = 0; k < deg(I2LiftedFull[1], vy); k++)
     2127    {
     2128      IOut[k+1] = intBas[2]*var(2)^(k);
     2129    }
     2130    for(k = 1; k <= size(intBas[1]); k++)
     2131    {
     2132      IOut[size(IOut)+1] = intBas[1][k]  * I2Lifted[1];
    20952133    }
    20962134  }
    … …  
    25762614          if(dMP > 1)
    25772615          {
     2616            dbprint(dbg, "Case where the base field is an algebraic extension of Q.");
    25782617            poly mp = composePolys(minPolys);
    25792618            poly mpX = subst(mp, var(2), var(1));
    … …  
    25912630              fFr = fFr + subst(rel[kk], var(1), par(1)) * CC[1, kk];
    25922631            }
    2593             dbprint(dbg, "Case where the base field is an algebraic extension of Q.");
     2632            polyGround = buildPolyGround(fFr, den);
    25942633          } else
    25952634          {
    25962635            // No extension in the base ring
    2597             poly fFB = buildPolyFrac(fF);
     2636            list polyGroundExt = buildPolyGroundXRoot(fF, den, d);
     2637            poly fFB = polyGroundExt[1];
     2638            kill polyGroundExt;
    25982639            setring R;
    25992640            fFr = imap(S, fFB);
     2641            polyGround = list(fFr, 1);
    26002642          }
    26012643
    … …  
    26072649          //polyGround = buildPolyGroundXRoot(fFr, den);
    26082650          //polyGround = list(fFr, 1);
    2609           polyGround = buildPolyGround(fFr, den);
     2651          //polyGround = buildPolyGround(fFr, den);
    26102652        } else
    26112653        {
    … …  
    29122954    if(deg(f2, (1,0,0)) > 0)
    29132955    {
    2914       //"f2", f2;
     2956      "Polynomial is not over the ground field";
     2957      "f2", f2;
     2958      ~;
    29152959      ERROR("Polynomial is not over the ground field");
    29162960      gfCheck = 0;
    … …  
    29683012  if(deg(f2, (1,0,0)) > 0)
    29693013  {
    2970     //"f2", f2;
     3014    "f2", f2;
     3015    ~;
    29713016    ERROR("Polynomial is not over the ground field");
    29723017    gfCheck = 0;
    … …  
    31633208  poly differ;
    31643209  number ordDiff;
     3210 
     3211  list l;
     3212 
    31653213  for(i = 1; i <= size(classes); i++)
    31663214  {
    … …  
    31853233          {
    31863234            // Finally, if they have the same order we use ordAtPol.
    3187             "Build fac";
     3235            //"Build fac";
    31883236            //~;
    3189             list l = buildPolyGroundXRoot(f(j), den(j), -1);
     3237            l = buildPolyGroundXRoot(f(j), den(j), -1);
    31903238            M[i, j] = ordAtPol(l[1], list(f(i), den(i)));
    3191             "ord of class i at poly corresponding to class j:", M[i,j];
     3239            //"ord of class i at poly corresponding to class j:", M[i,j];
    31923240            //~;
    31933241          }
    … …  
    39864034// [Example 4]
    39874035
    3988 static proc maxDegMerge(int n, int cl, int locBasis, int d0, matrix M, list MSelf, intvec md, list classes, list I2Lifted, list bestElem, matrix ordsFull, list ordsBest, int optimize);
     4036static proc maxDegMerge(int n, int cl, int locBasis, int d0, matrix M, list MSelf, intvec md, list classes, list I2Lifted, list bestElem, matrix ordsFull, list ordsBest);
    39894037{
    39904038  //"START - newib.lib - maxDegMerge - 4";
    … …  
    40304078    // The first term containing the locloc unit
    40314079    expBaseFactor = expUnit;
    4032 
    4033     "Check here";
    4034     //~;
    4035        
    4036     if(optimize == 0)
    4037     {
    4038       // If the order of hi is not integer, we add a factor.
    4039       // This is used in the Chinese remainder strategy
    4040       if(int(denominator(expBaseFactor)) > 1)
    4041       {
    4042         // We add the factor corresponding to the full expansions
    4043         list fullPoly;
    4044         for(j = 1; j <= nClasses; j++)
     4080       
     4081    // If the order of hi is not integer, we add a factor.
     4082    // This is used in the Chinese remainder strategy
     4083    if(int(denominator(expBaseFactor)) > 1)
     4084    {
     4085      // We add the factor corresponding to the full expansions
     4086      list fullPoly;
     4087      for(j = 1; j <= nClasses; j++)
     4088      {
     4089        fullPoly[j] = I2Lifted[j+1];
     4090      }
     4091
     4092      int whichElem = 0;
     4093      int den1;
     4094      int den2;
     4095
     4096      intvec degsCl;
     4097      for(i = 1; i <= size(MSelf); i++)
     4098      {
     4099        degsCl[i] = size(MSelf[i]) + 1;
     4100      }
     4101      intvec chVec = 0:size(MSelf);
     4102
     4103      den1 = 1;
     4104      den2 = int(denominator(expBaseFactor));
     4105
     4106      number ordNewFactor;
     4107
     4108      number oAP;
     4109      number oAPT;
     4110
     4111      while(den1 mod den2 != 0)
     4112      {
     4113        chVec = nextSummand(degsCl, chVec, cl);
     4114        ordNewFactor = 0;
     4115        for(i = 1; i <= nClasses; i++)
    40454116        {
    4046           fullPoly[j] = I2Lifted[j+1];
    4047         }
    4048 
    4049         int whichElem = 0;
    4050         int den1;
    4051         int den2;
    4052 
    4053         intvec degsCl;
    4054         for(i = 1; i <= size(MSelf); i++)
    4055         {
    4056           degsCl[i] = size(MSelf[i]) + 1;
    4057         }
    4058         intvec chVec = 0:size(MSelf);
    4059 
    4060         den1 = 1;
    4061         den2 = int(denominator(expBaseFactor));
    4062 
    4063         number ordNewFactor;
    4064 
    4065         number oAP;
    4066         number oAPT;
    4067 
    4068         while(den1 mod den2 != 0)
    4069         {
    4070           chVec = nextSummand(degsCl, chVec, cl);
    4071           ordNewFactor = 0;
    4072           for(i = 1; i <= nClasses; i++)
     4117          if(i != cl)
    40734118          {
    4074             if(i != cl)
     4119            k = chVec[i];
     4120            if(k > 0)
    40754121            {
    4076               k = chVec[i];
    4077               if(k > 0)
     4122              if(k <= size(bestElem[i]))
    40784123              {
    4079                 if(k <= size(bestElem[i]))
    4080                 {
    4081                   oAP = ordsBest[cl][i][k];
    4082                   ordNewFactor = ordNewFactor + oAP;
    4083                 } else
    4084                 {
    4085                   oAP = number(ordsFull[i, cl]);
    4086                   ordNewFactor = ordNewFactor + oAP;
    4087                 }
     4124                oAP = ordsBest[cl][i][k];
     4125                ordNewFactor = ordNewFactor + oAP;
     4126              } else
     4127              {
     4128                oAP = number(ordsFull[i, cl]);
     4129                ordNewFactor = ordNewFactor + oAP;
    40884130              }
    40894131            }
    40904132          }
    4091           den1 = int(denominator(ordNewFactor));
    40924133        }
    4093 
    4094         poly newFactor = 1;
    4095         int denNewFactor = den1;
    4096         int powFactor = 0;
    4097 
    4098         int nnE = int(numerator(expBaseFactor));
    4099         int ddE = int(denominator(expBaseFactor));
    4100         int nnY = int(numerator(ordNewFactor));
    4101 
    4102         while(ddE > 1)
    4103         {
    4104           powFactor = powFactor + 1;
    4105           expBaseFactor = expUnit + powFactor * ordNewFactor;
    4106           ddE = int(denominator(expBaseFactor));
    4107         }
     4134        den1 = int(denominator(ordNewFactor));
     4135      }
     4136
     4137      poly newFactor = 1;
     4138      int denNewFactor = den1;
     4139      int powFactor = 0;
     4140
     4141      int nnE = int(numerator(expBaseFactor));
     4142      int ddE = int(denominator(expBaseFactor));
     4143      int nnY = int(numerator(ordNewFactor));
     4144
     4145      while(ddE > 1)
     4146      {
     4147        powFactor = powFactor + 1;
     4148        expBaseFactor = expUnit + powFactor * ordNewFactor;
     4149        ddE = int(denominator(expBaseFactor));
    41084150      }
    41094151    }
    … …  
    56805722
    56815723
    5682 //////////////////////////////////////////////////////////////////////////////
    5683 //
    5684 //    OPTIMIZATION APPROACH TO COMPUTE THE INTEGRAL BASIS AS EXPLAINED IN THE
    5685 //    PAPER by Boehm, Decker, Laplagne & Pfister
    5686 //
    5687 //    In most cases this is faster than locLocAlgorithm
    5688 //
    5689 /////////////////////////////////////////////////////////////////////////////
    5690 
    5691 // Input:  n is the degree of f
    5692 // Output: local contribution to the integral basis at the origin
    5693 // Algorithm: incorporates the local contributions to the integral basis
    5694 //            one by one in an optimized way
    5695 // [Example 4]
    5696 
    5697 proc optimAlgorithm(int n, int locBasis, int d0, int d1, intvec vxy, matrix M, list MSelf, intvec md, list classes, list I2Lifted, list bestElem, int optimize)
    5698 {
    5699   dbprint(dbg, "START - intbasoptim.lib - optimAlgorithm - 1"); 
    5700   int i, j, k;
    5701   int t = timer;
    5702   list b;
    5703 
    5704   number maxExp;
    5705   poly num;
    5706   list ibNum;
    5707   list ibDen;
    5708   poly p, p1, p2;
    5709   t = timer;
    5710   intvec elem;
    5711 
    5712   int dbg = printlevel - voice + 5;
    5713 
    5714   intvec vy = (0,1);
    5715   intvec vx = (1,0);
    5716   poly x = var(1);
    5717   poly y = var(2);
    5718   int slD = vxy[1];
    5719 
    5720   dbprint(dbg, "Computing all the elements of the integral basis...");
    5721  
    5722   list bestElems;
    5723   if(optimize == 1){
    5724     bestElems = optiBest(M, MSelf, md, n);
    5725     //for(k = 1; k < d0; k++)
    5726     //{
    5727       //"elem: ", bestElems[k+1][2], " - ", bestElems[k+1][1];
    5728     //}
    5729   } else {
    5730     //"optimize = 0";
    5731     bestElems[1] = list();
    5732     for(k = 1; k < d0; k++)
    5733     {
    5734       bestElems[k+1] = ibElement(M, MSelf, k, md);
    5735       //"elem: ", bestElems[k+1][2], " - ", bestElems[k+1][1];
    5736     }
    5737   }
    5738  
    5739   "optiBest";
    5740   //~;
    5741  
    5742   if(locBasis == 0)
    5743   {
    5744     for(k = 1; k <= d1+1; k++)
    5745     {
    5746       ibNum[k] = y^(k-1);
    5747       ibDen[k] = 1;
    5748     }
    5749   } else
    5750   {
    5751     ibNum[1] = 1;
    5752     ibDen[1] = 1;
    5753   }
    5754 
    5755   poly pTemp;
    5756   for(k = 1; k < d0; k++)
    5757   {
    5758    
    5759     if(size(classes) > 1)
    5760     {
    5761       maxExp = bestElems[k+1][1];
    5762       elem = bestElems[k+1][2];
    5763     } else
    5764     {
    5765       maxExp = MSelf[1][k];
    5766       elem[1] = k;
    5767     }
    5768     if(dbg >= 3)
    5769     {
    5770       "element ", k, ": ", elem;
    5771     }
    5772 
    5773     num = 1;
    5774     for(i = 1; i <= size(elem); i++)
    5775     {
    5776       if((elem[i] > 0) and (elem[i] < md[i]))
    5777       {
    5778         pTemp = jet(bestElem[i][elem[i]], int(maxExp), vx);
    5779         pTemp = jet(pTemp, int(maxExp*slD), vxy);
    5780         num = num * pTemp;
    5781         num = jet(num, int(maxExp), vx);
    5782         num = jet(num, int(maxExp*slD), vxy);
    5783       }
    5784     }
    5785     for(i = 1; i <= size(elem); i++)
    5786     {
    5787       if(elem[i] == md[i])
    5788       {
    5789         pTemp = jet(I2Lifted[i + 1], int(maxExp), vx);
    5790         pTemp = jet(pTemp, int(maxExp*slD), vxy);
    5791         num = prodMod(num, pTemp, int(maxExp) + 1);
    5792         num = jet(num, int(maxExp), vx);
    5793         num = jet(num, int(maxExp*slD), vxy);
    5794       }
    5795     }
    5796 
    5797     p1 = jet(num, int(maxExp*slD), vxy);
    5798     if((locBasis == 1) or (d1 == 0))
    5799     {
    5800       ibNum[k+1] = p1;
    5801       ibDen[k+1] = x^int(maxExp);
    5802     } else
    5803     {
    5804       p2 = jet(I2Lifted[1], int(maxExp-1), vx);
    5805       p = jet(p1*p2, int(maxExp-1), vx);
    5806       ibNum[k+d1+1] = p;
    5807       ibDen[k+d1+1] = x^int(maxExp);
    5808     }
    5809   }
    5810   ///"Construction finished.";
    5811  
    5812   ideal IOut;
    5813   if(locBasis == 1)
    5814   {
    5815     for(i = 1; i<=d0; i++)
    5816     {
    5817       IOut[i] = (ibNum[i] * ibDen[d0]) / ibDen[i];
    5818     }
    5819   } else
    5820   {
    5821     for(i = 1; i<=n; i++)
    5822     {
    5823       IOut[i] = (ibNum[i] * ibDen[n]) / ibDen[i];
    5824     }
    5825   }
    5826   if(dbg >= 2)
    5827   {
    5828     "Time basis construction: ", timer - t;
    5829   }
    5830   return(IOut);
    5831 }
    5832 
    5833 // Efficient implementation of the optimization problem
    5834 // Output: A list, with the maximal valuation for each degree, and how
    5835 //         many elements to take in each conjugacy class of expansions
    5836 // [Example 1]
    5837 proc optiBest(matrix M, list MSelf, intvec md, int k)
    5838 {
    5839   //"START - newib.lib - optiBest - 1";
    5840   int i, j;
    5841   int nClasses = size(md);
    5842  
    5843   list beNew, beOld;
    5844   intvec choice;
    5845   number val;
    5846   int totNew = md[nClasses];
    5847   int totOld;
    5848   for(i = 1; i <= totNew + 1; i++)
    5849   {
    5850     choice[1] = i - 1;
    5851     if(i == 1){
    5852       val = 0;
    5853     } else {
    5854       if(i == md[nClasses] + 1){
    5855         val = -1;
    5856       } else {
    5857         val = MSelf[nClasses][i-1];
    5858       }
    5859     }
    5860     beOld[i] = list(val, choice);
    5861   }
    5862   //"M: "; M;
    5863   //"beOld"; beOld;
    5864   int best1, best2;
    5865   number bestV;
    5866   number val1, val2;
    5867   for(j = size(md) - 1; j>=1; j--)
    5868   {
    5869     totOld = totNew;
    5870     totNew = totNew + md[j];
    5871     for(i = 1; i <= totNew + 1; i++)
    5872     {
    5873       bestV = -1;
    5874       for(k = 1; k <= md[j] + 1; k++)
    5875       {
    5876         if(((k - 1) <= (i-1)) and (i - k <= totOld))
    5877         {
    5878           //"MSelf, i, j, k"; MSelf, i, j, k;
    5879           //"a, b: ", k-1, i-k;
    5880           if(k > 1)
    5881           {
    5882             if(k-1 == md[j])
    5883             {
    5884               val1 = -1;
    5885             } else
    5886             {
    5887               val1 = number(MSelf[j][k-1] + (i-k)*M[j, j+1]);
    5888             }
    5889           } else
    5890           {
    5891             val1 = number((i-k)*M[j, j+1]);
    5892           }
    5893           if(beOld[i-k+1][2] > -1)
    5894           {
    5895             val2 = number((k-1)*M[j, j+1] + beOld[i-k+1][1]);
    5896           } else
    5897           {
    5898             val2 = -1;
    5899           }
    5900           if(((val1 < val2) and (val1 > -1)) or (val2 == -1))
    5901           {
    5902             val = val1;
    5903           } else {
    5904             val = val2;
    5905           }
    5906           //"val1, val2: ", val1, val2;
    5907           //"val: ", val;
    5908           if((val > bestV) or (val == -1))
    5909           {
    5910             //"i-k", i-k;
    5911             choice = k-1, beOld[i-k+1][2];
    5912             beNew[i] = list(val, choice);
    5913           }
    5914         }
    5915       }
    5916       //"elem: ", beNew[i][2], " - ", beNew[i][1];
    5917       //"i: ", i;
    5918       //"beNew[i]", beNew[i];
    5919     }
    5920     beOld = beNew;
    5921   }
    5922   return(beNew);
    5923 }
    59245724
    59255725/////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.